Reference documentation for deal.II version Git 7e76c9f1da 2019-09-19 09:07:25 +0200
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
trilinos_sparse_matrix.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_trilinos_sparse_matrix_h
17 # define dealii_trilinos_sparse_matrix_h
18 
19 
20 # include <deal.II/base/config.h>
21 
22 # ifdef DEAL_II_WITH_TRILINOS
23 
24 # include <deal.II/base/index_set.h>
25 # include <deal.II/base/subscriptor.h>
26 
27 # include <deal.II/lac/exceptions.h>
28 # include <deal.II/lac/full_matrix.h>
29 # include <deal.II/lac/trilinos_epetra_vector.h>
30 # include <deal.II/lac/trilinos_index_access.h>
31 # include <deal.II/lac/trilinos_tpetra_vector.h>
32 # include <deal.II/lac/trilinos_vector.h>
33 # include <deal.II/lac/vector_memory.h>
34 # include <deal.II/lac/vector_operation.h>
35 
36 # include <Epetra_Comm.h>
37 # include <Epetra_CrsGraph.h>
38 # include <Epetra_Export.h>
39 # include <Epetra_FECrsMatrix.h>
40 # include <Epetra_Map.h>
41 # include <Epetra_MultiVector.h>
42 # include <Epetra_Operator.h>
43 
44 # include <cmath>
45 # include <memory>
46 # include <type_traits>
47 # include <vector>
48 # ifdef DEAL_II_WITH_MPI
49 # include <Epetra_MpiComm.h>
50 # include <mpi.h>
51 # else
52 # include <Epetra_SerialComm.h>
53 # endif
54 
55 DEAL_II_NAMESPACE_OPEN
56 
57 // forward declarations
58 # ifndef DOXYGEN
59 template <typename MatrixType>
60 class BlockMatrixBase;
61 
62 template <typename number>
63 class SparseMatrix;
64 class SparsityPattern;
66 
67 namespace TrilinosWrappers
68 {
69  class SparseMatrix;
70  class SparsityPattern;
71 
72  namespace SparseMatrixIterators
73  {
74  template <bool Constness>
75  class Iterator;
76  }
77 } // namespace TrilinosWrappers
78 # endif
79 
80 namespace TrilinosWrappers
81 {
86  {
91 
96  std::size_t,
97  std::size_t,
98  std::size_t,
99  << "You tried to access row " << arg1
100  << " of a distributed sparsity pattern, "
101  << " but only rows " << arg2 << " through " << arg3
102  << " are stored locally and can be accessed.");
103 
115  {
116  public:
121 
125  AccessorBase(SparseMatrix * matrix,
126  const size_type row,
127  const size_type index);
128 
132  size_type
133  row() const;
134 
138  size_type
139  index() const;
140 
144  size_type
145  column() const;
146 
147  protected:
159 
164 
170  void
171  visit_present_row();
172 
185  std::shared_ptr<std::vector<size_type>> colnum_cache;
186 
190  std::shared_ptr<std::vector<TrilinosScalar>> value_cache;
191  };
192 
203  template <bool Constess>
204  class Accessor : public AccessorBase
205  {
209  TrilinosScalar
210  value() const;
211 
215  TrilinosScalar &
216  value();
217  };
218 
222  template <>
223  class Accessor<true> : public AccessorBase
224  {
225  public:
230  using MatrixType = const SparseMatrix;
231 
236  Accessor(MatrixType *matrix, const size_type row, const size_type index);
237 
242  template <bool Other>
243  Accessor(const Accessor<Other> &a);
244 
248  TrilinosScalar
249  value() const;
250 
251  private:
252  // Make iterator class a friend.
253  template <bool>
254  friend class Iterator;
255  };
256 
260  template <>
261  class Accessor<false> : public AccessorBase
262  {
263  class Reference
264  {
265  public:
269  Reference(const Accessor<false> &accessor);
270 
274  operator TrilinosScalar() const;
275 
279  const Reference &
280  operator=(const TrilinosScalar n) const;
281 
285  const Reference &
286  operator+=(const TrilinosScalar n) const;
287 
291  const Reference &
292  operator-=(const TrilinosScalar n) const;
293 
297  const Reference &
298  operator*=(const TrilinosScalar n) const;
299 
303  const Reference &
304  operator/=(const TrilinosScalar n) const;
305 
306  private:
311  Accessor &accessor;
312  };
313 
314  public:
320 
325  Accessor(MatrixType *matrix, const size_type row, const size_type index);
326 
330  Reference
331  value() const;
332 
333  private:
334  // Make iterator class a friend.
335  template <bool>
336  friend class Iterator;
337 
338  // Make Reference object a friend.
339  friend class Reference;
340  };
341 
356  template <bool Constness>
357  class Iterator
358  {
359  public:
364 
370 
375  Iterator(MatrixType *matrix, const size_type row, const size_type index);
376 
380  template <bool Other>
381  Iterator(const Iterator<Other> &other);
382 
387  operator++();
388 
393  operator++(int);
394 
398  const Accessor<Constness> &operator*() const;
399 
403  const Accessor<Constness> *operator->() const;
404 
409  bool
410  operator==(const Iterator<Constness> &) const;
411 
415  bool
416  operator!=(const Iterator<Constness> &) const;
417 
423  bool
424  operator<(const Iterator<Constness> &) const;
425 
429  bool
430  operator>(const Iterator<Constness> &) const;
431 
435  DeclException2(ExcInvalidIndexWithinRow,
436  size_type,
437  size_type,
438  << "Attempt to access element " << arg2 << " of row "
439  << arg1 << " which doesn't have that many elements.");
440 
441  private:
446 
447  template <bool Other>
448  friend class Iterator;
449  };
450 
451  } // namespace SparseMatrixIterators
452 
453 
515  class SparseMatrix : public Subscriptor
516  {
517  public:
522 
527  std::size_t,
528  << "You tried to access row " << arg1
529  << " of a non-contiguous locally owned row set."
530  << " The row " << arg1
531  << " is not stored locally and can't be accessed.");
532 
540  struct Traits
541  {
546  static const bool zero_addition_can_be_elided = true;
547  };
548 
553 
558 
562  using value_type = TrilinosScalar;
563 
571  SparseMatrix();
572 
580  SparseMatrix(const size_type m,
581  const size_type n,
582  const unsigned int n_max_entries_per_row);
583 
591  SparseMatrix(const size_type m,
592  const size_type n,
593  const std::vector<unsigned int> &n_entries_per_row);
594 
598  SparseMatrix(const SparsityPattern &InputSparsityPattern);
599 
604  SparseMatrix(SparseMatrix &&other) noexcept;
605 
609  SparseMatrix(const SparseMatrix &) = delete;
610 
614  SparseMatrix &
615  operator=(const SparseMatrix &) = delete;
616 
620  virtual ~SparseMatrix() override = default;
621 
637  template <typename SparsityPatternType>
638  void
639  reinit(const SparsityPatternType &sparsity_pattern);
640 
653  void
654  reinit(const SparsityPattern &sparsity_pattern);
655 
664  void
665  reinit(const SparseMatrix &sparse_matrix);
666 
687  template <typename number>
688  void
689  reinit(const ::SparseMatrix<number> &dealii_sparse_matrix,
690  const double drop_tolerance = 1e-13,
691  const bool copy_values = true,
692  const ::SparsityPattern * use_this_sparsity = nullptr);
693 
699  void
700  reinit(const Epetra_CrsMatrix &input_matrix, const bool copy_values = true);
702 
720  DEAL_II_DEPRECATED
721  SparseMatrix(const Epetra_Map &parallel_partitioning,
722  const size_type n_max_entries_per_row = 0);
723 
733  DEAL_II_DEPRECATED
734  SparseMatrix(const Epetra_Map & parallel_partitioning,
735  const std::vector<unsigned int> &n_entries_per_row);
736 
755  DEAL_II_DEPRECATED
756  SparseMatrix(const Epetra_Map &row_parallel_partitioning,
757  const Epetra_Map &col_parallel_partitioning,
758  const size_type n_max_entries_per_row = 0);
759 
776  DEAL_II_DEPRECATED
777  SparseMatrix(const Epetra_Map & row_parallel_partitioning,
778  const Epetra_Map & col_parallel_partitioning,
779  const std::vector<unsigned int> &n_entries_per_row);
780 
807  template <typename SparsityPatternType>
808  DEAL_II_DEPRECATED void
809  reinit(const Epetra_Map & parallel_partitioning,
810  const SparsityPatternType &sparsity_pattern,
811  const bool exchange_data = false);
812 
827  template <typename SparsityPatternType>
828  DEAL_II_DEPRECATED void
829  reinit(const Epetra_Map & row_parallel_partitioning,
830  const Epetra_Map & col_parallel_partitioning,
831  const SparsityPatternType &sparsity_pattern,
832  const bool exchange_data = false);
833 
852  template <typename number>
853  DEAL_II_DEPRECATED void
854  reinit(const Epetra_Map & parallel_partitioning,
855  const ::SparseMatrix<number> &dealii_sparse_matrix,
856  const double drop_tolerance = 1e-13,
857  const bool copy_values = true,
858  const ::SparsityPattern * use_this_sparsity = nullptr);
859 
875  template <typename number>
876  DEAL_II_DEPRECATED void
877  reinit(const Epetra_Map & row_parallel_partitioning,
878  const Epetra_Map & col_parallel_partitioning,
879  const ::SparseMatrix<number> &dealii_sparse_matrix,
880  const double drop_tolerance = 1e-13,
881  const bool copy_values = true,
882  const ::SparsityPattern * use_this_sparsity = nullptr);
884 
900  SparseMatrix(const IndexSet & parallel_partitioning,
901  const MPI_Comm & communicator = MPI_COMM_WORLD,
902  const unsigned int n_max_entries_per_row = 0);
903 
911  SparseMatrix(const IndexSet & parallel_partitioning,
912  const MPI_Comm & communicator,
913  const std::vector<unsigned int> &n_entries_per_row);
914 
929  SparseMatrix(const IndexSet &row_parallel_partitioning,
930  const IndexSet &col_parallel_partitioning,
931  const MPI_Comm &communicator = MPI_COMM_WORLD,
932  const size_type n_max_entries_per_row = 0);
933 
948  SparseMatrix(const IndexSet & row_parallel_partitioning,
949  const IndexSet & col_parallel_partitioning,
950  const MPI_Comm & communicator,
951  const std::vector<unsigned int> &n_entries_per_row);
952 
973  template <typename SparsityPatternType>
974  void
975  reinit(const IndexSet & parallel_partitioning,
976  const SparsityPatternType &sparsity_pattern,
977  const MPI_Comm & communicator = MPI_COMM_WORLD,
978  const bool exchange_data = false);
979 
992  template <typename SparsityPatternType>
993  void
994  reinit(const IndexSet & row_parallel_partitioning,
995  const IndexSet & col_parallel_partitioning,
996  const SparsityPatternType &sparsity_pattern,
997  const MPI_Comm & communicator = MPI_COMM_WORLD,
998  const bool exchange_data = false);
999 
1016  template <typename number>
1017  void
1018  reinit(const IndexSet & parallel_partitioning,
1019  const ::SparseMatrix<number> &dealii_sparse_matrix,
1020  const MPI_Comm & communicator = MPI_COMM_WORLD,
1021  const double drop_tolerance = 1e-13,
1022  const bool copy_values = true,
1023  const ::SparsityPattern * use_this_sparsity = nullptr);
1024 
1038  template <typename number>
1039  void
1040  reinit(const IndexSet & row_parallel_partitioning,
1041  const IndexSet & col_parallel_partitioning,
1042  const ::SparseMatrix<number> &dealii_sparse_matrix,
1043  const MPI_Comm & communicator = MPI_COMM_WORLD,
1044  const double drop_tolerance = 1e-13,
1045  const bool copy_values = true,
1046  const ::SparsityPattern * use_this_sparsity = nullptr);
1048 
1052 
1056  size_type
1057  m() const;
1058 
1062  size_type
1063  n() const;
1064 
1073  unsigned int
1074  local_size() const;
1075 
1084  std::pair<size_type, size_type>
1085  local_range() const;
1086 
1091  bool
1092  in_local_range(const size_type index) const;
1093 
1098  size_type
1099  n_nonzero_elements() const;
1100 
1104  unsigned int
1105  row_length(const size_type row) const;
1106 
1113  bool
1114  is_compressed() const;
1115 
1121  size_type
1122  memory_consumption() const;
1123 
1127  MPI_Comm
1128  get_mpi_communicator() const;
1129 
1131 
1135 
1145  SparseMatrix &
1146  operator=(const double d);
1147 
1155  void
1156  clear();
1157 
1185  void
1186  compress(::VectorOperation::values operation);
1187 
1209  void
1210  set(const size_type i, const size_type j, const TrilinosScalar value);
1211 
1244  void
1245  set(const std::vector<size_type> & indices,
1246  const FullMatrix<TrilinosScalar> &full_matrix,
1247  const bool elide_zero_values = false);
1248 
1254  void
1255  set(const std::vector<size_type> & row_indices,
1256  const std::vector<size_type> & col_indices,
1257  const FullMatrix<TrilinosScalar> &full_matrix,
1258  const bool elide_zero_values = false);
1259 
1287  void
1288  set(const size_type row,
1289  const std::vector<size_type> & col_indices,
1290  const std::vector<TrilinosScalar> &values,
1291  const bool elide_zero_values = false);
1292 
1320  template <typename Number>
1321  void
1322  set(const size_type row,
1323  const size_type n_cols,
1324  const size_type *col_indices,
1325  const Number * values,
1326  const bool elide_zero_values = false);
1327 
1337  void
1338  add(const size_type i, const size_type j, const TrilinosScalar value);
1339 
1358  void
1359  add(const std::vector<size_type> & indices,
1360  const FullMatrix<TrilinosScalar> &full_matrix,
1361  const bool elide_zero_values = true);
1362 
1368  void
1369  add(const std::vector<size_type> & row_indices,
1370  const std::vector<size_type> & col_indices,
1371  const FullMatrix<TrilinosScalar> &full_matrix,
1372  const bool elide_zero_values = true);
1373 
1387  void
1388  add(const size_type row,
1389  const std::vector<size_type> & col_indices,
1390  const std::vector<TrilinosScalar> &values,
1391  const bool elide_zero_values = true);
1392 
1406  void
1407  add(const size_type row,
1408  const size_type n_cols,
1409  const size_type * col_indices,
1410  const TrilinosScalar *values,
1411  const bool elide_zero_values = true,
1412  const bool col_indices_are_sorted = false);
1413 
1417  SparseMatrix &
1418  operator*=(const TrilinosScalar factor);
1419 
1423  SparseMatrix &
1424  operator/=(const TrilinosScalar factor);
1425 
1429  void
1430  copy_from(const SparseMatrix &source);
1431 
1439  void
1440  add(const TrilinosScalar factor, const SparseMatrix &matrix);
1441 
1468  void
1469  clear_row(const size_type row, const TrilinosScalar new_diag_value = 0);
1470 
1491  void
1492  clear_rows(const std::vector<size_type> &rows,
1493  const TrilinosScalar new_diag_value = 0);
1494 
1504  void
1505  transpose();
1506 
1508 
1512 
1521  TrilinosScalar
1522  operator()(const size_type i, const size_type j) const;
1523 
1540  TrilinosScalar
1541  el(const size_type i, const size_type j) const;
1542 
1549  TrilinosScalar
1550  diag_element(const size_type i) const;
1551 
1553 
1557 
1588  template <typename VectorType>
1589  typename std::enable_if<std::is_same<typename VectorType::value_type,
1590  TrilinosScalar>::value>::type
1591  vmult(VectorType &dst, const VectorType &src) const;
1592 
1599  template <typename VectorType>
1600  typename std::enable_if<!std::is_same<typename VectorType::value_type,
1601  TrilinosScalar>::value>::type
1602  vmult(VectorType &dst, const VectorType &src) const;
1603 
1618  template <typename VectorType>
1619  typename std::enable_if<std::is_same<typename VectorType::value_type,
1620  TrilinosScalar>::value>::type
1621  Tvmult(VectorType &dst, const VectorType &src) const;
1622 
1629  template <typename VectorType>
1630  typename std::enable_if<!std::is_same<typename VectorType::value_type,
1631  TrilinosScalar>::value>::type
1632  Tvmult(VectorType &dst, const VectorType &src) const;
1633 
1643  template <typename VectorType>
1644  void
1645  vmult_add(VectorType &dst, const VectorType &src) const;
1646 
1657  template <typename VectorType>
1658  void
1659  Tvmult_add(VectorType &dst, const VectorType &src) const;
1660 
1682  TrilinosScalar
1683  matrix_norm_square(const MPI::Vector &v) const;
1684 
1704  TrilinosScalar
1705  matrix_scalar_product(const MPI::Vector &u, const MPI::Vector &v) const;
1706 
1723  TrilinosScalar
1724  residual(MPI::Vector & dst,
1725  const MPI::Vector &x,
1726  const MPI::Vector &b) const;
1727 
1742  void
1743  mmult(SparseMatrix & C,
1744  const SparseMatrix &B,
1745  const MPI::Vector & V = MPI::Vector()) const;
1746 
1747 
1764  void
1765  Tmmult(SparseMatrix & C,
1766  const SparseMatrix &B,
1767  const MPI::Vector & V = MPI::Vector()) const;
1768 
1770 
1774 
1782  TrilinosScalar
1783  l1_norm() const;
1784 
1793  TrilinosScalar
1794  linfty_norm() const;
1795 
1800  TrilinosScalar
1801  frobenius_norm() const;
1802 
1804 
1808 
1813  const Epetra_CrsMatrix &
1814  trilinos_matrix() const;
1815 
1820  const Epetra_CrsGraph &
1821  trilinos_sparsity_pattern() const;
1822 
1830  DEAL_II_DEPRECATED
1831  const Epetra_Map &
1832  domain_partitioner() const;
1833 
1842  DEAL_II_DEPRECATED
1843  const Epetra_Map &
1844  range_partitioner() const;
1845 
1853  DEAL_II_DEPRECATED
1854  const Epetra_Map &
1855  row_partitioner() const;
1856 
1866  DEAL_II_DEPRECATED
1867  const Epetra_Map &
1868  col_partitioner() const;
1870 
1875 
1880  IndexSet
1881  locally_owned_domain_indices() const;
1882 
1888  IndexSet
1889  locally_owned_range_indices() const;
1890 
1892 
1897 
1917  begin() const;
1918 
1922  iterator
1923  begin();
1924 
1930  end() const;
1931 
1935  iterator
1936  end();
1937 
1967  begin(const size_type r) const;
1968 
1972  iterator
1973  begin(const size_type r);
1974 
1985  end(const size_type r) const;
1986 
1990  iterator
1991  end(const size_type r);
1992 
1994 
1998 
2004  void
2005  write_ascii();
2006 
2014  void
2015  print(std::ostream &out,
2016  const bool write_extended_trilinos_info = false) const;
2017 
2019 
2028  int,
2029  << "An error with error number " << arg1
2030  << " occurred while calling a Trilinos function");
2031 
2035  DeclException2(ExcInvalidIndex,
2036  size_type,
2037  size_type,
2038  << "The entry with index <" << arg1 << ',' << arg2
2039  << "> does not exist.");
2040 
2044  DeclExceptionMsg(ExcSourceEqualsDestination,
2045  "You are attempting an operation on two matrices that "
2046  "are the same object, but the operation requires that the "
2047  "two objects are in fact different.");
2048 
2052  DeclException0(ExcMatrixNotCompressed);
2053 
2057  DeclException4(ExcAccessToNonLocalElement,
2058  size_type,
2059  size_type,
2060  size_type,
2061  size_type,
2062  << "You tried to access element (" << arg1 << "/" << arg2
2063  << ")"
2064  << " of a distributed matrix, but only rows " << arg3
2065  << " through " << arg4
2066  << " are stored locally and can be accessed.");
2067 
2071  DeclException2(ExcAccessToNonPresentElement,
2072  size_type,
2073  size_type,
2074  << "You tried to access element (" << arg1 << "/" << arg2
2075  << ")"
2076  << " of a sparse matrix, but it appears to not"
2077  << " exist in the Trilinos sparsity pattern.");
2079 
2080 
2081 
2082  protected:
2093  void
2094  prepare_add();
2095 
2101  void
2102  prepare_set();
2103 
2104 
2105 
2106  private:
2111  std::unique_ptr<Epetra_Map> column_space_map;
2112 
2118  std::unique_ptr<Epetra_FECrsMatrix> matrix;
2119 
2125  std::unique_ptr<Epetra_CrsMatrix> nonlocal_matrix;
2126 
2130  std::unique_ptr<Epetra_Export> nonlocal_matrix_exporter;
2131 
2143  Epetra_CombineMode last_action;
2144 
2150 
2151  // To allow calling protected prepare_add() and prepare_set().
2152  friend class BlockMatrixBase<SparseMatrix>;
2153  };
2154 
2155 
2156 
2157  // forwards declarations
2158  class SolverBase;
2159  class PreconditionBase;
2160 
2161  namespace internal
2162  {
2163  inline void
2164  check_vector_map_equality(const Epetra_CrsMatrix & mtrx,
2165  const Epetra_MultiVector &src,
2166  const Epetra_MultiVector &dst,
2167  const bool transpose)
2168  {
2169  if (transpose == false)
2170  {
2171  Assert(src.Map().SameAs(mtrx.DomainMap()) == true,
2172  ExcMessage(
2173  "Column map of matrix does not fit with vector map!"));
2174  Assert(dst.Map().SameAs(mtrx.RangeMap()) == true,
2175  ExcMessage("Row map of matrix does not fit with vector map!"));
2176  }
2177  else
2178  {
2179  Assert(src.Map().SameAs(mtrx.RangeMap()) == true,
2180  ExcMessage(
2181  "Column map of matrix does not fit with vector map!"));
2182  Assert(dst.Map().SameAs(mtrx.DomainMap()) == true,
2183  ExcMessage("Row map of matrix does not fit with vector map!"));
2184  }
2185  (void)mtrx; // removes -Wunused-variable in optimized mode
2186  (void)src;
2187  (void)dst;
2188  }
2189 
2190  inline void
2191  check_vector_map_equality(const Epetra_Operator & op,
2192  const Epetra_MultiVector &src,
2193  const Epetra_MultiVector &dst,
2194  const bool transpose)
2195  {
2196  if (transpose == false)
2197  {
2198  Assert(src.Map().SameAs(op.OperatorDomainMap()) == true,
2199  ExcMessage(
2200  "Column map of operator does not fit with vector map!"));
2201  Assert(dst.Map().SameAs(op.OperatorRangeMap()) == true,
2202  ExcMessage(
2203  "Row map of operator does not fit with vector map!"));
2204  }
2205  else
2206  {
2207  Assert(src.Map().SameAs(op.OperatorRangeMap()) == true,
2208  ExcMessage(
2209  "Column map of operator does not fit with vector map!"));
2210  Assert(dst.Map().SameAs(op.OperatorDomainMap()) == true,
2211  ExcMessage(
2212  "Row map of operator does not fit with vector map!"));
2213  }
2214  (void)op; // removes -Wunused-variable in optimized mode
2215  (void)src;
2216  (void)dst;
2217  }
2218 
2219  namespace LinearOperatorImplementation
2220  {
2241  class TrilinosPayload : public Epetra_Operator
2242  {
2243  public:
2247  using VectorType = Epetra_MultiVector;
2248 
2253 
2258 
2263 
2271  TrilinosPayload();
2272 
2276  TrilinosPayload(const TrilinosWrappers::SparseMatrix &matrix_exemplar,
2277  const TrilinosWrappers::SparseMatrix &matrix);
2278 
2283  const TrilinosWrappers::SparseMatrix & matrix_exemplar,
2284  const TrilinosWrappers::PreconditionBase &preconditioner);
2285 
2290  const TrilinosWrappers::PreconditionBase &preconditioner_exemplar,
2291  const TrilinosWrappers::PreconditionBase &preconditioner);
2292 
2296  TrilinosPayload(const TrilinosPayload &payload);
2297 
2305  TrilinosPayload(const TrilinosPayload &first_op,
2306  const TrilinosPayload &second_op);
2307 
2311  virtual ~TrilinosPayload() override = default;
2312 
2317  identity_payload() const;
2318 
2323  null_payload() const;
2324 
2329  transpose_payload() const;
2330 
2347  template <typename Solver, typename Preconditioner>
2348  typename std::enable_if<
2349  std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
2350  std::is_base_of<TrilinosWrappers::PreconditionBase,
2351  Preconditioner>::value,
2352  TrilinosPayload>::type
2353  inverse_payload(Solver &, const Preconditioner &) const;
2354 
2372  template <typename Solver, typename Preconditioner>
2373  typename std::enable_if<
2374  !(std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
2375  std::is_base_of<TrilinosWrappers::PreconditionBase,
2376  Preconditioner>::value),
2377  TrilinosPayload>::type
2378  inverse_payload(Solver &, const Preconditioner &) const;
2379 
2381 
2386 
2392  IndexSet
2393  locally_owned_domain_indices() const;
2394 
2400  IndexSet
2401  locally_owned_range_indices() const;
2402 
2406  MPI_Comm
2407  get_mpi_communicator() const;
2408 
2415  void
2416  transpose();
2417 
2425  std::function<void(VectorType &, const VectorType &)> vmult;
2426 
2434  std::function<void(VectorType &, const VectorType &)> Tvmult;
2435 
2444  std::function<void(VectorType &, const VectorType &)> inv_vmult;
2445 
2454  std::function<void(VectorType &, const VectorType &)> inv_Tvmult;
2455 
2457 
2462 
2469  virtual bool
2470  UseTranspose() const override;
2471 
2487  virtual int
2488  SetUseTranspose(bool UseTranspose) override;
2489 
2501  virtual int
2502  Apply(const VectorType &X, VectorType &Y) const override;
2503 
2522  virtual int
2523  ApplyInverse(const VectorType &Y, VectorType &X) const override;
2525 
2530 
2537  virtual const char *
2538  Label() const override;
2539 
2547  virtual const Epetra_Comm &
2548  Comm() const override;
2549 
2557  virtual const Epetra_Map &
2558  OperatorDomainMap() const override;
2559 
2568  virtual const Epetra_Map &
2569  OperatorRangeMap() const override;
2571 
2572  private:
2578 
2583 # ifdef DEAL_II_WITH_MPI
2584  Epetra_MpiComm communicator;
2585 # else
2586  Epetra_SerialComm communicator;
2587 # endif
2588 
2593  Epetra_Map domain_map;
2594 
2599  Epetra_Map range_map;
2600 
2609  virtual bool
2610  HasNormInf() const override;
2611 
2619  virtual double
2620  NormInf() const override;
2621  };
2622 
2628  operator+(const TrilinosPayload &first_op,
2629  const TrilinosPayload &second_op);
2630 
2635  TrilinosPayload operator*(const TrilinosPayload &first_op,
2636  const TrilinosPayload &second_op);
2637 
2638  } // namespace LinearOperatorImplementation
2639  } /* namespace internal */
2640 
2641 
2642 
2643  // ----------------------- inline and template functions --------------------
2644 
2645 # ifndef DOXYGEN
2646 
2647  namespace SparseMatrixIterators
2648  {
2649  inline AccessorBase::AccessorBase(SparseMatrix *matrix,
2650  size_type row,
2651  size_type index)
2652  : matrix(matrix)
2653  , a_row(row)
2654  , a_index(index)
2655  {
2656  visit_present_row();
2657  }
2658 
2659 
2660  inline AccessorBase::size_type
2661  AccessorBase::row() const
2662  {
2663  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2664  return a_row;
2665  }
2666 
2667 
2668  inline AccessorBase::size_type
2669  AccessorBase::column() const
2670  {
2671  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2672  return (*colnum_cache)[a_index];
2673  }
2674 
2675 
2676  inline AccessorBase::size_type
2677  AccessorBase::index() const
2678  {
2679  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2680  return a_index;
2681  }
2682 
2683 
2684  inline Accessor<true>::Accessor(MatrixType * matrix,
2685  const size_type row,
2686  const size_type index)
2687  : AccessorBase(const_cast<SparseMatrix *>(matrix), row, index)
2688  {}
2689 
2690 
2691  template <bool Other>
2692  inline Accessor<true>::Accessor(const Accessor<Other> &other)
2693  : AccessorBase(other)
2694  {}
2695 
2696 
2697  inline TrilinosScalar
2698  Accessor<true>::value() const
2699  {
2700  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2701  return (*value_cache)[a_index];
2702  }
2703 
2704 
2705  inline Accessor<false>::Reference::Reference(const Accessor<false> &acc)
2706  : accessor(const_cast<Accessor<false> &>(acc))
2707  {}
2708 
2709 
2710  inline Accessor<false>::Reference::operator TrilinosScalar() const
2711  {
2712  return (*accessor.value_cache)[accessor.a_index];
2713  }
2714 
2715  inline const Accessor<false>::Reference &
2716  Accessor<false>::Reference::operator=(const TrilinosScalar n) const
2717  {
2718  (*accessor.value_cache)[accessor.a_index] = n;
2719  accessor.matrix->set(accessor.row(),
2720  accessor.column(),
2721  static_cast<TrilinosScalar>(*this));
2722  return *this;
2723  }
2724 
2725 
2726  inline const Accessor<false>::Reference &
2727  Accessor<false>::Reference::operator+=(const TrilinosScalar n) const
2728  {
2729  (*accessor.value_cache)[accessor.a_index] += n;
2730  accessor.matrix->set(accessor.row(),
2731  accessor.column(),
2732  static_cast<TrilinosScalar>(*this));
2733  return *this;
2734  }
2735 
2736 
2737  inline const Accessor<false>::Reference &
2738  Accessor<false>::Reference::operator-=(const TrilinosScalar n) const
2739  {
2740  (*accessor.value_cache)[accessor.a_index] -= n;
2741  accessor.matrix->set(accessor.row(),
2742  accessor.column(),
2743  static_cast<TrilinosScalar>(*this));
2744  return *this;
2745  }
2746 
2747 
2748  inline const Accessor<false>::Reference &
2749  Accessor<false>::Reference::operator*=(const TrilinosScalar n) const
2750  {
2751  (*accessor.value_cache)[accessor.a_index] *= n;
2752  accessor.matrix->set(accessor.row(),
2753  accessor.column(),
2754  static_cast<TrilinosScalar>(*this));
2755  return *this;
2756  }
2757 
2758 
2759  inline const Accessor<false>::Reference &
2760  Accessor<false>::Reference::operator/=(const TrilinosScalar n) const
2761  {
2762  (*accessor.value_cache)[accessor.a_index] /= n;
2763  accessor.matrix->set(accessor.row(),
2764  accessor.column(),
2765  static_cast<TrilinosScalar>(*this));
2766  return *this;
2767  }
2768 
2769 
2770  inline Accessor<false>::Accessor(MatrixType * matrix,
2771  const size_type row,
2772  const size_type index)
2773  : AccessorBase(matrix, row, index)
2774  {}
2775 
2776 
2777  inline Accessor<false>::Reference
2778  Accessor<false>::value() const
2779  {
2780  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2781  return {*this};
2782  }
2783 
2784 
2785 
2786  template <bool Constness>
2787  inline Iterator<Constness>::Iterator(MatrixType * matrix,
2788  const size_type row,
2789  const size_type index)
2790  : accessor(matrix, row, index)
2791  {}
2792 
2793 
2794  template <bool Constness>
2795  template <bool Other>
2796  inline Iterator<Constness>::Iterator(const Iterator<Other> &other)
2797  : accessor(other.accessor)
2798  {}
2799 
2800 
2801  template <bool Constness>
2802  inline Iterator<Constness> &
2803  Iterator<Constness>::operator++()
2804  {
2805  Assert(accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
2806 
2807  ++accessor.a_index;
2808 
2809  // If at end of line: do one
2810  // step, then cycle until we
2811  // find a row with a nonzero
2812  // number of entries.
2813  if (accessor.a_index >= accessor.colnum_cache->size())
2814  {
2815  accessor.a_index = 0;
2816  ++accessor.a_row;
2817 
2818  while ((accessor.a_row < accessor.matrix->m()) &&
2819  ((accessor.matrix->in_local_range(accessor.a_row) == false) ||
2820  (accessor.matrix->row_length(accessor.a_row) == 0)))
2821  ++accessor.a_row;
2822 
2823  accessor.visit_present_row();
2824  }
2825  return *this;
2826  }
2827 
2828 
2829  template <bool Constness>
2830  inline Iterator<Constness>
2831  Iterator<Constness>::operator++(int)
2832  {
2833  const Iterator<Constness> old_state = *this;
2834  ++(*this);
2835  return old_state;
2836  }
2837 
2838 
2839 
2840  template <bool Constness>
2841  inline const Accessor<Constness> &Iterator<Constness>::operator*() const
2842  {
2843  return accessor;
2844  }
2845 
2846 
2847 
2848  template <bool Constness>
2849  inline const Accessor<Constness> *Iterator<Constness>::operator->() const
2850  {
2851  return &accessor;
2852  }
2853 
2854 
2855 
2856  template <bool Constness>
2857  inline bool
2858  Iterator<Constness>::operator==(const Iterator<Constness> &other) const
2859  {
2860  return (accessor.a_row == other.accessor.a_row &&
2861  accessor.a_index == other.accessor.a_index);
2862  }
2863 
2864 
2865 
2866  template <bool Constness>
2867  inline bool
2868  Iterator<Constness>::operator!=(const Iterator<Constness> &other) const
2869  {
2870  return !(*this == other);
2871  }
2872 
2873 
2874 
2875  template <bool Constness>
2876  inline bool
2877  Iterator<Constness>::operator<(const Iterator<Constness> &other) const
2878  {
2879  return (accessor.row() < other.accessor.row() ||
2880  (accessor.row() == other.accessor.row() &&
2881  accessor.index() < other.accessor.index()));
2882  }
2883 
2884 
2885  template <bool Constness>
2886  inline bool
2887  Iterator<Constness>::operator>(const Iterator<Constness> &other) const
2888  {
2889  return (other < *this);
2890  }
2891 
2892  } // namespace SparseMatrixIterators
2893 
2894 
2895 
2897  SparseMatrix::begin() const
2898  {
2899  return begin(0);
2900  }
2901 
2902 
2903 
2905  SparseMatrix::end() const
2906  {
2907  return const_iterator(this, m(), 0);
2908  }
2909 
2910 
2911 
2913  SparseMatrix::begin(const size_type r) const
2914  {
2915  Assert(r < m(), ExcIndexRange(r, 0, m()));
2916  if (in_local_range(r) && (row_length(r) > 0))
2917  return const_iterator(this, r, 0);
2918  else
2919  return end(r);
2920  }
2921 
2922 
2923 
2925  SparseMatrix::end(const size_type r) const
2926  {
2927  Assert(r < m(), ExcIndexRange(r, 0, m()));
2928 
2929  // place the iterator on the first entry
2930  // past this line, or at the end of the
2931  // matrix
2932  for (size_type i = r + 1; i < m(); ++i)
2933  if (in_local_range(i) && (row_length(i) > 0))
2934  return const_iterator(this, i, 0);
2935 
2936  // if there is no such line, then take the
2937  // end iterator of the matrix
2938  return end();
2939  }
2940 
2941 
2942 
2943  inline SparseMatrix::iterator
2945  {
2946  return begin(0);
2947  }
2948 
2949 
2950 
2951  inline SparseMatrix::iterator
2953  {
2954  return iterator(this, m(), 0);
2955  }
2956 
2957 
2958 
2959  inline SparseMatrix::iterator
2960  SparseMatrix::begin(const size_type r)
2961  {
2962  Assert(r < m(), ExcIndexRange(r, 0, m()));
2963  if (in_local_range(r) && (row_length(r) > 0))
2964  return iterator(this, r, 0);
2965  else
2966  return end(r);
2967  }
2968 
2969 
2970 
2971  inline SparseMatrix::iterator
2972  SparseMatrix::end(const size_type r)
2973  {
2974  Assert(r < m(), ExcIndexRange(r, 0, m()));
2975 
2976  // place the iterator on the first entry
2977  // past this line, or at the end of the
2978  // matrix
2979  for (size_type i = r + 1; i < m(); ++i)
2980  if (in_local_range(i) && (row_length(i) > 0))
2981  return iterator(this, i, 0);
2982 
2983  // if there is no such line, then take the
2984  // end iterator of the matrix
2985  return end();
2986  }
2987 
2988 
2989 
2990  inline bool
2991  SparseMatrix::in_local_range(const size_type index) const
2992  {
2993  TrilinosWrappers::types::int_type begin, end;
2994 # ifndef DEAL_II_WITH_64BIT_INDICES
2995  begin = matrix->RowMap().MinMyGID();
2996  end = matrix->RowMap().MaxMyGID() + 1;
2997 # else
2998  begin = matrix->RowMap().MinMyGID64();
2999  end = matrix->RowMap().MaxMyGID64() + 1;
3000 # endif
3001 
3002  return ((index >= static_cast<size_type>(begin)) &&
3003  (index < static_cast<size_type>(end)));
3004  }
3005 
3006 
3007 
3008  inline bool
3010  {
3011  return compressed;
3012  }
3013 
3014 
3015 
3016  // Inline the set() and add() functions, since they will be called
3017  // frequently, and the compiler can optimize away some unnecessary loops
3018  // when the sizes are given at compile time.
3019  template <>
3020  void
3021  SparseMatrix::set<TrilinosScalar>(const size_type row,
3022  const size_type n_cols,
3023  const size_type * col_indices,
3024  const TrilinosScalar *values,
3025  const bool elide_zero_values);
3026 
3027 
3028 
3029  template <typename Number>
3030  void
3031  SparseMatrix::set(const size_type row,
3032  const size_type n_cols,
3033  const size_type *col_indices,
3034  const Number * values,
3035  const bool elide_zero_values)
3036  {
3037  std::vector<TrilinosScalar> trilinos_values(n_cols);
3038  std::copy(values, values + n_cols, trilinos_values.begin());
3039  this->set(
3040  row, n_cols, col_indices, trilinos_values.data(), elide_zero_values);
3041  }
3042 
3043 
3044 
3045  inline void
3046  SparseMatrix::set(const size_type i,
3047  const size_type j,
3048  const TrilinosScalar value)
3049  {
3050  AssertIsFinite(value);
3051 
3052  set(i, 1, &j, &value, false);
3053  }
3054 
3055 
3056 
3057  inline void
3058  SparseMatrix::set(const std::vector<size_type> & indices,
3059  const FullMatrix<TrilinosScalar> &values,
3060  const bool elide_zero_values)
3061  {
3062  Assert(indices.size() == values.m(),
3063  ExcDimensionMismatch(indices.size(), values.m()));
3064  Assert(values.m() == values.n(), ExcNotQuadratic());
3065 
3066  for (size_type i = 0; i < indices.size(); ++i)
3067  set(indices[i],
3068  indices.size(),
3069  indices.data(),
3070  &values(i, 0),
3071  elide_zero_values);
3072  }
3073 
3074 
3075 
3076  inline void
3077  SparseMatrix::add(const size_type i,
3078  const size_type j,
3079  const TrilinosScalar value)
3080  {
3081  AssertIsFinite(value);
3082 
3083  if (value == 0)
3084  {
3085  // we have to check after Insert/Add in any case to be consistent
3086  // with the MPI communication model, but we can save some
3087  // work if the addend is zero. However, these actions are done in case
3088  // we pass on to the other function.
3089 
3090  // TODO: fix this (do not run compress here, but fail)
3091  if (last_action == Insert)
3092  {
3093  int ierr;
3094  ierr = matrix->GlobalAssemble(*column_space_map,
3095  matrix->RowMap(),
3096  false);
3097 
3098  Assert(ierr == 0, ExcTrilinosError(ierr));
3099  (void)ierr; // removes -Wunused-but-set-variable in optimized mode
3100  }
3101 
3102  last_action = Add;
3103 
3104  return;
3105  }
3106  else
3107  add(i, 1, &j, &value, false);
3108  }
3109 
3110 
3111 
3112  // inline "simple" functions that are called frequently and do only involve
3113  // a call to some Trilinos function.
3115  SparseMatrix::m() const
3116  {
3117 # ifndef DEAL_II_WITH_64BIT_INDICES
3118  return matrix->NumGlobalRows();
3119 # else
3120  return matrix->NumGlobalRows64();
3121 # endif
3122  }
3123 
3124 
3125 
3127  SparseMatrix::n() const
3128  {
3129  // If the matrix structure has not been fixed (i.e., we did not have a
3130  // sparsity pattern), it does not know about the number of columns so we
3131  // must always take this from the additional column space map
3132  Assert(column_space_map.get() != nullptr, ExcInternalError());
3133  return n_global_elements(*column_space_map);
3134  }
3135 
3136 
3137 
3138  inline unsigned int
3139  SparseMatrix::local_size() const
3140  {
3141  return matrix->NumMyRows();
3142  }
3143 
3144 
3145 
3146  inline std::pair<SparseMatrix::size_type, SparseMatrix::size_type>
3148  {
3149  size_type begin, end;
3150 # ifndef DEAL_II_WITH_64BIT_INDICES
3151  begin = matrix->RowMap().MinMyGID();
3152  end = matrix->RowMap().MaxMyGID() + 1;
3153 # else
3154  begin = matrix->RowMap().MinMyGID64();
3155  end = matrix->RowMap().MaxMyGID64() + 1;
3156 # endif
3157 
3158  return std::make_pair(begin, end);
3159  }
3160 
3161 
3162 
3165  {
3166 # ifndef DEAL_II_WITH_64BIT_INDICES
3167  return matrix->NumGlobalNonzeros();
3168 # else
3169  return matrix->NumGlobalNonzeros64();
3170 # endif
3171  }
3172 
3173 
3174 
3175  template <typename SparsityPatternType>
3176  inline void
3177  SparseMatrix::reinit(const IndexSet & parallel_partitioning,
3178  const SparsityPatternType &sparsity_pattern,
3179  const MPI_Comm & communicator,
3180  const bool exchange_data)
3181  {
3182  reinit(parallel_partitioning,
3183  parallel_partitioning,
3184  sparsity_pattern,
3185  communicator,
3186  exchange_data);
3187  }
3188 
3189 
3190 
3191  template <typename number>
3192  inline void
3193  SparseMatrix::reinit(const IndexSet &parallel_partitioning,
3194  const ::SparseMatrix<number> &sparse_matrix,
3195  const MPI_Comm & communicator,
3196  const double drop_tolerance,
3197  const bool copy_values,
3198  const ::SparsityPattern * use_this_sparsity)
3199  {
3200  Epetra_Map map =
3201  parallel_partitioning.make_trilinos_map(communicator, false);
3202  reinit(parallel_partitioning,
3203  parallel_partitioning,
3204  sparse_matrix,
3205  drop_tolerance,
3206  copy_values,
3207  use_this_sparsity);
3208  }
3209 
3210 
3211 
3212  inline const Epetra_CrsMatrix &
3214  {
3215  return static_cast<const Epetra_CrsMatrix &>(*matrix);
3216  }
3217 
3218 
3219 
3220  inline const Epetra_CrsGraph &
3222  {
3223  return matrix->Graph();
3224  }
3225 
3226 
3227 
3228  inline IndexSet
3230  {
3231  return IndexSet(matrix->DomainMap());
3232  }
3233 
3234 
3235 
3236  inline IndexSet
3238  {
3239  return IndexSet(matrix->RangeMap());
3240  }
3241 
3242 
3243 
3244  inline void
3246  {
3247  // nothing to do here
3248  }
3249 
3250 
3251 
3252  inline void
3254  {
3255  // nothing to do here
3256  }
3257 
3258 
3259  namespace internal
3260  {
3261  namespace LinearOperatorImplementation
3262  {
3263  template <typename Solver, typename Preconditioner>
3264  typename std::enable_if<
3265  std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
3266  std::is_base_of<TrilinosWrappers::PreconditionBase,
3267  Preconditioner>::value,
3268  TrilinosPayload>::type
3270  Solver & solver,
3271  const Preconditioner &preconditioner) const
3272  {
3273  const auto &payload = *this;
3274 
3275  TrilinosPayload return_op(payload);
3276 
3277  // Capture by copy so the payloads are always valid
3278 
3279  return_op.inv_vmult = [payload, &solver, &preconditioner](
3280  TrilinosPayload::Domain & tril_dst,
3281  const TrilinosPayload::Range &tril_src) {
3282  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3283  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3284  Assert(&tril_src != &tril_dst,
3286  internal::check_vector_map_equality(payload,
3287  tril_src,
3288  tril_dst,
3289  !payload.UseTranspose());
3290  solver.solve(payload, tril_dst, tril_src, preconditioner);
3291  };
3292 
3293  return_op.inv_Tvmult = [payload, &solver, &preconditioner](
3294  TrilinosPayload::Range & tril_dst,
3295  const TrilinosPayload::Domain &tril_src) {
3296  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3297  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3298  Assert(&tril_src != &tril_dst,
3300  internal::check_vector_map_equality(payload,
3301  tril_src,
3302  tril_dst,
3303  payload.UseTranspose());
3304 
3305  const_cast<TrilinosPayload &>(payload).transpose();
3306  solver.solve(payload, tril_dst, tril_src, preconditioner);
3307  const_cast<TrilinosPayload &>(payload).transpose();
3308  };
3309 
3310  // If the input operator is already setup for transpose operations, then
3311  // we must do similar with its inverse.
3312  if (return_op.UseTranspose() == true)
3313  std::swap(return_op.inv_vmult, return_op.inv_Tvmult);
3314 
3315  return return_op;
3316  }
3317 
3318  template <typename Solver, typename Preconditioner>
3319  typename std::enable_if<
3320  !(std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
3321  std::is_base_of<TrilinosWrappers::PreconditionBase,
3322  Preconditioner>::value),
3323  TrilinosPayload>::type
3324  TrilinosPayload::inverse_payload(Solver &, const Preconditioner &) const
3325  {
3326  TrilinosPayload return_op(*this);
3327 
3328  return_op.inv_vmult = [](TrilinosPayload::Domain &,
3329  const TrilinosPayload::Range &) {
3330  AssertThrow(false,
3331  ExcMessage("Payload inv_vmult disabled because of "
3332  "incompatible solver/preconditioner choice."));
3333  };
3334 
3335  return_op.inv_Tvmult = [](TrilinosPayload::Range &,
3336  const TrilinosPayload::Domain &) {
3337  AssertThrow(false,
3338  ExcMessage("Payload inv_vmult disabled because of "
3339  "incompatible solver/preconditioner choice."));
3340  };
3341 
3342  return return_op;
3343  }
3344  } // namespace LinearOperatorImplementation
3345  } // namespace internal
3346 
3347  template <>
3348  void
3349  SparseMatrix::set<TrilinosScalar>(const size_type row,
3350  const size_type n_cols,
3351  const size_type * col_indices,
3352  const TrilinosScalar *values,
3353  const bool elide_zero_values);
3354 # endif // DOXYGEN
3355 
3356 } /* namespace TrilinosWrappers */
3357 
3358 
3359 DEAL_II_NAMESPACE_CLOSE
3360 
3361 
3362 # endif // DEAL_II_WITH_TRILINOS
3363 
3364 
3365 /*----------------------- trilinos_sparse_matrix.h --------------------*/
3366 
3367 #endif
3368 /*----------------------- trilinos_sparse_matrix.h --------------------*/
DEAL_II_CONSTEXPR SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
static ::ExceptionBase & ExcTrilinosError(int arg1)
size_type m() const
const_iterator end() const
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:541
DEAL_II_CONSTEXPR SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
static ::ExceptionBase & ExcSourceEqualsDestination()
static ::ExceptionBase & ExcBeyondEndOfMatrix()
std::function< void(VectorType &, const VectorType &)> inv_Tvmult
static ::ExceptionBase & ExcAccessToNonlocalRow(std::size_t arg1, std::size_t arg2, std::size_t arg3)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
constexpr SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:594
std::shared_ptr< std::vector< size_type > > colnum_cache
std::unique_ptr< Epetra_FECrsMatrix > matrix
::types::global_dof_index size_type
bool in_local_range(const size_type index) const
static ::ExceptionBase & ExcMessage(std::string arg1)
size_type n() const
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:518
void add(const size_type i, const size_type j, const TrilinosScalar value)
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::unique_ptr< Epetra_Export > nonlocal_matrix_exporter
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:496
#define DeclException0(Exception0)
Definition: exceptions.h:473
typename Accessor< Constness >::MatrixType MatrixType
std::function< void(VectorType &, const VectorType &)> inv_vmult
std::enable_if< std::is_base_of< TrilinosWrappers::SolverBase, Solver >::value &&std::is_base_of< TrilinosWrappers::PreconditionBase, Preconditioner >::value, TrilinosPayload >::type inverse_payload(Solver &, const Preconditioner &) const
const Epetra_CrsMatrix & trilinos_matrix() const
AccessorBase(SparseMatrix *matrix, const size_type row, const size_type index)
IndexSet locally_owned_range_indices() const
size_type n_nonzero_elements() const
static ::ExceptionBase & ExcIteratorPastEnd()
const_iterator begin() const
static ::ExceptionBase & ExcNotQuadratic()
TrilinosWrappers::types::int_type n_global_elements(const Epetra_BlockMap &map)
void swap(Vector< Number > &u, Vector< Number > &v)
Definition: vector.h:1376
unsigned int global_dof_index
Definition: types.h:89
void reinit(const SparsityPatternType &sparsity_pattern)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:587
std::pair< size_type, size_type > local_range() const
IndexSet locally_owned_domain_indices() const
const Epetra_CrsGraph & trilinos_sparsity_pattern() const
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:564
std::unique_ptr< Epetra_Map > column_space_map
std::unique_ptr< Epetra_CrsMatrix > nonlocal_matrix
void set(const size_type i, const size_type j, const TrilinosScalar value)
#define AssertIsFinite(number)
Definition: exceptions.h:1669
unsigned int local_size() const
static ::ExceptionBase & ExcInternalError()
std::shared_ptr< std::vector< TrilinosScalar > > value_cache