Reference documentation for deal.II version GIT 8e09676776 2023-03-27 21:15: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 - 2022 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  bool
420 
424  bool
426 
432  bool
433  operator<(const Iterator<Constness> &) const;
434 
438  bool
440 
445  size_type,
446  size_type,
447  << "Attempt to access element " << arg2 << " of row "
448  << arg1 << " which doesn't have that many elements.");
449 
450  private:
455 
456  template <bool Other>
457  friend class Iterator;
458  };
459 
460  } // namespace SparseMatrixIterators
461 } // namespace TrilinosWrappers
462 
464 
465 namespace std
466 {
467  template <bool Constness>
468  struct iterator_traits<
470  {
471  using iterator_category = forward_iterator_tag;
472  using value_type =
473  typename ::TrilinosWrappers::SparseMatrixIterators::Iterator<
474  Constness>::value_type;
476  typename ::TrilinosWrappers::SparseMatrixIterators::Iterator<
477  Constness>::difference_type;
478  };
479 } // namespace std
480 
482 
483 
484 namespace TrilinosWrappers
485 {
546  class SparseMatrix : public Subscriptor
547  {
548  public:
553 
558  std::size_t,
559  << "You tried to access row " << arg1
560  << " of a non-contiguous locally owned row set."
561  << " The row " << arg1
562  << " is not stored locally and can't be accessed.");
563 
571  struct Traits
572  {
577  static const bool zero_addition_can_be_elided = true;
578  };
579 
584 
589 
594 
602  SparseMatrix();
603 
611  SparseMatrix(const size_type m,
612  const size_type n,
613  const unsigned int n_max_entries_per_row);
614 
622  SparseMatrix(const size_type m,
623  const size_type n,
624  const std::vector<unsigned int> &n_entries_per_row);
625 
629  SparseMatrix(const SparsityPattern &InputSparsityPattern);
630 
635  SparseMatrix(SparseMatrix &&other) noexcept;
636 
640  SparseMatrix(const SparseMatrix &) = delete;
641 
645  SparseMatrix &
646  operator=(const SparseMatrix &) = delete;
647 
651  virtual ~SparseMatrix() override = default;
652 
668  template <typename SparsityPatternType>
669  void
670  reinit(const SparsityPatternType &sparsity_pattern);
671 
684  void
685  reinit(const SparsityPattern &sparsity_pattern);
686 
695  void
696  reinit(const SparseMatrix &sparse_matrix);
697 
718  template <typename number>
719  void
720  reinit(const ::SparseMatrix<number> &dealii_sparse_matrix,
721  const double drop_tolerance = 1e-13,
722  const bool copy_values = true,
723  const ::SparsityPattern * use_this_sparsity = nullptr);
724 
730  void
731  reinit(const Epetra_CrsMatrix &input_matrix, const bool copy_values = true);
750  SparseMatrix(const IndexSet & parallel_partitioning,
751  const MPI_Comm & communicator = MPI_COMM_WORLD,
752  const unsigned int n_max_entries_per_row = 0);
753 
761  SparseMatrix(const IndexSet & parallel_partitioning,
762  const MPI_Comm & communicator,
763  const std::vector<unsigned int> &n_entries_per_row);
764 
779  SparseMatrix(const IndexSet &row_parallel_partitioning,
780  const IndexSet &col_parallel_partitioning,
781  const MPI_Comm &communicator = MPI_COMM_WORLD,
782  const size_type n_max_entries_per_row = 0);
783 
798  SparseMatrix(const IndexSet & row_parallel_partitioning,
799  const IndexSet & col_parallel_partitioning,
800  const MPI_Comm & communicator,
801  const std::vector<unsigned int> &n_entries_per_row);
802 
823  template <typename SparsityPatternType>
824  void
825  reinit(const IndexSet & parallel_partitioning,
826  const SparsityPatternType &sparsity_pattern,
827  const MPI_Comm & communicator = MPI_COMM_WORLD,
828  const bool exchange_data = false);
829 
842  template <typename SparsityPatternType>
843  std::enable_if_t<
844  !std::is_same<SparsityPatternType, ::SparseMatrix<double>>::value>
845  reinit(const IndexSet & row_parallel_partitioning,
846  const IndexSet & col_parallel_partitioning,
847  const SparsityPatternType &sparsity_pattern,
848  const MPI_Comm & communicator = MPI_COMM_WORLD,
849  const bool exchange_data = false);
850 
867  template <typename number>
868  void
869  reinit(const IndexSet & parallel_partitioning,
870  const ::SparseMatrix<number> &dealii_sparse_matrix,
871  const MPI_Comm & communicator = MPI_COMM_WORLD,
872  const double drop_tolerance = 1e-13,
873  const bool copy_values = true,
874  const ::SparsityPattern * use_this_sparsity = nullptr);
875 
889  template <typename number>
890  void
891  reinit(const IndexSet & row_parallel_partitioning,
892  const IndexSet & col_parallel_partitioning,
893  const ::SparseMatrix<number> &dealii_sparse_matrix,
894  const MPI_Comm & communicator = MPI_COMM_WORLD,
895  const double drop_tolerance = 1e-13,
896  const bool copy_values = true,
897  const ::SparsityPattern * use_this_sparsity = nullptr);
907  size_type
908  m() const;
909 
913  size_type
914  n() const;
915 
924  unsigned int
925  local_size() const;
926 
935  std::pair<size_type, size_type>
936  local_range() const;
937 
942  bool
943  in_local_range(const size_type index) const;
944 
949  std::uint64_t
951 
955  unsigned int
956  row_length(const size_type row) const;
957 
964  bool
965  is_compressed() const;
966 
972  size_type
973  memory_consumption() const;
974 
978  MPI_Comm
979  get_mpi_communicator() const;
980 
996  SparseMatrix &
997  operator=(const double d);
998 
1006  void
1007  clear();
1008 
1036  void
1037  compress(VectorOperation::values operation);
1038 
1060  void
1061  set(const size_type i, const size_type j, const TrilinosScalar value);
1062 
1095  void
1096  set(const std::vector<size_type> & indices,
1097  const FullMatrix<TrilinosScalar> &full_matrix,
1098  const bool elide_zero_values = false);
1099 
1105  void
1106  set(const std::vector<size_type> & row_indices,
1107  const std::vector<size_type> & col_indices,
1108  const FullMatrix<TrilinosScalar> &full_matrix,
1109  const bool elide_zero_values = false);
1110 
1138  void
1139  set(const size_type row,
1140  const std::vector<size_type> & col_indices,
1141  const std::vector<TrilinosScalar> &values,
1142  const bool elide_zero_values = false);
1143 
1171  template <typename Number>
1172  void
1173  set(const size_type row,
1174  const size_type n_cols,
1175  const size_type *col_indices,
1176  const Number * values,
1177  const bool elide_zero_values = false);
1178 
1188  void
1189  add(const size_type i, const size_type j, const TrilinosScalar value);
1190 
1209  void
1210  add(const std::vector<size_type> & indices,
1211  const FullMatrix<TrilinosScalar> &full_matrix,
1212  const bool elide_zero_values = true);
1213 
1219  void
1220  add(const std::vector<size_type> & row_indices,
1221  const std::vector<size_type> & col_indices,
1222  const FullMatrix<TrilinosScalar> &full_matrix,
1223  const bool elide_zero_values = true);
1224 
1238  void
1239  add(const size_type row,
1240  const std::vector<size_type> & col_indices,
1241  const std::vector<TrilinosScalar> &values,
1242  const bool elide_zero_values = true);
1243 
1257  void
1258  add(const size_type row,
1259  const size_type n_cols,
1260  const size_type * col_indices,
1261  const TrilinosScalar *values,
1262  const bool elide_zero_values = true,
1263  const bool col_indices_are_sorted = false);
1264 
1268  SparseMatrix &
1269  operator*=(const TrilinosScalar factor);
1270 
1274  SparseMatrix &
1275  operator/=(const TrilinosScalar factor);
1276 
1280  void
1281  copy_from(const SparseMatrix &source);
1282 
1290  void
1291  add(const TrilinosScalar factor, const SparseMatrix &matrix);
1292 
1319  void
1320  clear_row(const size_type row, const TrilinosScalar new_diag_value = 0);
1321 
1342  void
1343  clear_rows(const std::vector<size_type> &rows,
1344  const TrilinosScalar new_diag_value = 0);
1345 
1355  void
1356  transpose();
1357 
1373  operator()(const size_type i, const size_type j) const;
1374 
1392  el(const size_type i, const size_type j) const;
1393 
1401  diag_element(const size_type i) const;
1402 
1439  template <typename VectorType>
1440  std::enable_if_t<
1441  std::is_same<typename VectorType::value_type, TrilinosScalar>::value>
1442  vmult(VectorType &dst, const VectorType &src) const;
1443 
1450  template <typename VectorType>
1451  std::enable_if_t<
1452  !std::is_same<typename VectorType::value_type, TrilinosScalar>::value>
1453  vmult(VectorType &dst, const VectorType &src) const;
1454 
1469  template <typename VectorType>
1470  std::enable_if_t<
1471  std::is_same<typename VectorType::value_type, TrilinosScalar>::value>
1472  Tvmult(VectorType &dst, const VectorType &src) const;
1473 
1480  template <typename VectorType>
1481  std::enable_if_t<
1482  !std::is_same<typename VectorType::value_type, TrilinosScalar>::value>
1483  Tvmult(VectorType &dst, const VectorType &src) const;
1484 
1494  template <typename VectorType>
1495  void
1496  vmult_add(VectorType &dst, const VectorType &src) const;
1497 
1508  template <typename VectorType>
1509  void
1510  Tvmult_add(VectorType &dst, const VectorType &src) const;
1511 
1534  matrix_norm_square(const MPI::Vector &v) const;
1535 
1556  matrix_scalar_product(const MPI::Vector &u, const MPI::Vector &v) const;
1557 
1575  residual(MPI::Vector & dst,
1576  const MPI::Vector &x,
1577  const MPI::Vector &b) const;
1578 
1593  void
1594  mmult(SparseMatrix & C,
1595  const SparseMatrix &B,
1596  const MPI::Vector & V = MPI::Vector()) const;
1597 
1598 
1615  void
1616  Tmmult(SparseMatrix & C,
1617  const SparseMatrix &B,
1618  const MPI::Vector & V = MPI::Vector()) const;
1619 
1634  l1_norm() const;
1635 
1645  linfty_norm() const;
1646 
1652  frobenius_norm() const;
1653 
1664  const Epetra_CrsMatrix &
1666 
1671  const Epetra_CrsGraph &
1673 
1685  IndexSet
1687 
1693  IndexSet
1695 
1722  begin() const;
1723 
1727  iterator
1729 
1735  end() const;
1736 
1740  iterator
1741  end();
1742 
1772  begin(const size_type r) const;
1773 
1777  iterator
1778  begin(const size_type r);
1779 
1790  end(const size_type r) const;
1791 
1795  iterator
1796  end(const size_type r);
1797 
1809  void
1810  write_ascii();
1811 
1819  void
1820  print(std::ostream &out,
1821  const bool write_extended_trilinos_info = false) const;
1822 
1833  int,
1834  << "An error with error number " << arg1
1835  << " occurred while calling a Trilinos function. "
1836  "\n\n"
1837  "For historical reasons, many Trilinos functions express "
1838  "errors by returning specific integer values to indicate "
1839  "certain errors. Unfortunately, different Trilinos functions "
1840  "often use the same integer values for different kinds of "
1841  "errors, and in most cases it is also not documented what "
1842  "each error code actually means. As a consequence, it is often "
1843  "difficult to say what a particular error (in this case, "
1844  "the error with integer code '"
1845  << arg1
1846  << "') represents and how one should fix a code to avoid it. "
1847  "The best one can often do is to look up the call stack to "
1848  "see which deal.II function generated the error, and which "
1849  "Trilinos function the error code had originated from; "
1850  "then look up the Trilinos source code of that function (for "
1851  "example on github) to see what code path set that error "
1852  "code. Short of going through all of that, the only other "
1853  "option is to guess the cause of the error from "
1854  "the context in which the error appeared.");
1855 
1856 
1861  size_type,
1862  size_type,
1863  << "The entry with index <" << arg1 << ',' << arg2
1864  << "> does not exist.");
1865 
1870  "You are attempting an operation on two vectors that "
1871  "are the same object, but the operation requires that the "
1872  "two objects are in fact different.");
1873 
1878 
1883  size_type,
1884  size_type,
1885  size_type,
1886  size_type,
1887  << "You tried to access element (" << arg1 << '/' << arg2
1888  << ')'
1889  << " of a distributed matrix, but only rows in range ["
1890  << arg3 << ',' << arg4
1891  << "] are stored locally and can be accessed.");
1892 
1897  size_type,
1898  size_type,
1899  << "You tried to access element (" << arg1 << '/' << arg2
1900  << ')' << " of a sparse matrix, but it appears to not"
1901  << " exist in the Trilinos sparsity pattern.");
1906  protected:
1917  void
1919 
1925  void
1927 
1928 
1929 
1930  private:
1935  std::unique_ptr<Epetra_Map> column_space_map;
1936 
1942  std::unique_ptr<Epetra_FECrsMatrix> matrix;
1943 
1949  std::unique_ptr<Epetra_CrsMatrix> nonlocal_matrix;
1950 
1954  std::unique_ptr<Epetra_Export> nonlocal_matrix_exporter;
1955 
1967  Epetra_CombineMode last_action;
1968 
1974 
1975  // To allow calling protected prepare_add() and prepare_set().
1976  friend class BlockMatrixBase<SparseMatrix>;
1977  };
1978 
1979 
1980 
1981  // forwards declarations
1982  class SolverBase;
1983  class PreconditionBase;
1984 
1985  namespace internal
1986  {
1987  inline void
1988  check_vector_map_equality(const Epetra_CrsMatrix & mtrx,
1989  const Epetra_MultiVector &src,
1990  const Epetra_MultiVector &dst,
1991  const bool transpose)
1992  {
1993  if (transpose == false)
1994  {
1995  Assert(src.Map().SameAs(mtrx.DomainMap()) == true,
1996  ExcMessage(
1997  "Column map of matrix does not fit with vector map!"));
1998  Assert(dst.Map().SameAs(mtrx.RangeMap()) == true,
1999  ExcMessage("Row map of matrix does not fit with vector map!"));
2000  }
2001  else
2002  {
2003  Assert(src.Map().SameAs(mtrx.RangeMap()) == true,
2004  ExcMessage(
2005  "Column map of matrix does not fit with vector map!"));
2006  Assert(dst.Map().SameAs(mtrx.DomainMap()) == true,
2007  ExcMessage("Row map of matrix does not fit with vector map!"));
2008  }
2009  (void)mtrx; // removes -Wunused-variable in optimized mode
2010  (void)src;
2011  (void)dst;
2012  }
2013 
2014  inline void
2016  const Epetra_MultiVector &src,
2017  const Epetra_MultiVector &dst,
2018  const bool transpose)
2019  {
2020  if (transpose == false)
2021  {
2022  Assert(src.Map().SameAs(op.OperatorDomainMap()) == true,
2023  ExcMessage(
2024  "Column map of operator does not fit with vector map!"));
2025  Assert(dst.Map().SameAs(op.OperatorRangeMap()) == true,
2026  ExcMessage(
2027  "Row map of operator does not fit with vector map!"));
2028  }
2029  else
2030  {
2031  Assert(src.Map().SameAs(op.OperatorRangeMap()) == true,
2032  ExcMessage(
2033  "Column map of operator does not fit with vector map!"));
2034  Assert(dst.Map().SameAs(op.OperatorDomainMap()) == true,
2035  ExcMessage(
2036  "Row map of operator does not fit with vector map!"));
2037  }
2038  (void)op; // removes -Wunused-variable in optimized mode
2039  (void)src;
2040  (void)dst;
2041  }
2042 
2043  namespace LinearOperatorImplementation
2044  {
2065  {
2066  public:
2070  using VectorType = Epetra_MultiVector;
2071 
2076 
2081 
2094  TrilinosPayload();
2095 
2099  TrilinosPayload(const TrilinosWrappers::SparseMatrix &matrix_exemplar,
2101 
2105  TrilinosPayload(const TrilinosPayload & payload_exemplar,
2107 
2112  const TrilinosWrappers::SparseMatrix & matrix_exemplar,
2113  const TrilinosWrappers::PreconditionBase &preconditioner);
2114 
2119  const TrilinosWrappers::PreconditionBase &preconditioner_exemplar,
2120  const TrilinosWrappers::PreconditionBase &preconditioner);
2121 
2126  const TrilinosPayload & payload_exemplar,
2127  const TrilinosWrappers::PreconditionBase &preconditioner);
2128 
2132  TrilinosPayload(const TrilinosPayload &payload);
2133 
2141  TrilinosPayload(const TrilinosPayload &first_op,
2142  const TrilinosPayload &second_op);
2143 
2147  virtual ~TrilinosPayload() override = default;
2148 
2153  identity_payload() const;
2154 
2159  null_payload() const;
2160 
2165  transpose_payload() const;
2166 
2183  template <typename Solver, typename Preconditioner>
2184  std::enable_if_t<
2185  std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
2186  std::is_base_of<TrilinosWrappers::PreconditionBase,
2187  Preconditioner>::value,
2189  inverse_payload(Solver &, const Preconditioner &) const;
2190 
2208  template <typename Solver, typename Preconditioner>
2209  std::enable_if_t<
2210  !(std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
2211  std::is_base_of<TrilinosWrappers::PreconditionBase,
2212  Preconditioner>::value),
2214  inverse_payload(Solver &, const Preconditioner &) const;
2215 
2228  IndexSet
2230 
2236  IndexSet
2238 
2242  MPI_Comm
2243  get_mpi_communicator() const;
2244 
2251  void
2252  transpose();
2253 
2261  std::function<void(VectorType &, const VectorType &)> vmult;
2262 
2270  std::function<void(VectorType &, const VectorType &)> Tvmult;
2271 
2280  std::function<void(VectorType &, const VectorType &)> inv_vmult;
2281 
2290  std::function<void(VectorType &, const VectorType &)> inv_Tvmult;
2291 
2305  virtual bool
2306  UseTranspose() const override;
2307 
2323  virtual int
2324  SetUseTranspose(bool UseTranspose) override;
2325 
2337  virtual int
2338  Apply(const VectorType &X, VectorType &Y) const override;
2339 
2358  virtual int
2359  ApplyInverse(const VectorType &Y, VectorType &X) const override;
2373  virtual const char *
2374  Label() const override;
2375 
2383  virtual const Epetra_Comm &
2384  Comm() const override;
2385 
2393  virtual const Epetra_Map &
2394  OperatorDomainMap() const override;
2395 
2404  virtual const Epetra_Map &
2405  OperatorRangeMap() const override;
2408  private:
2418  template <typename EpetraOpType>
2419  TrilinosPayload(EpetraOpType & op,
2420  const bool supports_inverse_operations,
2421  const bool use_transpose,
2422  const MPI_Comm &mpi_communicator,
2425 
2431 
2436  Epetra_MpiComm communicator;
2437 
2442  Epetra_Map domain_map;
2443 
2448  Epetra_Map range_map;
2449 
2458  virtual bool
2459  HasNormInf() const override;
2460 
2468  virtual double
2469  NormInf() const override;
2470  };
2471 
2477  operator+(const TrilinosPayload &first_op,
2478  const TrilinosPayload &second_op);
2479 
2485  operator*(const TrilinosPayload &first_op,
2486  const TrilinosPayload &second_op);
2487 
2488  } // namespace LinearOperatorImplementation
2489  } /* namespace internal */
2490 
2491 
2492 
2493  // ----------------------- inline and template functions --------------------
2494 
2495 # ifndef DOXYGEN
2496 
2497  namespace SparseMatrixIterators
2498  {
2500  size_type row,
2501  size_type index)
2502  : matrix(matrix)
2503  , a_row(row)
2504  , a_index(index)
2505  {
2506  visit_present_row();
2507  }
2508 
2509 
2511  AccessorBase::row() const
2512  {
2513  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2514  return a_row;
2515  }
2516 
2517 
2519  AccessorBase::column() const
2520  {
2521  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2522  return (*colnum_cache)[a_index];
2523  }
2524 
2525 
2527  AccessorBase::index() const
2528  {
2529  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2530  return a_index;
2531  }
2532 
2533 
2534  inline Accessor<true>::Accessor(MatrixType * matrix,
2535  const size_type row,
2536  const size_type index)
2537  : AccessorBase(const_cast<SparseMatrix *>(matrix), row, index)
2538  {}
2539 
2540 
2541  template <bool Other>
2542  inline Accessor<true>::Accessor(const Accessor<Other> &other)
2543  : AccessorBase(other)
2544  {}
2545 
2546 
2547  inline TrilinosScalar
2548  Accessor<true>::value() const
2549  {
2550  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2551  return (*value_cache)[a_index];
2552  }
2553 
2554 
2555  inline Accessor<false>::Reference::Reference(const Accessor<false> &acc)
2556  : accessor(const_cast<Accessor<false> &>(acc))
2557  {}
2558 
2559 
2560  inline Accessor<false>::Reference::operator TrilinosScalar() const
2561  {
2562  return (*accessor.value_cache)[accessor.a_index];
2563  }
2564 
2565  inline const Accessor<false>::Reference &
2566  Accessor<false>::Reference::operator=(const TrilinosScalar n) const
2567  {
2568  (*accessor.value_cache)[accessor.a_index] = n;
2569  accessor.matrix->set(accessor.row(),
2570  accessor.column(),
2571  static_cast<TrilinosScalar>(*this));
2572  return *this;
2573  }
2574 
2575 
2576  inline const Accessor<false>::Reference &
2577  Accessor<false>::Reference::operator+=(const TrilinosScalar n) const
2578  {
2579  (*accessor.value_cache)[accessor.a_index] += n;
2580  accessor.matrix->set(accessor.row(),
2581  accessor.column(),
2582  static_cast<TrilinosScalar>(*this));
2583  return *this;
2584  }
2585 
2586 
2587  inline const Accessor<false>::Reference &
2588  Accessor<false>::Reference::operator-=(const TrilinosScalar n) const
2589  {
2590  (*accessor.value_cache)[accessor.a_index] -= n;
2591  accessor.matrix->set(accessor.row(),
2592  accessor.column(),
2593  static_cast<TrilinosScalar>(*this));
2594  return *this;
2595  }
2596 
2597 
2598  inline const Accessor<false>::Reference &
2599  Accessor<false>::Reference::operator*=(const TrilinosScalar n) const
2600  {
2601  (*accessor.value_cache)[accessor.a_index] *= n;
2602  accessor.matrix->set(accessor.row(),
2603  accessor.column(),
2604  static_cast<TrilinosScalar>(*this));
2605  return *this;
2606  }
2607 
2608 
2609  inline const Accessor<false>::Reference &
2610  Accessor<false>::Reference::operator/=(const TrilinosScalar n) const
2611  {
2612  (*accessor.value_cache)[accessor.a_index] /= n;
2613  accessor.matrix->set(accessor.row(),
2614  accessor.column(),
2615  static_cast<TrilinosScalar>(*this));
2616  return *this;
2617  }
2618 
2619 
2620  inline Accessor<false>::Accessor(MatrixType * matrix,
2621  const size_type row,
2622  const size_type index)
2623  : AccessorBase(matrix, row, index)
2624  {}
2625 
2626 
2627  inline Accessor<false>::Reference
2628  Accessor<false>::value() const
2629  {
2630  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2631  return {*this};
2632  }
2633 
2634 
2635 
2636  template <bool Constness>
2637  inline Iterator<Constness>::Iterator(MatrixType * matrix,
2638  const size_type row,
2639  const size_type index)
2640  : accessor(matrix, row, index)
2641  {}
2642 
2643 
2644  template <bool Constness>
2645  template <bool Other>
2646  inline Iterator<Constness>::Iterator(const Iterator<Other> &other)
2647  : accessor(other.accessor)
2648  {}
2649 
2650 
2651  template <bool Constness>
2652  inline Iterator<Constness> &
2654  {
2655  Assert(accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
2656 
2657  ++accessor.a_index;
2658 
2659  // If at end of line: do one
2660  // step, then cycle until we
2661  // find a row with a nonzero
2662  // number of entries.
2663  if (accessor.a_index >= accessor.colnum_cache->size())
2664  {
2665  accessor.a_index = 0;
2666  ++accessor.a_row;
2667 
2668  while ((accessor.a_row < accessor.matrix->m()) &&
2669  ((accessor.matrix->in_local_range(accessor.a_row) == false) ||
2670  (accessor.matrix->row_length(accessor.a_row) == 0)))
2671  ++accessor.a_row;
2672 
2673  accessor.visit_present_row();
2674  }
2675  return *this;
2676  }
2677 
2678 
2679  template <bool Constness>
2680  inline Iterator<Constness>
2682  {
2683  const Iterator<Constness> old_state = *this;
2684  ++(*this);
2685  return old_state;
2686  }
2687 
2688 
2689 
2690  template <bool Constness>
2691  inline const Accessor<Constness> &
2693  {
2694  return accessor;
2695  }
2696 
2697 
2698 
2699  template <bool Constness>
2700  inline const Accessor<Constness> *
2701  Iterator<Constness>::operator->() const
2702  {
2703  return &accessor;
2704  }
2705 
2706 
2707 
2708  template <bool Constness>
2709  inline bool
2710  Iterator<Constness>::operator==(const Iterator<Constness> &other) const
2711  {
2712  return (accessor.a_row == other.accessor.a_row &&
2713  accessor.a_index == other.accessor.a_index);
2714  }
2715 
2716 
2717 
2718  template <bool Constness>
2719  inline bool
2720  Iterator<Constness>::operator!=(const Iterator<Constness> &other) const
2721  {
2722  return !(*this == other);
2723  }
2724 
2725 
2726 
2727  template <bool Constness>
2728  inline bool
2729  Iterator<Constness>::operator<(const Iterator<Constness> &other) const
2730  {
2731  return (accessor.row() < other.accessor.row() ||
2732  (accessor.row() == other.accessor.row() &&
2733  accessor.index() < other.accessor.index()));
2734  }
2735 
2736 
2737  template <bool Constness>
2738  inline bool
2739  Iterator<Constness>::operator>(const Iterator<Constness> &other) const
2740  {
2741  return (other < *this);
2742  }
2743 
2744  } // namespace SparseMatrixIterators
2745 
2746 
2747 
2749  SparseMatrix::begin() const
2750  {
2751  return begin(0);
2752  }
2753 
2754 
2755 
2757  SparseMatrix::end() const
2758  {
2759  return const_iterator(this, m(), 0);
2760  }
2761 
2762 
2763 
2765  SparseMatrix::begin(const size_type r) const
2766  {
2767  AssertIndexRange(r, m());
2768  if (in_local_range(r) && (row_length(r) > 0))
2769  return const_iterator(this, r, 0);
2770  else
2771  return end(r);
2772  }
2773 
2774 
2775 
2777  SparseMatrix::end(const size_type r) const
2778  {
2779  AssertIndexRange(r, m());
2780 
2781  // place the iterator on the first entry
2782  // past this line, or at the end of the
2783  // matrix
2784  for (size_type i = r + 1; i < m(); ++i)
2785  if (in_local_range(i) && (row_length(i) > 0))
2786  return const_iterator(this, i, 0);
2787 
2788  // if there is no such line, then take the
2789  // end iterator of the matrix
2790  return end();
2791  }
2792 
2793 
2794 
2795  inline SparseMatrix::iterator
2797  {
2798  return begin(0);
2799  }
2800 
2801 
2802 
2803  inline SparseMatrix::iterator
2805  {
2806  return iterator(this, m(), 0);
2807  }
2808 
2809 
2810 
2811  inline SparseMatrix::iterator
2813  {
2814  AssertIndexRange(r, m());
2815  if (in_local_range(r) && (row_length(r) > 0))
2816  return iterator(this, r, 0);
2817  else
2818  return end(r);
2819  }
2820 
2821 
2822 
2823  inline SparseMatrix::iterator
2824  SparseMatrix::end(const size_type r)
2825  {
2826  AssertIndexRange(r, m());
2827 
2828  // place the iterator on the first entry
2829  // past this line, or at the end of the
2830  // matrix
2831  for (size_type i = r + 1; i < m(); ++i)
2832  if (in_local_range(i) && (row_length(i) > 0))
2833  return iterator(this, i, 0);
2834 
2835  // if there is no such line, then take the
2836  // end iterator of the matrix
2837  return end();
2838  }
2839 
2840 
2841 
2842  inline bool
2843  SparseMatrix::in_local_range(const size_type index) const
2844  {
2846 # ifndef DEAL_II_WITH_64BIT_INDICES
2847  begin = matrix->RowMap().MinMyGID();
2848  end = matrix->RowMap().MaxMyGID() + 1;
2849 # else
2850  begin = matrix->RowMap().MinMyGID64();
2851  end = matrix->RowMap().MaxMyGID64() + 1;
2852 # endif
2853 
2854  return ((index >= static_cast<size_type>(begin)) &&
2855  (index < static_cast<size_type>(end)));
2856  }
2857 
2858 
2859 
2860  inline bool
2862  {
2863  return compressed;
2864  }
2865 
2866 
2867 
2868  // Inline the set() and add() functions, since they will be called
2869  // frequently, and the compiler can optimize away some unnecessary loops
2870  // when the sizes are given at compile time.
2871  template <>
2872  void
2873  SparseMatrix::set<TrilinosScalar>(const size_type row,
2874  const size_type n_cols,
2875  const size_type * col_indices,
2876  const TrilinosScalar *values,
2877  const bool elide_zero_values);
2878 
2879 
2880 
2881  template <typename Number>
2882  void
2883  SparseMatrix::set(const size_type row,
2884  const size_type n_cols,
2885  const size_type *col_indices,
2886  const Number * values,
2887  const bool elide_zero_values)
2888  {
2889  std::vector<TrilinosScalar> trilinos_values(n_cols);
2890  std::copy(values, values + n_cols, trilinos_values.begin());
2891  this->set(
2892  row, n_cols, col_indices, trilinos_values.data(), elide_zero_values);
2893  }
2894 
2895 
2896 
2897  inline void
2898  SparseMatrix::set(const size_type i,
2899  const size_type j,
2900  const TrilinosScalar value)
2901  {
2902  AssertIsFinite(value);
2903 
2904  set(i, 1, &j, &value, false);
2905  }
2906 
2907 
2908 
2909  inline void
2910  SparseMatrix::set(const std::vector<size_type> & indices,
2912  const bool elide_zero_values)
2913  {
2914  Assert(indices.size() == values.m(),
2915  ExcDimensionMismatch(indices.size(), values.m()));
2916  Assert(values.m() == values.n(), ExcNotQuadratic());
2917 
2918  for (size_type i = 0; i < indices.size(); ++i)
2919  set(indices[i],
2920  indices.size(),
2921  indices.data(),
2922  &values(i, 0),
2923  elide_zero_values);
2924  }
2925 
2926 
2927 
2928  inline void
2929  SparseMatrix::add(const size_type i,
2930  const size_type j,
2931  const TrilinosScalar value)
2932  {
2933  AssertIsFinite(value);
2934 
2935  if (value == 0)
2936  {
2937  // we have to check after Insert/Add in any case to be consistent
2938  // with the MPI communication model, but we can save some
2939  // work if the addend is zero. However, these actions are done in case
2940  // we pass on to the other function.
2941 
2942  // TODO: fix this (do not run compress here, but fail)
2943  if (last_action == Insert)
2944  {
2945  const int ierr = matrix->GlobalAssemble(*column_space_map,
2946  matrix->RowMap(),
2947  false);
2948 
2949  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2950  }
2951 
2952  last_action = Add;
2953 
2954  return;
2955  }
2956  else
2957  add(i, 1, &j, &value, false);
2958  }
2959 
2960 
2961 
2962  // inline "simple" functions that are called frequently and do only involve
2963  // a call to some Trilinos function.
2965  SparseMatrix::m() const
2966  {
2967 # ifndef DEAL_II_WITH_64BIT_INDICES
2968  return matrix->NumGlobalRows();
2969 # else
2970  return matrix->NumGlobalRows64();
2971 # endif
2972  }
2973 
2974 
2975 
2977  SparseMatrix::n() const
2978  {
2979  // If the matrix structure has not been fixed (i.e., we did not have a
2980  // sparsity pattern), it does not know about the number of columns so we
2981  // must always take this from the additional column space map
2982  Assert(column_space_map.get() != nullptr, ExcInternalError());
2983  return n_global_elements(*column_space_map);
2984  }
2985 
2986 
2987 
2988  inline unsigned int
2989  SparseMatrix::local_size() const
2990  {
2991  return matrix->NumMyRows();
2992  }
2993 
2994 
2995 
2996  inline std::pair<SparseMatrix::size_type, SparseMatrix::size_type>
2998  {
2999  size_type begin, end;
3000 # ifndef DEAL_II_WITH_64BIT_INDICES
3001  begin = matrix->RowMap().MinMyGID();
3002  end = matrix->RowMap().MaxMyGID() + 1;
3003 # else
3004  begin = matrix->RowMap().MinMyGID64();
3005  end = matrix->RowMap().MaxMyGID64() + 1;
3006 # endif
3007 
3008  return std::make_pair(begin, end);
3009  }
3010 
3011 
3012 
3013  inline std::uint64_t
3015  {
3016  // Trilinos uses 64bit functions internally for attribute access, which
3017  // return `long long`. They also offer 32bit variants that return `int`,
3018  // however those call the 64bit version and convert the values to 32bit.
3019  // There is no necessity in using the 32bit versions at all.
3020  return static_cast<std::uint64_t>(matrix->NumGlobalNonzeros64());
3021  }
3022 
3023 
3024 
3025  template <typename SparsityPatternType>
3026  inline void
3027  SparseMatrix::reinit(const IndexSet & parallel_partitioning,
3028  const SparsityPatternType &sparsity_pattern,
3029  const MPI_Comm & communicator,
3030  const bool exchange_data)
3031  {
3032  reinit(parallel_partitioning,
3033  parallel_partitioning,
3034  sparsity_pattern,
3035  communicator,
3036  exchange_data);
3037  }
3038 
3039 
3040 
3041  template <typename number>
3042  inline void
3043  SparseMatrix::reinit(const IndexSet &parallel_partitioning,
3044  const ::SparseMatrix<number> &sparse_matrix,
3045  const MPI_Comm & communicator,
3046  const double drop_tolerance,
3047  const bool copy_values,
3048  const ::SparsityPattern * use_this_sparsity)
3049  {
3050  Epetra_Map map =
3051  parallel_partitioning.make_trilinos_map(communicator, false);
3052  reinit(parallel_partitioning,
3053  parallel_partitioning,
3054  sparse_matrix,
3055  drop_tolerance,
3056  copy_values,
3057  use_this_sparsity);
3058  }
3059 
3060 
3061 
3062  inline const Epetra_CrsMatrix &
3064  {
3065  return static_cast<const Epetra_CrsMatrix &>(*matrix);
3066  }
3067 
3068 
3069 
3070  inline const Epetra_CrsGraph &
3072  {
3073  return matrix->Graph();
3074  }
3075 
3076 
3077 
3078  inline IndexSet
3080  {
3081  return IndexSet(matrix->DomainMap());
3082  }
3083 
3084 
3085 
3086  inline IndexSet
3088  {
3089  return IndexSet(matrix->RangeMap());
3090  }
3091 
3092 
3093 
3094  inline void
3096  {
3097  // nothing to do here
3098  }
3099 
3100 
3101 
3102  inline void
3104  {
3105  // nothing to do here
3106  }
3107 
3108 
3109  namespace internal
3110  {
3111  namespace LinearOperatorImplementation
3112  {
3113  template <typename EpetraOpType>
3114  TrilinosPayload::TrilinosPayload(
3115  EpetraOpType & op,
3116  const bool supports_inverse_operations,
3117  const bool use_transpose,
3118  const MPI_Comm &mpi_communicator,
3119  const IndexSet &locally_owned_domain_indices,
3120  const IndexSet &locally_owned_range_indices)
3121  : use_transpose(use_transpose)
3122  , communicator(mpi_communicator)
3123  , domain_map(
3124  locally_owned_domain_indices.make_trilinos_map(communicator.Comm()))
3125  , range_map(
3126  locally_owned_range_indices.make_trilinos_map(communicator.Comm()))
3127  {
3128  vmult = [&op](Range &tril_dst, const Domain &tril_src) {
3129  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3130  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3131  Assert(&tril_src != &tril_dst,
3134  tril_src,
3135  tril_dst,
3136  op.UseTranspose());
3137 
3138  const int ierr = op.Apply(tril_src, tril_dst);
3139  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3140  };
3141 
3142  Tvmult = [&op](Domain &tril_dst, const Range &tril_src) {
3143  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3144  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3145  Assert(&tril_src != &tril_dst,
3148  tril_src,
3149  tril_dst,
3150  !op.UseTranspose());
3151 
3152  op.SetUseTranspose(!op.UseTranspose());
3153  const int ierr = op.Apply(tril_src, tril_dst);
3154  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3155  op.SetUseTranspose(!op.UseTranspose());
3156  };
3157 
3158  if (supports_inverse_operations)
3159  {
3160  inv_vmult = [&op](Domain &tril_dst, const Range &tril_src) {
3161  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3162  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3163  Assert(
3164  &tril_src != &tril_dst,
3167  tril_src,
3168  tril_dst,
3169  !op.UseTranspose());
3170 
3171  const int ierr = op.ApplyInverse(tril_src, tril_dst);
3172  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3173  };
3174 
3175  inv_Tvmult = [&op](Range &tril_dst, const Domain &tril_src) {
3176  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3177  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3178  Assert(
3179  &tril_src != &tril_dst,
3182  tril_src,
3183  tril_dst,
3184  op.UseTranspose());
3185 
3186  op.SetUseTranspose(!op.UseTranspose());
3187  const int ierr = op.ApplyInverse(tril_src, tril_dst);
3188  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3189  op.SetUseTranspose(!op.UseTranspose());
3190  };
3191  }
3192  else
3193  {
3194  inv_vmult = [](Domain &, const Range &) {
3195  Assert(false,
3196  ExcMessage(
3197  "Uninitialized TrilinosPayload::inv_vmult called. "
3198  "The operator does not support inverse operations."));
3199  };
3200 
3201  inv_Tvmult = [](Range &, const Domain &) {
3202  Assert(false,
3203  ExcMessage(
3204  "Uninitialized TrilinosPayload::inv_Tvmult called. "
3205  "The operator does not support inverse operations."));
3206  };
3207  }
3208  }
3209 
3210 
3211  template <typename Solver, typename Preconditioner>
3212  std::enable_if_t<
3213  std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
3214  std::is_base_of<TrilinosWrappers::PreconditionBase,
3215  Preconditioner>::value,
3216  TrilinosPayload>
3218  Solver & solver,
3219  const Preconditioner &preconditioner) const
3220  {
3221  const auto &payload = *this;
3222 
3223  TrilinosPayload return_op(payload);
3224 
3225  // Capture by copy so the payloads are always valid
3226 
3227  return_op.inv_vmult = [payload, &solver, &preconditioner](
3228  TrilinosPayload::Domain & tril_dst,
3229  const TrilinosPayload::Range &tril_src) {
3230  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3231  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3232  Assert(&tril_src != &tril_dst,
3235  tril_src,
3236  tril_dst,
3237  !payload.UseTranspose());
3238  solver.solve(payload, tril_dst, tril_src, preconditioner);
3239  };
3240 
3241  return_op.inv_Tvmult = [payload, &solver, &preconditioner](
3242  TrilinosPayload::Range & tril_dst,
3243  const TrilinosPayload::Domain &tril_src) {
3244  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3245  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3246  Assert(&tril_src != &tril_dst,
3249  tril_src,
3250  tril_dst,
3251  payload.UseTranspose());
3252 
3253  const_cast<TrilinosPayload &>(payload).transpose();
3254  solver.solve(payload, tril_dst, tril_src, preconditioner);
3255  const_cast<TrilinosPayload &>(payload).transpose();
3256  };
3257 
3258  // If the input operator is already setup for transpose operations, then
3259  // we must do similar with its inverse.
3260  if (return_op.UseTranspose() == true)
3261  std::swap(return_op.inv_vmult, return_op.inv_Tvmult);
3262 
3263  return return_op;
3264  }
3265 
3266  template <typename Solver, typename Preconditioner>
3267  std::enable_if_t<
3268  !(std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
3269  std::is_base_of<TrilinosWrappers::PreconditionBase,
3270  Preconditioner>::value),
3271  TrilinosPayload>
3272  TrilinosPayload::inverse_payload(Solver &, const Preconditioner &) const
3273  {
3274  TrilinosPayload return_op(*this);
3275 
3276  return_op.inv_vmult = [](TrilinosPayload::Domain &,
3277  const TrilinosPayload::Range &) {
3278  AssertThrow(false,
3279  ExcMessage("Payload inv_vmult disabled because of "
3280  "incompatible solver/preconditioner choice."));
3281  };
3282 
3283  return_op.inv_Tvmult = [](TrilinosPayload::Range &,
3284  const TrilinosPayload::Domain &) {
3285  AssertThrow(false,
3286  ExcMessage("Payload inv_vmult disabled because of "
3287  "incompatible solver/preconditioner choice."));
3288  };
3289 
3290  return return_op;
3291  }
3292  } // namespace LinearOperatorImplementation
3293  } // namespace internal
3294 
3295  template <>
3296  void
3297  SparseMatrix::set<TrilinosScalar>(const size_type row,
3298  const size_type n_cols,
3299  const size_type * col_indices,
3300  const TrilinosScalar *values,
3301  const bool elide_zero_values);
3302 # endif // DOXYGEN
3303 
3304 } /* namespace TrilinosWrappers */
3305 
3306 
3308 
3309 
3310 # endif // DEAL_II_WITH_TRILINOS
3311 
3312 
3313 /*----------------------- trilinos_sparse_matrix.h --------------------*/
3314 
3315 #endif
3316 /*----------------------- 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:789
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)
bool operator>(const Iterator< Constness > &) const
typename Accessor< Constness >::MatrixType MatrixType
Iterator(const Iterator< Other > &other)
bool operator==(const Iterator< Constness > &) const
bool operator!=(const Iterator< Constness > &) const
Iterator(MatrixType *matrix, const size_type row, const size_type index)
const Accessor< Constness > * operator->() const
bool operator<(const Iterator< Constness > &) 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
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)
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< SparsityPatternType, ::SparseMatrix< double > >::value > 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 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< typename VectorType::value_type, TrilinosScalar >::value > Tvmult(VectorType &dst, const VectorType &src) const
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
const_iterator begin(const size_type r) const
SparseMatrix(const SparseMatrix &)=delete
iterator end(const size_type r)
void reinit(const IndexSet &parallel_partitioning, const SparsityPatternType &sparsity_pattern, const MPI_Comm &communicator=MPI_COMM_WORLD, const bool exchange_data=false)
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
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)
std::enable_if_t< std::is_same< typename VectorType::value_type, TrilinosScalar >::value > vmult(VectorType &dst, const VectorType &src) const
SparseMatrix & operator=(const SparseMatrix &)=delete
std::unique_ptr< Epetra_CrsMatrix > nonlocal_matrix
std::function< void(VectorType &, const VectorType &)> inv_Tvmult
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)
virtual int ApplyInverse(const VectorType &Y, VectorType &X) const override
virtual int Apply(const VectorType &X, VectorType &Y) const override
std::enable_if_t< !(std::is_base_of< TrilinosWrappers::SolverBase, Solver >::value &&std::is_base_of< TrilinosWrappers::PreconditionBase, Preconditioner >::value), TrilinosPayload > inverse_payload(Solver &, const Preconditioner &) const
std::function< void(VectorType &, const VectorType &)> inv_vmult
std::enable_if_t< std::is_base_of< TrilinosWrappers::SolverBase, Solver >::value &&std::is_base_of< TrilinosWrappers::PreconditionBase, Preconditioner >::value, TrilinosPayload > inverse_payload(Solver &, const Preconditioner &) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
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:465
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:579
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:1586
#define AssertIsFinite(number)
Definition: exceptions.h:1854
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:533
static ::ExceptionBase & ExcAccessToNonLocalElement(size_type arg1, size_type arg2, size_type arg3, size_type arg4)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1827
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:488
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:556
static ::ExceptionBase & ExcMatrixNotCompressed()
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:510
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:1675
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