Reference documentation for deal.II version Git f0919993dd 2020-09-21 18:25:06 -0600
\(\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 - 2020 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_MultiVector.h>
42 # include <Epetra_Operator.h>
43 
44 # include <cmath>
45 # include <memory>
46 # include <type_traits>
47 # include <vector>
48 # ifdef DEAL_II_WITH_MPI
49 # include <Epetra_MpiComm.h>
50 # include <mpi.h>
51 # else
52 # include <Epetra_SerialComm.h>
53 # endif
54 
56 
57 // forward declarations
58 # ifndef DOXYGEN
59 template <typename MatrixType>
60 class BlockMatrixBase;
61 
62 template <typename number>
63 class SparseMatrix;
64 class SparsityPattern;
66 
67 namespace TrilinosWrappers
68 {
69  class SparseMatrix;
70  class SparsityPattern;
71 
72  namespace SparseMatrixIterators
73  {
74  template <bool Constness>
75  class Iterator;
76  }
77 } // namespace TrilinosWrappers
78 # endif
79 
80 namespace TrilinosWrappers
81 {
86  {
91 
96  std::size_t,
97  std::size_t,
98  std::size_t,
99  << "You tried to access row " << arg1
100  << " of a distributed sparsity pattern, "
101  << " but only rows " << arg2 << " through " << arg3
102  << " are stored locally and can be accessed.");
103 
112  {
113  public:
118 
123  const size_type row,
124  const size_type index);
125 
129  size_type
130  row() const;
131 
135  size_type
136  index() const;
137 
141  size_type
142  column() const;
143 
144  protected:
156 
161 
167  void
168  visit_present_row();
169 
182  std::shared_ptr<std::vector<size_type>> colnum_cache;
183 
187  std::shared_ptr<std::vector<TrilinosScalar>> value_cache;
188  };
189 
200  template <bool Constess>
201  class Accessor : public AccessorBase
202  {
207  value() const;
208 
213  value();
214  };
215 
219  template <>
220  class Accessor<true> : public AccessorBase
221  {
222  public:
227  using MatrixType = const SparseMatrix;
228 
233  Accessor(MatrixType *matrix, const size_type row, const size_type index);
234 
239  template <bool Other>
240  Accessor(const Accessor<Other> &a);
241 
246  value() const;
247 
248  private:
249  // Make iterator class a friend.
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:
309  };
310 
311  public:
317 
322  Accessor(MatrixType *matrix, const size_type row, const size_type index);
323 
327  Reference
328  value() const;
329 
330  private:
331  // Make iterator class a friend.
332  template <bool>
333  friend class Iterator;
334 
335  // Make Reference object a friend.
336  friend class Reference;
337  };
338 
352  template <bool Constness>
353  class Iterator
354  {
355  public:
360 
366 
371  Iterator(MatrixType *matrix, const size_type row, const size_type index);
372 
376  template <bool Other>
377  Iterator(const Iterator<Other> &other);
378 
383  operator++();
384 
389  operator++(int);
390 
394  const Accessor<Constness> &operator*() const;
395 
399  const Accessor<Constness> *operator->() const;
400 
405  bool
406  operator==(const Iterator<Constness> &) const;
407 
411  bool
412  operator!=(const Iterator<Constness> &) const;
413 
419  bool
420  operator<(const Iterator<Constness> &) const;
421 
425  bool
426  operator>(const Iterator<Constness> &) const;
427 
431  DeclException2(ExcInvalidIndexWithinRow,
432  size_type,
433  size_type,
434  << "Attempt to access element " << arg2 << " of row "
435  << arg1 << " which doesn't have that many elements.");
436 
437  private:
442 
443  template <bool Other>
444  friend class Iterator;
445  };
446 
447  } // namespace SparseMatrixIterators
448 
449 
510  class SparseMatrix : public Subscriptor
511  {
512  public:
517 
522  std::size_t,
523  << "You tried to access row " << arg1
524  << " of a non-contiguous locally owned row set."
525  << " The row " << arg1
526  << " is not stored locally and can't be accessed.");
527 
535  struct Traits
536  {
541  static const bool zero_addition_can_be_elided = true;
542  };
543 
548 
553 
558 
566  SparseMatrix();
567 
575  SparseMatrix(const size_type m,
576  const size_type n,
577  const unsigned int n_max_entries_per_row);
578 
586  SparseMatrix(const size_type m,
587  const size_type n,
588  const std::vector<unsigned int> &n_entries_per_row);
589 
593  SparseMatrix(const SparsityPattern &InputSparsityPattern);
594 
599  SparseMatrix(SparseMatrix &&other) noexcept;
600 
604  SparseMatrix(const SparseMatrix &) = delete;
605 
609  SparseMatrix &
610  operator=(const SparseMatrix &) = delete;
611 
615  virtual ~SparseMatrix() override = default;
616 
632  template <typename SparsityPatternType>
633  void
634  reinit(const SparsityPatternType &sparsity_pattern);
635 
648  void
649  reinit(const SparsityPattern &sparsity_pattern);
650 
659  void
660  reinit(const SparseMatrix &sparse_matrix);
661 
682  template <typename number>
683  void
684  reinit(const ::SparseMatrix<number> &dealii_sparse_matrix,
685  const double drop_tolerance = 1e-13,
686  const bool copy_values = true,
687  const ::SparsityPattern * use_this_sparsity = nullptr);
688 
694  void
695  reinit(const Epetra_CrsMatrix &input_matrix, const bool copy_values = true);
697 
714  SparseMatrix(const IndexSet & parallel_partitioning,
715  const MPI_Comm & communicator = MPI_COMM_WORLD,
716  const unsigned int n_max_entries_per_row = 0);
717 
725  SparseMatrix(const IndexSet & parallel_partitioning,
726  const MPI_Comm & communicator,
727  const std::vector<unsigned int> &n_entries_per_row);
728 
743  SparseMatrix(const IndexSet &row_parallel_partitioning,
744  const IndexSet &col_parallel_partitioning,
745  const MPI_Comm &communicator = MPI_COMM_WORLD,
746  const size_type n_max_entries_per_row = 0);
747 
762  SparseMatrix(const IndexSet & row_parallel_partitioning,
763  const IndexSet & col_parallel_partitioning,
764  const MPI_Comm & communicator,
765  const std::vector<unsigned int> &n_entries_per_row);
766 
787  template <typename SparsityPatternType>
788  void
789  reinit(const IndexSet & parallel_partitioning,
790  const SparsityPatternType &sparsity_pattern,
791  const MPI_Comm & communicator = MPI_COMM_WORLD,
792  const bool exchange_data = false);
793 
806  template <typename SparsityPatternType>
807  typename std::enable_if<
808  !std::is_same<SparsityPatternType,
809  ::SparseMatrix<double>>::value>::type
810  reinit(const IndexSet & row_parallel_partitioning,
811  const IndexSet & col_parallel_partitioning,
812  const SparsityPatternType &sparsity_pattern,
813  const MPI_Comm & communicator = MPI_COMM_WORLD,
814  const bool exchange_data = false);
815 
832  template <typename number>
833  void
834  reinit(const IndexSet & parallel_partitioning,
835  const ::SparseMatrix<number> &dealii_sparse_matrix,
836  const MPI_Comm & communicator = MPI_COMM_WORLD,
837  const double drop_tolerance = 1e-13,
838  const bool copy_values = true,
839  const ::SparsityPattern * use_this_sparsity = nullptr);
840 
854  template <typename number>
855  void
856  reinit(const IndexSet & row_parallel_partitioning,
857  const IndexSet & col_parallel_partitioning,
858  const ::SparseMatrix<number> &dealii_sparse_matrix,
859  const MPI_Comm & communicator = MPI_COMM_WORLD,
860  const double drop_tolerance = 1e-13,
861  const bool copy_values = true,
862  const ::SparsityPattern * use_this_sparsity = nullptr);
864 
868 
872  size_type
873  m() const;
874 
878  size_type
879  n() const;
880 
889  unsigned int
890  local_size() const;
891 
900  std::pair<size_type, size_type>
901  local_range() const;
902 
907  bool
908  in_local_range(const size_type index) const;
909 
914  size_type
915  n_nonzero_elements() const;
916 
920  unsigned int
921  row_length(const size_type row) const;
922 
929  bool
930  is_compressed() const;
931 
937  size_type
938  memory_consumption() const;
939 
943  MPI_Comm
944  get_mpi_communicator() const;
945 
947 
951 
961  SparseMatrix &
962  operator=(const double d);
963 
971  void
972  clear();
973 
1001  void
1002  compress(::VectorOperation::values operation);
1003 
1025  void
1026  set(const size_type i, const size_type j, const TrilinosScalar value);
1027 
1060  void
1061  set(const std::vector<size_type> & indices,
1062  const FullMatrix<TrilinosScalar> &full_matrix,
1063  const bool elide_zero_values = false);
1064 
1070  void
1071  set(const std::vector<size_type> & row_indices,
1072  const std::vector<size_type> & col_indices,
1073  const FullMatrix<TrilinosScalar> &full_matrix,
1074  const bool elide_zero_values = false);
1075 
1103  void
1104  set(const size_type row,
1105  const std::vector<size_type> & col_indices,
1106  const std::vector<TrilinosScalar> &values,
1107  const bool elide_zero_values = false);
1108 
1136  template <typename Number>
1137  void
1138  set(const size_type row,
1139  const size_type n_cols,
1140  const size_type *col_indices,
1141  const Number * values,
1142  const bool elide_zero_values = false);
1143 
1153  void
1154  add(const size_type i, const size_type j, const TrilinosScalar value);
1155 
1174  void
1175  add(const std::vector<size_type> & indices,
1176  const FullMatrix<TrilinosScalar> &full_matrix,
1177  const bool elide_zero_values = true);
1178 
1184  void
1185  add(const std::vector<size_type> & row_indices,
1186  const std::vector<size_type> & col_indices,
1187  const FullMatrix<TrilinosScalar> &full_matrix,
1188  const bool elide_zero_values = true);
1189 
1203  void
1204  add(const size_type row,
1205  const std::vector<size_type> & col_indices,
1206  const std::vector<TrilinosScalar> &values,
1207  const bool elide_zero_values = true);
1208 
1222  void
1223  add(const size_type row,
1224  const size_type n_cols,
1225  const size_type * col_indices,
1226  const TrilinosScalar *values,
1227  const bool elide_zero_values = true,
1228  const bool col_indices_are_sorted = false);
1229 
1233  SparseMatrix &
1234  operator*=(const TrilinosScalar factor);
1235 
1239  SparseMatrix &
1240  operator/=(const TrilinosScalar factor);
1241 
1245  void
1246  copy_from(const SparseMatrix &source);
1247 
1255  void
1256  add(const TrilinosScalar factor, const SparseMatrix &matrix);
1257 
1284  void
1285  clear_row(const size_type row, const TrilinosScalar new_diag_value = 0);
1286 
1307  void
1308  clear_rows(const std::vector<size_type> &rows,
1309  const TrilinosScalar new_diag_value = 0);
1310 
1320  void
1321  transpose();
1322 
1324 
1328 
1338  operator()(const size_type i, const size_type j) const;
1339 
1357  el(const size_type i, const size_type j) const;
1358 
1366  diag_element(const size_type i) const;
1367 
1369 
1373 
1404  template <typename VectorType>
1405  typename std::enable_if<std::is_same<typename VectorType::value_type,
1406  TrilinosScalar>::value>::type
1407  vmult(VectorType &dst, const VectorType &src) const;
1408 
1415  template <typename VectorType>
1416  typename std::enable_if<!std::is_same<typename VectorType::value_type,
1417  TrilinosScalar>::value>::type
1418  vmult(VectorType &dst, const VectorType &src) const;
1419 
1434  template <typename VectorType>
1435  typename std::enable_if<std::is_same<typename VectorType::value_type,
1436  TrilinosScalar>::value>::type
1437  Tvmult(VectorType &dst, const VectorType &src) const;
1438 
1445  template <typename VectorType>
1446  typename std::enable_if<!std::is_same<typename VectorType::value_type,
1447  TrilinosScalar>::value>::type
1448  Tvmult(VectorType &dst, const VectorType &src) const;
1449 
1459  template <typename VectorType>
1460  void
1461  vmult_add(VectorType &dst, const VectorType &src) const;
1462 
1473  template <typename VectorType>
1474  void
1475  Tvmult_add(VectorType &dst, const VectorType &src) const;
1476 
1499  matrix_norm_square(const MPI::Vector &v) const;
1500 
1521  matrix_scalar_product(const MPI::Vector &u, const MPI::Vector &v) const;
1522 
1540  residual(MPI::Vector & dst,
1541  const MPI::Vector &x,
1542  const MPI::Vector &b) const;
1543 
1558  void
1559  mmult(SparseMatrix & C,
1560  const SparseMatrix &B,
1561  const MPI::Vector & V = MPI::Vector()) const;
1562 
1563 
1580  void
1581  Tmmult(SparseMatrix & C,
1582  const SparseMatrix &B,
1583  const MPI::Vector & V = MPI::Vector()) const;
1584 
1586 
1590 
1599  l1_norm() const;
1600 
1610  linfty_norm() const;
1611 
1617  frobenius_norm() const;
1618 
1620 
1624 
1629  const Epetra_CrsMatrix &
1630  trilinos_matrix() const;
1631 
1636  const Epetra_CrsGraph &
1637  trilinos_sparsity_pattern() const;
1638 
1640 
1645 
1650  IndexSet
1651  locally_owned_domain_indices() const;
1652 
1658  IndexSet
1659  locally_owned_range_indices() const;
1660 
1662 
1667 
1687  begin() const;
1688 
1692  iterator
1693  begin();
1694 
1700  end() const;
1701 
1705  iterator
1706  end();
1707 
1737  begin(const size_type r) const;
1738 
1742  iterator
1743  begin(const size_type r);
1744 
1755  end(const size_type r) const;
1756 
1760  iterator
1761  end(const size_type r);
1762 
1764 
1768 
1774  void
1775  write_ascii();
1776 
1784  void
1785  print(std::ostream &out,
1786  const bool write_extended_trilinos_info = false) const;
1787 
1789 
1797  int,
1798  << "An error with error number " << arg1
1799  << " occurred while calling a Trilinos function");
1800 
1804  DeclException2(ExcInvalidIndex,
1805  size_type,
1806  size_type,
1807  << "The entry with index <" << arg1 << ',' << arg2
1808  << "> does not exist.");
1809 
1813  DeclExceptionMsg(ExcSourceEqualsDestination,
1814  "You are attempting an operation on two matrices that "
1815  "are the same object, but the operation requires that the "
1816  "two objects are in fact different.");
1817 
1821  DeclException0(ExcMatrixNotCompressed);
1822 
1826  DeclException4(ExcAccessToNonLocalElement,
1827  size_type,
1828  size_type,
1829  size_type,
1830  size_type,
1831  << "You tried to access element (" << arg1 << "/" << arg2
1832  << ")"
1833  << " of a distributed matrix, but only rows " << arg3
1834  << " through " << arg4
1835  << " are stored locally and can be accessed.");
1836 
1840  DeclException2(ExcAccessToNonPresentElement,
1841  size_type,
1842  size_type,
1843  << "You tried to access element (" << arg1 << "/" << arg2
1844  << ")"
1845  << " of a sparse matrix, but it appears to not"
1846  << " exist in the Trilinos sparsity pattern.");
1848 
1849 
1850 
1851  protected:
1862  void
1863  prepare_add();
1864 
1870  void
1871  prepare_set();
1872 
1873 
1874 
1875  private:
1880  std::unique_ptr<Epetra_Map> column_space_map;
1881 
1887  std::unique_ptr<Epetra_FECrsMatrix> matrix;
1888 
1894  std::unique_ptr<Epetra_CrsMatrix> nonlocal_matrix;
1895 
1899  std::unique_ptr<Epetra_Export> nonlocal_matrix_exporter;
1900 
1912  Epetra_CombineMode last_action;
1913 
1919 
1920  // To allow calling protected prepare_add() and prepare_set().
1922  };
1923 
1924 
1925 
1926  // forwards declarations
1927  class SolverBase;
1928  class PreconditionBase;
1929 
1930  namespace internal
1931  {
1932  inline void
1933  check_vector_map_equality(const Epetra_CrsMatrix & mtrx,
1934  const Epetra_MultiVector &src,
1935  const Epetra_MultiVector &dst,
1936  const bool transpose)
1937  {
1938  if (transpose == false)
1939  {
1940  Assert(src.Map().SameAs(mtrx.DomainMap()) == true,
1941  ExcMessage(
1942  "Column map of matrix does not fit with vector map!"));
1943  Assert(dst.Map().SameAs(mtrx.RangeMap()) == true,
1944  ExcMessage("Row map of matrix does not fit with vector map!"));
1945  }
1946  else
1947  {
1948  Assert(src.Map().SameAs(mtrx.RangeMap()) == true,
1949  ExcMessage(
1950  "Column map of matrix does not fit with vector map!"));
1951  Assert(dst.Map().SameAs(mtrx.DomainMap()) == true,
1952  ExcMessage("Row map of matrix does not fit with vector map!"));
1953  }
1954  (void)mtrx; // removes -Wunused-variable in optimized mode
1955  (void)src;
1956  (void)dst;
1957  }
1958 
1959  inline void
1961  const Epetra_MultiVector &src,
1962  const Epetra_MultiVector &dst,
1963  const bool transpose)
1964  {
1965  if (transpose == false)
1966  {
1967  Assert(src.Map().SameAs(op.OperatorDomainMap()) == true,
1968  ExcMessage(
1969  "Column map of operator does not fit with vector map!"));
1970  Assert(dst.Map().SameAs(op.OperatorRangeMap()) == true,
1971  ExcMessage(
1972  "Row map of operator does not fit with vector map!"));
1973  }
1974  else
1975  {
1976  Assert(src.Map().SameAs(op.OperatorRangeMap()) == true,
1977  ExcMessage(
1978  "Column map of operator does not fit with vector map!"));
1979  Assert(dst.Map().SameAs(op.OperatorDomainMap()) == true,
1980  ExcMessage(
1981  "Row map of operator does not fit with vector map!"));
1982  }
1983  (void)op; // removes -Wunused-variable in optimized mode
1984  (void)src;
1985  (void)dst;
1986  }
1987 
1988  namespace LinearOperatorImplementation
1989  {
2010  {
2011  public:
2015  using VectorType = Epetra_MultiVector;
2016 
2021 
2026 
2031 
2039  TrilinosPayload();
2040 
2044  TrilinosPayload(const TrilinosWrappers::SparseMatrix &matrix_exemplar,
2046 
2051  const TrilinosWrappers::SparseMatrix & matrix_exemplar,
2052  const TrilinosWrappers::PreconditionBase &preconditioner);
2053 
2058  const TrilinosWrappers::PreconditionBase &preconditioner_exemplar,
2059  const TrilinosWrappers::PreconditionBase &preconditioner);
2060 
2064  TrilinosPayload(const TrilinosPayload &payload);
2065 
2073  TrilinosPayload(const TrilinosPayload &first_op,
2074  const TrilinosPayload &second_op);
2075 
2079  virtual ~TrilinosPayload() override = default;
2080 
2085  identity_payload() const;
2086 
2091  null_payload() const;
2092 
2097  transpose_payload() const;
2098 
2115  template <typename Solver, typename Preconditioner>
2116  typename std::enable_if<
2117  std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
2118  std::is_base_of<TrilinosWrappers::PreconditionBase,
2119  Preconditioner>::value,
2120  TrilinosPayload>::type
2121  inverse_payload(Solver &, const Preconditioner &) const;
2122 
2140  template <typename Solver, typename Preconditioner>
2141  typename std::enable_if<
2142  !(std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
2143  std::is_base_of<TrilinosWrappers::PreconditionBase,
2144  Preconditioner>::value),
2145  TrilinosPayload>::type
2146  inverse_payload(Solver &, const Preconditioner &) const;
2147 
2149 
2154 
2160  IndexSet
2161  locally_owned_domain_indices() const;
2162 
2168  IndexSet
2169  locally_owned_range_indices() const;
2170 
2174  MPI_Comm
2175  get_mpi_communicator() const;
2176 
2183  void
2184  transpose();
2185 
2193  std::function<void(VectorType &, const VectorType &)> vmult;
2194 
2202  std::function<void(VectorType &, const VectorType &)> Tvmult;
2203 
2212  std::function<void(VectorType &, const VectorType &)> inv_vmult;
2213 
2222  std::function<void(VectorType &, const VectorType &)> inv_Tvmult;
2223 
2225 
2230 
2237  virtual bool
2238  UseTranspose() const override;
2239 
2255  virtual int
2256  SetUseTranspose(bool UseTranspose) override;
2257 
2269  virtual int
2270  Apply(const VectorType &X, VectorType &Y) const override;
2271 
2290  virtual int
2291  ApplyInverse(const VectorType &Y, VectorType &X) const override;
2293 
2298 
2305  virtual const char *
2306  Label() const override;
2307 
2315  virtual const Epetra_Comm &
2316  Comm() const override;
2317 
2325  virtual const Epetra_Map &
2326  OperatorDomainMap() const override;
2327 
2336  virtual const Epetra_Map &
2337  OperatorRangeMap() const override;
2339 
2340  private:
2346 
2351 # ifdef DEAL_II_WITH_MPI
2352  Epetra_MpiComm communicator;
2353 # else
2354  Epetra_SerialComm communicator;
2355 # endif
2356 
2361  Epetra_Map domain_map;
2362 
2367  Epetra_Map range_map;
2368 
2377  virtual bool
2378  HasNormInf() const override;
2379 
2387  virtual double
2388  NormInf() const override;
2389  };
2390 
2396  operator+(const TrilinosPayload &first_op,
2397  const TrilinosPayload &second_op);
2398 
2403  TrilinosPayload operator*(const TrilinosPayload &first_op,
2404  const TrilinosPayload &second_op);
2405 
2406  } // namespace LinearOperatorImplementation
2407  } /* namespace internal */
2408 
2409 
2410 
2411  // ----------------------- inline and template functions --------------------
2412 
2413 # ifndef DOXYGEN
2414 
2415  namespace SparseMatrixIterators
2416  {
2417  inline AccessorBase::AccessorBase(SparseMatrix *matrix,
2418  size_type row,
2419  size_type index)
2420  : matrix(matrix)
2421  , a_row(row)
2422  , a_index(index)
2423  {
2424  visit_present_row();
2425  }
2426 
2427 
2429  AccessorBase::row() const
2430  {
2431  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2432  return a_row;
2433  }
2434 
2435 
2437  AccessorBase::column() const
2438  {
2439  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2440  return (*colnum_cache)[a_index];
2441  }
2442 
2443 
2445  AccessorBase::index() const
2446  {
2447  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2448  return a_index;
2449  }
2450 
2451 
2452  inline Accessor<true>::Accessor(MatrixType * matrix,
2453  const size_type row,
2454  const size_type index)
2455  : AccessorBase(const_cast<SparseMatrix *>(matrix), row, index)
2456  {}
2457 
2458 
2459  template <bool Other>
2460  inline Accessor<true>::Accessor(const Accessor<Other> &other)
2461  : AccessorBase(other)
2462  {}
2463 
2464 
2465  inline TrilinosScalar
2466  Accessor<true>::value() const
2467  {
2468  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2469  return (*value_cache)[a_index];
2470  }
2471 
2472 
2473  inline Accessor<false>::Reference::Reference(const Accessor<false> &acc)
2474  : accessor(const_cast<Accessor<false> &>(acc))
2475  {}
2476 
2477 
2478  inline Accessor<false>::Reference::operator TrilinosScalar() const
2479  {
2480  return (*accessor.value_cache)[accessor.a_index];
2481  }
2482 
2483  inline const Accessor<false>::Reference &
2484  Accessor<false>::Reference::operator=(const TrilinosScalar n) const
2485  {
2486  (*accessor.value_cache)[accessor.a_index] = n;
2487  accessor.matrix->set(accessor.row(),
2488  accessor.column(),
2489  static_cast<TrilinosScalar>(*this));
2490  return *this;
2491  }
2492 
2493 
2494  inline const Accessor<false>::Reference &
2495  Accessor<false>::Reference::operator+=(const TrilinosScalar n) const
2496  {
2497  (*accessor.value_cache)[accessor.a_index] += n;
2498  accessor.matrix->set(accessor.row(),
2499  accessor.column(),
2500  static_cast<TrilinosScalar>(*this));
2501  return *this;
2502  }
2503 
2504 
2505  inline const Accessor<false>::Reference &
2506  Accessor<false>::Reference::operator-=(const TrilinosScalar n) const
2507  {
2508  (*accessor.value_cache)[accessor.a_index] -= n;
2509  accessor.matrix->set(accessor.row(),
2510  accessor.column(),
2511  static_cast<TrilinosScalar>(*this));
2512  return *this;
2513  }
2514 
2515 
2516  inline const Accessor<false>::Reference &
2517  Accessor<false>::Reference::operator*=(const TrilinosScalar n) const
2518  {
2519  (*accessor.value_cache)[accessor.a_index] *= n;
2520  accessor.matrix->set(accessor.row(),
2521  accessor.column(),
2522  static_cast<TrilinosScalar>(*this));
2523  return *this;
2524  }
2525 
2526 
2527  inline const Accessor<false>::Reference &
2528  Accessor<false>::Reference::operator/=(const TrilinosScalar n) const
2529  {
2530  (*accessor.value_cache)[accessor.a_index] /= n;
2531  accessor.matrix->set(accessor.row(),
2532  accessor.column(),
2533  static_cast<TrilinosScalar>(*this));
2534  return *this;
2535  }
2536 
2537 
2538  inline Accessor<false>::Accessor(MatrixType * matrix,
2539  const size_type row,
2540  const size_type index)
2541  : AccessorBase(matrix, row, index)
2542  {}
2543 
2544 
2545  inline Accessor<false>::Reference
2546  Accessor<false>::value() const
2547  {
2548  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2549  return {*this};
2550  }
2551 
2552 
2553 
2554  template <bool Constness>
2555  inline Iterator<Constness>::Iterator(MatrixType * matrix,
2556  const size_type row,
2557  const size_type index)
2558  : accessor(matrix, row, index)
2559  {}
2560 
2561 
2562  template <bool Constness>
2563  template <bool Other>
2564  inline Iterator<Constness>::Iterator(const Iterator<Other> &other)
2565  : accessor(other.accessor)
2566  {}
2567 
2568 
2569  template <bool Constness>
2570  inline Iterator<Constness> &
2572  {
2573  Assert(accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
2574 
2575  ++accessor.a_index;
2576 
2577  // If at end of line: do one
2578  // step, then cycle until we
2579  // find a row with a nonzero
2580  // number of entries.
2581  if (accessor.a_index >= accessor.colnum_cache->size())
2582  {
2583  accessor.a_index = 0;
2584  ++accessor.a_row;
2585 
2586  while ((accessor.a_row < accessor.matrix->m()) &&
2587  ((accessor.matrix->in_local_range(accessor.a_row) == false) ||
2588  (accessor.matrix->row_length(accessor.a_row) == 0)))
2589  ++accessor.a_row;
2590 
2591  accessor.visit_present_row();
2592  }
2593  return *this;
2594  }
2595 
2596 
2597  template <bool Constness>
2598  inline Iterator<Constness>
2600  {
2601  const Iterator<Constness> old_state = *this;
2602  ++(*this);
2603  return old_state;
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 const Accessor<Constness> *Iterator<Constness>::operator->() const
2618  {
2619  return &accessor;
2620  }
2621 
2622 
2623 
2624  template <bool Constness>
2625  inline bool
2626  Iterator<Constness>::operator==(const Iterator<Constness> &other) const
2627  {
2628  return (accessor.a_row == other.accessor.a_row &&
2629  accessor.a_index == other.accessor.a_index);
2630  }
2631 
2632 
2633 
2634  template <bool Constness>
2635  inline bool
2636  Iterator<Constness>::operator!=(const Iterator<Constness> &other) const
2637  {
2638  return !(*this == other);
2639  }
2640 
2641 
2642 
2643  template <bool Constness>
2644  inline bool
2645  Iterator<Constness>::operator<(const Iterator<Constness> &other) const
2646  {
2647  return (accessor.row() < other.accessor.row() ||
2648  (accessor.row() == other.accessor.row() &&
2649  accessor.index() < other.accessor.index()));
2650  }
2651 
2652 
2653  template <bool Constness>
2654  inline bool
2655  Iterator<Constness>::operator>(const Iterator<Constness> &other) const
2656  {
2657  return (other < *this);
2658  }
2659 
2660  } // namespace SparseMatrixIterators
2661 
2662 
2663 
2665  SparseMatrix::begin() const
2666  {
2667  return begin(0);
2668  }
2669 
2670 
2671 
2673  SparseMatrix::end() const
2674  {
2675  return const_iterator(this, m(), 0);
2676  }
2677 
2678 
2679 
2681  SparseMatrix::begin(const size_type r) const
2682  {
2683  AssertIndexRange(r, m());
2684  if (in_local_range(r) && (row_length(r) > 0))
2685  return const_iterator(this, r, 0);
2686  else
2687  return end(r);
2688  }
2689 
2690 
2691 
2693  SparseMatrix::end(const size_type r) const
2694  {
2695  AssertIndexRange(r, m());
2696 
2697  // place the iterator on the first entry
2698  // past this line, or at the end of the
2699  // matrix
2700  for (size_type i = r + 1; i < m(); ++i)
2701  if (in_local_range(i) && (row_length(i) > 0))
2702  return const_iterator(this, i, 0);
2703 
2704  // if there is no such line, then take the
2705  // end iterator of the matrix
2706  return end();
2707  }
2708 
2709 
2710 
2711  inline SparseMatrix::iterator
2713  {
2714  return begin(0);
2715  }
2716 
2717 
2718 
2719  inline SparseMatrix::iterator
2721  {
2722  return iterator(this, m(), 0);
2723  }
2724 
2725 
2726 
2727  inline SparseMatrix::iterator
2729  {
2730  AssertIndexRange(r, m());
2731  if (in_local_range(r) && (row_length(r) > 0))
2732  return iterator(this, r, 0);
2733  else
2734  return end(r);
2735  }
2736 
2737 
2738 
2739  inline SparseMatrix::iterator
2740  SparseMatrix::end(const size_type r)
2741  {
2742  AssertIndexRange(r, m());
2743 
2744  // place the iterator on the first entry
2745  // past this line, or at the end of the
2746  // matrix
2747  for (size_type i = r + 1; i < m(); ++i)
2748  if (in_local_range(i) && (row_length(i) > 0))
2749  return iterator(this, i, 0);
2750 
2751  // if there is no such line, then take the
2752  // end iterator of the matrix
2753  return end();
2754  }
2755 
2756 
2757 
2758  inline bool
2759  SparseMatrix::in_local_range(const size_type index) const
2760  {
2762 # ifndef DEAL_II_WITH_64BIT_INDICES
2763  begin = matrix->RowMap().MinMyGID();
2764  end = matrix->RowMap().MaxMyGID() + 1;
2765 # else
2766  begin = matrix->RowMap().MinMyGID64();
2767  end = matrix->RowMap().MaxMyGID64() + 1;
2768 # endif
2769 
2770  return ((index >= static_cast<size_type>(begin)) &&
2771  (index < static_cast<size_type>(end)));
2772  }
2773 
2774 
2775 
2776  inline bool
2778  {
2779  return compressed;
2780  }
2781 
2782 
2783 
2784  // Inline the set() and add() functions, since they will be called
2785  // frequently, and the compiler can optimize away some unnecessary loops
2786  // when the sizes are given at compile time.
2787  template <>
2788  void
2789  SparseMatrix::set<TrilinosScalar>(const size_type row,
2790  const size_type n_cols,
2791  const size_type * col_indices,
2792  const TrilinosScalar *values,
2793  const bool elide_zero_values);
2794 
2795 
2796 
2797  template <typename Number>
2798  void
2799  SparseMatrix::set(const size_type row,
2800  const size_type n_cols,
2801  const size_type *col_indices,
2802  const Number * values,
2803  const bool elide_zero_values)
2804  {
2805  std::vector<TrilinosScalar> trilinos_values(n_cols);
2806  std::copy(values, values + n_cols, trilinos_values.begin());
2807  this->set(
2808  row, n_cols, col_indices, trilinos_values.data(), elide_zero_values);
2809  }
2810 
2811 
2812 
2813  inline void
2814  SparseMatrix::set(const size_type i,
2815  const size_type j,
2816  const TrilinosScalar value)
2817  {
2818  AssertIsFinite(value);
2819 
2820  set(i, 1, &j, &value, false);
2821  }
2822 
2823 
2824 
2825  inline void
2826  SparseMatrix::set(const std::vector<size_type> & indices,
2827  const FullMatrix<TrilinosScalar> &values,
2828  const bool elide_zero_values)
2829  {
2830  Assert(indices.size() == values.m(),
2831  ExcDimensionMismatch(indices.size(), values.m()));
2832  Assert(values.m() == values.n(), ExcNotQuadratic());
2833 
2834  for (size_type i = 0; i < indices.size(); ++i)
2835  set(indices[i],
2836  indices.size(),
2837  indices.data(),
2838  &values(i, 0),
2839  elide_zero_values);
2840  }
2841 
2842 
2843 
2844  inline void
2845  SparseMatrix::add(const size_type i,
2846  const size_type j,
2847  const TrilinosScalar value)
2848  {
2849  AssertIsFinite(value);
2850 
2851  if (value == 0)
2852  {
2853  // we have to check after Insert/Add in any case to be consistent
2854  // with the MPI communication model, but we can save some
2855  // work if the addend is zero. However, these actions are done in case
2856  // we pass on to the other function.
2857 
2858  // TODO: fix this (do not run compress here, but fail)
2859  if (last_action == Insert)
2860  {
2861  int ierr;
2862  ierr = matrix->GlobalAssemble(*column_space_map,
2863  matrix->RowMap(),
2864  false);
2865 
2866  Assert(ierr == 0, ExcTrilinosError(ierr));
2867  (void)ierr; // removes -Wunused-but-set-variable in optimized mode
2868  }
2869 
2870  last_action = Add;
2871 
2872  return;
2873  }
2874  else
2875  add(i, 1, &j, &value, false);
2876  }
2877 
2878 
2879 
2880  // inline "simple" functions that are called frequently and do only involve
2881  // a call to some Trilinos function.
2883  SparseMatrix::m() const
2884  {
2885 # ifndef DEAL_II_WITH_64BIT_INDICES
2886  return matrix->NumGlobalRows();
2887 # else
2888  return matrix->NumGlobalRows64();
2889 # endif
2890  }
2891 
2892 
2893 
2895  SparseMatrix::n() const
2896  {
2897  // If the matrix structure has not been fixed (i.e., we did not have a
2898  // sparsity pattern), it does not know about the number of columns so we
2899  // must always take this from the additional column space map
2900  Assert(column_space_map.get() != nullptr, ExcInternalError());
2901  return n_global_elements(*column_space_map);
2902  }
2903 
2904 
2905 
2906  inline unsigned int
2907  SparseMatrix::local_size() const
2908  {
2909  return matrix->NumMyRows();
2910  }
2911 
2912 
2913 
2914  inline std::pair<SparseMatrix::size_type, SparseMatrix::size_type>
2916  {
2917  size_type begin, end;
2918 # ifndef DEAL_II_WITH_64BIT_INDICES
2919  begin = matrix->RowMap().MinMyGID();
2920  end = matrix->RowMap().MaxMyGID() + 1;
2921 # else
2922  begin = matrix->RowMap().MinMyGID64();
2923  end = matrix->RowMap().MaxMyGID64() + 1;
2924 # endif
2925 
2926  return std::make_pair(begin, end);
2927  }
2928 
2929 
2930 
2933  {
2934 # ifndef DEAL_II_WITH_64BIT_INDICES
2935  return matrix->NumGlobalNonzeros();
2936 # else
2937  return matrix->NumGlobalNonzeros64();
2938 # endif
2939  }
2940 
2941 
2942 
2943  template <typename SparsityPatternType>
2944  inline void
2945  SparseMatrix::reinit(const IndexSet & parallel_partitioning,
2946  const SparsityPatternType &sparsity_pattern,
2947  const MPI_Comm & communicator,
2948  const bool exchange_data)
2949  {
2950  reinit(parallel_partitioning,
2951  parallel_partitioning,
2952  sparsity_pattern,
2953  communicator,
2954  exchange_data);
2955  }
2956 
2957 
2958 
2959  template <typename number>
2960  inline void
2961  SparseMatrix::reinit(const IndexSet &parallel_partitioning,
2962  const ::SparseMatrix<number> &sparse_matrix,
2963  const MPI_Comm & communicator,
2964  const double drop_tolerance,
2965  const bool copy_values,
2966  const ::SparsityPattern * use_this_sparsity)
2967  {
2968  Epetra_Map map =
2969  parallel_partitioning.make_trilinos_map(communicator, false);
2970  reinit(parallel_partitioning,
2971  parallel_partitioning,
2972  sparse_matrix,
2973  drop_tolerance,
2974  copy_values,
2975  use_this_sparsity);
2976  }
2977 
2978 
2979 
2980  inline const Epetra_CrsMatrix &
2982  {
2983  return static_cast<const Epetra_CrsMatrix &>(*matrix);
2984  }
2985 
2986 
2987 
2988  inline const Epetra_CrsGraph &
2990  {
2991  return matrix->Graph();
2992  }
2993 
2994 
2995 
2996  inline IndexSet
2998  {
2999  return IndexSet(matrix->DomainMap());
3000  }
3001 
3002 
3003 
3004  inline IndexSet
3006  {
3007  return IndexSet(matrix->RangeMap());
3008  }
3009 
3010 
3011 
3012  inline void
3014  {
3015  // nothing to do here
3016  }
3017 
3018 
3019 
3020  inline void
3022  {
3023  // nothing to do here
3024  }
3025 
3026 
3027  namespace internal
3028  {
3029  namespace LinearOperatorImplementation
3030  {
3031  template <typename Solver, typename Preconditioner>
3032  typename std::enable_if<
3033  std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
3034  std::is_base_of<TrilinosWrappers::PreconditionBase,
3035  Preconditioner>::value,
3036  TrilinosPayload>::type
3038  Solver & solver,
3039  const Preconditioner &preconditioner) const
3040  {
3041  const auto &payload = *this;
3042 
3043  TrilinosPayload return_op(payload);
3044 
3045  // Capture by copy so the payloads are always valid
3046 
3047  return_op.inv_vmult = [payload, &solver, &preconditioner](
3048  TrilinosPayload::Domain & tril_dst,
3049  const TrilinosPayload::Range &tril_src) {
3050  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3051  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3052  Assert(&tril_src != &tril_dst,
3055  tril_src,
3056  tril_dst,
3057  !payload.UseTranspose());
3058  solver.solve(payload, tril_dst, tril_src, preconditioner);
3059  };
3060 
3061  return_op.inv_Tvmult = [payload, &solver, &preconditioner](
3062  TrilinosPayload::Range & tril_dst,
3063  const TrilinosPayload::Domain &tril_src) {
3064  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3065  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3066  Assert(&tril_src != &tril_dst,
3069  tril_src,
3070  tril_dst,
3071  payload.UseTranspose());
3072 
3073  const_cast<TrilinosPayload &>(payload).transpose();
3074  solver.solve(payload, tril_dst, tril_src, preconditioner);
3075  const_cast<TrilinosPayload &>(payload).transpose();
3076  };
3077 
3078  // If the input operator is already setup for transpose operations, then
3079  // we must do similar with its inverse.
3080  if (return_op.UseTranspose() == true)
3081  std::swap(return_op.inv_vmult, return_op.inv_Tvmult);
3082 
3083  return return_op;
3084  }
3085 
3086  template <typename Solver, typename Preconditioner>
3087  typename std::enable_if<
3088  !(std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
3089  std::is_base_of<TrilinosWrappers::PreconditionBase,
3090  Preconditioner>::value),
3091  TrilinosPayload>::type
3092  TrilinosPayload::inverse_payload(Solver &, const Preconditioner &) const
3093  {
3094  TrilinosPayload return_op(*this);
3095 
3096  return_op.inv_vmult = [](TrilinosPayload::Domain &,
3097  const TrilinosPayload::Range &) {
3098  AssertThrow(false,
3099  ExcMessage("Payload inv_vmult disabled because of "
3100  "incompatible solver/preconditioner choice."));
3101  };
3102 
3103  return_op.inv_Tvmult = [](TrilinosPayload::Range &,
3104  const TrilinosPayload::Domain &) {
3105  AssertThrow(false,
3106  ExcMessage("Payload inv_vmult disabled because of "
3107  "incompatible solver/preconditioner choice."));
3108  };
3109 
3110  return return_op;
3111  }
3112  } // namespace LinearOperatorImplementation
3113  } // namespace internal
3114 
3115  template <>
3116  void
3117  SparseMatrix::set<TrilinosScalar>(const size_type row,
3118  const size_type n_cols,
3119  const size_type * col_indices,
3120  const TrilinosScalar *values,
3121  const bool elide_zero_values);
3122 # endif // DOXYGEN
3123 
3124 } /* namespace TrilinosWrappers */
3125 
3126 
3128 
3129 
3130 # endif // DEAL_II_WITH_TRILINOS
3131 
3132 
3133 /*----------------------- trilinos_sparse_matrix.h --------------------*/
3134 
3135 #endif
3136 /*----------------------- 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:1636
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:1521
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:596
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:1411
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:2811
#define DeclException0(Exception0)
Definition: exceptions.h:470
typename Accessor< Constness >::MatrixType MatrixType
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
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 > &)
Definition: memory_space.h:103
void reinit(const SparsityPatternType &sparsity_pattern)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
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:2785
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:1667
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