Reference documentation for deal.II version Git b3b9f019ec 2021-07-26 17:39:47 -0400
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
trilinos_sparse_matrix.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 2021 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>
26 
27 # include <deal.II/lac/exceptions.h>
28 # include <deal.II/lac/full_matrix.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_MpiComm.h>
42 # include <Epetra_MultiVector.h>
43 # include <Epetra_Operator.h>
44 # include <mpi.h>
45 
46 # include <cmath>
47 # include <memory>
48 # include <type_traits>
49 # include <vector>
50 
52 
53 // forward declarations
54 # ifndef DOXYGEN
55 template <typename MatrixType>
56 class BlockMatrixBase;
57 
58 template <typename number>
59 class SparseMatrix;
60 class SparsityPattern;
62 
63 namespace TrilinosWrappers
64 {
65  class SparseMatrix;
66  class SparsityPattern;
67 
68  namespace SparseMatrixIterators
69  {
70  template <bool Constness>
71  class Iterator;
72  }
73 } // namespace TrilinosWrappers
74 # endif
75 
76 namespace TrilinosWrappers
77 {
82  {
87 
92  std::size_t,
93  std::size_t,
94  std::size_t,
95  << "You tried to access row " << arg1
96  << " of a distributed sparsity pattern, "
97  << " but only rows " << arg2 << " through " << arg3
98  << " are stored locally and can be accessed.");
99 
108  {
109  public:
114 
119  const size_type row,
120  const size_type index);
121 
125  size_type
126  row() const;
127 
131  size_type
132  index() const;
133 
137  size_type
138  column() const;
139 
140  protected:
152 
157 
163  void
164  visit_present_row();
165 
178  std::shared_ptr<std::vector<size_type>> colnum_cache;
179 
183  std::shared_ptr<std::vector<TrilinosScalar>> value_cache;
184  };
185 
196  template <bool Constess>
197  class Accessor : public AccessorBase
198  {
203  value() const;
204 
209  value();
210  };
211 
215  template <>
216  class Accessor<true> : public AccessorBase
217  {
218  public:
223  using MatrixType = const SparseMatrix;
224 
229  Accessor(MatrixType *matrix, const size_type row, const size_type index);
230 
235  template <bool Other>
236  Accessor(const Accessor<Other> &a);
237 
242  value() const;
243 
244  private:
245  // Make iterator class a friend.
246  template <bool>
247  friend class Iterator;
248  };
249 
253  template <>
254  class Accessor<false> : public AccessorBase
255  {
256  class Reference
257  {
258  public:
262  Reference(const Accessor<false> &accessor);
263 
267  operator TrilinosScalar() const;
268 
272  const Reference &
273  operator=(const TrilinosScalar n) const;
274 
278  const Reference &
279  operator+=(const TrilinosScalar n) const;
280 
284  const Reference &
285  operator-=(const TrilinosScalar n) const;
286 
290  const Reference &
291  operator*=(const TrilinosScalar n) const;
292 
296  const Reference &
297  operator/=(const TrilinosScalar n) const;
298 
299  private:
305  };
306 
307  public:
313 
318  Accessor(MatrixType *matrix, const size_type row, const size_type index);
319 
323  Reference
324  value() const;
325 
326  private:
327  // Make iterator class a friend.
328  template <bool>
329  friend class Iterator;
330 
331  // Make Reference object a friend.
332  friend class Reference;
333  };
334 
348  template <bool Constness>
349  class Iterator
350  {
351  public:
356 
362 
367  Iterator(MatrixType *matrix, const size_type row, const size_type index);
368 
372  template <bool Other>
373  Iterator(const Iterator<Other> &other);
374 
379  operator++();
380 
385  operator++(int);
386 
390  const Accessor<Constness> &operator*() const;
391 
395  const Accessor<Constness> *operator->() const;
396 
401  bool
402  operator==(const Iterator<Constness> &) const;
403 
407  bool
408  operator!=(const Iterator<Constness> &) const;
409 
415  bool
416  operator<(const Iterator<Constness> &) const;
417 
421  bool
422  operator>(const Iterator<Constness> &) const;
423 
427  DeclException2(ExcInvalidIndexWithinRow,
428  size_type,
429  size_type,
430  << "Attempt to access element " << arg2 << " of row "
431  << arg1 << " which doesn't have that many elements.");
432 
433  private:
438 
439  template <bool Other>
440  friend class Iterator;
441  };
442 
443  } // namespace SparseMatrixIterators
444 
445 
506  class SparseMatrix : public Subscriptor
507  {
508  public:
513 
518  std::size_t,
519  << "You tried to access row " << arg1
520  << " of a non-contiguous locally owned row set."
521  << " The row " << arg1
522  << " is not stored locally and can't be accessed.");
523 
531  struct Traits
532  {
537  static const bool zero_addition_can_be_elided = true;
538  };
539 
544 
549 
554 
562  SparseMatrix();
563 
571  SparseMatrix(const size_type m,
572  const size_type n,
573  const unsigned int n_max_entries_per_row);
574 
582  SparseMatrix(const size_type m,
583  const size_type n,
584  const std::vector<unsigned int> &n_entries_per_row);
585 
589  SparseMatrix(const SparsityPattern &InputSparsityPattern);
590 
595  SparseMatrix(SparseMatrix &&other) noexcept;
596 
600  SparseMatrix(const SparseMatrix &) = delete;
601 
605  SparseMatrix &
606  operator=(const SparseMatrix &) = delete;
607 
611  virtual ~SparseMatrix() override = default;
612 
628  template <typename SparsityPatternType>
629  void
630  reinit(const SparsityPatternType &sparsity_pattern);
631 
644  void
645  reinit(const SparsityPattern &sparsity_pattern);
646 
655  void
656  reinit(const SparseMatrix &sparse_matrix);
657 
678  template <typename number>
679  void
680  reinit(const ::SparseMatrix<number> &dealii_sparse_matrix,
681  const double drop_tolerance = 1e-13,
682  const bool copy_values = true,
683  const ::SparsityPattern * use_this_sparsity = nullptr);
684 
690  void
691  reinit(const Epetra_CrsMatrix &input_matrix, const bool copy_values = true);
693 
710  SparseMatrix(const IndexSet & parallel_partitioning,
711  const MPI_Comm & communicator = MPI_COMM_WORLD,
712  const unsigned int n_max_entries_per_row = 0);
713 
721  SparseMatrix(const IndexSet & parallel_partitioning,
722  const MPI_Comm & communicator,
723  const std::vector<unsigned int> &n_entries_per_row);
724 
739  SparseMatrix(const IndexSet &row_parallel_partitioning,
740  const IndexSet &col_parallel_partitioning,
741  const MPI_Comm &communicator = MPI_COMM_WORLD,
742  const size_type n_max_entries_per_row = 0);
743 
758  SparseMatrix(const IndexSet & row_parallel_partitioning,
759  const IndexSet & col_parallel_partitioning,
760  const MPI_Comm & communicator,
761  const std::vector<unsigned int> &n_entries_per_row);
762 
783  template <typename SparsityPatternType>
784  void
785  reinit(const IndexSet & parallel_partitioning,
786  const SparsityPatternType &sparsity_pattern,
787  const MPI_Comm & communicator = MPI_COMM_WORLD,
788  const bool exchange_data = false);
789 
802  template <typename SparsityPatternType>
803  typename std::enable_if<
804  !std::is_same<SparsityPatternType,
805  ::SparseMatrix<double>>::value>::type
806  reinit(const IndexSet & row_parallel_partitioning,
807  const IndexSet & col_parallel_partitioning,
808  const SparsityPatternType &sparsity_pattern,
809  const MPI_Comm & communicator = MPI_COMM_WORLD,
810  const bool exchange_data = false);
811 
828  template <typename number>
829  void
830  reinit(const IndexSet & parallel_partitioning,
831  const ::SparseMatrix<number> &dealii_sparse_matrix,
832  const MPI_Comm & communicator = MPI_COMM_WORLD,
833  const double drop_tolerance = 1e-13,
834  const bool copy_values = true,
835  const ::SparsityPattern * use_this_sparsity = nullptr);
836 
850  template <typename number>
851  void
852  reinit(const IndexSet & row_parallel_partitioning,
853  const IndexSet & col_parallel_partitioning,
854  const ::SparseMatrix<number> &dealii_sparse_matrix,
855  const MPI_Comm & communicator = MPI_COMM_WORLD,
856  const double drop_tolerance = 1e-13,
857  const bool copy_values = true,
858  const ::SparsityPattern * use_this_sparsity = nullptr);
860 
864 
868  size_type
869  m() const;
870 
874  size_type
875  n() const;
876 
885  unsigned int
886  local_size() const;
887 
896  std::pair<size_type, size_type>
897  local_range() const;
898 
903  bool
904  in_local_range(const size_type index) const;
905 
910  size_type
911  n_nonzero_elements() const;
912 
916  unsigned int
917  row_length(const size_type row) const;
918 
925  bool
926  is_compressed() const;
927 
933  size_type
934  memory_consumption() const;
935 
939  MPI_Comm
940  get_mpi_communicator() const;
941 
943 
947 
957  SparseMatrix &
958  operator=(const double d);
959 
967  void
968  clear();
969 
997  void
998  compress(::VectorOperation::values operation);
999 
1021  void
1022  set(const size_type i, const size_type j, const TrilinosScalar value);
1023 
1056  void
1057  set(const std::vector<size_type> & indices,
1058  const FullMatrix<TrilinosScalar> &full_matrix,
1059  const bool elide_zero_values = false);
1060 
1066  void
1067  set(const std::vector<size_type> & row_indices,
1068  const std::vector<size_type> & col_indices,
1069  const FullMatrix<TrilinosScalar> &full_matrix,
1070  const bool elide_zero_values = false);
1071 
1099  void
1100  set(const size_type row,
1101  const std::vector<size_type> & col_indices,
1102  const std::vector<TrilinosScalar> &values,
1103  const bool elide_zero_values = false);
1104 
1132  template <typename Number>
1133  void
1134  set(const size_type row,
1135  const size_type n_cols,
1136  const size_type *col_indices,
1137  const Number * values,
1138  const bool elide_zero_values = false);
1139 
1149  void
1150  add(const size_type i, const size_type j, const TrilinosScalar value);
1151 
1170  void
1171  add(const std::vector<size_type> & indices,
1172  const FullMatrix<TrilinosScalar> &full_matrix,
1173  const bool elide_zero_values = true);
1174 
1180  void
1181  add(const std::vector<size_type> & row_indices,
1182  const std::vector<size_type> & col_indices,
1183  const FullMatrix<TrilinosScalar> &full_matrix,
1184  const bool elide_zero_values = true);
1185 
1199  void
1200  add(const size_type row,
1201  const std::vector<size_type> & col_indices,
1202  const std::vector<TrilinosScalar> &values,
1203  const bool elide_zero_values = true);
1204 
1218  void
1219  add(const size_type row,
1220  const size_type n_cols,
1221  const size_type * col_indices,
1222  const TrilinosScalar *values,
1223  const bool elide_zero_values = true,
1224  const bool col_indices_are_sorted = false);
1225 
1229  SparseMatrix &
1230  operator*=(const TrilinosScalar factor);
1231 
1235  SparseMatrix &
1236  operator/=(const TrilinosScalar factor);
1237 
1241  void
1242  copy_from(const SparseMatrix &source);
1243 
1251  void
1252  add(const TrilinosScalar factor, const SparseMatrix &matrix);
1253 
1280  void
1281  clear_row(const size_type row, const TrilinosScalar new_diag_value = 0);
1282 
1303  void
1304  clear_rows(const std::vector<size_type> &rows,
1305  const TrilinosScalar new_diag_value = 0);
1306 
1316  void
1317  transpose();
1318 
1320 
1324 
1334  operator()(const size_type i, const size_type j) const;
1335 
1353  el(const size_type i, const size_type j) const;
1354 
1362  diag_element(const size_type i) const;
1363 
1365 
1369 
1400  template <typename VectorType>
1401  typename std::enable_if<std::is_same<typename VectorType::value_type,
1402  TrilinosScalar>::value>::type
1403  vmult(VectorType &dst, const VectorType &src) const;
1404 
1411  template <typename VectorType>
1412  typename std::enable_if<!std::is_same<typename VectorType::value_type,
1413  TrilinosScalar>::value>::type
1414  vmult(VectorType &dst, const VectorType &src) const;
1415 
1430  template <typename VectorType>
1431  typename std::enable_if<std::is_same<typename VectorType::value_type,
1432  TrilinosScalar>::value>::type
1433  Tvmult(VectorType &dst, const VectorType &src) const;
1434 
1441  template <typename VectorType>
1442  typename std::enable_if<!std::is_same<typename VectorType::value_type,
1443  TrilinosScalar>::value>::type
1444  Tvmult(VectorType &dst, const VectorType &src) const;
1445 
1455  template <typename VectorType>
1456  void
1457  vmult_add(VectorType &dst, const VectorType &src) const;
1458 
1469  template <typename VectorType>
1470  void
1471  Tvmult_add(VectorType &dst, const VectorType &src) const;
1472 
1495  matrix_norm_square(const MPI::Vector &v) const;
1496 
1517  matrix_scalar_product(const MPI::Vector &u, const MPI::Vector &v) const;
1518 
1536  residual(MPI::Vector & dst,
1537  const MPI::Vector &x,
1538  const MPI::Vector &b) const;
1539 
1554  void
1555  mmult(SparseMatrix & C,
1556  const SparseMatrix &B,
1557  const MPI::Vector & V = MPI::Vector()) const;
1558 
1559 
1576  void
1577  Tmmult(SparseMatrix & C,
1578  const SparseMatrix &B,
1579  const MPI::Vector & V = MPI::Vector()) const;
1580 
1582 
1586 
1595  l1_norm() const;
1596 
1606  linfty_norm() const;
1607 
1613  frobenius_norm() const;
1614 
1616 
1620 
1625  const Epetra_CrsMatrix &
1626  trilinos_matrix() const;
1627 
1632  const Epetra_CrsGraph &
1633  trilinos_sparsity_pattern() const;
1634 
1636 
1641 
1646  IndexSet
1647  locally_owned_domain_indices() const;
1648 
1654  IndexSet
1655  locally_owned_range_indices() const;
1656 
1658 
1663 
1683  begin() const;
1684 
1688  iterator
1689  begin();
1690 
1696  end() const;
1697 
1701  iterator
1702  end();
1703 
1733  begin(const size_type r) const;
1734 
1738  iterator
1739  begin(const size_type r);
1740 
1751  end(const size_type r) const;
1752 
1756  iterator
1757  end(const size_type r);
1758 
1760 
1764 
1770  void
1771  write_ascii();
1772 
1780  void
1781  print(std::ostream &out,
1782  const bool write_extended_trilinos_info = false) const;
1783 
1785 
1793  int,
1794  << "An error with error number " << arg1
1795  << " occurred while calling a Trilinos function");
1796 
1800  DeclException2(ExcInvalidIndex,
1801  size_type,
1802  size_type,
1803  << "The entry with index <" << arg1 << ',' << arg2
1804  << "> does not exist.");
1805 
1809  DeclExceptionMsg(ExcSourceEqualsDestination,
1810  "You are attempting an operation on two matrices that "
1811  "are the same object, but the operation requires that the "
1812  "two objects are in fact different.");
1813 
1817  DeclException0(ExcMatrixNotCompressed);
1818 
1822  DeclException4(ExcAccessToNonLocalElement,
1823  size_type,
1824  size_type,
1825  size_type,
1826  size_type,
1827  << "You tried to access element (" << arg1 << "/" << arg2
1828  << ")"
1829  << " of a distributed matrix, but only rows in range ["
1830  << arg3 << "," << arg4
1831  << "] are stored locally and can be accessed.");
1832 
1836  DeclException2(ExcAccessToNonPresentElement,
1837  size_type,
1838  size_type,
1839  << "You tried to access element (" << arg1 << "/" << arg2
1840  << ")"
1841  << " of a sparse matrix, but it appears to not"
1842  << " exist in the Trilinos sparsity pattern.");
1844 
1845 
1846 
1847  protected:
1858  void
1859  prepare_add();
1860 
1866  void
1867  prepare_set();
1868 
1869 
1870 
1871  private:
1876  std::unique_ptr<Epetra_Map> column_space_map;
1877 
1883  std::unique_ptr<Epetra_FECrsMatrix> matrix;
1884 
1890  std::unique_ptr<Epetra_CrsMatrix> nonlocal_matrix;
1891 
1895  std::unique_ptr<Epetra_Export> nonlocal_matrix_exporter;
1896 
1908  Epetra_CombineMode last_action;
1909 
1915 
1916  // To allow calling protected prepare_add() and prepare_set().
1918  };
1919 
1920 
1921 
1922  // forwards declarations
1923  class SolverBase;
1924  class PreconditionBase;
1925 
1926  namespace internal
1927  {
1928  inline void
1929  check_vector_map_equality(const Epetra_CrsMatrix & mtrx,
1930  const Epetra_MultiVector &src,
1931  const Epetra_MultiVector &dst,
1932  const bool transpose)
1933  {
1934  if (transpose == false)
1935  {
1936  Assert(src.Map().SameAs(mtrx.DomainMap()) == true,
1937  ExcMessage(
1938  "Column map of matrix does not fit with vector map!"));
1939  Assert(dst.Map().SameAs(mtrx.RangeMap()) == true,
1940  ExcMessage("Row map of matrix does not fit with vector map!"));
1941  }
1942  else
1943  {
1944  Assert(src.Map().SameAs(mtrx.RangeMap()) == true,
1945  ExcMessage(
1946  "Column map of matrix does not fit with vector map!"));
1947  Assert(dst.Map().SameAs(mtrx.DomainMap()) == true,
1948  ExcMessage("Row map of matrix does not fit with vector map!"));
1949  }
1950  (void)mtrx; // removes -Wunused-variable in optimized mode
1951  (void)src;
1952  (void)dst;
1953  }
1954 
1955  inline void
1957  const Epetra_MultiVector &src,
1958  const Epetra_MultiVector &dst,
1959  const bool transpose)
1960  {
1961  if (transpose == false)
1962  {
1963  Assert(src.Map().SameAs(op.OperatorDomainMap()) == true,
1964  ExcMessage(
1965  "Column map of operator does not fit with vector map!"));
1966  Assert(dst.Map().SameAs(op.OperatorRangeMap()) == true,
1967  ExcMessage(
1968  "Row map of operator does not fit with vector map!"));
1969  }
1970  else
1971  {
1972  Assert(src.Map().SameAs(op.OperatorRangeMap()) == true,
1973  ExcMessage(
1974  "Column map of operator does not fit with vector map!"));
1975  Assert(dst.Map().SameAs(op.OperatorDomainMap()) == true,
1976  ExcMessage(
1977  "Row map of operator does not fit with vector map!"));
1978  }
1979  (void)op; // removes -Wunused-variable in optimized mode
1980  (void)src;
1981  (void)dst;
1982  }
1983 
1984  namespace LinearOperatorImplementation
1985  {
2006  {
2007  public:
2011  using VectorType = Epetra_MultiVector;
2012 
2017 
2022 
2027 
2035  TrilinosPayload();
2036 
2040  TrilinosPayload(const TrilinosWrappers::SparseMatrix &matrix_exemplar,
2042 
2047  const TrilinosWrappers::SparseMatrix & matrix_exemplar,
2048  const TrilinosWrappers::PreconditionBase &preconditioner);
2049 
2054  const TrilinosWrappers::PreconditionBase &preconditioner_exemplar,
2055  const TrilinosWrappers::PreconditionBase &preconditioner);
2056 
2060  TrilinosPayload(const TrilinosPayload &payload);
2061 
2069  TrilinosPayload(const TrilinosPayload &first_op,
2070  const TrilinosPayload &second_op);
2071 
2075  virtual ~TrilinosPayload() override = default;
2076 
2081  identity_payload() const;
2082 
2087  null_payload() const;
2088 
2093  transpose_payload() const;
2094 
2111  template <typename Solver, typename Preconditioner>
2112  typename std::enable_if<
2113  std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
2114  std::is_base_of<TrilinosWrappers::PreconditionBase,
2115  Preconditioner>::value,
2116  TrilinosPayload>::type
2117  inverse_payload(Solver &, const Preconditioner &) const;
2118 
2136  template <typename Solver, typename Preconditioner>
2137  typename std::enable_if<
2138  !(std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
2139  std::is_base_of<TrilinosWrappers::PreconditionBase,
2140  Preconditioner>::value),
2141  TrilinosPayload>::type
2142  inverse_payload(Solver &, const Preconditioner &) const;
2143 
2145 
2150 
2156  IndexSet
2157  locally_owned_domain_indices() const;
2158 
2164  IndexSet
2165  locally_owned_range_indices() const;
2166 
2170  MPI_Comm
2171  get_mpi_communicator() const;
2172 
2179  void
2180  transpose();
2181 
2189  std::function<void(VectorType &, const VectorType &)> vmult;
2190 
2198  std::function<void(VectorType &, const VectorType &)> Tvmult;
2199 
2208  std::function<void(VectorType &, const VectorType &)> inv_vmult;
2209 
2218  std::function<void(VectorType &, const VectorType &)> inv_Tvmult;
2219 
2221 
2226 
2233  virtual bool
2234  UseTranspose() const override;
2235 
2251  virtual int
2252  SetUseTranspose(bool UseTranspose) override;
2253 
2265  virtual int
2266  Apply(const VectorType &X, VectorType &Y) const override;
2267 
2286  virtual int
2287  ApplyInverse(const VectorType &Y, VectorType &X) const override;
2289 
2294 
2301  virtual const char *
2302  Label() const override;
2303 
2311  virtual const Epetra_Comm &
2312  Comm() const override;
2313 
2321  virtual const Epetra_Map &
2322  OperatorDomainMap() const override;
2323 
2332  virtual const Epetra_Map &
2333  OperatorRangeMap() const override;
2335 
2336  private:
2342 
2347  Epetra_MpiComm communicator;
2348 
2353  Epetra_Map domain_map;
2354 
2359  Epetra_Map range_map;
2360 
2369  virtual bool
2370  HasNormInf() const override;
2371 
2379  virtual double
2380  NormInf() const override;
2381  };
2382 
2388  operator+(const TrilinosPayload &first_op,
2389  const TrilinosPayload &second_op);
2390 
2395  TrilinosPayload operator*(const TrilinosPayload &first_op,
2396  const TrilinosPayload &second_op);
2397 
2398  } // namespace LinearOperatorImplementation
2399  } /* namespace internal */
2400 
2401 
2402 
2403  // ----------------------- inline and template functions --------------------
2404 
2405 # ifndef DOXYGEN
2406 
2407  namespace SparseMatrixIterators
2408  {
2409  inline AccessorBase::AccessorBase(SparseMatrix *matrix,
2410  size_type row,
2411  size_type index)
2412  : matrix(matrix)
2413  , a_row(row)
2414  , a_index(index)
2415  {
2416  visit_present_row();
2417  }
2418 
2419 
2421  AccessorBase::row() const
2422  {
2423  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2424  return a_row;
2425  }
2426 
2427 
2429  AccessorBase::column() const
2430  {
2431  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2432  return (*colnum_cache)[a_index];
2433  }
2434 
2435 
2437  AccessorBase::index() const
2438  {
2439  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2440  return a_index;
2441  }
2442 
2443 
2444  inline Accessor<true>::Accessor(MatrixType * matrix,
2445  const size_type row,
2446  const size_type index)
2447  : AccessorBase(const_cast<SparseMatrix *>(matrix), row, index)
2448  {}
2449 
2450 
2451  template <bool Other>
2452  inline Accessor<true>::Accessor(const Accessor<Other> &other)
2453  : AccessorBase(other)
2454  {}
2455 
2456 
2457  inline TrilinosScalar
2458  Accessor<true>::value() const
2459  {
2460  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2461  return (*value_cache)[a_index];
2462  }
2463 
2464 
2465  inline Accessor<false>::Reference::Reference(const Accessor<false> &acc)
2466  : accessor(const_cast<Accessor<false> &>(acc))
2467  {}
2468 
2469 
2470  inline Accessor<false>::Reference::operator TrilinosScalar() const
2471  {
2472  return (*accessor.value_cache)[accessor.a_index];
2473  }
2474 
2475  inline const Accessor<false>::Reference &
2476  Accessor<false>::Reference::operator=(const TrilinosScalar n) const
2477  {
2478  (*accessor.value_cache)[accessor.a_index] = n;
2479  accessor.matrix->set(accessor.row(),
2480  accessor.column(),
2481  static_cast<TrilinosScalar>(*this));
2482  return *this;
2483  }
2484 
2485 
2486  inline const Accessor<false>::Reference &
2487  Accessor<false>::Reference::operator+=(const TrilinosScalar n) const
2488  {
2489  (*accessor.value_cache)[accessor.a_index] += n;
2490  accessor.matrix->set(accessor.row(),
2491  accessor.column(),
2492  static_cast<TrilinosScalar>(*this));
2493  return *this;
2494  }
2495 
2496 
2497  inline const Accessor<false>::Reference &
2498  Accessor<false>::Reference::operator-=(const TrilinosScalar n) const
2499  {
2500  (*accessor.value_cache)[accessor.a_index] -= n;
2501  accessor.matrix->set(accessor.row(),
2502  accessor.column(),
2503  static_cast<TrilinosScalar>(*this));
2504  return *this;
2505  }
2506 
2507 
2508  inline const Accessor<false>::Reference &
2509  Accessor<false>::Reference::operator*=(const TrilinosScalar n) const
2510  {
2511  (*accessor.value_cache)[accessor.a_index] *= n;
2512  accessor.matrix->set(accessor.row(),
2513  accessor.column(),
2514  static_cast<TrilinosScalar>(*this));
2515  return *this;
2516  }
2517 
2518 
2519  inline const Accessor<false>::Reference &
2520  Accessor<false>::Reference::operator/=(const TrilinosScalar n) const
2521  {
2522  (*accessor.value_cache)[accessor.a_index] /= n;
2523  accessor.matrix->set(accessor.row(),
2524  accessor.column(),
2525  static_cast<TrilinosScalar>(*this));
2526  return *this;
2527  }
2528 
2529 
2530  inline Accessor<false>::Accessor(MatrixType * matrix,
2531  const size_type row,
2532  const size_type index)
2533  : AccessorBase(matrix, row, index)
2534  {}
2535 
2536 
2537  inline Accessor<false>::Reference
2538  Accessor<false>::value() const
2539  {
2540  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2541  return {*this};
2542  }
2543 
2544 
2545 
2546  template <bool Constness>
2547  inline Iterator<Constness>::Iterator(MatrixType * matrix,
2548  const size_type row,
2549  const size_type index)
2550  : accessor(matrix, row, index)
2551  {}
2552 
2553 
2554  template <bool Constness>
2555  template <bool Other>
2556  inline Iterator<Constness>::Iterator(const Iterator<Other> &other)
2557  : accessor(other.accessor)
2558  {}
2559 
2560 
2561  template <bool Constness>
2562  inline Iterator<Constness> &
2564  {
2565  Assert(accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
2566 
2567  ++accessor.a_index;
2568 
2569  // If at end of line: do one
2570  // step, then cycle until we
2571  // find a row with a nonzero
2572  // number of entries.
2573  if (accessor.a_index >= accessor.colnum_cache->size())
2574  {
2575  accessor.a_index = 0;
2576  ++accessor.a_row;
2577 
2578  while ((accessor.a_row < accessor.matrix->m()) &&
2579  ((accessor.matrix->in_local_range(accessor.a_row) == false) ||
2580  (accessor.matrix->row_length(accessor.a_row) == 0)))
2581  ++accessor.a_row;
2582 
2583  accessor.visit_present_row();
2584  }
2585  return *this;
2586  }
2587 
2588 
2589  template <bool Constness>
2590  inline Iterator<Constness>
2592  {
2593  const Iterator<Constness> old_state = *this;
2594  ++(*this);
2595  return old_state;
2596  }
2597 
2598 
2599 
2600  template <bool Constness>
2601  inline const Accessor<Constness> &Iterator<Constness>::operator*() const
2602  {
2603  return accessor;
2604  }
2605 
2606 
2607 
2608  template <bool Constness>
2609  inline const Accessor<Constness> *Iterator<Constness>::operator->() const
2610  {
2611  return &accessor;
2612  }
2613 
2614 
2615 
2616  template <bool Constness>
2617  inline bool
2618  Iterator<Constness>::operator==(const Iterator<Constness> &other) const
2619  {
2620  return (accessor.a_row == other.accessor.a_row &&
2621  accessor.a_index == other.accessor.a_index);
2622  }
2623 
2624 
2625 
2626  template <bool Constness>
2627  inline bool
2628  Iterator<Constness>::operator!=(const Iterator<Constness> &other) const
2629  {
2630  return !(*this == other);
2631  }
2632 
2633 
2634 
2635  template <bool Constness>
2636  inline bool
2637  Iterator<Constness>::operator<(const Iterator<Constness> &other) const
2638  {
2639  return (accessor.row() < other.accessor.row() ||
2640  (accessor.row() == other.accessor.row() &&
2641  accessor.index() < other.accessor.index()));
2642  }
2643 
2644 
2645  template <bool Constness>
2646  inline bool
2647  Iterator<Constness>::operator>(const Iterator<Constness> &other) const
2648  {
2649  return (other < *this);
2650  }
2651 
2652  } // namespace SparseMatrixIterators
2653 
2654 
2655 
2657  SparseMatrix::begin() const
2658  {
2659  return begin(0);
2660  }
2661 
2662 
2663 
2665  SparseMatrix::end() const
2666  {
2667  return const_iterator(this, m(), 0);
2668  }
2669 
2670 
2671 
2673  SparseMatrix::begin(const size_type r) const
2674  {
2675  AssertIndexRange(r, m());
2676  if (in_local_range(r) && (row_length(r) > 0))
2677  return const_iterator(this, r, 0);
2678  else
2679  return end(r);
2680  }
2681 
2682 
2683 
2685  SparseMatrix::end(const size_type r) const
2686  {
2687  AssertIndexRange(r, m());
2688 
2689  // place the iterator on the first entry
2690  // past this line, or at the end of the
2691  // matrix
2692  for (size_type i = r + 1; i < m(); ++i)
2693  if (in_local_range(i) && (row_length(i) > 0))
2694  return const_iterator(this, i, 0);
2695 
2696  // if there is no such line, then take the
2697  // end iterator of the matrix
2698  return end();
2699  }
2700 
2701 
2702 
2703  inline SparseMatrix::iterator
2705  {
2706  return begin(0);
2707  }
2708 
2709 
2710 
2711  inline SparseMatrix::iterator
2713  {
2714  return iterator(this, m(), 0);
2715  }
2716 
2717 
2718 
2719  inline SparseMatrix::iterator
2721  {
2722  AssertIndexRange(r, m());
2723  if (in_local_range(r) && (row_length(r) > 0))
2724  return iterator(this, r, 0);
2725  else
2726  return end(r);
2727  }
2728 
2729 
2730 
2731  inline SparseMatrix::iterator
2732  SparseMatrix::end(const size_type r)
2733  {
2734  AssertIndexRange(r, m());
2735 
2736  // place the iterator on the first entry
2737  // past this line, or at the end of the
2738  // matrix
2739  for (size_type i = r + 1; i < m(); ++i)
2740  if (in_local_range(i) && (row_length(i) > 0))
2741  return iterator(this, i, 0);
2742 
2743  // if there is no such line, then take the
2744  // end iterator of the matrix
2745  return end();
2746  }
2747 
2748 
2749 
2750  inline bool
2751  SparseMatrix::in_local_range(const size_type index) const
2752  {
2754 # ifndef DEAL_II_WITH_64BIT_INDICES
2755  begin = matrix->RowMap().MinMyGID();
2756  end = matrix->RowMap().MaxMyGID() + 1;
2757 # else
2758  begin = matrix->RowMap().MinMyGID64();
2759  end = matrix->RowMap().MaxMyGID64() + 1;
2760 # endif
2761 
2762  return ((index >= static_cast<size_type>(begin)) &&
2763  (index < static_cast<size_type>(end)));
2764  }
2765 
2766 
2767 
2768  inline bool
2770  {
2771  return compressed;
2772  }
2773 
2774 
2775 
2776  // Inline the set() and add() functions, since they will be called
2777  // frequently, and the compiler can optimize away some unnecessary loops
2778  // when the sizes are given at compile time.
2779  template <>
2780  void
2781  SparseMatrix::set<TrilinosScalar>(const size_type row,
2782  const size_type n_cols,
2783  const size_type * col_indices,
2784  const TrilinosScalar *values,
2785  const bool elide_zero_values);
2786 
2787 
2788 
2789  template <typename Number>
2790  void
2791  SparseMatrix::set(const size_type row,
2792  const size_type n_cols,
2793  const size_type *col_indices,
2794  const Number * values,
2795  const bool elide_zero_values)
2796  {
2797  std::vector<TrilinosScalar> trilinos_values(n_cols);
2798  std::copy(values, values + n_cols, trilinos_values.begin());
2799  this->set(
2800  row, n_cols, col_indices, trilinos_values.data(), elide_zero_values);
2801  }
2802 
2803 
2804 
2805  inline void
2806  SparseMatrix::set(const size_type i,
2807  const size_type j,
2808  const TrilinosScalar value)
2809  {
2810  AssertIsFinite(value);
2811 
2812  set(i, 1, &j, &value, false);
2813  }
2814 
2815 
2816 
2817  inline void
2818  SparseMatrix::set(const std::vector<size_type> & indices,
2819  const FullMatrix<TrilinosScalar> &values,
2820  const bool elide_zero_values)
2821  {
2822  Assert(indices.size() == values.m(),
2823  ExcDimensionMismatch(indices.size(), values.m()));
2824  Assert(values.m() == values.n(), ExcNotQuadratic());
2825 
2826  for (size_type i = 0; i < indices.size(); ++i)
2827  set(indices[i],
2828  indices.size(),
2829  indices.data(),
2830  &values(i, 0),
2831  elide_zero_values);
2832  }
2833 
2834 
2835 
2836  inline void
2837  SparseMatrix::add(const size_type i,
2838  const size_type j,
2839  const TrilinosScalar value)
2840  {
2841  AssertIsFinite(value);
2842 
2843  if (value == 0)
2844  {
2845  // we have to check after Insert/Add in any case to be consistent
2846  // with the MPI communication model, but we can save some
2847  // work if the addend is zero. However, these actions are done in case
2848  // we pass on to the other function.
2849 
2850  // TODO: fix this (do not run compress here, but fail)
2851  if (last_action == Insert)
2852  {
2853  int ierr;
2854  ierr = matrix->GlobalAssemble(*column_space_map,
2855  matrix->RowMap(),
2856  false);
2857 
2858  Assert(ierr == 0, ExcTrilinosError(ierr));
2859  (void)ierr; // removes -Wunused-but-set-variable in optimized mode
2860  }
2861 
2862  last_action = Add;
2863 
2864  return;
2865  }
2866  else
2867  add(i, 1, &j, &value, false);
2868  }
2869 
2870 
2871 
2872  // inline "simple" functions that are called frequently and do only involve
2873  // a call to some Trilinos function.
2875  SparseMatrix::m() const
2876  {
2877 # ifndef DEAL_II_WITH_64BIT_INDICES
2878  return matrix->NumGlobalRows();
2879 # else
2880  return matrix->NumGlobalRows64();
2881 # endif
2882  }
2883 
2884 
2885 
2887  SparseMatrix::n() const
2888  {
2889  // If the matrix structure has not been fixed (i.e., we did not have a
2890  // sparsity pattern), it does not know about the number of columns so we
2891  // must always take this from the additional column space map
2892  Assert(column_space_map.get() != nullptr, ExcInternalError());
2893  return n_global_elements(*column_space_map);
2894  }
2895 
2896 
2897 
2898  inline unsigned int
2899  SparseMatrix::local_size() const
2900  {
2901  return matrix->NumMyRows();
2902  }
2903 
2904 
2905 
2906  inline std::pair<SparseMatrix::size_type, SparseMatrix::size_type>
2908  {
2909  size_type begin, end;
2910 # ifndef DEAL_II_WITH_64BIT_INDICES
2911  begin = matrix->RowMap().MinMyGID();
2912  end = matrix->RowMap().MaxMyGID() + 1;
2913 # else
2914  begin = matrix->RowMap().MinMyGID64();
2915  end = matrix->RowMap().MaxMyGID64() + 1;
2916 # endif
2917 
2918  return std::make_pair(begin, end);
2919  }
2920 
2921 
2922 
2925  {
2926 # ifndef DEAL_II_WITH_64BIT_INDICES
2927  return matrix->NumGlobalNonzeros();
2928 # else
2929  return matrix->NumGlobalNonzeros64();
2930 # endif
2931  }
2932 
2933 
2934 
2935  template <typename SparsityPatternType>
2936  inline void
2937  SparseMatrix::reinit(const IndexSet & parallel_partitioning,
2938  const SparsityPatternType &sparsity_pattern,
2939  const MPI_Comm & communicator,
2940  const bool exchange_data)
2941  {
2942  reinit(parallel_partitioning,
2943  parallel_partitioning,
2944  sparsity_pattern,
2945  communicator,
2946  exchange_data);
2947  }
2948 
2949 
2950 
2951  template <typename number>
2952  inline void
2953  SparseMatrix::reinit(const IndexSet &parallel_partitioning,
2954  const ::SparseMatrix<number> &sparse_matrix,
2955  const MPI_Comm & communicator,
2956  const double drop_tolerance,
2957  const bool copy_values,
2958  const ::SparsityPattern * use_this_sparsity)
2959  {
2960  Epetra_Map map =
2961  parallel_partitioning.make_trilinos_map(communicator, false);
2962  reinit(parallel_partitioning,
2963  parallel_partitioning,
2964  sparse_matrix,
2965  drop_tolerance,
2966  copy_values,
2967  use_this_sparsity);
2968  }
2969 
2970 
2971 
2972  inline const Epetra_CrsMatrix &
2974  {
2975  return static_cast<const Epetra_CrsMatrix &>(*matrix);
2976  }
2977 
2978 
2979 
2980  inline const Epetra_CrsGraph &
2982  {
2983  return matrix->Graph();
2984  }
2985 
2986 
2987 
2988  inline IndexSet
2990  {
2991  return IndexSet(matrix->DomainMap());
2992  }
2993 
2994 
2995 
2996  inline IndexSet
2998  {
2999  return IndexSet(matrix->RangeMap());
3000  }
3001 
3002 
3003 
3004  inline void
3006  {
3007  // nothing to do here
3008  }
3009 
3010 
3011 
3012  inline void
3014  {
3015  // nothing to do here
3016  }
3017 
3018 
3019  namespace internal
3020  {
3021  namespace LinearOperatorImplementation
3022  {
3023  template <typename Solver, typename Preconditioner>
3024  typename std::enable_if<
3025  std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
3026  std::is_base_of<TrilinosWrappers::PreconditionBase,
3027  Preconditioner>::value,
3028  TrilinosPayload>::type
3030  Solver & solver,
3031  const Preconditioner &preconditioner) const
3032  {
3033  const auto &payload = *this;
3034 
3035  TrilinosPayload return_op(payload);
3036 
3037  // Capture by copy so the payloads are always valid
3038 
3039  return_op.inv_vmult = [payload, &solver, &preconditioner](
3040  TrilinosPayload::Domain & tril_dst,
3041  const TrilinosPayload::Range &tril_src) {
3042  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3043  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3044  Assert(&tril_src != &tril_dst,
3047  tril_src,
3048  tril_dst,
3049  !payload.UseTranspose());
3050  solver.solve(payload, tril_dst, tril_src, preconditioner);
3051  };
3052 
3053  return_op.inv_Tvmult = [payload, &solver, &preconditioner](
3054  TrilinosPayload::Range & tril_dst,
3055  const TrilinosPayload::Domain &tril_src) {
3056  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3057  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3058  Assert(&tril_src != &tril_dst,
3061  tril_src,
3062  tril_dst,
3063  payload.UseTranspose());
3064 
3065  const_cast<TrilinosPayload &>(payload).transpose();
3066  solver.solve(payload, tril_dst, tril_src, preconditioner);
3067  const_cast<TrilinosPayload &>(payload).transpose();
3068  };
3069 
3070  // If the input operator is already setup for transpose operations, then
3071  // we must do similar with its inverse.
3072  if (return_op.UseTranspose() == true)
3073  std::swap(return_op.inv_vmult, return_op.inv_Tvmult);
3074 
3075  return return_op;
3076  }
3077 
3078  template <typename Solver, typename Preconditioner>
3079  typename std::enable_if<
3080  !(std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
3081  std::is_base_of<TrilinosWrappers::PreconditionBase,
3082  Preconditioner>::value),
3083  TrilinosPayload>::type
3084  TrilinosPayload::inverse_payload(Solver &, const Preconditioner &) const
3085  {
3086  TrilinosPayload return_op(*this);
3087 
3088  return_op.inv_vmult = [](TrilinosPayload::Domain &,
3089  const TrilinosPayload::Range &) {
3090  AssertThrow(false,
3091  ExcMessage("Payload inv_vmult disabled because of "
3092  "incompatible solver/preconditioner choice."));
3093  };
3094 
3095  return_op.inv_Tvmult = [](TrilinosPayload::Range &,
3096  const TrilinosPayload::Domain &) {
3097  AssertThrow(false,
3098  ExcMessage("Payload inv_vmult disabled because of "
3099  "incompatible solver/preconditioner choice."));
3100  };
3101 
3102  return return_op;
3103  }
3104  } // namespace LinearOperatorImplementation
3105  } // namespace internal
3106 
3107  template <>
3108  void
3109  SparseMatrix::set<TrilinosScalar>(const size_type row,
3110  const size_type n_cols,
3111  const size_type * col_indices,
3112  const TrilinosScalar *values,
3113  const bool elide_zero_values);
3114 # endif // DOXYGEN
3115 
3116 } /* namespace TrilinosWrappers */
3117 
3118 
3120 
3121 
3122 # endif // DEAL_II_WITH_TRILINOS
3123 
3124 
3125 /*----------------------- trilinos_sparse_matrix.h --------------------*/
3126 
3127 #endif
3128 /*----------------------- trilinos_sparse_matrix.h --------------------*/
static ::ExceptionBase & ExcTrilinosError(int arg1)
size_type m() const
const_iterator end() const
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:538
Contents is actually a matrix.
static ::ExceptionBase & ExcSourceEqualsDestination()
TrilinosPayload operator*(const TrilinosPayload &first_op, const TrilinosPayload &second_op)
types::global_dof_index size_type
Definition: cuda_kernels.h:45
static ::ExceptionBase & ExcBeyondEndOfMatrix()
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
std::function< void(VectorType &, const VectorType &)> inv_Tvmult
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1722
static const char V
TrilinosPayload operator+(const TrilinosPayload &first_op, const TrilinosPayload &second_op)
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)
static ::ExceptionBase & ExcAccessToNonlocalRow(std::size_t arg1, std::size_t arg2, std::size_t arg3)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:603
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)
std::string compress(const std::string &input)
Definition: utilities.cc:392
size_type n() const
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:515
void add(const size_type i, const size_type j, const TrilinosScalar value)
#define Assert(cond, exc)
Definition: exceptions.h:1465
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
std::unique_ptr< Epetra_Export > nonlocal_matrix_exporter
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:493
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2939
#define DeclException0(Exception0)
Definition: exceptions.h:470
typename Accessor< Constness >::MatrixType MatrixType
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
std::function< void(VectorType &, const VectorType &)> inv_vmult
VectorType::value_type * end(VectorType &V)
double TrilinosScalar
Definition: types.h:163
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
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
Expression operator>(const Expression &lhs, const Expression &rhs)
const Epetra_CrsMatrix & trilinos_matrix() const
AccessorBase(SparseMatrix *matrix, const size_type row, const size_type index)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SparseMatrix< double > SparseMatrix
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
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)
unsigned int global_dof_index
Definition: types.h:76
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
void reinit(const SparsityPatternType &sparsity_pattern)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
VectorType::value_type * begin(VectorType &V)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:584
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2913
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:561
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)
void check_vector_map_equality(const Epetra_Operator &op, const Epetra_MultiVector &src, const Epetra_MultiVector &dst, const bool transpose)
void copy(const T *begin, const T *end, U *dest)
#define AssertIsFinite(number)
Definition: exceptions.h:1748
unsigned int local_size() const
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
void check_vector_map_equality(const Epetra_CrsMatrix &mtrx, const Epetra_MultiVector &src, const Epetra_MultiVector &dst, const bool transpose)
static ::ExceptionBase & ExcInternalError()
std::shared_ptr< std::vector< TrilinosScalar > > value_cache