Reference documentation for deal.II version 9.1.1
\(\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_tpetra_vector.h>
31 # include <deal.II/lac/trilinos_vector.h>
32 # include <deal.II/lac/vector_memory.h>
33 # include <deal.II/lac/vector_operation.h>
34 
35 # include <Epetra_Comm.h>
36 # include <Epetra_CrsGraph.h>
37 # include <Epetra_Export.h>
38 # include <Epetra_FECrsMatrix.h>
39 # include <Epetra_Map.h>
40 # include <Epetra_MultiVector.h>
41 # include <Epetra_Operator.h>
42 
43 # include <cmath>
44 # include <memory>
45 # include <type_traits>
46 # include <vector>
47 # ifdef DEAL_II_WITH_MPI
48 # include <Epetra_MpiComm.h>
49 # include <mpi.h>
50 # else
51 # include <Epetra_SerialComm.h>
52 # endif
53 
54 DEAL_II_NAMESPACE_OPEN
55 
56 // forward declarations
57 template <typename MatrixType>
58 class BlockMatrixBase;
59 
60 template <typename number>
61 class SparseMatrix;
62 class SparsityPattern;
64 
65 
66 
67 namespace TrilinosWrappers
68 {
69  // forward declarations
70  class SparseMatrix;
71  class SparsityPattern;
72 
77  {
78  // forward declaration
79  template <bool Constness>
80  class Iterator;
81 
86 
91  std::size_t,
92  std::size_t,
93  std::size_t,
94  << "You tried to access row " << arg1
95  << " of a distributed sparsity pattern, "
96  << " but only rows " << arg2 << " through " << arg3
97  << " are stored locally and can be accessed.");
98 
110  {
111  public:
116 
121  const size_type row,
122  const size_type index);
123 
127  size_type
128  row() const;
129 
133  size_type
134  index() const;
135 
139  size_type
140  column() const;
141 
142  protected:
154 
159 
165  void
167 
180  std::shared_ptr<std::vector<size_type>> colnum_cache;
181 
185  std::shared_ptr<std::vector<TrilinosScalar>> value_cache;
186  };
187 
198  template <bool Constess>
199  class Accessor : public AccessorBase
200  {
204  TrilinosScalar
205  value() const;
206 
210  TrilinosScalar &
211  value();
212  };
213 
217  template <>
218  class Accessor<true> : public AccessorBase
219  {
220  public:
225  using MatrixType = const SparseMatrix;
226 
232 
237  template <bool Other>
238  Accessor(const Accessor<Other> &a);
239 
243  TrilinosScalar
244  value() const;
245 
246  private:
250  template <bool>
251  friend class Iterator;
252  };
253 
257  template <>
258  class Accessor<false> : public AccessorBase
259  {
260  class Reference
261  {
262  public:
266  Reference(const Accessor<false> &accessor);
267 
271  operator TrilinosScalar() const;
272 
276  const Reference &
277  operator=(const TrilinosScalar n) const;
278 
282  const Reference &
283  operator+=(const TrilinosScalar n) const;
284 
288  const Reference &
289  operator-=(const TrilinosScalar n) const;
290 
294  const Reference &
295  operator*=(const TrilinosScalar n) const;
296 
300  const Reference &
301  operator/=(const TrilinosScalar n) const;
302 
303  private:
308  Accessor &accessor;
309  };
310 
311  public:
317 
323 
327  Reference
328  value() const;
329 
330  private:
334  template <bool>
335  friend class Iterator;
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 
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
1882 
1888  IndexSet
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 
2036  size_type,
2037  size_type,
2038  << "The entry with index <" << arg1 << ',' << arg2
2039  << "> does not exist.");
2040 
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 
2053 
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 
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 
2155  };
2156 
2157 
2158 
2159  // forwards declarations
2160  class SolverBase;
2161  class PreconditionBase;
2162 
2163  namespace internal
2164  {
2165  inline void
2166  check_vector_map_equality(const Epetra_CrsMatrix & mtrx,
2167  const Epetra_MultiVector &src,
2168  const Epetra_MultiVector &dst,
2169  const bool transpose)
2170  {
2171  if (transpose == false)
2172  {
2173  Assert(src.Map().SameAs(mtrx.DomainMap()) == true,
2174  ExcMessage(
2175  "Column map of matrix does not fit with vector map!"));
2176  Assert(dst.Map().SameAs(mtrx.RangeMap()) == true,
2177  ExcMessage("Row map of matrix does not fit with vector map!"));
2178  }
2179  else
2180  {
2181  Assert(src.Map().SameAs(mtrx.RangeMap()) == true,
2182  ExcMessage(
2183  "Column map of matrix does not fit with vector map!"));
2184  Assert(dst.Map().SameAs(mtrx.DomainMap()) == true,
2185  ExcMessage("Row map of matrix does not fit with vector map!"));
2186  }
2187  (void)mtrx; // removes -Wunused-variable in optimized mode
2188  (void)src;
2189  (void)dst;
2190  }
2191 
2192  inline void
2193  check_vector_map_equality(const Epetra_Operator & op,
2194  const Epetra_MultiVector &src,
2195  const Epetra_MultiVector &dst,
2196  const bool transpose)
2197  {
2198  if (transpose == false)
2199  {
2200  Assert(src.Map().SameAs(op.OperatorDomainMap()) == true,
2201  ExcMessage(
2202  "Column map of operator does not fit with vector map!"));
2203  Assert(dst.Map().SameAs(op.OperatorRangeMap()) == true,
2204  ExcMessage(
2205  "Row map of operator does not fit with vector map!"));
2206  }
2207  else
2208  {
2209  Assert(src.Map().SameAs(op.OperatorRangeMap()) == true,
2210  ExcMessage(
2211  "Column map of operator does not fit with vector map!"));
2212  Assert(dst.Map().SameAs(op.OperatorDomainMap()) == true,
2213  ExcMessage(
2214  "Row map of operator does not fit with vector map!"));
2215  }
2216  (void)op; // removes -Wunused-variable in optimized mode
2217  (void)src;
2218  (void)dst;
2219  }
2220 
2221  namespace LinearOperatorImplementation
2222  {
2243  class TrilinosPayload : public Epetra_Operator
2244  {
2245  public:
2249  using VectorType = Epetra_MultiVector;
2250 
2255 
2260 
2265 
2273  TrilinosPayload();
2274 
2278  TrilinosPayload(const TrilinosWrappers::SparseMatrix &matrix_exemplar,
2279  const TrilinosWrappers::SparseMatrix &matrix);
2280 
2285  const TrilinosWrappers::SparseMatrix & matrix_exemplar,
2286  const TrilinosWrappers::PreconditionBase &preconditioner);
2287 
2292  const TrilinosWrappers::PreconditionBase &preconditioner_exemplar,
2293  const TrilinosWrappers::PreconditionBase &preconditioner);
2294 
2298  TrilinosPayload(const TrilinosPayload &payload);
2299 
2307  TrilinosPayload(const TrilinosPayload &first_op,
2308  const TrilinosPayload &second_op);
2309 
2313  virtual ~TrilinosPayload() override = default;
2314 
2319  identity_payload() const;
2320 
2325  null_payload() const;
2326 
2331  transpose_payload() const;
2332 
2349  template <typename Solver, typename Preconditioner>
2350  typename std::enable_if<
2351  std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
2352  std::is_base_of<TrilinosWrappers::PreconditionBase,
2353  Preconditioner>::value,
2354  TrilinosPayload>::type
2355  inverse_payload(Solver &, const Preconditioner &) const;
2356 
2374  template <typename Solver, typename Preconditioner>
2375  typename std::enable_if<
2376  !(std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
2377  std::is_base_of<TrilinosWrappers::PreconditionBase,
2378  Preconditioner>::value),
2379  TrilinosPayload>::type
2380  inverse_payload(Solver &, const Preconditioner &) const;
2381 
2383 
2388 
2394  IndexSet
2396 
2402  IndexSet
2404 
2408  MPI_Comm
2409  get_mpi_communicator() const;
2410 
2417  void
2418  transpose();
2419 
2427  std::function<void(VectorType &, const VectorType &)> vmult;
2428 
2436  std::function<void(VectorType &, const VectorType &)> Tvmult;
2437 
2446  std::function<void(VectorType &, const VectorType &)> inv_vmult;
2447 
2456  std::function<void(VectorType &, const VectorType &)> inv_Tvmult;
2457 
2459 
2464 
2471  virtual bool
2472  UseTranspose() const override;
2473 
2489  virtual int
2490  SetUseTranspose(bool UseTranspose) override;
2491 
2503  virtual int
2504  Apply(const VectorType &X, VectorType &Y) const override;
2505 
2524  virtual int
2525  ApplyInverse(const VectorType &Y, VectorType &X) const override;
2527 
2532 
2539  virtual const char *
2540  Label() const override;
2541 
2549  virtual const Epetra_Comm &
2550  Comm() const override;
2551 
2559  virtual const Epetra_Map &
2560  OperatorDomainMap() const override;
2561 
2570  virtual const Epetra_Map &
2571  OperatorRangeMap() const override;
2573 
2574  private:
2580 
2585 # ifdef DEAL_II_WITH_MPI
2586  Epetra_MpiComm communicator;
2587 # else
2588  Epetra_SerialComm communicator;
2589 # endif
2590 
2595  Epetra_Map domain_map;
2596 
2601  Epetra_Map range_map;
2602 
2611  virtual bool
2612  HasNormInf() const override;
2613 
2621  virtual double
2622  NormInf() const override;
2623  };
2624 
2630  operator+(const TrilinosPayload &first_op,
2631  const TrilinosPayload &second_op);
2632 
2637  TrilinosPayload operator*(const TrilinosPayload &first_op,
2638  const TrilinosPayload &second_op);
2639 
2640  } // namespace LinearOperatorImplementation
2641  } /* namespace internal */
2642 
2643 
2644 
2645  // ----------------------- inline and template functions --------------------
2646 
2647 # ifndef DOXYGEN
2648 
2649  namespace SparseMatrixIterators
2650  {
2652  size_type row,
2653  size_type index)
2654  : matrix(matrix)
2655  , a_row(row)
2656  , a_index(index)
2657  {
2658  visit_present_row();
2659  }
2660 
2661 
2662  inline AccessorBase::size_type
2663  AccessorBase::row() const
2664  {
2665  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2666  return a_row;
2667  }
2668 
2669 
2670  inline AccessorBase::size_type
2671  AccessorBase::column() const
2672  {
2673  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2674  return (*colnum_cache)[a_index];
2675  }
2676 
2677 
2678  inline AccessorBase::size_type
2679  AccessorBase::index() const
2680  {
2681  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2682  return a_index;
2683  }
2684 
2685 
2686  inline Accessor<true>::Accessor(MatrixType * matrix,
2687  const size_type row,
2688  const size_type index)
2689  : AccessorBase(const_cast<SparseMatrix *>(matrix), row, index)
2690  {}
2691 
2692 
2693  template <bool Other>
2694  inline Accessor<true>::Accessor(const Accessor<Other> &other)
2695  : AccessorBase(other)
2696  {}
2697 
2698 
2699  inline TrilinosScalar
2700  Accessor<true>::value() const
2701  {
2702  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2703  return (*value_cache)[a_index];
2704  }
2705 
2706 
2707  inline Accessor<false>::Reference::Reference(const Accessor<false> &acc)
2708  : accessor(const_cast<Accessor<false> &>(acc))
2709  {}
2710 
2711 
2712  inline Accessor<false>::Reference::operator TrilinosScalar() const
2713  {
2714  return (*accessor.value_cache)[accessor.a_index];
2715  }
2716 
2717  inline const Accessor<false>::Reference &
2718  Accessor<false>::Reference::operator=(const TrilinosScalar n) const
2719  {
2720  (*accessor.value_cache)[accessor.a_index] = n;
2721  accessor.matrix->set(accessor.row(),
2722  accessor.column(),
2723  static_cast<TrilinosScalar>(*this));
2724  return *this;
2725  }
2726 
2727 
2728  inline const Accessor<false>::Reference &
2729  Accessor<false>::Reference::operator+=(const TrilinosScalar n) const
2730  {
2731  (*accessor.value_cache)[accessor.a_index] += n;
2732  accessor.matrix->set(accessor.row(),
2733  accessor.column(),
2734  static_cast<TrilinosScalar>(*this));
2735  return *this;
2736  }
2737 
2738 
2739  inline const Accessor<false>::Reference &
2740  Accessor<false>::Reference::operator-=(const TrilinosScalar n) const
2741  {
2742  (*accessor.value_cache)[accessor.a_index] -= n;
2743  accessor.matrix->set(accessor.row(),
2744  accessor.column(),
2745  static_cast<TrilinosScalar>(*this));
2746  return *this;
2747  }
2748 
2749 
2750  inline const Accessor<false>::Reference &
2751  Accessor<false>::Reference::operator*=(const TrilinosScalar n) const
2752  {
2753  (*accessor.value_cache)[accessor.a_index] *= n;
2754  accessor.matrix->set(accessor.row(),
2755  accessor.column(),
2756  static_cast<TrilinosScalar>(*this));
2757  return *this;
2758  }
2759 
2760 
2761  inline const Accessor<false>::Reference &
2762  Accessor<false>::Reference::operator/=(const TrilinosScalar n) const
2763  {
2764  (*accessor.value_cache)[accessor.a_index] /= n;
2765  accessor.matrix->set(accessor.row(),
2766  accessor.column(),
2767  static_cast<TrilinosScalar>(*this));
2768  return *this;
2769  }
2770 
2771 
2772  inline Accessor<false>::Accessor(MatrixType * matrix,
2773  const size_type row,
2774  const size_type index)
2775  : AccessorBase(matrix, row, index)
2776  {}
2777 
2778 
2779  inline Accessor<false>::Reference
2780  Accessor<false>::value() const
2781  {
2782  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2783  return {*this};
2784  }
2785 
2786 
2787 
2788  template <bool Constness>
2789  inline Iterator<Constness>::Iterator(MatrixType * matrix,
2790  const size_type row,
2791  const size_type index)
2792  : accessor(matrix, row, index)
2793  {}
2794 
2795 
2796  template <bool Constness>
2797  template <bool Other>
2798  inline Iterator<Constness>::Iterator(const Iterator<Other> &other)
2799  : accessor(other.accessor)
2800  {}
2801 
2802 
2803  template <bool Constness>
2804  inline Iterator<Constness> &
2805  Iterator<Constness>::operator++()
2806  {
2807  Assert(accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
2808 
2809  ++accessor.a_index;
2810 
2811  // If at end of line: do one
2812  // step, then cycle until we
2813  // find a row with a nonzero
2814  // number of entries.
2815  if (accessor.a_index >= accessor.colnum_cache->size())
2816  {
2817  accessor.a_index = 0;
2818  ++accessor.a_row;
2819 
2820  while ((accessor.a_row < accessor.matrix->m()) &&
2821  ((accessor.matrix->in_local_range(accessor.a_row) == false) ||
2822  (accessor.matrix->row_length(accessor.a_row) == 0)))
2823  ++accessor.a_row;
2824 
2825  accessor.visit_present_row();
2826  }
2827  return *this;
2828  }
2829 
2830 
2831  template <bool Constness>
2832  inline Iterator<Constness>
2833  Iterator<Constness>::operator++(int)
2834  {
2835  const Iterator<Constness> old_state = *this;
2836  ++(*this);
2837  return old_state;
2838  }
2839 
2840 
2841 
2842  template <bool Constness>
2843  inline const Accessor<Constness> &Iterator<Constness>::operator*() const
2844  {
2845  return accessor;
2846  }
2847 
2848 
2849 
2850  template <bool Constness>
2851  inline const Accessor<Constness> *Iterator<Constness>::operator->() const
2852  {
2853  return &accessor;
2854  }
2855 
2856 
2857 
2858  template <bool Constness>
2859  inline bool
2860  Iterator<Constness>::operator==(const Iterator<Constness> &other) const
2861  {
2862  return (accessor.a_row == other.accessor.a_row &&
2863  accessor.a_index == other.accessor.a_index);
2864  }
2865 
2866 
2867 
2868  template <bool Constness>
2869  inline bool
2870  Iterator<Constness>::operator!=(const Iterator<Constness> &other) const
2871  {
2872  return !(*this == other);
2873  }
2874 
2875 
2876 
2877  template <bool Constness>
2878  inline bool
2879  Iterator<Constness>::operator<(const Iterator<Constness> &other) const
2880  {
2881  return (accessor.row() < other.accessor.row() ||
2882  (accessor.row() == other.accessor.row() &&
2883  accessor.index() < other.accessor.index()));
2884  }
2885 
2886 
2887  template <bool Constness>
2888  inline bool
2889  Iterator<Constness>::operator>(const Iterator<Constness> &other) const
2890  {
2891  return (other < *this);
2892  }
2893 
2894  } // namespace SparseMatrixIterators
2895 
2896 
2897 
2899  SparseMatrix::begin() const
2900  {
2901  return begin(0);
2902  }
2903 
2904 
2905 
2907  SparseMatrix::end() const
2908  {
2909  return const_iterator(this, m(), 0);
2910  }
2911 
2912 
2913 
2915  SparseMatrix::begin(const size_type r) const
2916  {
2917  Assert(r < m(), ExcIndexRange(r, 0, m()));
2918  if (in_local_range(r) && (row_length(r) > 0))
2919  return const_iterator(this, r, 0);
2920  else
2921  return end(r);
2922  }
2923 
2924 
2925 
2927  SparseMatrix::end(const size_type r) const
2928  {
2929  Assert(r < m(), ExcIndexRange(r, 0, m()));
2930 
2931  // place the iterator on the first entry
2932  // past this line, or at the end of the
2933  // matrix
2934  for (size_type i = r + 1; i < m(); ++i)
2935  if (in_local_range(i) && (row_length(i) > 0))
2936  return const_iterator(this, i, 0);
2937 
2938  // if there is no such line, then take the
2939  // end iterator of the matrix
2940  return end();
2941  }
2942 
2943 
2944 
2945  inline SparseMatrix::iterator
2947  {
2948  return begin(0);
2949  }
2950 
2951 
2952 
2953  inline SparseMatrix::iterator
2955  {
2956  return iterator(this, m(), 0);
2957  }
2958 
2959 
2960 
2961  inline SparseMatrix::iterator
2962  SparseMatrix::begin(const size_type r)
2963  {
2964  Assert(r < m(), ExcIndexRange(r, 0, m()));
2965  if (in_local_range(r) && (row_length(r) > 0))
2966  return iterator(this, r, 0);
2967  else
2968  return end(r);
2969  }
2970 
2971 
2972 
2973  inline SparseMatrix::iterator
2974  SparseMatrix::end(const size_type r)
2975  {
2976  Assert(r < m(), ExcIndexRange(r, 0, m()));
2977 
2978  // place the iterator on the first entry
2979  // past this line, or at the end of the
2980  // matrix
2981  for (size_type i = r + 1; i < m(); ++i)
2982  if (in_local_range(i) && (row_length(i) > 0))
2983  return iterator(this, i, 0);
2984 
2985  // if there is no such line, then take the
2986  // end iterator of the matrix
2987  return end();
2988  }
2989 
2990 
2991 
2992  inline bool
2993  SparseMatrix::in_local_range(const size_type index) const
2994  {
2995  TrilinosWrappers::types::int_type begin, end;
2996 # ifndef DEAL_II_WITH_64BIT_INDICES
2997  begin = matrix->RowMap().MinMyGID();
2998  end = matrix->RowMap().MaxMyGID() + 1;
2999 # else
3000  begin = matrix->RowMap().MinMyGID64();
3001  end = matrix->RowMap().MaxMyGID64() + 1;
3002 # endif
3003 
3004  return ((index >= static_cast<size_type>(begin)) &&
3005  (index < static_cast<size_type>(end)));
3006  }
3007 
3008 
3009 
3010  inline bool
3012  {
3013  return compressed;
3014  }
3015 
3016 
3017 
3018  // Inline the set() and add() functions, since they will be called
3019  // frequently, and the compiler can optimize away some unnecessary loops
3020  // when the sizes are given at compile time.
3021  template <>
3022  void
3023  SparseMatrix::set<TrilinosScalar>(const size_type row,
3024  const size_type n_cols,
3025  const size_type * col_indices,
3026  const TrilinosScalar *values,
3027  const bool elide_zero_values);
3028 
3029 
3030 
3031  template <typename Number>
3032  void
3033  SparseMatrix::set(const size_type row,
3034  const size_type n_cols,
3035  const size_type *col_indices,
3036  const Number * values,
3037  const bool elide_zero_values)
3038  {
3039  std::vector<TrilinosScalar> trilinos_values(n_cols);
3040  std::copy(values, values + n_cols, trilinos_values.begin());
3041  this->set(
3042  row, n_cols, col_indices, trilinos_values.data(), elide_zero_values);
3043  }
3044 
3045 
3046 
3047  inline void
3048  SparseMatrix::set(const size_type i,
3049  const size_type j,
3050  const TrilinosScalar value)
3051  {
3052  AssertIsFinite(value);
3053 
3054  set(i, 1, &j, &value, false);
3055  }
3056 
3057 
3058 
3059  inline void
3060  SparseMatrix::set(const std::vector<size_type> & indices,
3061  const FullMatrix<TrilinosScalar> &values,
3062  const bool elide_zero_values)
3063  {
3064  Assert(indices.size() == values.m(),
3065  ExcDimensionMismatch(indices.size(), values.m()));
3066  Assert(values.m() == values.n(), ExcNotQuadratic());
3067 
3068  for (size_type i = 0; i < indices.size(); ++i)
3069  set(indices[i],
3070  indices.size(),
3071  indices.data(),
3072  &values(i, 0),
3073  elide_zero_values);
3074  }
3075 
3076 
3077 
3078  inline void
3079  SparseMatrix::add(const size_type i,
3080  const size_type j,
3081  const TrilinosScalar value)
3082  {
3083  AssertIsFinite(value);
3084 
3085  if (value == 0)
3086  {
3087  // we have to check after Insert/Add in any case to be consistent
3088  // with the MPI communication model, but we can save some
3089  // work if the addend is zero. However, these actions are done in case
3090  // we pass on to the other function.
3091 
3092  // TODO: fix this (do not run compress here, but fail)
3093  if (last_action == Insert)
3094  {
3095  int ierr;
3096  ierr = matrix->GlobalAssemble(*column_space_map,
3097  matrix->RowMap(),
3098  false);
3099 
3100  Assert(ierr == 0, ExcTrilinosError(ierr));
3101  (void)ierr; // removes -Wunused-but-set-variable in optimized mode
3102  }
3103 
3104  last_action = Add;
3105 
3106  return;
3107  }
3108  else
3109  add(i, 1, &j, &value, false);
3110  }
3111 
3112 
3113 
3114  // inline "simple" functions that are called frequently and do only involve
3115  // a call to some Trilinos function.
3117  SparseMatrix::m() const
3118  {
3119 # ifndef DEAL_II_WITH_64BIT_INDICES
3120  return matrix->NumGlobalRows();
3121 # else
3122  return matrix->NumGlobalRows64();
3123 # endif
3124  }
3125 
3126 
3127 
3129  SparseMatrix::n() const
3130  {
3131  // If the matrix structure has not been fixed (i.e., we did not have a
3132  // sparsity pattern), it does not know about the number of columns so we
3133  // must always take this from the additional column space map
3134  Assert(column_space_map.get() != nullptr, ExcInternalError());
3135 # ifndef DEAL_II_WITH_64BIT_INDICES
3136  return column_space_map->NumGlobalElements();
3137 # else
3138  return column_space_map->NumGlobalElements64();
3139 # endif
3140  }
3141 
3142 
3143 
3144  inline unsigned int
3145  SparseMatrix::local_size() const
3146  {
3147  return matrix->NumMyRows();
3148  }
3149 
3150 
3151 
3152  inline std::pair<SparseMatrix::size_type, SparseMatrix::size_type>
3154  {
3155  size_type begin, end;
3156 # ifndef DEAL_II_WITH_64BIT_INDICES
3157  begin = matrix->RowMap().MinMyGID();
3158  end = matrix->RowMap().MaxMyGID() + 1;
3159 # else
3160  begin = matrix->RowMap().MinMyGID64();
3161  end = matrix->RowMap().MaxMyGID64() + 1;
3162 # endif
3163 
3164  return std::make_pair(begin, end);
3165  }
3166 
3167 
3168 
3171  {
3172 # ifndef DEAL_II_WITH_64BIT_INDICES
3173  return matrix->NumGlobalNonzeros();
3174 # else
3175  return matrix->NumGlobalNonzeros64();
3176 # endif
3177  }
3178 
3179 
3180 
3181  template <typename SparsityPatternType>
3182  inline void
3183  SparseMatrix::reinit(const IndexSet & parallel_partitioning,
3184  const SparsityPatternType &sparsity_pattern,
3185  const MPI_Comm & communicator,
3186  const bool exchange_data)
3187  {
3188  reinit(parallel_partitioning,
3189  parallel_partitioning,
3190  sparsity_pattern,
3191  communicator,
3192  exchange_data);
3193  }
3194 
3195 
3196 
3197  template <typename number>
3198  inline void
3199  SparseMatrix::reinit(const IndexSet &parallel_partitioning,
3200  const ::SparseMatrix<number> &sparse_matrix,
3201  const MPI_Comm & communicator,
3202  const double drop_tolerance,
3203  const bool copy_values,
3204  const ::SparsityPattern * use_this_sparsity)
3205  {
3206  Epetra_Map map =
3207  parallel_partitioning.make_trilinos_map(communicator, false);
3208  reinit(parallel_partitioning,
3209  parallel_partitioning,
3210  sparse_matrix,
3211  drop_tolerance,
3212  copy_values,
3213  use_this_sparsity);
3214  }
3215 
3216 
3217 
3218  inline const Epetra_CrsMatrix &
3220  {
3221  return static_cast<const Epetra_CrsMatrix &>(*matrix);
3222  }
3223 
3224 
3225 
3226  inline const Epetra_CrsGraph &
3228  {
3229  return matrix->Graph();
3230  }
3231 
3232 
3233 
3234  inline IndexSet
3236  {
3237  return IndexSet(matrix->DomainMap());
3238  }
3239 
3240 
3241 
3242  inline IndexSet
3244  {
3245  return IndexSet(matrix->RangeMap());
3246  }
3247 
3248 
3249 
3250  inline void
3252  {
3253  // nothing to do here
3254  }
3255 
3256 
3257 
3258  inline void
3260  {
3261  // nothing to do here
3262  }
3263 
3264 
3265  namespace internal
3266  {
3267  namespace LinearOperatorImplementation
3268  {
3269  template <typename Solver, typename Preconditioner>
3270  typename std::enable_if<
3271  std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
3272  std::is_base_of<TrilinosWrappers::PreconditionBase,
3273  Preconditioner>::value,
3274  TrilinosPayload>::type
3275  TrilinosPayload::inverse_payload(
3276  Solver & solver,
3277  const Preconditioner &preconditioner) const
3278  {
3279  const auto &payload = *this;
3280 
3281  TrilinosPayload return_op(payload);
3282 
3283  // Capture by copy so the payloads are always valid
3284 
3285  return_op.inv_vmult = [payload, &solver, &preconditioner](
3286  TrilinosPayload::Domain & tril_dst,
3287  const TrilinosPayload::Range &tril_src) {
3288  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3289  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3290  Assert(&tril_src != &tril_dst,
3292  internal::check_vector_map_equality(payload,
3293  tril_src,
3294  tril_dst,
3295  !payload.UseTranspose());
3296  solver.solve(payload, tril_dst, tril_src, preconditioner);
3297  };
3298 
3299  return_op.inv_Tvmult = [payload, &solver, &preconditioner](
3300  TrilinosPayload::Range & tril_dst,
3301  const TrilinosPayload::Domain &tril_src) {
3302  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3303  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3304  Assert(&tril_src != &tril_dst,
3306  internal::check_vector_map_equality(payload,
3307  tril_src,
3308  tril_dst,
3309  payload.UseTranspose());
3310 
3311  const_cast<TrilinosPayload &>(payload).transpose();
3312  solver.solve(payload, tril_dst, tril_src, preconditioner);
3313  const_cast<TrilinosPayload &>(payload).transpose();
3314  };
3315 
3316  // If the input operator is already setup for transpose operations, then
3317  // we must do similar with its inverse.
3318  if (return_op.UseTranspose() == true)
3319  std::swap(return_op.inv_vmult, return_op.inv_Tvmult);
3320 
3321  return return_op;
3322  }
3323 
3324  template <typename Solver, typename Preconditioner>
3325  typename std::enable_if<
3326  !(std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
3327  std::is_base_of<TrilinosWrappers::PreconditionBase,
3328  Preconditioner>::value),
3329  TrilinosPayload>::type
3330  TrilinosPayload::inverse_payload(Solver &, const Preconditioner &) const
3331  {
3332  TrilinosPayload return_op(*this);
3333 
3334  return_op.inv_vmult = [](TrilinosPayload::Domain &,
3335  const TrilinosPayload::Range &) {
3336  AssertThrow(false,
3337  ExcMessage("Payload inv_vmult disabled because of "
3338  "incompatible solver/preconditioner choice."));
3339  };
3340 
3341  return_op.inv_Tvmult = [](TrilinosPayload::Range &,
3342  const TrilinosPayload::Domain &) {
3343  AssertThrow(false,
3344  ExcMessage("Payload inv_vmult disabled because of "
3345  "incompatible solver/preconditioner choice."));
3346  };
3347 
3348  return return_op;
3349  }
3350  } // namespace LinearOperatorImplementation
3351  } // namespace internal
3352 
3353  template <>
3354  void
3355  SparseMatrix::set<TrilinosScalar>(const size_type row,
3356  const size_type n_cols,
3357  const size_type * col_indices,
3358  const TrilinosScalar *values,
3359  const bool elide_zero_values);
3360 # endif // DOXYGEN
3361 
3362 } /* namespace TrilinosWrappers */
3363 
3364 
3365 DEAL_II_NAMESPACE_CLOSE
3366 
3367 
3368 # endif // DEAL_II_WITH_TRILINOS
3369 
3370 
3371 /*----------------------- trilinos_sparse_matrix.h --------------------*/
3372 
3373 #endif
3374 /*----------------------- trilinos_sparse_matrix.h --------------------*/
TrilinosScalar diag_element(const size_type i) const
static ::ExceptionBase & ExcTrilinosError(int arg1)
TrilinosScalar el(const size_type i, const size_type j) const
size_type m() const
const_iterator end() const
const Epetra_Map & col_partitioner() const
bool operator>(const Iterator< Constness > &) const
void mmult(SparseMatrix &C, const SparseMatrix &B, const MPI::Vector &V=MPI::Vector()) const
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:541
Contents is actually a matrix.
static ::ExceptionBase & ExcSourceEqualsDestination()
static ::ExceptionBase & ExcBeyondEndOfMatrix()
bool operator==(const Iterator< Constness > &) const
void clear_row(const size_type row, const TrilinosScalar new_diag_value=0)
const Epetra_Map & range_partitioner() const
void vmult_add(VectorType &dst, const VectorType &src) const
std::function< void(VectorType &, const VectorType &)> inv_Tvmult
TrilinosScalar operator()(const size_type i, const size_type j) const
static ::ExceptionBase & ExcInvalidIndex(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcAccessToNonlocalRow(std::size_t arg1, std::size_t arg2, std::size_t arg3)
const Epetra_Map & row_partitioner() const
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
static ::ExceptionBase & ExcAccessToNonPresentElement(size_type arg1, size_type arg2)
SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
TrilinosScalar residual(MPI::Vector &dst, const MPI::Vector &x, const MPI::Vector &b) const
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
static ::ExceptionBase & ExcTrilinosError(int arg1)
SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
TrilinosScalar matrix_scalar_product(const MPI::Vector &u, const MPI::Vector &v) const
void clear_rows(const std::vector< size_type > &rows, const TrilinosScalar new_diag_value=0)
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
SparseMatrix & operator*=(const TrilinosScalar factor)
std::unique_ptr< Epetra_FECrsMatrix > matrix
virtual int Apply(const VectorType &X, VectorType &Y) const override
SparseMatrix & operator=(const SparseMatrix &)=delete
static ::ExceptionBase & ExcAccessToNonlocalRow(std::size_t arg1)
::types::global_dof_index size_type
bool in_local_range(const size_type index) const
static ::ExceptionBase & ExcMessage(std::string arg1)
const Accessor< Constness > & operator*() const
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
void print(std::ostream &out, const bool write_extended_trilinos_info=false) const
bool operator!=(const Iterator< Constness > &) const
#define DeclException0(Exception0)
Definition: exceptions.h:473
typename Accessor< Constness >::MatrixType MatrixType
std::function< void(VectorType &, const VectorType &)> inv_vmult
TrilinosScalar matrix_norm_square(const MPI::Vector &v) const
const Accessor< Constness > * operator->() const
virtual int ApplyInverse(const VectorType &Y, VectorType &X) const override
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
SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
AccessorBase(SparseMatrix *matrix, const size_type row, const size_type index)
static ::ExceptionBase & ExcInvalidIndexWithinRow(size_type arg1, size_type arg2)
virtual ~SparseMatrix() override=default
IndexSet locally_owned_range_indices() const
size_type n_nonzero_elements() const
std::enable_if< std::is_same< typename VectorType::value_type, TrilinosScalar >::value >::type Tvmult(VectorType &dst, const VectorType &src) const
static ::ExceptionBase & ExcIteratorPastEnd()
types::global_dof_index size_type
const_iterator begin() const
static ::ExceptionBase & ExcNotQuadratic()
void swap(Vector< Number > &u, Vector< Number > &v)
Definition: vector.h:1376
unsigned int row_length(const size_type row) const
void Tvmult_add(VectorType &dst, const VectorType &src) const
unsigned int global_dof_index
Definition: types.h:89
std::enable_if< std::is_same< typename VectorType::value_type, TrilinosScalar >::value >::type vmult(VectorType &dst, const VectorType &src) const
void reinit(const SparsityPatternType &sparsity_pattern)
void copy_from(const SparseMatrix &source)
const Epetra_Map & domain_partitioner() const
SparseMatrix & operator/=(const TrilinosScalar factor)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:587
std::pair< size_type, size_type > local_range() const
static ::ExceptionBase & ExcMatrixNotCompressed()
IndexSet locally_owned_domain_indices() const
const Epetra_CrsGraph & trilinos_sparsity_pattern() const
void Tmmult(SparseMatrix &C, const SparseMatrix &B, const MPI::Vector &V=MPI::Vector()) const
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:564
std::unique_ptr< Epetra_Map > column_space_map
void compress(::VectorOperation::values operation)
static ::ExceptionBase & ExcAccessToNonLocalElement(size_type arg1, size_type arg2, size_type arg3, size_type arg4)
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