Reference documentation for deal.II version Git 50c3491829 2021-08-01 13:40:40 +0200
\(\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> &
391  operator*() const;
392 
396  const Accessor<Constness> *
397  operator->() const;
398 
403  bool
404  operator==(const Iterator<Constness> &) const;
405 
409  bool
410  operator!=(const Iterator<Constness> &) const;
411 
417  bool
418  operator<(const Iterator<Constness> &) const;
419 
423  bool
424  operator>(const Iterator<Constness> &) const;
425 
429  DeclException2(ExcInvalidIndexWithinRow,
430  size_type,
431  size_type,
432  << "Attempt to access element " << arg2 << " of row "
433  << arg1 << " which doesn't have that many elements.");
434 
435  private:
440 
441  template <bool Other>
442  friend class Iterator;
443  };
444 
445  } // namespace SparseMatrixIterators
446 
447 
508  class SparseMatrix : public Subscriptor
509  {
510  public:
515 
520  std::size_t,
521  << "You tried to access row " << arg1
522  << " of a non-contiguous locally owned row set."
523  << " The row " << arg1
524  << " is not stored locally and can't be accessed.");
525 
533  struct Traits
534  {
539  static const bool zero_addition_can_be_elided = true;
540  };
541 
546 
551 
556 
564  SparseMatrix();
565 
573  SparseMatrix(const size_type m,
574  const size_type n,
575  const unsigned int n_max_entries_per_row);
576 
584  SparseMatrix(const size_type m,
585  const size_type n,
586  const std::vector<unsigned int> &n_entries_per_row);
587 
591  SparseMatrix(const SparsityPattern &InputSparsityPattern);
592 
597  SparseMatrix(SparseMatrix &&other) noexcept;
598 
602  SparseMatrix(const SparseMatrix &) = delete;
603 
607  SparseMatrix &
608  operator=(const SparseMatrix &) = delete;
609 
613  virtual ~SparseMatrix() override = default;
614 
630  template <typename SparsityPatternType>
631  void
632  reinit(const SparsityPatternType &sparsity_pattern);
633 
646  void
647  reinit(const SparsityPattern &sparsity_pattern);
648 
657  void
658  reinit(const SparseMatrix &sparse_matrix);
659 
680  template <typename number>
681  void
682  reinit(const ::SparseMatrix<number> &dealii_sparse_matrix,
683  const double drop_tolerance = 1e-13,
684  const bool copy_values = true,
685  const ::SparsityPattern * use_this_sparsity = nullptr);
686 
692  void
693  reinit(const Epetra_CrsMatrix &input_matrix, const bool copy_values = true);
695 
712  SparseMatrix(const IndexSet & parallel_partitioning,
713  const MPI_Comm & communicator = MPI_COMM_WORLD,
714  const unsigned int n_max_entries_per_row = 0);
715 
723  SparseMatrix(const IndexSet & parallel_partitioning,
724  const MPI_Comm & communicator,
725  const std::vector<unsigned int> &n_entries_per_row);
726 
741  SparseMatrix(const IndexSet &row_parallel_partitioning,
742  const IndexSet &col_parallel_partitioning,
743  const MPI_Comm &communicator = MPI_COMM_WORLD,
744  const size_type n_max_entries_per_row = 0);
745 
760  SparseMatrix(const IndexSet & row_parallel_partitioning,
761  const IndexSet & col_parallel_partitioning,
762  const MPI_Comm & communicator,
763  const std::vector<unsigned int> &n_entries_per_row);
764 
785  template <typename SparsityPatternType>
786  void
787  reinit(const IndexSet & parallel_partitioning,
788  const SparsityPatternType &sparsity_pattern,
789  const MPI_Comm & communicator = MPI_COMM_WORLD,
790  const bool exchange_data = false);
791 
804  template <typename SparsityPatternType>
805  typename std::enable_if<
806  !std::is_same<SparsityPatternType,
807  ::SparseMatrix<double>>::value>::type
808  reinit(const IndexSet & row_parallel_partitioning,
809  const IndexSet & col_parallel_partitioning,
810  const SparsityPatternType &sparsity_pattern,
811  const MPI_Comm & communicator = MPI_COMM_WORLD,
812  const bool exchange_data = false);
813 
830  template <typename number>
831  void
832  reinit(const IndexSet & parallel_partitioning,
833  const ::SparseMatrix<number> &dealii_sparse_matrix,
834  const MPI_Comm & communicator = MPI_COMM_WORLD,
835  const double drop_tolerance = 1e-13,
836  const bool copy_values = true,
837  const ::SparsityPattern * use_this_sparsity = nullptr);
838 
852  template <typename number>
853  void
854  reinit(const IndexSet & row_parallel_partitioning,
855  const IndexSet & col_parallel_partitioning,
856  const ::SparseMatrix<number> &dealii_sparse_matrix,
857  const MPI_Comm & communicator = MPI_COMM_WORLD,
858  const double drop_tolerance = 1e-13,
859  const bool copy_values = true,
860  const ::SparsityPattern * use_this_sparsity = nullptr);
862 
866 
870  size_type
871  m() const;
872 
876  size_type
877  n() const;
878 
887  unsigned int
888  local_size() const;
889 
898  std::pair<size_type, size_type>
899  local_range() const;
900 
905  bool
906  in_local_range(const size_type index) const;
907 
912  size_type
913  n_nonzero_elements() const;
914 
918  unsigned int
919  row_length(const size_type row) const;
920 
927  bool
928  is_compressed() const;
929 
935  size_type
936  memory_consumption() const;
937 
941  MPI_Comm
942  get_mpi_communicator() const;
943 
945 
949 
959  SparseMatrix &
960  operator=(const double d);
961 
969  void
970  clear();
971 
999  void
1000  compress(::VectorOperation::values operation);
1001 
1023  void
1024  set(const size_type i, const size_type j, const TrilinosScalar value);
1025 
1058  void
1059  set(const std::vector<size_type> & indices,
1060  const FullMatrix<TrilinosScalar> &full_matrix,
1061  const bool elide_zero_values = false);
1062 
1068  void
1069  set(const std::vector<size_type> & row_indices,
1070  const std::vector<size_type> & col_indices,
1071  const FullMatrix<TrilinosScalar> &full_matrix,
1072  const bool elide_zero_values = false);
1073 
1101  void
1102  set(const size_type row,
1103  const std::vector<size_type> & col_indices,
1104  const std::vector<TrilinosScalar> &values,
1105  const bool elide_zero_values = false);
1106 
1134  template <typename Number>
1135  void
1136  set(const size_type row,
1137  const size_type n_cols,
1138  const size_type *col_indices,
1139  const Number * values,
1140  const bool elide_zero_values = false);
1141 
1151  void
1152  add(const size_type i, const size_type j, const TrilinosScalar value);
1153 
1172  void
1173  add(const std::vector<size_type> & indices,
1174  const FullMatrix<TrilinosScalar> &full_matrix,
1175  const bool elide_zero_values = true);
1176 
1182  void
1183  add(const std::vector<size_type> & row_indices,
1184  const std::vector<size_type> & col_indices,
1185  const FullMatrix<TrilinosScalar> &full_matrix,
1186  const bool elide_zero_values = true);
1187 
1201  void
1202  add(const size_type row,
1203  const std::vector<size_type> & col_indices,
1204  const std::vector<TrilinosScalar> &values,
1205  const bool elide_zero_values = true);
1206 
1220  void
1221  add(const size_type row,
1222  const size_type n_cols,
1223  const size_type * col_indices,
1224  const TrilinosScalar *values,
1225  const bool elide_zero_values = true,
1226  const bool col_indices_are_sorted = false);
1227 
1231  SparseMatrix &
1232  operator*=(const TrilinosScalar factor);
1233 
1237  SparseMatrix &
1238  operator/=(const TrilinosScalar factor);
1239 
1243  void
1244  copy_from(const SparseMatrix &source);
1245 
1253  void
1254  add(const TrilinosScalar factor, const SparseMatrix &matrix);
1255 
1282  void
1283  clear_row(const size_type row, const TrilinosScalar new_diag_value = 0);
1284 
1305  void
1306  clear_rows(const std::vector<size_type> &rows,
1307  const TrilinosScalar new_diag_value = 0);
1308 
1318  void
1319  transpose();
1320 
1322 
1326 
1336  operator()(const size_type i, const size_type j) const;
1337 
1355  el(const size_type i, const size_type j) const;
1356 
1364  diag_element(const size_type i) const;
1365 
1367 
1371 
1402  template <typename VectorType>
1403  typename std::enable_if<std::is_same<typename VectorType::value_type,
1404  TrilinosScalar>::value>::type
1405  vmult(VectorType &dst, const VectorType &src) const;
1406 
1413  template <typename VectorType>
1414  typename std::enable_if<!std::is_same<typename VectorType::value_type,
1415  TrilinosScalar>::value>::type
1416  vmult(VectorType &dst, const VectorType &src) const;
1417 
1432  template <typename VectorType>
1433  typename std::enable_if<std::is_same<typename VectorType::value_type,
1434  TrilinosScalar>::value>::type
1435  Tvmult(VectorType &dst, const VectorType &src) const;
1436 
1443  template <typename VectorType>
1444  typename std::enable_if<!std::is_same<typename VectorType::value_type,
1445  TrilinosScalar>::value>::type
1446  Tvmult(VectorType &dst, const VectorType &src) const;
1447 
1457  template <typename VectorType>
1458  void
1459  vmult_add(VectorType &dst, const VectorType &src) const;
1460 
1471  template <typename VectorType>
1472  void
1473  Tvmult_add(VectorType &dst, const VectorType &src) const;
1474 
1497  matrix_norm_square(const MPI::Vector &v) const;
1498 
1519  matrix_scalar_product(const MPI::Vector &u, const MPI::Vector &v) const;
1520 
1538  residual(MPI::Vector & dst,
1539  const MPI::Vector &x,
1540  const MPI::Vector &b) const;
1541 
1556  void
1557  mmult(SparseMatrix & C,
1558  const SparseMatrix &B,
1559  const MPI::Vector & V = MPI::Vector()) const;
1560 
1561 
1578  void
1579  Tmmult(SparseMatrix & C,
1580  const SparseMatrix &B,
1581  const MPI::Vector & V = MPI::Vector()) const;
1582 
1584 
1588 
1597  l1_norm() const;
1598 
1608  linfty_norm() const;
1609 
1615  frobenius_norm() const;
1616 
1618 
1622 
1627  const Epetra_CrsMatrix &
1628  trilinos_matrix() const;
1629 
1634  const Epetra_CrsGraph &
1635  trilinos_sparsity_pattern() const;
1636 
1638 
1643 
1648  IndexSet
1649  locally_owned_domain_indices() const;
1650 
1656  IndexSet
1657  locally_owned_range_indices() const;
1658 
1660 
1665 
1685  begin() const;
1686 
1690  iterator
1691  begin();
1692 
1698  end() const;
1699 
1703  iterator
1704  end();
1705 
1735  begin(const size_type r) const;
1736 
1740  iterator
1741  begin(const size_type r);
1742 
1753  end(const size_type r) const;
1754 
1758  iterator
1759  end(const size_type r);
1760 
1762 
1766 
1772  void
1773  write_ascii();
1774 
1782  void
1783  print(std::ostream &out,
1784  const bool write_extended_trilinos_info = false) const;
1785 
1787 
1795  int,
1796  << "An error with error number " << arg1
1797  << " occurred while calling a Trilinos function");
1798 
1802  DeclException2(ExcInvalidIndex,
1803  size_type,
1804  size_type,
1805  << "The entry with index <" << arg1 << ',' << arg2
1806  << "> does not exist.");
1807 
1811  DeclExceptionMsg(ExcSourceEqualsDestination,
1812  "You are attempting an operation on two matrices that "
1813  "are the same object, but the operation requires that the "
1814  "two objects are in fact different.");
1815 
1819  DeclException0(ExcMatrixNotCompressed);
1820 
1824  DeclException4(ExcAccessToNonLocalElement,
1825  size_type,
1826  size_type,
1827  size_type,
1828  size_type,
1829  << "You tried to access element (" << arg1 << "/" << arg2
1830  << ")"
1831  << " of a distributed matrix, but only rows in range ["
1832  << arg3 << "," << arg4
1833  << "] are stored locally and can be accessed.");
1834 
1838  DeclException2(ExcAccessToNonPresentElement,
1839  size_type,
1840  size_type,
1841  << "You tried to access element (" << arg1 << "/" << arg2
1842  << ")"
1843  << " of a sparse matrix, but it appears to not"
1844  << " exist in the Trilinos sparsity pattern.");
1846 
1847 
1848 
1849  protected:
1860  void
1861  prepare_add();
1862 
1868  void
1869  prepare_set();
1870 
1871 
1872 
1873  private:
1878  std::unique_ptr<Epetra_Map> column_space_map;
1879 
1885  std::unique_ptr<Epetra_FECrsMatrix> matrix;
1886 
1892  std::unique_ptr<Epetra_CrsMatrix> nonlocal_matrix;
1893 
1897  std::unique_ptr<Epetra_Export> nonlocal_matrix_exporter;
1898 
1910  Epetra_CombineMode last_action;
1911 
1917 
1918  // To allow calling protected prepare_add() and prepare_set().
1920  };
1921 
1922 
1923 
1924  // forwards declarations
1925  class SolverBase;
1926  class PreconditionBase;
1927 
1928  namespace internal
1929  {
1930  inline void
1931  check_vector_map_equality(const Epetra_CrsMatrix & mtrx,
1932  const Epetra_MultiVector &src,
1933  const Epetra_MultiVector &dst,
1934  const bool transpose)
1935  {
1936  if (transpose == false)
1937  {
1938  Assert(src.Map().SameAs(mtrx.DomainMap()) == true,
1939  ExcMessage(
1940  "Column map of matrix does not fit with vector map!"));
1941  Assert(dst.Map().SameAs(mtrx.RangeMap()) == true,
1942  ExcMessage("Row map of matrix does not fit with vector map!"));
1943  }
1944  else
1945  {
1946  Assert(src.Map().SameAs(mtrx.RangeMap()) == true,
1947  ExcMessage(
1948  "Column map of matrix does not fit with vector map!"));
1949  Assert(dst.Map().SameAs(mtrx.DomainMap()) == true,
1950  ExcMessage("Row map of matrix does not fit with vector map!"));
1951  }
1952  (void)mtrx; // removes -Wunused-variable in optimized mode
1953  (void)src;
1954  (void)dst;
1955  }
1956 
1957  inline void
1959  const Epetra_MultiVector &src,
1960  const Epetra_MultiVector &dst,
1961  const bool transpose)
1962  {
1963  if (transpose == false)
1964  {
1965  Assert(src.Map().SameAs(op.OperatorDomainMap()) == true,
1966  ExcMessage(
1967  "Column map of operator does not fit with vector map!"));
1968  Assert(dst.Map().SameAs(op.OperatorRangeMap()) == true,
1969  ExcMessage(
1970  "Row map of operator does not fit with vector map!"));
1971  }
1972  else
1973  {
1974  Assert(src.Map().SameAs(op.OperatorRangeMap()) == true,
1975  ExcMessage(
1976  "Column map of operator does not fit with vector map!"));
1977  Assert(dst.Map().SameAs(op.OperatorDomainMap()) == true,
1978  ExcMessage(
1979  "Row map of operator does not fit with vector map!"));
1980  }
1981  (void)op; // removes -Wunused-variable in optimized mode
1982  (void)src;
1983  (void)dst;
1984  }
1985 
1986  namespace LinearOperatorImplementation
1987  {
2008  {
2009  public:
2013  using VectorType = Epetra_MultiVector;
2014 
2019 
2024 
2029 
2037  TrilinosPayload();
2038 
2042  TrilinosPayload(const TrilinosWrappers::SparseMatrix &matrix_exemplar,
2044 
2049  const TrilinosWrappers::SparseMatrix & matrix_exemplar,
2050  const TrilinosWrappers::PreconditionBase &preconditioner);
2051 
2056  const TrilinosWrappers::PreconditionBase &preconditioner_exemplar,
2057  const TrilinosWrappers::PreconditionBase &preconditioner);
2058 
2062  TrilinosPayload(const TrilinosPayload &payload);
2063 
2071  TrilinosPayload(const TrilinosPayload &first_op,
2072  const TrilinosPayload &second_op);
2073 
2077  virtual ~TrilinosPayload() override = default;
2078 
2083  identity_payload() const;
2084 
2089  null_payload() const;
2090 
2095  transpose_payload() const;
2096 
2113  template <typename Solver, typename Preconditioner>
2114  typename std::enable_if<
2115  std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
2116  std::is_base_of<TrilinosWrappers::PreconditionBase,
2117  Preconditioner>::value,
2118  TrilinosPayload>::type
2119  inverse_payload(Solver &, const Preconditioner &) const;
2120 
2138  template <typename Solver, typename Preconditioner>
2139  typename std::enable_if<
2140  !(std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
2141  std::is_base_of<TrilinosWrappers::PreconditionBase,
2142  Preconditioner>::value),
2143  TrilinosPayload>::type
2144  inverse_payload(Solver &, const Preconditioner &) const;
2145 
2147 
2152 
2158  IndexSet
2159  locally_owned_domain_indices() const;
2160 
2166  IndexSet
2167  locally_owned_range_indices() const;
2168 
2172  MPI_Comm
2173  get_mpi_communicator() const;
2174 
2181  void
2182  transpose();
2183 
2191  std::function<void(VectorType &, const VectorType &)> vmult;
2192 
2200  std::function<void(VectorType &, const VectorType &)> Tvmult;
2201 
2210  std::function<void(VectorType &, const VectorType &)> inv_vmult;
2211 
2220  std::function<void(VectorType &, const VectorType &)> inv_Tvmult;
2221 
2223 
2228 
2235  virtual bool
2236  UseTranspose() const override;
2237 
2253  virtual int
2254  SetUseTranspose(bool UseTranspose) override;
2255 
2267  virtual int
2268  Apply(const VectorType &X, VectorType &Y) const override;
2269 
2288  virtual int
2289  ApplyInverse(const VectorType &Y, VectorType &X) const override;
2291 
2296 
2303  virtual const char *
2304  Label() const override;
2305 
2313  virtual const Epetra_Comm &
2314  Comm() const override;
2315 
2323  virtual const Epetra_Map &
2324  OperatorDomainMap() const override;
2325 
2334  virtual const Epetra_Map &
2335  OperatorRangeMap() const override;
2337 
2338  private:
2344 
2349  Epetra_MpiComm communicator;
2350 
2355  Epetra_Map domain_map;
2356 
2361  Epetra_Map range_map;
2362 
2371  virtual bool
2372  HasNormInf() const override;
2373 
2381  virtual double
2382  NormInf() const override;
2383  };
2384 
2390  operator+(const TrilinosPayload &first_op,
2391  const TrilinosPayload &second_op);
2392 
2398  operator*(const TrilinosPayload &first_op,
2399  const TrilinosPayload &second_op);
2400 
2401  } // namespace LinearOperatorImplementation
2402  } /* namespace internal */
2403 
2404 
2405 
2406  // ----------------------- inline and template functions --------------------
2407 
2408 # ifndef DOXYGEN
2409 
2410  namespace SparseMatrixIterators
2411  {
2412  inline AccessorBase::AccessorBase(SparseMatrix *matrix,
2413  size_type row,
2414  size_type index)
2415  : matrix(matrix)
2416  , a_row(row)
2417  , a_index(index)
2418  {
2419  visit_present_row();
2420  }
2421 
2422 
2424  AccessorBase::row() const
2425  {
2426  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2427  return a_row;
2428  }
2429 
2430 
2432  AccessorBase::column() const
2433  {
2434  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2435  return (*colnum_cache)[a_index];
2436  }
2437 
2438 
2440  AccessorBase::index() const
2441  {
2442  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2443  return a_index;
2444  }
2445 
2446 
2447  inline Accessor<true>::Accessor(MatrixType * matrix,
2448  const size_type row,
2449  const size_type index)
2450  : AccessorBase(const_cast<SparseMatrix *>(matrix), row, index)
2451  {}
2452 
2453 
2454  template <bool Other>
2455  inline Accessor<true>::Accessor(const Accessor<Other> &other)
2456  : AccessorBase(other)
2457  {}
2458 
2459 
2460  inline TrilinosScalar
2461  Accessor<true>::value() const
2462  {
2463  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2464  return (*value_cache)[a_index];
2465  }
2466 
2467 
2468  inline Accessor<false>::Reference::Reference(const Accessor<false> &acc)
2469  : accessor(const_cast<Accessor<false> &>(acc))
2470  {}
2471 
2472 
2473  inline Accessor<false>::Reference::operator TrilinosScalar() const
2474  {
2475  return (*accessor.value_cache)[accessor.a_index];
2476  }
2477 
2478  inline const Accessor<false>::Reference &
2479  Accessor<false>::Reference::operator=(const TrilinosScalar n) const
2480  {
2481  (*accessor.value_cache)[accessor.a_index] = n;
2482  accessor.matrix->set(accessor.row(),
2483  accessor.column(),
2484  static_cast<TrilinosScalar>(*this));
2485  return *this;
2486  }
2487 
2488 
2489  inline const Accessor<false>::Reference &
2490  Accessor<false>::Reference::operator+=(const TrilinosScalar n) const
2491  {
2492  (*accessor.value_cache)[accessor.a_index] += n;
2493  accessor.matrix->set(accessor.row(),
2494  accessor.column(),
2495  static_cast<TrilinosScalar>(*this));
2496  return *this;
2497  }
2498 
2499 
2500  inline const Accessor<false>::Reference &
2501  Accessor<false>::Reference::operator-=(const TrilinosScalar n) const
2502  {
2503  (*accessor.value_cache)[accessor.a_index] -= n;
2504  accessor.matrix->set(accessor.row(),
2505  accessor.column(),
2506  static_cast<TrilinosScalar>(*this));
2507  return *this;
2508  }
2509 
2510 
2511  inline const Accessor<false>::Reference &
2512  Accessor<false>::Reference::operator*=(const TrilinosScalar n) const
2513  {
2514  (*accessor.value_cache)[accessor.a_index] *= n;
2515  accessor.matrix->set(accessor.row(),
2516  accessor.column(),
2517  static_cast<TrilinosScalar>(*this));
2518  return *this;
2519  }
2520 
2521 
2522  inline const Accessor<false>::Reference &
2523  Accessor<false>::Reference::operator/=(const TrilinosScalar n) const
2524  {
2525  (*accessor.value_cache)[accessor.a_index] /= n;
2526  accessor.matrix->set(accessor.row(),
2527  accessor.column(),
2528  static_cast<TrilinosScalar>(*this));
2529  return *this;
2530  }
2531 
2532 
2533  inline Accessor<false>::Accessor(MatrixType * matrix,
2534  const size_type row,
2535  const size_type index)
2536  : AccessorBase(matrix, row, index)
2537  {}
2538 
2539 
2540  inline Accessor<false>::Reference
2541  Accessor<false>::value() const
2542  {
2543  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2544  return {*this};
2545  }
2546 
2547 
2548 
2549  template <bool Constness>
2550  inline Iterator<Constness>::Iterator(MatrixType * matrix,
2551  const size_type row,
2552  const size_type index)
2553  : accessor(matrix, row, index)
2554  {}
2555 
2556 
2557  template <bool Constness>
2558  template <bool Other>
2559  inline Iterator<Constness>::Iterator(const Iterator<Other> &other)
2560  : accessor(other.accessor)
2561  {}
2562 
2563 
2564  template <bool Constness>
2565  inline Iterator<Constness> &
2567  {
2568  Assert(accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
2569 
2570  ++accessor.a_index;
2571 
2572  // If at end of line: do one
2573  // step, then cycle until we
2574  // find a row with a nonzero
2575  // number of entries.
2576  if (accessor.a_index >= accessor.colnum_cache->size())
2577  {
2578  accessor.a_index = 0;
2579  ++accessor.a_row;
2580 
2581  while ((accessor.a_row < accessor.matrix->m()) &&
2582  ((accessor.matrix->in_local_range(accessor.a_row) == false) ||
2583  (accessor.matrix->row_length(accessor.a_row) == 0)))
2584  ++accessor.a_row;
2585 
2586  accessor.visit_present_row();
2587  }
2588  return *this;
2589  }
2590 
2591 
2592  template <bool Constness>
2593  inline Iterator<Constness>
2595  {
2596  const Iterator<Constness> old_state = *this;
2597  ++(*this);
2598  return old_state;
2599  }
2600 
2601 
2602 
2603  template <bool Constness>
2604  inline const Accessor<Constness> &
2606  {
2607  return accessor;
2608  }
2609 
2610 
2611 
2612  template <bool Constness>
2613  inline const Accessor<Constness> *
2614  Iterator<Constness>::operator->() const
2615  {
2616  return &accessor;
2617  }
2618 
2619 
2620 
2621  template <bool Constness>
2622  inline bool
2623  Iterator<Constness>::operator==(const Iterator<Constness> &other) const
2624  {
2625  return (accessor.a_row == other.accessor.a_row &&
2626  accessor.a_index == other.accessor.a_index);
2627  }
2628 
2629 
2630 
2631  template <bool Constness>
2632  inline bool
2633  Iterator<Constness>::operator!=(const Iterator<Constness> &other) const
2634  {
2635  return !(*this == other);
2636  }
2637 
2638 
2639 
2640  template <bool Constness>
2641  inline bool
2642  Iterator<Constness>::operator<(const Iterator<Constness> &other) const
2643  {
2644  return (accessor.row() < other.accessor.row() ||
2645  (accessor.row() == other.accessor.row() &&
2646  accessor.index() < other.accessor.index()));
2647  }
2648 
2649 
2650  template <bool Constness>
2651  inline bool
2652  Iterator<Constness>::operator>(const Iterator<Constness> &other) const
2653  {
2654  return (other < *this);
2655  }
2656 
2657  } // namespace SparseMatrixIterators
2658 
2659 
2660 
2662  SparseMatrix::begin() const
2663  {
2664  return begin(0);
2665  }
2666 
2667 
2668 
2670  SparseMatrix::end() const
2671  {
2672  return const_iterator(this, m(), 0);
2673  }
2674 
2675 
2676 
2678  SparseMatrix::begin(const size_type r) const
2679  {
2680  AssertIndexRange(r, m());
2681  if (in_local_range(r) && (row_length(r) > 0))
2682  return const_iterator(this, r, 0);
2683  else
2684  return end(r);
2685  }
2686 
2687 
2688 
2690  SparseMatrix::end(const size_type r) const
2691  {
2692  AssertIndexRange(r, m());
2693 
2694  // place the iterator on the first entry
2695  // past this line, or at the end of the
2696  // matrix
2697  for (size_type i = r + 1; i < m(); ++i)
2698  if (in_local_range(i) && (row_length(i) > 0))
2699  return const_iterator(this, i, 0);
2700 
2701  // if there is no such line, then take the
2702  // end iterator of the matrix
2703  return end();
2704  }
2705 
2706 
2707 
2708  inline SparseMatrix::iterator
2710  {
2711  return begin(0);
2712  }
2713 
2714 
2715 
2716  inline SparseMatrix::iterator
2718  {
2719  return iterator(this, m(), 0);
2720  }
2721 
2722 
2723 
2724  inline SparseMatrix::iterator
2726  {
2727  AssertIndexRange(r, m());
2728  if (in_local_range(r) && (row_length(r) > 0))
2729  return iterator(this, r, 0);
2730  else
2731  return end(r);
2732  }
2733 
2734 
2735 
2736  inline SparseMatrix::iterator
2737  SparseMatrix::end(const size_type r)
2738  {
2739  AssertIndexRange(r, m());
2740 
2741  // place the iterator on the first entry
2742  // past this line, or at the end of the
2743  // matrix
2744  for (size_type i = r + 1; i < m(); ++i)
2745  if (in_local_range(i) && (row_length(i) > 0))
2746  return iterator(this, i, 0);
2747 
2748  // if there is no such line, then take the
2749  // end iterator of the matrix
2750  return end();
2751  }
2752 
2753 
2754 
2755  inline bool
2756  SparseMatrix::in_local_range(const size_type index) const
2757  {
2759 # ifndef DEAL_II_WITH_64BIT_INDICES
2760  begin = matrix->RowMap().MinMyGID();
2761  end = matrix->RowMap().MaxMyGID() + 1;
2762 # else
2763  begin = matrix->RowMap().MinMyGID64();
2764  end = matrix->RowMap().MaxMyGID64() + 1;
2765 # endif
2766 
2767  return ((index >= static_cast<size_type>(begin)) &&
2768  (index < static_cast<size_type>(end)));
2769  }
2770 
2771 
2772 
2773  inline bool
2775  {
2776  return compressed;
2777  }
2778 
2779 
2780 
2781  // Inline the set() and add() functions, since they will be called
2782  // frequently, and the compiler can optimize away some unnecessary loops
2783  // when the sizes are given at compile time.
2784  template <>
2785  void
2786  SparseMatrix::set<TrilinosScalar>(const size_type row,
2787  const size_type n_cols,
2788  const size_type * col_indices,
2789  const TrilinosScalar *values,
2790  const bool elide_zero_values);
2791 
2792 
2793 
2794  template <typename Number>
2795  void
2796  SparseMatrix::set(const size_type row,
2797  const size_type n_cols,
2798  const size_type *col_indices,
2799  const Number * values,
2800  const bool elide_zero_values)
2801  {
2802  std::vector<TrilinosScalar> trilinos_values(n_cols);
2803  std::copy(values, values + n_cols, trilinos_values.begin());
2804  this->set(
2805  row, n_cols, col_indices, trilinos_values.data(), elide_zero_values);
2806  }
2807 
2808 
2809 
2810  inline void
2811  SparseMatrix::set(const size_type i,
2812  const size_type j,
2813  const TrilinosScalar value)
2814  {
2815  AssertIsFinite(value);
2816 
2817  set(i, 1, &j, &value, false);
2818  }
2819 
2820 
2821 
2822  inline void
2823  SparseMatrix::set(const std::vector<size_type> & indices,
2824  const FullMatrix<TrilinosScalar> &values,
2825  const bool elide_zero_values)
2826  {
2827  Assert(indices.size() == values.m(),
2828  ExcDimensionMismatch(indices.size(), values.m()));
2829  Assert(values.m() == values.n(), ExcNotQuadratic());
2830 
2831  for (size_type i = 0; i < indices.size(); ++i)
2832  set(indices[i],
2833  indices.size(),
2834  indices.data(),
2835  &values(i, 0),
2836  elide_zero_values);
2837  }
2838 
2839 
2840 
2841  inline void
2842  SparseMatrix::add(const size_type i,
2843  const size_type j,
2844  const TrilinosScalar value)
2845  {
2846  AssertIsFinite(value);
2847 
2848  if (value == 0)
2849  {
2850  // we have to check after Insert/Add in any case to be consistent
2851  // with the MPI communication model, but we can save some
2852  // work if the addend is zero. However, these actions are done in case
2853  // we pass on to the other function.
2854 
2855  // TODO: fix this (do not run compress here, but fail)
2856  if (last_action == Insert)
2857  {
2858  int ierr;
2859  ierr = matrix->GlobalAssemble(*column_space_map,
2860  matrix->RowMap(),
2861  false);
2862 
2863  Assert(ierr == 0, ExcTrilinosError(ierr));
2864  (void)ierr; // removes -Wunused-but-set-variable in optimized mode
2865  }
2866 
2867  last_action = Add;
2868 
2869  return;
2870  }
2871  else
2872  add(i, 1, &j, &value, false);
2873  }
2874 
2875 
2876 
2877  // inline "simple" functions that are called frequently and do only involve
2878  // a call to some Trilinos function.
2880  SparseMatrix::m() const
2881  {
2882 # ifndef DEAL_II_WITH_64BIT_INDICES
2883  return matrix->NumGlobalRows();
2884 # else
2885  return matrix->NumGlobalRows64();
2886 # endif
2887  }
2888 
2889 
2890 
2892  SparseMatrix::n() const
2893  {
2894  // If the matrix structure has not been fixed (i.e., we did not have a
2895  // sparsity pattern), it does not know about the number of columns so we
2896  // must always take this from the additional column space map
2897  Assert(column_space_map.get() != nullptr, ExcInternalError());
2898  return n_global_elements(*column_space_map);
2899  }
2900 
2901 
2902 
2903  inline unsigned int
2904  SparseMatrix::local_size() const
2905  {
2906  return matrix->NumMyRows();
2907  }
2908 
2909 
2910 
2911  inline std::pair<SparseMatrix::size_type, SparseMatrix::size_type>
2913  {
2914  size_type begin, end;
2915 # ifndef DEAL_II_WITH_64BIT_INDICES
2916  begin = matrix->RowMap().MinMyGID();
2917  end = matrix->RowMap().MaxMyGID() + 1;
2918 # else
2919  begin = matrix->RowMap().MinMyGID64();
2920  end = matrix->RowMap().MaxMyGID64() + 1;
2921 # endif
2922 
2923  return std::make_pair(begin, end);
2924  }
2925 
2926 
2927 
2930  {
2931 # ifndef DEAL_II_WITH_64BIT_INDICES
2932  return matrix->NumGlobalNonzeros();
2933 # else
2934  return matrix->NumGlobalNonzeros64();
2935 # endif
2936  }
2937 
2938 
2939 
2940  template <typename SparsityPatternType>
2941  inline void
2942  SparseMatrix::reinit(const IndexSet & parallel_partitioning,
2943  const SparsityPatternType &sparsity_pattern,
2944  const MPI_Comm & communicator,
2945  const bool exchange_data)
2946  {
2947  reinit(parallel_partitioning,
2948  parallel_partitioning,
2949  sparsity_pattern,
2950  communicator,
2951  exchange_data);
2952  }
2953 
2954 
2955 
2956  template <typename number>
2957  inline void
2958  SparseMatrix::reinit(const IndexSet &parallel_partitioning,
2959  const ::SparseMatrix<number> &sparse_matrix,
2960  const MPI_Comm & communicator,
2961  const double drop_tolerance,
2962  const bool copy_values,
2963  const ::SparsityPattern * use_this_sparsity)
2964  {
2965  Epetra_Map map =
2966  parallel_partitioning.make_trilinos_map(communicator, false);
2967  reinit(parallel_partitioning,
2968  parallel_partitioning,
2969  sparse_matrix,
2970  drop_tolerance,
2971  copy_values,
2972  use_this_sparsity);
2973  }
2974 
2975 
2976 
2977  inline const Epetra_CrsMatrix &
2979  {
2980  return static_cast<const Epetra_CrsMatrix &>(*matrix);
2981  }
2982 
2983 
2984 
2985  inline const Epetra_CrsGraph &
2987  {
2988  return matrix->Graph();
2989  }
2990 
2991 
2992 
2993  inline IndexSet
2995  {
2996  return IndexSet(matrix->DomainMap());
2997  }
2998 
2999 
3000 
3001  inline IndexSet
3003  {
3004  return IndexSet(matrix->RangeMap());
3005  }
3006 
3007 
3008 
3009  inline void
3011  {
3012  // nothing to do here
3013  }
3014 
3015 
3016 
3017  inline void
3019  {
3020  // nothing to do here
3021  }
3022 
3023 
3024  namespace internal
3025  {
3026  namespace LinearOperatorImplementation
3027  {
3028  template <typename Solver, typename Preconditioner>
3029  typename std::enable_if<
3030  std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
3031  std::is_base_of<TrilinosWrappers::PreconditionBase,
3032  Preconditioner>::value,
3033  TrilinosPayload>::type
3035  Solver & solver,
3036  const Preconditioner &preconditioner) const
3037  {
3038  const auto &payload = *this;
3039 
3040  TrilinosPayload return_op(payload);
3041 
3042  // Capture by copy so the payloads are always valid
3043 
3044  return_op.inv_vmult = [payload, &solver, &preconditioner](
3045  TrilinosPayload::Domain & tril_dst,
3046  const TrilinosPayload::Range &tril_src) {
3047  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3048  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3049  Assert(&tril_src != &tril_dst,
3052  tril_src,
3053  tril_dst,
3054  !payload.UseTranspose());
3055  solver.solve(payload, tril_dst, tril_src, preconditioner);
3056  };
3057 
3058  return_op.inv_Tvmult = [payload, &solver, &preconditioner](
3059  TrilinosPayload::Range & tril_dst,
3060  const TrilinosPayload::Domain &tril_src) {
3061  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3062  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3063  Assert(&tril_src != &tril_dst,
3066  tril_src,
3067  tril_dst,
3068  payload.UseTranspose());
3069 
3070  const_cast<TrilinosPayload &>(payload).transpose();
3071  solver.solve(payload, tril_dst, tril_src, preconditioner);
3072  const_cast<TrilinosPayload &>(payload).transpose();
3073  };
3074 
3075  // If the input operator is already setup for transpose operations, then
3076  // we must do similar with its inverse.
3077  if (return_op.UseTranspose() == true)
3078  std::swap(return_op.inv_vmult, return_op.inv_Tvmult);
3079 
3080  return return_op;
3081  }
3082 
3083  template <typename Solver, typename Preconditioner>
3084  typename std::enable_if<
3085  !(std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
3086  std::is_base_of<TrilinosWrappers::PreconditionBase,
3087  Preconditioner>::value),
3088  TrilinosPayload>::type
3089  TrilinosPayload::inverse_payload(Solver &, const Preconditioner &) const
3090  {
3091  TrilinosPayload return_op(*this);
3092 
3093  return_op.inv_vmult = [](TrilinosPayload::Domain &,
3094  const TrilinosPayload::Range &) {
3095  AssertThrow(false,
3096  ExcMessage("Payload inv_vmult disabled because of "
3097  "incompatible solver/preconditioner choice."));
3098  };
3099 
3100  return_op.inv_Tvmult = [](TrilinosPayload::Range &,
3101  const TrilinosPayload::Domain &) {
3102  AssertThrow(false,
3103  ExcMessage("Payload inv_vmult disabled because of "
3104  "incompatible solver/preconditioner choice."));
3105  };
3106 
3107  return return_op;
3108  }
3109  } // namespace LinearOperatorImplementation
3110  } // namespace internal
3111 
3112  template <>
3113  void
3114  SparseMatrix::set<TrilinosScalar>(const size_type row,
3115  const size_type n_cols,
3116  const size_type * col_indices,
3117  const TrilinosScalar *values,
3118  const bool elide_zero_values);
3119 # endif // DOXYGEN
3120 
3121 } /* namespace TrilinosWrappers */
3122 
3123 
3125 
3126 
3127 # endif // DEAL_II_WITH_TRILINOS
3128 
3129 
3130 /*----------------------- trilinos_sparse_matrix.h --------------------*/
3131 
3132 #endif
3133 /*----------------------- 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:1724
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:1577
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:604
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:1467
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:2944
#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:2918
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:1750
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