Reference documentation for deal.II version GIT 58febcd5cf 2023-09-30 20:00:01+00:00
\(\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 - 2023 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_trilinos_sparse_matrix_h
17 # define dealii_trilinos_sparse_matrix_h
18 
19 
20 # include <deal.II/base/config.h>
21 
22 # ifdef DEAL_II_WITH_TRILINOS
23 
24 # include <deal.II/base/index_set.h>
25 # include <deal.II/base/mpi_stub.h>
27 
28 # include <deal.II/lac/exceptions.h>
29 # include <deal.II/lac/full_matrix.h>
36 
37 # include <Epetra_Comm.h>
38 # include <Epetra_CrsGraph.h>
39 # include <Epetra_Export.h>
40 # include <Epetra_FECrsMatrix.h>
41 # include <Epetra_Map.h>
42 # include <Epetra_MpiComm.h>
43 # include <Epetra_MultiVector.h>
44 # include <Epetra_Operator.h>
45 
46 # include <cmath>
47 # include <iterator>
48 # include <memory>
49 # include <type_traits>
50 # include <vector>
51 
53 
54 // forward declarations
55 # ifndef DOXYGEN
56 template <typename MatrixType>
57 class BlockMatrixBase;
58 
59 template <typename number>
60 class SparseMatrix;
61 class SparsityPattern;
63 
64 namespace TrilinosWrappers
65 {
66  class SparseMatrix;
67  class SparsityPattern;
68 
69  namespace SparseMatrixIterators
70  {
71  template <bool Constness>
72  class Iterator;
73  }
74 } // namespace TrilinosWrappers
75 # endif
76 
77 namespace TrilinosWrappers
78 {
83  {
88 
93  std::size_t,
94  std::size_t,
95  std::size_t,
96  << "You tried to access row " << arg1
97  << " of a distributed sparsity pattern, "
98  << " but only rows " << arg2 << " through " << arg3
99  << " are stored locally and can be accessed.");
100 
109  {
110  public:
115 
120  const size_type row,
121  const size_type index);
122 
126  size_type
127  row() const;
128 
132  size_type
133  index() const;
134 
138  size_type
139  column() const;
140 
141  protected:
153 
158 
164  void
166 
179  std::shared_ptr<std::vector<size_type>> colnum_cache;
180 
184  std::shared_ptr<std::vector<TrilinosScalar>> value_cache;
185  };
186 
197  template <bool Constess>
198  class Accessor : public AccessorBase
199  {
204  value() const;
205 
210  value();
211  };
212 
213 
214 
218  template <>
219  class Accessor<true> : public AccessorBase
220  {
221  public:
226  using MatrixType = const SparseMatrix;
227 
233 
238  template <bool Other>
240 
245  value() const;
246 
247  private:
248  // Make iterator class a friend.
249  template <bool>
250  friend class Iterator;
251  };
252 
256  template <>
257  class Accessor<false> : public AccessorBase
258  {
259  class Reference
260  {
261  public:
265  Reference(const Accessor<false> &accessor);
266 
270  operator TrilinosScalar() const;
271 
275  const Reference &
276  operator=(const TrilinosScalar n) const;
277 
281  const Reference &
282  operator+=(const TrilinosScalar n) const;
283 
287  const Reference &
288  operator-=(const TrilinosScalar n) const;
289 
293  const Reference &
294  operator*=(const TrilinosScalar n) const;
295 
299  const Reference &
300  operator/=(const TrilinosScalar n) const;
301 
302  private:
308  };
309 
310  public:
316 
322 
326  Reference
327  value() const;
328 
329  private:
330  // Make iterator class a friend.
331  template <bool>
332  friend class Iterator;
333 
334  // Make Reference object a friend.
335  friend class Reference;
336  };
337 
351  template <bool Constness>
352  class Iterator
353  {
354  public:
359 
365 
371 
377 
382  Iterator(MatrixType *matrix, const size_type row, const size_type index);
383 
387  template <bool Other>
388  Iterator(const Iterator<Other> &other);
389 
395 
401 
405  const Accessor<Constness> &
406  operator*() const;
407 
411  const Accessor<Constness> *
412  operator->() const;
413 
418  template <bool OtherConstness>
419  bool
421 
425  template <bool OtherConstness>
426  bool
428 
434  template <bool OtherConstness>
435  bool
436  operator<(const Iterator<OtherConstness> &) const;
437 
441  template <bool OtherConstness>
442  bool
444 
449  size_type,
450  size_type,
451  << "Attempt to access element " << arg2 << " of row "
452  << arg1 << " which doesn't have that many elements.");
453 
454  private:
459 
460  template <bool Other>
461  friend class Iterator;
462  };
463 
464  } // namespace SparseMatrixIterators
465 } // namespace TrilinosWrappers
466 
468 
469 namespace std
470 {
471  template <bool Constness>
474  {
475  using iterator_category = forward_iterator_tag;
476  using value_type =
477  typename ::TrilinosWrappers::SparseMatrixIterators::Iterator<
478  Constness>::value_type;
480  typename ::TrilinosWrappers::SparseMatrixIterators::Iterator<
481  Constness>::difference_type;
482  };
483 } // namespace std
484 
486 
487 
488 namespace TrilinosWrappers
489 {
550  class SparseMatrix : public Subscriptor
551  {
552  public:
557 
562  std::size_t,
563  << "You tried to access row " << arg1
564  << " of a non-contiguous locally owned row set."
565  << " The row " << arg1
566  << " is not stored locally and can't be accessed.");
567 
575  struct Traits
576  {
581  static const bool zero_addition_can_be_elided = true;
582  };
583 
588 
593 
598 
606  SparseMatrix();
607 
615  SparseMatrix(const size_type m,
616  const size_type n,
617  const unsigned int n_max_entries_per_row);
618 
626  SparseMatrix(const size_type m,
627  const size_type n,
628  const std::vector<unsigned int> &n_entries_per_row);
629 
633  SparseMatrix(const SparsityPattern &InputSparsityPattern);
634 
639  SparseMatrix(SparseMatrix &&other) noexcept;
640 
644  SparseMatrix(const SparseMatrix &) = delete;
645 
649  SparseMatrix &
650  operator=(const SparseMatrix &) = delete;
651 
655  virtual ~SparseMatrix() override = default;
656 
672  template <typename SparsityPatternType>
673  void
674  reinit(const SparsityPatternType &sparsity_pattern);
675 
688  void
689  reinit(const SparsityPattern &sparsity_pattern);
690 
699  void
700  reinit(const SparseMatrix &sparse_matrix);
701 
722  template <typename number>
723  void
724  reinit(const ::SparseMatrix<number> &dealii_sparse_matrix,
725  const double drop_tolerance = 1e-13,
726  const bool copy_values = true,
727  const ::SparsityPattern *use_this_sparsity = nullptr);
728 
734  void
735  reinit(const Epetra_CrsMatrix &input_matrix, const bool copy_values = true);
754  SparseMatrix(const IndexSet &parallel_partitioning,
755  const MPI_Comm communicator = MPI_COMM_WORLD,
756  const unsigned int n_max_entries_per_row = 0);
757 
765  SparseMatrix(const IndexSet &parallel_partitioning,
766  const MPI_Comm communicator,
767  const std::vector<unsigned int> &n_entries_per_row);
768 
783  SparseMatrix(const IndexSet &row_parallel_partitioning,
784  const IndexSet &col_parallel_partitioning,
785  const MPI_Comm communicator = MPI_COMM_WORLD,
786  const size_type n_max_entries_per_row = 0);
787 
802  SparseMatrix(const IndexSet &row_parallel_partitioning,
803  const IndexSet &col_parallel_partitioning,
804  const MPI_Comm communicator,
805  const std::vector<unsigned int> &n_entries_per_row);
806 
827  template <typename SparsityPatternType>
828  void
829  reinit(const IndexSet &parallel_partitioning,
830  const SparsityPatternType &sparsity_pattern,
831  const MPI_Comm communicator = MPI_COMM_WORLD,
832  const bool exchange_data = false);
833 
846  template <typename SparsityPatternType>
847  std::enable_if_t<
848  !std::is_same_v<SparsityPatternType, ::SparseMatrix<double>>>
849  reinit(const IndexSet &row_parallel_partitioning,
850  const IndexSet &col_parallel_partitioning,
851  const SparsityPatternType &sparsity_pattern,
852  const MPI_Comm communicator = MPI_COMM_WORLD,
853  const bool exchange_data = false);
854 
871  template <typename number>
872  void
873  reinit(const IndexSet &parallel_partitioning,
874  const ::SparseMatrix<number> &dealii_sparse_matrix,
875  const MPI_Comm communicator = MPI_COMM_WORLD,
876  const double drop_tolerance = 1e-13,
877  const bool copy_values = true,
878  const ::SparsityPattern *use_this_sparsity = nullptr);
879 
893  template <typename number>
894  void
895  reinit(const IndexSet &row_parallel_partitioning,
896  const IndexSet &col_parallel_partitioning,
897  const ::SparseMatrix<number> &dealii_sparse_matrix,
898  const MPI_Comm communicator = MPI_COMM_WORLD,
899  const double drop_tolerance = 1e-13,
900  const bool copy_values = true,
901  const ::SparsityPattern *use_this_sparsity = nullptr);
911  size_type
912  m() const;
913 
917  size_type
918  n() const;
919 
928  unsigned int
929  local_size() const;
930 
939  std::pair<size_type, size_type>
940  local_range() const;
941 
946  bool
947  in_local_range(const size_type index) const;
948 
953  std::uint64_t
955 
959  unsigned int
960  row_length(const size_type row) const;
961 
968  bool
969  is_compressed() const;
970 
976  size_type
977  memory_consumption() const;
978 
982  MPI_Comm
983  get_mpi_communicator() const;
984 
1000  SparseMatrix &
1001  operator=(const double d);
1002 
1010  void
1011  clear();
1012 
1040  void
1041  compress(VectorOperation::values operation);
1042 
1064  void
1065  set(const size_type i, const size_type j, const TrilinosScalar value);
1066 
1099  void
1100  set(const std::vector<size_type> &indices,
1101  const FullMatrix<TrilinosScalar> &full_matrix,
1102  const bool elide_zero_values = false);
1103 
1109  void
1110  set(const std::vector<size_type> &row_indices,
1111  const std::vector<size_type> &col_indices,
1112  const FullMatrix<TrilinosScalar> &full_matrix,
1113  const bool elide_zero_values = false);
1114 
1142  void
1143  set(const size_type row,
1144  const std::vector<size_type> &col_indices,
1145  const std::vector<TrilinosScalar> &values,
1146  const bool elide_zero_values = false);
1147 
1175  template <typename Number>
1176  void
1177  set(const size_type row,
1178  const size_type n_cols,
1179  const size_type *col_indices,
1180  const Number *values,
1181  const bool elide_zero_values = false);
1182 
1192  void
1193  add(const size_type i, const size_type j, const TrilinosScalar value);
1194 
1213  void
1214  add(const std::vector<size_type> &indices,
1215  const FullMatrix<TrilinosScalar> &full_matrix,
1216  const bool elide_zero_values = true);
1217 
1223  void
1224  add(const std::vector<size_type> &row_indices,
1225  const std::vector<size_type> &col_indices,
1226  const FullMatrix<TrilinosScalar> &full_matrix,
1227  const bool elide_zero_values = true);
1228 
1242  void
1243  add(const size_type row,
1244  const std::vector<size_type> &col_indices,
1245  const std::vector<TrilinosScalar> &values,
1246  const bool elide_zero_values = true);
1247 
1261  void
1262  add(const size_type row,
1263  const size_type n_cols,
1264  const size_type *col_indices,
1265  const TrilinosScalar *values,
1266  const bool elide_zero_values = true,
1267  const bool col_indices_are_sorted = false);
1268 
1272  SparseMatrix &
1273  operator*=(const TrilinosScalar factor);
1274 
1278  SparseMatrix &
1279  operator/=(const TrilinosScalar factor);
1280 
1284  void
1285  copy_from(const SparseMatrix &source);
1286 
1294  void
1295  add(const TrilinosScalar factor, const SparseMatrix &matrix);
1296 
1323  void
1324  clear_row(const size_type row, const TrilinosScalar new_diag_value = 0);
1325 
1346  void
1347  clear_rows(const std::vector<size_type> &rows,
1348  const TrilinosScalar new_diag_value = 0);
1349 
1359  void
1360  transpose();
1361 
1377  operator()(const size_type i, const size_type j) const;
1378 
1396  el(const size_type i, const size_type j) const;
1397 
1405  diag_element(const size_type i) const;
1406 
1443  template <typename VectorType>
1444  std::enable_if_t<
1445  std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
1446  vmult(VectorType &dst, const VectorType &src) const;
1447 
1454  template <typename VectorType>
1455  std::enable_if_t<
1456  !std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
1457  vmult(VectorType &dst, const VectorType &src) const;
1458 
1473  template <typename VectorType>
1474  std::enable_if_t<
1475  std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
1476  Tvmult(VectorType &dst, const VectorType &src) const;
1477 
1484  template <typename VectorType>
1485  std::enable_if_t<
1486  !std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
1487  Tvmult(VectorType &dst, const VectorType &src) const;
1488 
1498  template <typename VectorType>
1499  void
1500  vmult_add(VectorType &dst, const VectorType &src) const;
1501 
1512  template <typename VectorType>
1513  void
1514  Tvmult_add(VectorType &dst, const VectorType &src) const;
1515 
1538  matrix_norm_square(const MPI::Vector &v) const;
1539 
1560  matrix_scalar_product(const MPI::Vector &u, const MPI::Vector &v) const;
1561 
1579  residual(MPI::Vector &dst,
1580  const MPI::Vector &x,
1581  const MPI::Vector &b) const;
1582 
1597  void
1598  mmult(SparseMatrix &C,
1599  const SparseMatrix &B,
1600  const MPI::Vector &V = MPI::Vector()) const;
1601 
1602 
1619  void
1621  const SparseMatrix &B,
1622  const MPI::Vector &V = MPI::Vector()) const;
1623 
1638  l1_norm() const;
1639 
1649  linfty_norm() const;
1650 
1656  frobenius_norm() const;
1657 
1668  const Epetra_CrsMatrix &
1670 
1675  const Epetra_CrsGraph &
1677 
1689  IndexSet
1691 
1697  IndexSet
1699 
1726  begin() const;
1727 
1731  iterator
1733 
1739  end() const;
1740 
1744  iterator
1745  end();
1746 
1776  begin(const size_type r) const;
1777 
1781  iterator
1782  begin(const size_type r);
1783 
1794  end(const size_type r) const;
1795 
1799  iterator
1800  end(const size_type r);
1801 
1813  void
1814  write_ascii();
1815 
1823  void
1824  print(std::ostream &out,
1825  const bool write_extended_trilinos_info = false) const;
1826 
1837  int,
1838  << "An error with error number " << arg1
1839  << " occurred while calling a Trilinos function. "
1840  "\n\n"
1841  "For historical reasons, many Trilinos functions express "
1842  "errors by returning specific integer values to indicate "
1843  "certain errors. Unfortunately, different Trilinos functions "
1844  "often use the same integer values for different kinds of "
1845  "errors, and in most cases it is also not documented what "
1846  "each error code actually means. As a consequence, it is often "
1847  "difficult to say what a particular error (in this case, "
1848  "the error with integer code '"
1849  << arg1
1850  << "') represents and how one should fix a code to avoid it. "
1851  "The best one can often do is to look up the call stack to "
1852  "see which deal.II function generated the error, and which "
1853  "Trilinos function the error code had originated from; "
1854  "then look up the Trilinos source code of that function (for "
1855  "example on github) to see what code path set that error "
1856  "code. Short of going through all of that, the only other "
1857  "option is to guess the cause of the error from "
1858  "the context in which the error appeared.");
1859 
1860 
1865  size_type,
1866  size_type,
1867  << "The entry with index <" << arg1 << ',' << arg2
1868  << "> does not exist.");
1869 
1874  "You are attempting an operation on two vectors that "
1875  "are the same object, but the operation requires that the "
1876  "two objects are in fact different.");
1877 
1882 
1887  size_type,
1888  size_type,
1889  size_type,
1890  size_type,
1891  << "You tried to access element (" << arg1 << '/' << arg2
1892  << ')'
1893  << " of a distributed matrix, but only rows in range ["
1894  << arg3 << ',' << arg4
1895  << "] are stored locally and can be accessed.");
1896 
1901  size_type,
1902  size_type,
1903  << "You tried to access element (" << arg1 << '/' << arg2
1904  << ')' << " of a sparse matrix, but it appears to not"
1905  << " exist in the Trilinos sparsity pattern.");
1910  protected:
1921  void
1923 
1929  void
1931 
1932 
1933 
1934  private:
1939  std::unique_ptr<Epetra_Map> column_space_map;
1940 
1946  std::unique_ptr<Epetra_FECrsMatrix> matrix;
1947 
1953  std::unique_ptr<Epetra_CrsMatrix> nonlocal_matrix;
1954 
1958  std::unique_ptr<Epetra_Export> nonlocal_matrix_exporter;
1959 
1971  Epetra_CombineMode last_action;
1972 
1978 
1979  // To allow calling protected prepare_add() and prepare_set().
1980  friend class BlockMatrixBase<SparseMatrix>;
1981  };
1982 
1983 
1984 
1985  // forwards declarations
1986  class SolverBase;
1987  class PreconditionBase;
1988 
1989  namespace internal
1990  {
1991  inline void
1992  check_vector_map_equality(const Epetra_CrsMatrix &mtrx,
1993  const Epetra_MultiVector &src,
1994  const Epetra_MultiVector &dst,
1995  const bool transpose)
1996  {
1997  if (transpose == false)
1998  {
1999  Assert(src.Map().SameAs(mtrx.DomainMap()) == true,
2000  ExcMessage(
2001  "Column map of matrix does not fit with vector map!"));
2002  Assert(dst.Map().SameAs(mtrx.RangeMap()) == true,
2003  ExcMessage("Row map of matrix does not fit with vector map!"));
2004  }
2005  else
2006  {
2007  Assert(src.Map().SameAs(mtrx.RangeMap()) == true,
2008  ExcMessage(
2009  "Column map of matrix does not fit with vector map!"));
2010  Assert(dst.Map().SameAs(mtrx.DomainMap()) == true,
2011  ExcMessage("Row map of matrix does not fit with vector map!"));
2012  }
2013  (void)mtrx; // removes -Wunused-variable in optimized mode
2014  (void)src;
2015  (void)dst;
2016  }
2017 
2018  inline void
2020  const Epetra_MultiVector &src,
2021  const Epetra_MultiVector &dst,
2022  const bool transpose)
2023  {
2024  if (transpose == false)
2025  {
2026  Assert(src.Map().SameAs(op.OperatorDomainMap()) == true,
2027  ExcMessage(
2028  "Column map of operator does not fit with vector map!"));
2029  Assert(dst.Map().SameAs(op.OperatorRangeMap()) == true,
2030  ExcMessage(
2031  "Row map of operator does not fit with vector map!"));
2032  }
2033  else
2034  {
2035  Assert(src.Map().SameAs(op.OperatorRangeMap()) == true,
2036  ExcMessage(
2037  "Column map of operator does not fit with vector map!"));
2038  Assert(dst.Map().SameAs(op.OperatorDomainMap()) == true,
2039  ExcMessage(
2040  "Row map of operator does not fit with vector map!"));
2041  }
2042  (void)op; // removes -Wunused-variable in optimized mode
2043  (void)src;
2044  (void)dst;
2045  }
2046 
2047  namespace LinearOperatorImplementation
2048  {
2069  {
2070  public:
2074  using VectorType = Epetra_MultiVector;
2075 
2080 
2085 
2098  TrilinosPayload();
2099 
2103  TrilinosPayload(const TrilinosWrappers::SparseMatrix &matrix_exemplar,
2105 
2109  TrilinosPayload(const TrilinosPayload &payload_exemplar,
2111 
2116  const TrilinosWrappers::SparseMatrix &matrix_exemplar,
2117  const TrilinosWrappers::PreconditionBase &preconditioner);
2118 
2123  const TrilinosWrappers::PreconditionBase &preconditioner_exemplar,
2124  const TrilinosWrappers::PreconditionBase &preconditioner);
2125 
2130  const TrilinosPayload &payload_exemplar,
2131  const TrilinosWrappers::PreconditionBase &preconditioner);
2132 
2136  TrilinosPayload(const TrilinosPayload &payload);
2137 
2145  TrilinosPayload(const TrilinosPayload &first_op,
2146  const TrilinosPayload &second_op);
2147 
2151  virtual ~TrilinosPayload() override = default;
2152 
2157  identity_payload() const;
2158 
2163  null_payload() const;
2164 
2169  transpose_payload() const;
2170 
2187  template <typename Solver, typename Preconditioner>
2188  std::enable_if_t<
2189  std::is_base_of_v<TrilinosWrappers::SolverBase, Solver> &&
2190  std::is_base_of_v<TrilinosWrappers::PreconditionBase,
2191  Preconditioner>,
2193  inverse_payload(Solver &, const Preconditioner &) const;
2194 
2212  template <typename Solver, typename Preconditioner>
2213  std::enable_if_t<
2214  !(std::is_base_of_v<TrilinosWrappers::SolverBase, Solver> &&
2215  std::is_base_of_v<TrilinosWrappers::PreconditionBase,
2216  Preconditioner>),
2218  inverse_payload(Solver &, const Preconditioner &) const;
2219 
2232  IndexSet
2234 
2240  IndexSet
2242 
2246  MPI_Comm
2247  get_mpi_communicator() const;
2248 
2255  void
2256  transpose();
2257 
2265  std::function<void(VectorType &, const VectorType &)> vmult;
2266 
2274  std::function<void(VectorType &, const VectorType &)> Tvmult;
2275 
2284  std::function<void(VectorType &, const VectorType &)> inv_vmult;
2285 
2294  std::function<void(VectorType &, const VectorType &)> inv_Tvmult;
2295 
2309  virtual bool
2310  UseTranspose() const override;
2311 
2327  virtual int
2328  SetUseTranspose(bool UseTranspose) override;
2329 
2341  virtual int
2342  Apply(const VectorType &X, VectorType &Y) const override;
2343 
2362  virtual int
2363  ApplyInverse(const VectorType &Y, VectorType &X) const override;
2377  virtual const char *
2378  Label() const override;
2379 
2387  virtual const Epetra_Comm &
2388  Comm() const override;
2389 
2397  virtual const Epetra_Map &
2398  OperatorDomainMap() const override;
2399 
2408  virtual const Epetra_Map &
2409  OperatorRangeMap() const override;
2412  private:
2422  template <typename EpetraOpType>
2423  TrilinosPayload(EpetraOpType &op,
2424  const bool supports_inverse_operations,
2425  const bool use_transpose,
2426  const MPI_Comm mpi_communicator,
2429 
2435 
2440  Epetra_MpiComm communicator;
2441 
2446  Epetra_Map domain_map;
2447 
2452  Epetra_Map range_map;
2453 
2462  virtual bool
2463  HasNormInf() const override;
2464 
2472  virtual double
2473  NormInf() const override;
2474  };
2475 
2481  operator+(const TrilinosPayload &first_op,
2482  const TrilinosPayload &second_op);
2483 
2489  operator*(const TrilinosPayload &first_op,
2490  const TrilinosPayload &second_op);
2491 
2492  } // namespace LinearOperatorImplementation
2493  } /* namespace internal */
2494 
2495 
2496 
2497  // ----------------------- inline and template functions --------------------
2498 
2499 # ifndef DOXYGEN
2500 
2501  namespace SparseMatrixIterators
2502  {
2504  size_type row,
2505  size_type index)
2506  : matrix(matrix)
2507  , a_row(row)
2508  , a_index(index)
2509  {
2510  visit_present_row();
2511  }
2512 
2513 
2515  AccessorBase::row() const
2516  {
2517  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2518  return a_row;
2519  }
2520 
2521 
2523  AccessorBase::column() const
2524  {
2525  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2526  return (*colnum_cache)[a_index];
2527  }
2528 
2529 
2531  AccessorBase::index() const
2532  {
2533  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2534  return a_index;
2535  }
2536 
2537 
2538  inline Accessor<true>::Accessor(MatrixType *matrix,
2539  const size_type row,
2540  const size_type index)
2541  : AccessorBase(const_cast<SparseMatrix *>(matrix), row, index)
2542  {}
2543 
2544 
2545  template <bool Other>
2546  inline Accessor<true>::Accessor(const Accessor<Other> &other)
2547  : AccessorBase(other)
2548  {}
2549 
2550 
2551  inline TrilinosScalar
2552  Accessor<true>::value() const
2553  {
2554  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2555  return (*value_cache)[a_index];
2556  }
2557 
2558 
2559  inline Accessor<false>::Reference::Reference(const Accessor<false> &acc)
2560  : accessor(const_cast<Accessor<false> &>(acc))
2561  {}
2562 
2563 
2564  inline Accessor<false>::Reference::operator TrilinosScalar() const
2565  {
2566  return (*accessor.value_cache)[accessor.a_index];
2567  }
2568 
2569  inline const Accessor<false>::Reference &
2570  Accessor<false>::Reference::operator=(const TrilinosScalar n) const
2571  {
2572  (*accessor.value_cache)[accessor.a_index] = n;
2573  accessor.matrix->set(accessor.row(),
2574  accessor.column(),
2575  static_cast<TrilinosScalar>(*this));
2576  return *this;
2577  }
2578 
2579 
2580  inline const Accessor<false>::Reference &
2581  Accessor<false>::Reference::operator+=(const TrilinosScalar n) const
2582  {
2583  (*accessor.value_cache)[accessor.a_index] += n;
2584  accessor.matrix->set(accessor.row(),
2585  accessor.column(),
2586  static_cast<TrilinosScalar>(*this));
2587  return *this;
2588  }
2589 
2590 
2591  inline const Accessor<false>::Reference &
2592  Accessor<false>::Reference::operator-=(const TrilinosScalar n) const
2593  {
2594  (*accessor.value_cache)[accessor.a_index] -= n;
2595  accessor.matrix->set(accessor.row(),
2596  accessor.column(),
2597  static_cast<TrilinosScalar>(*this));
2598  return *this;
2599  }
2600 
2601 
2602  inline const Accessor<false>::Reference &
2603  Accessor<false>::Reference::operator*=(const TrilinosScalar n) const
2604  {
2605  (*accessor.value_cache)[accessor.a_index] *= n;
2606  accessor.matrix->set(accessor.row(),
2607  accessor.column(),
2608  static_cast<TrilinosScalar>(*this));
2609  return *this;
2610  }
2611 
2612 
2613  inline const Accessor<false>::Reference &
2614  Accessor<false>::Reference::operator/=(const TrilinosScalar n) const
2615  {
2616  (*accessor.value_cache)[accessor.a_index] /= n;
2617  accessor.matrix->set(accessor.row(),
2618  accessor.column(),
2619  static_cast<TrilinosScalar>(*this));
2620  return *this;
2621  }
2622 
2623 
2624  inline Accessor<false>::Accessor(MatrixType *matrix,
2625  const size_type row,
2626  const size_type index)
2627  : AccessorBase(matrix, row, index)
2628  {}
2629 
2630 
2631  inline Accessor<false>::Reference
2632  Accessor<false>::value() const
2633  {
2634  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2635  return {*this};
2636  }
2637 
2638 
2639 
2640  template <bool Constness>
2641  inline Iterator<Constness>::Iterator(MatrixType *matrix,
2642  const size_type row,
2643  const size_type index)
2644  : accessor(matrix, row, index)
2645  {}
2646 
2647 
2648  template <bool Constness>
2649  template <bool Other>
2650  inline Iterator<Constness>::Iterator(const Iterator<Other> &other)
2651  : accessor(other.accessor)
2652  {}
2653 
2654 
2655  template <bool Constness>
2656  inline Iterator<Constness> &
2658  {
2659  Assert(accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
2660 
2661  ++accessor.a_index;
2662 
2663  // If at end of line: do one
2664  // step, then cycle until we
2665  // find a row with a nonzero
2666  // number of entries.
2667  if (accessor.a_index >= accessor.colnum_cache->size())
2668  {
2669  accessor.a_index = 0;
2670  ++accessor.a_row;
2671 
2672  while ((accessor.a_row < accessor.matrix->m()) &&
2673  ((accessor.matrix->in_local_range(accessor.a_row) == false) ||
2674  (accessor.matrix->row_length(accessor.a_row) == 0)))
2675  ++accessor.a_row;
2676 
2677  accessor.visit_present_row();
2678  }
2679  return *this;
2680  }
2681 
2682 
2683  template <bool Constness>
2684  inline Iterator<Constness>
2686  {
2687  const Iterator<Constness> old_state = *this;
2688  ++(*this);
2689  return old_state;
2690  }
2691 
2692 
2693 
2694  template <bool Constness>
2695  inline const Accessor<Constness> &
2697  {
2698  return accessor;
2699  }
2700 
2701 
2702 
2703  template <bool Constness>
2704  inline const Accessor<Constness> *
2705  Iterator<Constness>::operator->() const
2706  {
2707  return &accessor;
2708  }
2709 
2710 
2711 
2712  template <bool Constness>
2713  template <bool OtherConstness>
2714  inline bool
2715  Iterator<Constness>::operator==(const Iterator<OtherConstness> &other) const
2716  {
2717  return (accessor.a_row == other.accessor.a_row &&
2718  accessor.a_index == other.accessor.a_index);
2719  }
2720 
2721 
2722 
2723  template <bool Constness>
2724  template <bool OtherConstness>
2725  inline bool
2726  Iterator<Constness>::operator!=(const Iterator<OtherConstness> &other) const
2727  {
2728  return !(*this == other);
2729  }
2730 
2731 
2732 
2733  template <bool Constness>
2734  template <bool OtherConstness>
2735  inline bool
2736  Iterator<Constness>::operator<(const Iterator<OtherConstness> &other) const
2737  {
2738  return (accessor.row() < other.accessor.row() ||
2739  (accessor.row() == other.accessor.row() &&
2740  accessor.index() < other.accessor.index()));
2741  }
2742 
2743 
2744  template <bool Constness>
2745  template <bool OtherConstness>
2746  inline bool
2747  Iterator<Constness>::operator>(const Iterator<OtherConstness> &other) const
2748  {
2749  return (other < *this);
2750  }
2751 
2752  } // namespace SparseMatrixIterators
2753 
2754 
2755 
2757  SparseMatrix::begin() const
2758  {
2759  return begin(0);
2760  }
2761 
2762 
2763 
2765  SparseMatrix::end() const
2766  {
2767  return const_iterator(this, m(), 0);
2768  }
2769 
2770 
2771 
2773  SparseMatrix::begin(const size_type r) const
2774  {
2775  AssertIndexRange(r, m());
2776  if (in_local_range(r) && (row_length(r) > 0))
2777  return const_iterator(this, r, 0);
2778  else
2779  return end(r);
2780  }
2781 
2782 
2783 
2785  SparseMatrix::end(const size_type r) const
2786  {
2787  AssertIndexRange(r, m());
2788 
2789  // place the iterator on the first entry
2790  // past this line, or at the end of the
2791  // matrix
2792  for (size_type i = r + 1; i < m(); ++i)
2793  if (in_local_range(i) && (row_length(i) > 0))
2794  return const_iterator(this, i, 0);
2795 
2796  // if there is no such line, then take the
2797  // end iterator of the matrix
2798  return end();
2799  }
2800 
2801 
2802 
2803  inline SparseMatrix::iterator
2805  {
2806  return begin(0);
2807  }
2808 
2809 
2810 
2811  inline SparseMatrix::iterator
2813  {
2814  return iterator(this, m(), 0);
2815  }
2816 
2817 
2818 
2819  inline SparseMatrix::iterator
2821  {
2822  AssertIndexRange(r, m());
2823  if (in_local_range(r) && (row_length(r) > 0))
2824  return iterator(this, r, 0);
2825  else
2826  return end(r);
2827  }
2828 
2829 
2830 
2831  inline SparseMatrix::iterator
2832  SparseMatrix::end(const size_type r)
2833  {
2834  AssertIndexRange(r, m());
2835 
2836  // place the iterator on the first entry
2837  // past this line, or at the end of the
2838  // matrix
2839  for (size_type i = r + 1; i < m(); ++i)
2840  if (in_local_range(i) && (row_length(i) > 0))
2841  return iterator(this, i, 0);
2842 
2843  // if there is no such line, then take the
2844  // end iterator of the matrix
2845  return end();
2846  }
2847 
2848 
2849 
2850  inline bool
2851  SparseMatrix::in_local_range(const size_type index) const
2852  {
2854 # ifndef DEAL_II_WITH_64BIT_INDICES
2855  begin = matrix->RowMap().MinMyGID();
2856  end = matrix->RowMap().MaxMyGID() + 1;
2857 # else
2858  begin = matrix->RowMap().MinMyGID64();
2859  end = matrix->RowMap().MaxMyGID64() + 1;
2860 # endif
2861 
2862  return ((index >= static_cast<size_type>(begin)) &&
2863  (index < static_cast<size_type>(end)));
2864  }
2865 
2866 
2867 
2868  inline bool
2870  {
2871  return compressed;
2872  }
2873 
2874 
2875 
2876  // Inline the set() and add() functions, since they will be called
2877  // frequently, and the compiler can optimize away some unnecessary loops
2878  // when the sizes are given at compile time.
2879  template <>
2880  void
2881  SparseMatrix::set<TrilinosScalar>(const size_type row,
2882  const size_type n_cols,
2883  const size_type *col_indices,
2884  const TrilinosScalar *values,
2885  const bool elide_zero_values);
2886 
2887 
2888 
2889  template <typename Number>
2890  void
2891  SparseMatrix::set(const size_type row,
2892  const size_type n_cols,
2893  const size_type *col_indices,
2894  const Number *values,
2895  const bool elide_zero_values)
2896  {
2897  std::vector<TrilinosScalar> trilinos_values(n_cols);
2898  std::copy(values, values + n_cols, trilinos_values.begin());
2899  this->set(
2900  row, n_cols, col_indices, trilinos_values.data(), elide_zero_values);
2901  }
2902 
2903 
2904 
2905  inline void
2906  SparseMatrix::set(const size_type i,
2907  const size_type j,
2908  const TrilinosScalar value)
2909  {
2910  AssertIsFinite(value);
2911 
2912  set(i, 1, &j, &value, false);
2913  }
2914 
2915 
2916 
2917  inline void
2918  SparseMatrix::set(const std::vector<size_type> &indices,
2920  const bool elide_zero_values)
2921  {
2922  Assert(indices.size() == values.m(),
2923  ExcDimensionMismatch(indices.size(), values.m()));
2924  Assert(values.m() == values.n(), ExcNotQuadratic());
2925 
2926  for (size_type i = 0; i < indices.size(); ++i)
2927  set(indices[i],
2928  indices.size(),
2929  indices.data(),
2930  &values(i, 0),
2931  elide_zero_values);
2932  }
2933 
2934 
2935 
2936  inline void
2937  SparseMatrix::add(const size_type i,
2938  const size_type j,
2939  const TrilinosScalar value)
2940  {
2941  AssertIsFinite(value);
2942 
2943  if (value == 0)
2944  {
2945  // we have to check after Insert/Add in any case to be consistent
2946  // with the MPI communication model, but we can save some
2947  // work if the addend is zero. However, these actions are done in case
2948  // we pass on to the other function.
2949 
2950  // TODO: fix this (do not run compress here, but fail)
2951  if (last_action == Insert)
2952  {
2953  const int ierr = matrix->GlobalAssemble(*column_space_map,
2954  matrix->RowMap(),
2955  false);
2956 
2957  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2958  }
2959 
2960  last_action = Add;
2961 
2962  return;
2963  }
2964  else
2965  add(i, 1, &j, &value, false);
2966  }
2967 
2968 
2969 
2970  // inline "simple" functions that are called frequently and do only involve
2971  // a call to some Trilinos function.
2973  SparseMatrix::m() const
2974  {
2975 # ifndef DEAL_II_WITH_64BIT_INDICES
2976  return matrix->NumGlobalRows();
2977 # else
2978  return matrix->NumGlobalRows64();
2979 # endif
2980  }
2981 
2982 
2983 
2985  SparseMatrix::n() const
2986  {
2987  // If the matrix structure has not been fixed (i.e., we did not have a
2988  // sparsity pattern), it does not know about the number of columns so we
2989  // must always take this from the additional column space map
2990  Assert(column_space_map.get() != nullptr, ExcInternalError());
2991  return n_global_elements(*column_space_map);
2992  }
2993 
2994 
2995 
2996  inline unsigned int
2997  SparseMatrix::local_size() const
2998  {
2999  return matrix->NumMyRows();
3000  }
3001 
3002 
3003 
3004  inline std::pair<SparseMatrix::size_type, SparseMatrix::size_type>
3006  {
3007  size_type begin, end;
3008 # ifndef DEAL_II_WITH_64BIT_INDICES
3009  begin = matrix->RowMap().MinMyGID();
3010  end = matrix->RowMap().MaxMyGID() + 1;
3011 # else
3012  begin = matrix->RowMap().MinMyGID64();
3013  end = matrix->RowMap().MaxMyGID64() + 1;
3014 # endif
3015 
3016  return std::make_pair(begin, end);
3017  }
3018 
3019 
3020 
3021  inline std::uint64_t
3023  {
3024  // Trilinos uses 64bit functions internally for attribute access, which
3025  // return `long long`. They also offer 32bit variants that return `int`,
3026  // however those call the 64bit version and convert the values to 32bit.
3027  // There is no necessity in using the 32bit versions at all.
3028  return static_cast<std::uint64_t>(matrix->NumGlobalNonzeros64());
3029  }
3030 
3031 
3032 
3033  template <typename SparsityPatternType>
3034  inline void
3035  SparseMatrix::reinit(const IndexSet &parallel_partitioning,
3036  const SparsityPatternType &sparsity_pattern,
3037  const MPI_Comm communicator,
3038  const bool exchange_data)
3039  {
3040  reinit(parallel_partitioning,
3041  parallel_partitioning,
3042  sparsity_pattern,
3043  communicator,
3044  exchange_data);
3045  }
3046 
3047 
3048 
3049  template <typename number>
3050  inline void
3051  SparseMatrix::reinit(const IndexSet &parallel_partitioning,
3052  const ::SparseMatrix<number> &sparse_matrix,
3053  const MPI_Comm communicator,
3054  const double drop_tolerance,
3055  const bool copy_values,
3056  const ::SparsityPattern *use_this_sparsity)
3057  {
3058  Epetra_Map map =
3059  parallel_partitioning.make_trilinos_map(communicator, false);
3060  reinit(parallel_partitioning,
3061  parallel_partitioning,
3062  sparse_matrix,
3063  drop_tolerance,
3064  copy_values,
3065  use_this_sparsity);
3066  }
3067 
3068 
3069 
3070  inline const Epetra_CrsMatrix &
3072  {
3073  return static_cast<const Epetra_CrsMatrix &>(*matrix);
3074  }
3075 
3076 
3077 
3078  inline const Epetra_CrsGraph &
3080  {
3081  return matrix->Graph();
3082  }
3083 
3084 
3085 
3086  inline IndexSet
3088  {
3089  return IndexSet(matrix->DomainMap());
3090  }
3091 
3092 
3093 
3094  inline IndexSet
3096  {
3097  return IndexSet(matrix->RangeMap());
3098  }
3099 
3100 
3101 
3102  inline void
3104  {
3105  // nothing to do here
3106  }
3107 
3108 
3109 
3110  inline void
3112  {
3113  // nothing to do here
3114  }
3115 
3116 
3117  namespace internal
3118  {
3119  namespace LinearOperatorImplementation
3120  {
3121  template <typename EpetraOpType>
3122  TrilinosPayload::TrilinosPayload(
3123  EpetraOpType &op,
3124  const bool supports_inverse_operations,
3125  const bool use_transpose,
3126  const MPI_Comm mpi_communicator,
3127  const IndexSet &locally_owned_domain_indices,
3128  const IndexSet &locally_owned_range_indices)
3129  : use_transpose(use_transpose)
3130  , communicator(mpi_communicator)
3131  , domain_map(
3132  locally_owned_domain_indices.make_trilinos_map(communicator.Comm()))
3133  , range_map(
3134  locally_owned_range_indices.make_trilinos_map(communicator.Comm()))
3135  {
3136  vmult = [&op](Range &tril_dst, const Domain &tril_src) {
3137  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3138  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3139  Assert(&tril_src != &tril_dst,
3142  tril_src,
3143  tril_dst,
3144  op.UseTranspose());
3145 
3146  const int ierr = op.Apply(tril_src, tril_dst);
3147  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3148  };
3149 
3150  Tvmult = [&op](Domain &tril_dst, const Range &tril_src) {
3151  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3152  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3153  Assert(&tril_src != &tril_dst,
3156  tril_src,
3157  tril_dst,
3158  !op.UseTranspose());
3159 
3160  op.SetUseTranspose(!op.UseTranspose());
3161  const int ierr = op.Apply(tril_src, tril_dst);
3162  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3163  op.SetUseTranspose(!op.UseTranspose());
3164  };
3165 
3166  if (supports_inverse_operations)
3167  {
3168  inv_vmult = [&op](Domain &tril_dst, const Range &tril_src) {
3169  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3170  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3171  Assert(
3172  &tril_src != &tril_dst,
3175  tril_src,
3176  tril_dst,
3177  !op.UseTranspose());
3178 
3179  const int ierr = op.ApplyInverse(tril_src, tril_dst);
3180  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3181  };
3182 
3183  inv_Tvmult = [&op](Range &tril_dst, const Domain &tril_src) {
3184  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3185  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3186  Assert(
3187  &tril_src != &tril_dst,
3190  tril_src,
3191  tril_dst,
3192  op.UseTranspose());
3193 
3194  op.SetUseTranspose(!op.UseTranspose());
3195  const int ierr = op.ApplyInverse(tril_src, tril_dst);
3196  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3197  op.SetUseTranspose(!op.UseTranspose());
3198  };
3199  }
3200  else
3201  {
3202  inv_vmult = [](Domain &, const Range &) {
3203  Assert(false,
3204  ExcMessage(
3205  "Uninitialized TrilinosPayload::inv_vmult called. "
3206  "The operator does not support inverse operations."));
3207  };
3208 
3209  inv_Tvmult = [](Range &, const Domain &) {
3210  Assert(false,
3211  ExcMessage(
3212  "Uninitialized TrilinosPayload::inv_Tvmult called. "
3213  "The operator does not support inverse operations."));
3214  };
3215  }
3216  }
3217 
3218 
3219  template <typename Solver, typename Preconditioner>
3220  std::enable_if_t<
3221  std::is_base_of_v<TrilinosWrappers::SolverBase, Solver> &&
3222  std::is_base_of_v<TrilinosWrappers::PreconditionBase, Preconditioner>,
3223  TrilinosPayload>
3225  Solver &solver,
3226  const Preconditioner &preconditioner) const
3227  {
3228  const auto &payload = *this;
3229 
3230  TrilinosPayload return_op(payload);
3231 
3232  // Capture by copy so the payloads are always valid
3233 
3234  return_op.inv_vmult = [payload, &solver, &preconditioner](
3235  TrilinosPayload::Domain &tril_dst,
3236  const TrilinosPayload::Range &tril_src) {
3237  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3238  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3239  Assert(&tril_src != &tril_dst,
3242  tril_src,
3243  tril_dst,
3244  !payload.UseTranspose());
3245  solver.solve(payload, tril_dst, tril_src, preconditioner);
3246  };
3247 
3248  return_op.inv_Tvmult = [payload, &solver, &preconditioner](
3249  TrilinosPayload::Range &tril_dst,
3250  const TrilinosPayload::Domain &tril_src) {
3251  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3252  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3253  Assert(&tril_src != &tril_dst,
3256  tril_src,
3257  tril_dst,
3258  payload.UseTranspose());
3259 
3260  const_cast<TrilinosPayload &>(payload).transpose();
3261  solver.solve(payload, tril_dst, tril_src, preconditioner);
3262  const_cast<TrilinosPayload &>(payload).transpose();
3263  };
3264 
3265  // If the input operator is already setup for transpose operations, then
3266  // we must do similar with its inverse.
3267  if (return_op.UseTranspose() == true)
3268  std::swap(return_op.inv_vmult, return_op.inv_Tvmult);
3269 
3270  return return_op;
3271  }
3272 
3273  template <typename Solver, typename Preconditioner>
3274  std::enable_if_t<
3275  !(std::is_base_of_v<TrilinosWrappers::SolverBase, Solver> &&
3276  std::is_base_of_v<TrilinosWrappers::PreconditionBase,
3277  Preconditioner>),
3278  TrilinosPayload>
3279  TrilinosPayload::inverse_payload(Solver &, const Preconditioner &) const
3280  {
3281  TrilinosPayload return_op(*this);
3282 
3283  return_op.inv_vmult = [](TrilinosPayload::Domain &,
3284  const TrilinosPayload::Range &) {
3285  AssertThrow(false,
3286  ExcMessage("Payload inv_vmult disabled because of "
3287  "incompatible solver/preconditioner choice."));
3288  };
3289 
3290  return_op.inv_Tvmult = [](TrilinosPayload::Range &,
3291  const TrilinosPayload::Domain &) {
3292  AssertThrow(false,
3293  ExcMessage("Payload inv_vmult disabled because of "
3294  "incompatible solver/preconditioner choice."));
3295  };
3296 
3297  return return_op;
3298  }
3299  } // namespace LinearOperatorImplementation
3300  } // namespace internal
3301 
3302  template <>
3303  void
3304  SparseMatrix::set<TrilinosScalar>(const size_type row,
3305  const size_type n_cols,
3306  const size_type *col_indices,
3307  const TrilinosScalar *values,
3308  const bool elide_zero_values);
3309 # endif // DOXYGEN
3310 
3311 } /* namespace TrilinosWrappers */
3312 
3313 
3315 
3316 
3317 # endif // DEAL_II_WITH_TRILINOS
3318 
3319 
3320 /*----------------------- trilinos_sparse_matrix.h --------------------*/
3321 
3322 #endif
3323 /*----------------------- trilinos_sparse_matrix.h --------------------*/
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
Epetra_Map make_trilinos_map(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:840
types::global_dof_index size_type
AccessorBase(SparseMatrix *matrix, const size_type row, const size_type index)
std::shared_ptr< std::vector< TrilinosScalar > > value_cache
std::shared_ptr< std::vector< size_type > > colnum_cache
const Reference & operator/=(const TrilinosScalar n) const
const Reference & operator=(const TrilinosScalar n) const
const Reference & operator*=(const TrilinosScalar n) const
const Reference & operator-=(const TrilinosScalar n) const
const Reference & operator+=(const TrilinosScalar n) const
Accessor(MatrixType *matrix, const size_type row, const size_type index)
Accessor(MatrixType *matrix, const size_type row, const size_type index)
typename Accessor< Constness >::MatrixType MatrixType
bool operator==(const Iterator< OtherConstness > &) const
bool operator!=(const Iterator< OtherConstness > &) const
bool operator<(const Iterator< OtherConstness > &) const
Iterator(const Iterator< Other > &other)
bool operator>(const Iterator< OtherConstness > &) const
Iterator(MatrixType *matrix, const size_type row, const size_type index)
const Accessor< Constness > * operator->() const
const Accessor< Constness > & operator*() const
void set(const size_type i, const size_type j, const TrilinosScalar value)
std::unique_ptr< Epetra_Map > column_space_map
SparseMatrix & operator*=(const TrilinosScalar factor)
void mmult(SparseMatrix &C, const SparseMatrix &B, const MPI::Vector &V=MPI::Vector()) const
std::unique_ptr< Epetra_Export > nonlocal_matrix_exporter
void compress(VectorOperation::values operation)
std::unique_ptr< Epetra_FECrsMatrix > matrix
const Epetra_CrsMatrix & trilinos_matrix() const
::types::global_dof_index size_type
TrilinosScalar matrix_norm_square(const MPI::Vector &v) const
IndexSet locally_owned_range_indices() const
void print(std::ostream &out, const bool write_extended_trilinos_info=false) const
void Tmmult(SparseMatrix &C, const SparseMatrix &B, const MPI::Vector &V=MPI::Vector()) const
std::enable_if_t< std::is_same_v< typename VectorType::value_type, TrilinosScalar > > Tvmult(VectorType &dst, const VectorType &src) const
void clear_row(const size_type row, const TrilinosScalar new_diag_value=0)
const_iterator begin() const
std::pair< size_type, size_type > local_range() const
const Epetra_CrsGraph & trilinos_sparsity_pattern() const
void reinit(const SparsityPatternType &sparsity_pattern)
std::enable_if_t< !std::is_same_v< SparsityPatternType, ::SparseMatrix< double > > > reinit(const IndexSet &row_parallel_partitioning, const IndexSet &col_parallel_partitioning, const SparsityPatternType &sparsity_pattern, const MPI_Comm communicator=MPI_COMM_WORLD, const bool exchange_data=false)
void vmult_add(VectorType &dst, const VectorType &src) const
const_iterator end(const size_type r) const
SparseMatrix & operator/=(const TrilinosScalar factor)
void Tvmult_add(VectorType &dst, const VectorType &src) const
void clear_rows(const std::vector< size_type > &rows, const TrilinosScalar new_diag_value=0)
TrilinosScalar el(const size_type i, const size_type j) const
IndexSet locally_owned_domain_indices() const
TrilinosScalar residual(MPI::Vector &dst, const MPI::Vector &x, const MPI::Vector &b) const
bool in_local_range(const size_type index) const
void copy_from(const SparseMatrix &source)
unsigned int row_length(const size_type row) const
std::uint64_t n_nonzero_elements() const
std::enable_if_t< std::is_same_v< typename VectorType::value_type, TrilinosScalar > > vmult(VectorType &dst, const VectorType &src) const
const_iterator begin(const size_type r) const
SparseMatrix(const SparseMatrix &)=delete
iterator end(const size_type r)
TrilinosScalar diag_element(const size_type i) const
void add(const size_type i, const size_type j, const TrilinosScalar value)
unsigned int local_size() const
iterator begin(const size_type r)
TrilinosScalar operator()(const size_type i, const size_type j) const
void reinit(const IndexSet &parallel_partitioning, const ::SparseMatrix< number > &dealii_sparse_matrix, const MPI_Comm communicator=MPI_COMM_WORLD, const double drop_tolerance=1e-13, const bool copy_values=true, const ::SparsityPattern *use_this_sparsity=nullptr)
virtual ~SparseMatrix() override=default
void set(const size_type row, const size_type n_cols, const size_type *col_indices, const Number *values, const bool elide_zero_values=false)
TrilinosScalar matrix_scalar_product(const MPI::Vector &u, const MPI::Vector &v) const
const_iterator end() const
void set(const std::vector< size_type > &indices, const FullMatrix< TrilinosScalar > &full_matrix, const bool elide_zero_values=false)
SparseMatrix & operator=(const SparseMatrix &)=delete
void reinit(const IndexSet &parallel_partitioning, const SparsityPatternType &sparsity_pattern, const MPI_Comm communicator=MPI_COMM_WORLD, const bool exchange_data=false)
std::unique_ptr< Epetra_CrsMatrix > nonlocal_matrix
std::function< void(VectorType &, const VectorType &)> inv_Tvmult
virtual int ApplyInverse(const VectorType &Y, VectorType &X) const override
virtual int Apply(const VectorType &X, VectorType &Y) const override
std::function< void(VectorType &, const VectorType &)> inv_vmult
TrilinosPayload(EpetraOpType &op, const bool supports_inverse_operations, const bool use_transpose, const MPI_Comm mpi_communicator, const IndexSet &locally_owned_domain_indices, const IndexSet &locally_owned_range_indices)
std::enable_if_t< std::is_base_of_v< TrilinosWrappers::SolverBase, Solver > &&std::is_base_of_v< TrilinosWrappers::PreconditionBase, Preconditioner >, TrilinosPayload > inverse_payload(Solver &, const Preconditioner &) const
std::enable_if_t< !(std::is_base_of_v< TrilinosWrappers::SolverBase, Solver > &&std::is_base_of_v< TrilinosWrappers::PreconditionBase, Preconditioner >), TrilinosPayload > inverse_payload(Solver &, const Preconditioner &) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
__global__ void set(Number *val, const Number s, const size_type N)
#define DeclException0(Exception0)
Definition: exceptions.h:467
static ::ExceptionBase & ExcInvalidIndexWithinRow(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcAccessToNonlocalRow(std::size_t arg1)
static ::ExceptionBase & ExcInternalError()
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:581
static ::ExceptionBase & ExcInvalidIndex(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcSourceEqualsDestination()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcIteratorPastEnd()
static ::ExceptionBase & ExcAccessToNonlocalRow(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcBeyondEndOfMatrix()
#define Assert(cond, exc)
Definition: exceptions.h:1616
#define AssertIsFinite(number)
Definition: exceptions.h:1884
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:535
static ::ExceptionBase & ExcAccessToNonLocalElement(size_type arg1, size_type arg2, size_type arg3, size_type arg4)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1857
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:490
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:558
static ::ExceptionBase & ExcMatrixNotCompressed()
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:512
static ::ExceptionBase & ExcAccessToNonPresentElement(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcNotQuadratic()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
Expression operator>(const Expression &lhs, const Expression &rhs)
@ matrix
Contents is actually a matrix.
static const char V
SparseMatrix< double > SparseMatrix
types::global_dof_index size_type
Definition: cuda_kernels.h:45
void swap(MemorySpaceData< T, MemorySpace > &u, MemorySpaceData< T, MemorySpace > &v)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
TrilinosPayload operator+(const TrilinosPayload &first_op, const TrilinosPayload &second_op)
TrilinosPayload operator*(const TrilinosPayload &first_op, const TrilinosPayload &second_op)
void check_vector_map_equality(const Epetra_CrsMatrix &mtrx, const Epetra_MultiVector &src, const Epetra_MultiVector &dst, const bool transpose)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
TrilinosWrappers::types::int64_type n_global_elements(const Epetra_BlockMap &map)
void copy(const T *begin, const T *end, U *dest)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
unsigned int global_dof_index
Definition: types.h:82
typename ::TrilinosWrappers::SparseMatrixIterators::Iterator< Constness >::value_type value_type
typename ::TrilinosWrappers::SparseMatrixIterators::Iterator< Constness >::difference_type difference_type
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)
double TrilinosScalar
Definition: types.h:175