Reference documentation for deal.II version GIT 6da2e5d553 2022-07-01 18:55:02+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>
26 
27 # include <deal.II/lac/exceptions.h>
28 # include <deal.II/lac/full_matrix.h>
35 
36 # include <Epetra_Comm.h>
37 # include <Epetra_CrsGraph.h>
38 # include <Epetra_Export.h>
39 # include <Epetra_FECrsMatrix.h>
40 # include <Epetra_Map.h>
41 # include <Epetra_MpiComm.h>
42 # include <Epetra_MultiVector.h>
43 # include <Epetra_Operator.h>
44 # include <mpi.h>
45 
46 # include <cmath>
47 # include <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);
733 
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  typename std::enable_if<
844  !std::is_same<SparsityPatternType,
845  ::SparseMatrix<double>>::value>::type
846  reinit(const IndexSet & row_parallel_partitioning,
847  const IndexSet & col_parallel_partitioning,
848  const SparsityPatternType &sparsity_pattern,
849  const MPI_Comm & communicator = MPI_COMM_WORLD,
850  const bool exchange_data = false);
851 
868  template <typename number>
869  void
870  reinit(const IndexSet & parallel_partitioning,
871  const ::SparseMatrix<number> &dealii_sparse_matrix,
872  const MPI_Comm & communicator = MPI_COMM_WORLD,
873  const double drop_tolerance = 1e-13,
874  const bool copy_values = true,
875  const ::SparsityPattern * use_this_sparsity = nullptr);
876 
890  template <typename number>
891  void
892  reinit(const IndexSet & row_parallel_partitioning,
893  const IndexSet & col_parallel_partitioning,
894  const ::SparseMatrix<number> &dealii_sparse_matrix,
895  const MPI_Comm & communicator = MPI_COMM_WORLD,
896  const double drop_tolerance = 1e-13,
897  const bool copy_values = true,
898  const ::SparsityPattern * use_this_sparsity = nullptr);
900 
904 
908  size_type
909  m() const;
910 
914  size_type
915  n() const;
916 
925  unsigned int
926  local_size() const;
927 
936  std::pair<size_type, size_type>
937  local_range() const;
938 
943  bool
944  in_local_range(const size_type index) const;
945 
950  std::uint64_t
952 
956  unsigned int
957  row_length(const size_type row) const;
958 
965  bool
966  is_compressed() const;
967 
973  size_type
974  memory_consumption() const;
975 
979  MPI_Comm
980  get_mpi_communicator() const;
981 
983 
987 
997  SparseMatrix &
998  operator=(const double d);
999 
1007  void
1008  clear();
1009 
1037  void
1038  compress(::VectorOperation::values operation);
1039 
1061  void
1062  set(const size_type i, const size_type j, const TrilinosScalar value);
1063 
1096  void
1097  set(const std::vector<size_type> & indices,
1098  const FullMatrix<TrilinosScalar> &full_matrix,
1099  const bool elide_zero_values = false);
1100 
1106  void
1107  set(const std::vector<size_type> & row_indices,
1108  const std::vector<size_type> & col_indices,
1109  const FullMatrix<TrilinosScalar> &full_matrix,
1110  const bool elide_zero_values = false);
1111 
1139  void
1140  set(const size_type row,
1141  const std::vector<size_type> & col_indices,
1142  const std::vector<TrilinosScalar> &values,
1143  const bool elide_zero_values = false);
1144 
1172  template <typename Number>
1173  void
1174  set(const size_type row,
1175  const size_type n_cols,
1176  const size_type *col_indices,
1177  const Number * values,
1178  const bool elide_zero_values = false);
1179 
1189  void
1190  add(const size_type i, const size_type j, const TrilinosScalar value);
1191 
1210  void
1211  add(const std::vector<size_type> & indices,
1212  const FullMatrix<TrilinosScalar> &full_matrix,
1213  const bool elide_zero_values = true);
1214 
1220  void
1221  add(const std::vector<size_type> & row_indices,
1222  const std::vector<size_type> & col_indices,
1223  const FullMatrix<TrilinosScalar> &full_matrix,
1224  const bool elide_zero_values = true);
1225 
1239  void
1240  add(const size_type row,
1241  const std::vector<size_type> & col_indices,
1242  const std::vector<TrilinosScalar> &values,
1243  const bool elide_zero_values = true);
1244 
1258  void
1259  add(const size_type row,
1260  const size_type n_cols,
1261  const size_type * col_indices,
1262  const TrilinosScalar *values,
1263  const bool elide_zero_values = true,
1264  const bool col_indices_are_sorted = false);
1265 
1269  SparseMatrix &
1270  operator*=(const TrilinosScalar factor);
1271 
1275  SparseMatrix &
1276  operator/=(const TrilinosScalar factor);
1277 
1281  void
1282  copy_from(const SparseMatrix &source);
1283 
1291  void
1292  add(const TrilinosScalar factor, const SparseMatrix &matrix);
1293 
1320  void
1321  clear_row(const size_type row, const TrilinosScalar new_diag_value = 0);
1322 
1343  void
1344  clear_rows(const std::vector<size_type> &rows,
1345  const TrilinosScalar new_diag_value = 0);
1346 
1356  void
1357  transpose();
1358 
1360 
1364 
1374  operator()(const size_type i, const size_type j) const;
1375 
1393  el(const size_type i, const size_type j) const;
1394 
1402  diag_element(const size_type i) const;
1403 
1405 
1409 
1440  template <typename VectorType>
1441  typename std::enable_if<std::is_same<typename VectorType::value_type,
1442  TrilinosScalar>::value>::type
1443  vmult(VectorType &dst, const VectorType &src) const;
1444 
1451  template <typename VectorType>
1452  typename std::enable_if<!std::is_same<typename VectorType::value_type,
1453  TrilinosScalar>::value>::type
1454  vmult(VectorType &dst, const VectorType &src) const;
1455 
1470  template <typename VectorType>
1471  typename std::enable_if<std::is_same<typename VectorType::value_type,
1472  TrilinosScalar>::value>::type
1473  Tvmult(VectorType &dst, const VectorType &src) const;
1474 
1481  template <typename VectorType>
1482  typename std::enable_if<!std::is_same<typename VectorType::value_type,
1483  TrilinosScalar>::value>::type
1484  Tvmult(VectorType &dst, const VectorType &src) const;
1485 
1495  template <typename VectorType>
1496  void
1497  vmult_add(VectorType &dst, const VectorType &src) const;
1498 
1509  template <typename VectorType>
1510  void
1511  Tvmult_add(VectorType &dst, const VectorType &src) const;
1512 
1535  matrix_norm_square(const MPI::Vector &v) const;
1536 
1557  matrix_scalar_product(const MPI::Vector &u, const MPI::Vector &v) const;
1558 
1576  residual(MPI::Vector & dst,
1577  const MPI::Vector &x,
1578  const MPI::Vector &b) const;
1579 
1594  void
1595  mmult(SparseMatrix & C,
1596  const SparseMatrix &B,
1597  const MPI::Vector & V = MPI::Vector()) const;
1598 
1599 
1616  void
1617  Tmmult(SparseMatrix & C,
1618  const SparseMatrix &B,
1619  const MPI::Vector & V = MPI::Vector()) const;
1620 
1622 
1626 
1635  l1_norm() const;
1636 
1646  linfty_norm() const;
1647 
1653  frobenius_norm() const;
1654 
1656 
1660 
1665  const Epetra_CrsMatrix &
1667 
1672  const Epetra_CrsGraph &
1674 
1676 
1681 
1686  IndexSet
1688 
1694  IndexSet
1696 
1698 
1703 
1723  begin() const;
1724 
1728  iterator
1730 
1736  end() const;
1737 
1741  iterator
1742  end();
1743 
1773  begin(const size_type r) const;
1774 
1778  iterator
1779  begin(const size_type r);
1780 
1791  end(const size_type r) const;
1792 
1796  iterator
1797  end(const size_type r);
1798 
1800 
1804 
1810  void
1811  write_ascii();
1812 
1820  void
1821  print(std::ostream &out,
1822  const bool write_extended_trilinos_info = false) const;
1823 
1825 
1834  int,
1835  << "An error with error number " << arg1
1836  << " occurred while calling a Trilinos function. "
1837  "\n\n"
1838  "For historical reasons, many Trilinos functions express "
1839  "errors by returning specific integer values to indicate "
1840  "certain errors. Unfortunately, different Trilinos functions "
1841  "often use the same integer values for different kinds of "
1842  "errors, and in most cases it is also not documented what "
1843  "each error code actually means. As a consequence, it is often "
1844  "difficult to say what a particular error (in this case, "
1845  "the error with integer code '"
1846  << arg1
1847  << "') represents and how one should fix a code to avoid it. "
1848  "The best one can often do is to look up the call stack to "
1849  "see which deal.II function generated the error, and which "
1850  "Trilinos function the error code had originated from; "
1851  "then look up the Trilinos source code of that function (for "
1852  "example on github) to see what code path set that error "
1853  "code. Short of going through all of that, the only other "
1854  "option is to guess the cause of the error from "
1855  "the context in which the error appeared.");
1856 
1857 
1862  size_type,
1863  size_type,
1864  << "The entry with index <" << arg1 << ',' << arg2
1865  << "> does not exist.");
1866 
1871  "You are attempting an operation on two matrices that "
1872  "are the same object, but the operation requires that the "
1873  "two objects are in fact different.");
1874 
1879 
1884  size_type,
1885  size_type,
1886  size_type,
1887  size_type,
1888  << "You tried to access element (" << arg1 << '/' << arg2
1889  << ')'
1890  << " of a distributed matrix, but only rows in range ["
1891  << arg3 << ',' << arg4
1892  << "] are stored locally and can be accessed.");
1893 
1898  size_type,
1899  size_type,
1900  << "You tried to access element (" << arg1 << '/' << arg2
1901  << ')' << " of a sparse matrix, but it appears to not"
1902  << " exist in the Trilinos sparsity pattern.");
1904 
1905 
1906 
1907  protected:
1918  void
1920 
1926  void
1928 
1929 
1930 
1931  private:
1936  std::unique_ptr<Epetra_Map> column_space_map;
1937 
1943  std::unique_ptr<Epetra_FECrsMatrix> matrix;
1944 
1950  std::unique_ptr<Epetra_CrsMatrix> nonlocal_matrix;
1951 
1955  std::unique_ptr<Epetra_Export> nonlocal_matrix_exporter;
1956 
1968  Epetra_CombineMode last_action;
1969 
1975 
1976  // To allow calling protected prepare_add() and prepare_set().
1977  friend class BlockMatrixBase<SparseMatrix>;
1978  };
1979 
1980 
1981 
1982  // forwards declarations
1983  class SolverBase;
1984  class PreconditionBase;
1985 
1986  namespace internal
1987  {
1988  inline void
1989  check_vector_map_equality(const Epetra_CrsMatrix & mtrx,
1990  const Epetra_MultiVector &src,
1991  const Epetra_MultiVector &dst,
1992  const bool transpose)
1993  {
1994  if (transpose == false)
1995  {
1996  Assert(src.Map().SameAs(mtrx.DomainMap()) == true,
1997  ExcMessage(
1998  "Column map of matrix does not fit with vector map!"));
1999  Assert(dst.Map().SameAs(mtrx.RangeMap()) == true,
2000  ExcMessage("Row map of matrix does not fit with vector map!"));
2001  }
2002  else
2003  {
2004  Assert(src.Map().SameAs(mtrx.RangeMap()) == true,
2005  ExcMessage(
2006  "Column map of matrix does not fit with vector map!"));
2007  Assert(dst.Map().SameAs(mtrx.DomainMap()) == true,
2008  ExcMessage("Row map of matrix does not fit with vector map!"));
2009  }
2010  (void)mtrx; // removes -Wunused-variable in optimized mode
2011  (void)src;
2012  (void)dst;
2013  }
2014 
2015  inline void
2017  const Epetra_MultiVector &src,
2018  const Epetra_MultiVector &dst,
2019  const bool transpose)
2020  {
2021  if (transpose == false)
2022  {
2023  Assert(src.Map().SameAs(op.OperatorDomainMap()) == true,
2024  ExcMessage(
2025  "Column map of operator does not fit with vector map!"));
2026  Assert(dst.Map().SameAs(op.OperatorRangeMap()) == true,
2027  ExcMessage(
2028  "Row map of operator does not fit with vector map!"));
2029  }
2030  else
2031  {
2032  Assert(src.Map().SameAs(op.OperatorRangeMap()) == true,
2033  ExcMessage(
2034  "Column map of operator does not fit with vector map!"));
2035  Assert(dst.Map().SameAs(op.OperatorDomainMap()) == true,
2036  ExcMessage(
2037  "Row map of operator does not fit with vector map!"));
2038  }
2039  (void)op; // removes -Wunused-variable in optimized mode
2040  (void)src;
2041  (void)dst;
2042  }
2043 
2044  namespace LinearOperatorImplementation
2045  {
2066  {
2067  public:
2071  using VectorType = Epetra_MultiVector;
2072 
2077 
2082 
2087 
2095  TrilinosPayload();
2096 
2100  TrilinosPayload(const TrilinosWrappers::SparseMatrix &matrix_exemplar,
2102 
2106  TrilinosPayload(const TrilinosPayload & payload_exemplar,
2108 
2113  const TrilinosWrappers::SparseMatrix & matrix_exemplar,
2114  const TrilinosWrappers::PreconditionBase &preconditioner);
2115 
2120  const TrilinosWrappers::PreconditionBase &preconditioner_exemplar,
2121  const TrilinosWrappers::PreconditionBase &preconditioner);
2122 
2127  const TrilinosPayload & payload_exemplar,
2128  const TrilinosWrappers::PreconditionBase &preconditioner);
2129 
2133  TrilinosPayload(const TrilinosPayload &payload);
2134 
2142  TrilinosPayload(const TrilinosPayload &first_op,
2143  const TrilinosPayload &second_op);
2144 
2148  virtual ~TrilinosPayload() override = default;
2149 
2154  identity_payload() const;
2155 
2160  null_payload() const;
2161 
2166  transpose_payload() const;
2167 
2184  template <typename Solver, typename Preconditioner>
2185  typename std::enable_if<
2186  std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
2187  std::is_base_of<TrilinosWrappers::PreconditionBase,
2188  Preconditioner>::value,
2189  TrilinosPayload>::type
2190  inverse_payload(Solver &, const Preconditioner &) const;
2191 
2209  template <typename Solver, typename Preconditioner>
2210  typename std::enable_if<
2211  !(std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
2212  std::is_base_of<TrilinosWrappers::PreconditionBase,
2213  Preconditioner>::value),
2214  TrilinosPayload>::type
2215  inverse_payload(Solver &, const Preconditioner &) const;
2216 
2218 
2223 
2229  IndexSet
2231 
2237  IndexSet
2239 
2243  MPI_Comm
2244  get_mpi_communicator() const;
2245 
2252  void
2253  transpose();
2254 
2262  std::function<void(VectorType &, const VectorType &)> vmult;
2263 
2271  std::function<void(VectorType &, const VectorType &)> Tvmult;
2272 
2281  std::function<void(VectorType &, const VectorType &)> inv_vmult;
2282 
2291  std::function<void(VectorType &, const VectorType &)> inv_Tvmult;
2292 
2294 
2299 
2306  virtual bool
2307  UseTranspose() const override;
2308 
2324  virtual int
2325  SetUseTranspose(bool UseTranspose) override;
2326 
2338  virtual int
2339  Apply(const VectorType &X, VectorType &Y) const override;
2340 
2359  virtual int
2360  ApplyInverse(const VectorType &Y, VectorType &X) const override;
2362 
2367 
2374  virtual const char *
2375  Label() const override;
2376 
2384  virtual const Epetra_Comm &
2385  Comm() const override;
2386 
2394  virtual const Epetra_Map &
2395  OperatorDomainMap() const override;
2396 
2405  virtual const Epetra_Map &
2406  OperatorRangeMap() const override;
2408 
2409  private:
2419  template <typename EpetraOpType>
2420  TrilinosPayload(EpetraOpType & op,
2421  const bool supports_inverse_operations,
2422  const bool use_transpose,
2423  const MPI_Comm &mpi_communicator,
2426 
2432 
2437  Epetra_MpiComm communicator;
2438 
2443  Epetra_Map domain_map;
2444 
2449  Epetra_Map range_map;
2450 
2459  virtual bool
2460  HasNormInf() const override;
2461 
2469  virtual double
2470  NormInf() const override;
2471  };
2472 
2478  operator+(const TrilinosPayload &first_op,
2479  const TrilinosPayload &second_op);
2480 
2486  operator*(const TrilinosPayload &first_op,
2487  const TrilinosPayload &second_op);
2488 
2489  } // namespace LinearOperatorImplementation
2490  } /* namespace internal */
2491 
2492 
2493 
2494  // ----------------------- inline and template functions --------------------
2495 
2496 # ifndef DOXYGEN
2497 
2498  namespace SparseMatrixIterators
2499  {
2501  size_type row,
2502  size_type index)
2503  : matrix(matrix)
2504  , a_row(row)
2505  , a_index(index)
2506  {
2507  visit_present_row();
2508  }
2509 
2510 
2512  AccessorBase::row() const
2513  {
2514  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2515  return a_row;
2516  }
2517 
2518 
2520  AccessorBase::column() const
2521  {
2522  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2523  return (*colnum_cache)[a_index];
2524  }
2525 
2526 
2528  AccessorBase::index() const
2529  {
2530  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2531  return a_index;
2532  }
2533 
2534 
2535  inline Accessor<true>::Accessor(MatrixType * matrix,
2536  const size_type row,
2537  const size_type index)
2538  : AccessorBase(const_cast<SparseMatrix *>(matrix), row, index)
2539  {}
2540 
2541 
2542  template <bool Other>
2543  inline Accessor<true>::Accessor(const Accessor<Other> &other)
2544  : AccessorBase(other)
2545  {}
2546 
2547 
2548  inline TrilinosScalar
2549  Accessor<true>::value() const
2550  {
2551  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2552  return (*value_cache)[a_index];
2553  }
2554 
2555 
2556  inline Accessor<false>::Reference::Reference(const Accessor<false> &acc)
2557  : accessor(const_cast<Accessor<false> &>(acc))
2558  {}
2559 
2560 
2561  inline Accessor<false>::Reference::operator TrilinosScalar() const
2562  {
2563  return (*accessor.value_cache)[accessor.a_index];
2564  }
2565 
2566  inline const Accessor<false>::Reference &
2567  Accessor<false>::Reference::operator=(const TrilinosScalar n) const
2568  {
2569  (*accessor.value_cache)[accessor.a_index] = n;
2570  accessor.matrix->set(accessor.row(),
2571  accessor.column(),
2572  static_cast<TrilinosScalar>(*this));
2573  return *this;
2574  }
2575 
2576 
2577  inline const Accessor<false>::Reference &
2578  Accessor<false>::Reference::operator+=(const TrilinosScalar n) const
2579  {
2580  (*accessor.value_cache)[accessor.a_index] += n;
2581  accessor.matrix->set(accessor.row(),
2582  accessor.column(),
2583  static_cast<TrilinosScalar>(*this));
2584  return *this;
2585  }
2586 
2587 
2588  inline const Accessor<false>::Reference &
2589  Accessor<false>::Reference::operator-=(const TrilinosScalar n) const
2590  {
2591  (*accessor.value_cache)[accessor.a_index] -= n;
2592  accessor.matrix->set(accessor.row(),
2593  accessor.column(),
2594  static_cast<TrilinosScalar>(*this));
2595  return *this;
2596  }
2597 
2598 
2599  inline const Accessor<false>::Reference &
2600  Accessor<false>::Reference::operator*=(const TrilinosScalar n) const
2601  {
2602  (*accessor.value_cache)[accessor.a_index] *= n;
2603  accessor.matrix->set(accessor.row(),
2604  accessor.column(),
2605  static_cast<TrilinosScalar>(*this));
2606  return *this;
2607  }
2608 
2609 
2610  inline const Accessor<false>::Reference &
2611  Accessor<false>::Reference::operator/=(const TrilinosScalar n) const
2612  {
2613  (*accessor.value_cache)[accessor.a_index] /= n;
2614  accessor.matrix->set(accessor.row(),
2615  accessor.column(),
2616  static_cast<TrilinosScalar>(*this));
2617  return *this;
2618  }
2619 
2620 
2621  inline Accessor<false>::Accessor(MatrixType * matrix,
2622  const size_type row,
2623  const size_type index)
2624  : AccessorBase(matrix, row, index)
2625  {}
2626 
2627 
2628  inline Accessor<false>::Reference
2629  Accessor<false>::value() const
2630  {
2631  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2632  return {*this};
2633  }
2634 
2635 
2636 
2637  template <bool Constness>
2638  inline Iterator<Constness>::Iterator(MatrixType * matrix,
2639  const size_type row,
2640  const size_type index)
2641  : accessor(matrix, row, index)
2642  {}
2643 
2644 
2645  template <bool Constness>
2646  template <bool Other>
2647  inline Iterator<Constness>::Iterator(const Iterator<Other> &other)
2648  : accessor(other.accessor)
2649  {}
2650 
2651 
2652  template <bool Constness>
2653  inline Iterator<Constness> &
2655  {
2656  Assert(accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
2657 
2658  ++accessor.a_index;
2659 
2660  // If at end of line: do one
2661  // step, then cycle until we
2662  // find a row with a nonzero
2663  // number of entries.
2664  if (accessor.a_index >= accessor.colnum_cache->size())
2665  {
2666  accessor.a_index = 0;
2667  ++accessor.a_row;
2668 
2669  while ((accessor.a_row < accessor.matrix->m()) &&
2670  ((accessor.matrix->in_local_range(accessor.a_row) == false) ||
2671  (accessor.matrix->row_length(accessor.a_row) == 0)))
2672  ++accessor.a_row;
2673 
2674  accessor.visit_present_row();
2675  }
2676  return *this;
2677  }
2678 
2679 
2680  template <bool Constness>
2681  inline Iterator<Constness>
2683  {
2684  const Iterator<Constness> old_state = *this;
2685  ++(*this);
2686  return old_state;
2687  }
2688 
2689 
2690 
2691  template <bool Constness>
2692  inline const Accessor<Constness> &
2694  {
2695  return accessor;
2696  }
2697 
2698 
2699 
2700  template <bool Constness>
2701  inline const Accessor<Constness> *
2702  Iterator<Constness>::operator->() const
2703  {
2704  return &accessor;
2705  }
2706 
2707 
2708 
2709  template <bool Constness>
2710  inline bool
2711  Iterator<Constness>::operator==(const Iterator<Constness> &other) const
2712  {
2713  return (accessor.a_row == other.accessor.a_row &&
2714  accessor.a_index == other.accessor.a_index);
2715  }
2716 
2717 
2718 
2719  template <bool Constness>
2720  inline bool
2721  Iterator<Constness>::operator!=(const Iterator<Constness> &other) const
2722  {
2723  return !(*this == other);
2724  }
2725 
2726 
2727 
2728  template <bool Constness>
2729  inline bool
2730  Iterator<Constness>::operator<(const Iterator<Constness> &other) const
2731  {
2732  return (accessor.row() < other.accessor.row() ||
2733  (accessor.row() == other.accessor.row() &&
2734  accessor.index() < other.accessor.index()));
2735  }
2736 
2737 
2738  template <bool Constness>
2739  inline bool
2740  Iterator<Constness>::operator>(const Iterator<Constness> &other) const
2741  {
2742  return (other < *this);
2743  }
2744 
2745  } // namespace SparseMatrixIterators
2746 
2747 
2748 
2750  SparseMatrix::begin() const
2751  {
2752  return begin(0);
2753  }
2754 
2755 
2756 
2758  SparseMatrix::end() const
2759  {
2760  return const_iterator(this, m(), 0);
2761  }
2762 
2763 
2764 
2766  SparseMatrix::begin(const size_type r) const
2767  {
2768  AssertIndexRange(r, m());
2769  if (in_local_range(r) && (row_length(r) > 0))
2770  return const_iterator(this, r, 0);
2771  else
2772  return end(r);
2773  }
2774 
2775 
2776 
2778  SparseMatrix::end(const size_type r) const
2779  {
2780  AssertIndexRange(r, m());
2781 
2782  // place the iterator on the first entry
2783  // past this line, or at the end of the
2784  // matrix
2785  for (size_type i = r + 1; i < m(); ++i)
2786  if (in_local_range(i) && (row_length(i) > 0))
2787  return const_iterator(this, i, 0);
2788 
2789  // if there is no such line, then take the
2790  // end iterator of the matrix
2791  return end();
2792  }
2793 
2794 
2795 
2796  inline SparseMatrix::iterator
2798  {
2799  return begin(0);
2800  }
2801 
2802 
2803 
2804  inline SparseMatrix::iterator
2806  {
2807  return iterator(this, m(), 0);
2808  }
2809 
2810 
2811 
2812  inline SparseMatrix::iterator
2814  {
2815  AssertIndexRange(r, m());
2816  if (in_local_range(r) && (row_length(r) > 0))
2817  return iterator(this, r, 0);
2818  else
2819  return end(r);
2820  }
2821 
2822 
2823 
2824  inline SparseMatrix::iterator
2825  SparseMatrix::end(const size_type r)
2826  {
2827  AssertIndexRange(r, m());
2828 
2829  // place the iterator on the first entry
2830  // past this line, or at the end of the
2831  // matrix
2832  for (size_type i = r + 1; i < m(); ++i)
2833  if (in_local_range(i) && (row_length(i) > 0))
2834  return iterator(this, i, 0);
2835 
2836  // if there is no such line, then take the
2837  // end iterator of the matrix
2838  return end();
2839  }
2840 
2841 
2842 
2843  inline bool
2844  SparseMatrix::in_local_range(const size_type index) const
2845  {
2847 # ifndef DEAL_II_WITH_64BIT_INDICES
2848  begin = matrix->RowMap().MinMyGID();
2849  end = matrix->RowMap().MaxMyGID() + 1;
2850 # else
2851  begin = matrix->RowMap().MinMyGID64();
2852  end = matrix->RowMap().MaxMyGID64() + 1;
2853 # endif
2854 
2855  return ((index >= static_cast<size_type>(begin)) &&
2856  (index < static_cast<size_type>(end)));
2857  }
2858 
2859 
2860 
2861  inline bool
2863  {
2864  return compressed;
2865  }
2866 
2867 
2868 
2869  // Inline the set() and add() functions, since they will be called
2870  // frequently, and the compiler can optimize away some unnecessary loops
2871  // when the sizes are given at compile time.
2872  template <>
2873  void
2874  SparseMatrix::set<TrilinosScalar>(const size_type row,
2875  const size_type n_cols,
2876  const size_type * col_indices,
2877  const TrilinosScalar *values,
2878  const bool elide_zero_values);
2879 
2880 
2881 
2882  template <typename Number>
2883  void
2884  SparseMatrix::set(const size_type row,
2885  const size_type n_cols,
2886  const size_type *col_indices,
2887  const Number * values,
2888  const bool elide_zero_values)
2889  {
2890  std::vector<TrilinosScalar> trilinos_values(n_cols);
2891  std::copy(values, values + n_cols, trilinos_values.begin());
2892  this->set(
2893  row, n_cols, col_indices, trilinos_values.data(), elide_zero_values);
2894  }
2895 
2896 
2897 
2898  inline void
2899  SparseMatrix::set(const size_type i,
2900  const size_type j,
2901  const TrilinosScalar value)
2902  {
2903  AssertIsFinite(value);
2904 
2905  set(i, 1, &j, &value, false);
2906  }
2907 
2908 
2909 
2910  inline void
2911  SparseMatrix::set(const std::vector<size_type> & indices,
2913  const bool elide_zero_values)
2914  {
2915  Assert(indices.size() == values.m(),
2916  ExcDimensionMismatch(indices.size(), values.m()));
2917  Assert(values.m() == values.n(), ExcNotQuadratic());
2918 
2919  for (size_type i = 0; i < indices.size(); ++i)
2920  set(indices[i],
2921  indices.size(),
2922  indices.data(),
2923  &values(i, 0),
2924  elide_zero_values);
2925  }
2926 
2927 
2928 
2929  inline void
2930  SparseMatrix::add(const size_type i,
2931  const size_type j,
2932  const TrilinosScalar value)
2933  {
2934  AssertIsFinite(value);
2935 
2936  if (value == 0)
2937  {
2938  // we have to check after Insert/Add in any case to be consistent
2939  // with the MPI communication model, but we can save some
2940  // work if the addend is zero. However, these actions are done in case
2941  // we pass on to the other function.
2942 
2943  // TODO: fix this (do not run compress here, but fail)
2944  if (last_action == Insert)
2945  {
2946  const int ierr = matrix->GlobalAssemble(*column_space_map,
2947  matrix->RowMap(),
2948  false);
2949 
2950  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2951  }
2952 
2953  last_action = Add;
2954 
2955  return;
2956  }
2957  else
2958  add(i, 1, &j, &value, false);
2959  }
2960 
2961 
2962 
2963  // inline "simple" functions that are called frequently and do only involve
2964  // a call to some Trilinos function.
2966  SparseMatrix::m() const
2967  {
2968 # ifndef DEAL_II_WITH_64BIT_INDICES
2969  return matrix->NumGlobalRows();
2970 # else
2971  return matrix->NumGlobalRows64();
2972 # endif
2973  }
2974 
2975 
2976 
2978  SparseMatrix::n() const
2979  {
2980  // If the matrix structure has not been fixed (i.e., we did not have a
2981  // sparsity pattern), it does not know about the number of columns so we
2982  // must always take this from the additional column space map
2983  Assert(column_space_map.get() != nullptr, ExcInternalError());
2984  return n_global_elements(*column_space_map);
2985  }
2986 
2987 
2988 
2989  inline unsigned int
2990  SparseMatrix::local_size() const
2991  {
2992  return matrix->NumMyRows();
2993  }
2994 
2995 
2996 
2997  inline std::pair<SparseMatrix::size_type, SparseMatrix::size_type>
2999  {
3000  size_type begin, end;
3001 # ifndef DEAL_II_WITH_64BIT_INDICES
3002  begin = matrix->RowMap().MinMyGID();
3003  end = matrix->RowMap().MaxMyGID() + 1;
3004 # else
3005  begin = matrix->RowMap().MinMyGID64();
3006  end = matrix->RowMap().MaxMyGID64() + 1;
3007 # endif
3008 
3009  return std::make_pair(begin, end);
3010  }
3011 
3012 
3013 
3014  inline std::uint64_t
3016  {
3017  // Trilinos uses 64bit functions internally for attribute access, which
3018  // return `long long`. They also offer 32bit variants that return `int`,
3019  // however those call the 64bit version and convert the values to 32bit.
3020  // There is no necessity in using the 32bit versions at all.
3021  return static_cast<std::uint64_t>(matrix->NumGlobalNonzeros64());
3022  }
3023 
3024 
3025 
3026  template <typename SparsityPatternType>
3027  inline void
3028  SparseMatrix::reinit(const IndexSet & parallel_partitioning,
3029  const SparsityPatternType &sparsity_pattern,
3030  const MPI_Comm & communicator,
3031  const bool exchange_data)
3032  {
3033  reinit(parallel_partitioning,
3034  parallel_partitioning,
3035  sparsity_pattern,
3036  communicator,
3037  exchange_data);
3038  }
3039 
3040 
3041 
3042  template <typename number>
3043  inline void
3044  SparseMatrix::reinit(const IndexSet &parallel_partitioning,
3045  const ::SparseMatrix<number> &sparse_matrix,
3046  const MPI_Comm & communicator,
3047  const double drop_tolerance,
3048  const bool copy_values,
3049  const ::SparsityPattern * use_this_sparsity)
3050  {
3051  Epetra_Map map =
3052  parallel_partitioning.make_trilinos_map(communicator, false);
3053  reinit(parallel_partitioning,
3054  parallel_partitioning,
3055  sparse_matrix,
3056  drop_tolerance,
3057  copy_values,
3058  use_this_sparsity);
3059  }
3060 
3061 
3062 
3063  inline const Epetra_CrsMatrix &
3065  {
3066  return static_cast<const Epetra_CrsMatrix &>(*matrix);
3067  }
3068 
3069 
3070 
3071  inline const Epetra_CrsGraph &
3073  {
3074  return matrix->Graph();
3075  }
3076 
3077 
3078 
3079  inline IndexSet
3081  {
3082  return IndexSet(matrix->DomainMap());
3083  }
3084 
3085 
3086 
3087  inline IndexSet
3089  {
3090  return IndexSet(matrix->RangeMap());
3091  }
3092 
3093 
3094 
3095  inline void
3097  {
3098  // nothing to do here
3099  }
3100 
3101 
3102 
3103  inline void
3105  {
3106  // nothing to do here
3107  }
3108 
3109 
3110  namespace internal
3111  {
3112  namespace LinearOperatorImplementation
3113  {
3114  template <typename EpetraOpType>
3115  TrilinosPayload::TrilinosPayload(
3116  EpetraOpType & op,
3117  const bool supports_inverse_operations,
3118  const bool use_transpose,
3119  const MPI_Comm &mpi_communicator,
3120  const IndexSet &locally_owned_domain_indices,
3121  const IndexSet &locally_owned_range_indices)
3122  : use_transpose(use_transpose)
3123  , communicator(mpi_communicator)
3124  , domain_map(
3125  locally_owned_domain_indices.make_trilinos_map(communicator.Comm()))
3126  , range_map(
3127  locally_owned_range_indices.make_trilinos_map(communicator.Comm()))
3128  {
3129  vmult = [&op](Range &tril_dst, const Domain &tril_src) {
3130  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3131  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3132  Assert(&tril_src != &tril_dst,
3135  tril_src,
3136  tril_dst,
3137  op.UseTranspose());
3138 
3139  const int ierr = op.Apply(tril_src, tril_dst);
3140  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3141  };
3142 
3143  Tvmult = [&op](Domain &tril_dst, const Range &tril_src) {
3144  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3145  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3146  Assert(&tril_src != &tril_dst,
3149  tril_src,
3150  tril_dst,
3151  !op.UseTranspose());
3152 
3153  op.SetUseTranspose(!op.UseTranspose());
3154  const int ierr = op.Apply(tril_src, tril_dst);
3155  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3156  op.SetUseTranspose(!op.UseTranspose());
3157  };
3158 
3159  if (supports_inverse_operations)
3160  {
3161  inv_vmult = [&op](Domain &tril_dst, const Range &tril_src) {
3162  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3163  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3164  Assert(
3165  &tril_src != &tril_dst,
3168  tril_src,
3169  tril_dst,
3170  !op.UseTranspose());
3171 
3172  const int ierr = op.ApplyInverse(tril_src, tril_dst);
3173  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3174  };
3175 
3176  inv_Tvmult = [&op](Range &tril_dst, const Domain &tril_src) {
3177  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3178  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3179  Assert(
3180  &tril_src != &tril_dst,
3183  tril_src,
3184  tril_dst,
3185  op.UseTranspose());
3186 
3187  op.SetUseTranspose(!op.UseTranspose());
3188  const int ierr = op.ApplyInverse(tril_src, tril_dst);
3189  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3190  op.SetUseTranspose(!op.UseTranspose());
3191  };
3192  }
3193  else
3194  {
3195  inv_vmult = [](Domain &, const Range &) {
3196  Assert(false,
3197  ExcMessage(
3198  "Uninitialized TrilinosPayload::inv_vmult called. "
3199  "The operator does not support inverse operations."));
3200  };
3201 
3202  inv_Tvmult = [](Range &, const Domain &) {
3203  Assert(false,
3204  ExcMessage(
3205  "Uninitialized TrilinosPayload::inv_Tvmult called. "
3206  "The operator does not support inverse operations."));
3207  };
3208  }
3209  }
3210 
3211 
3212  template <typename Solver, typename Preconditioner>
3213  typename std::enable_if<
3214  std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
3215  std::is_base_of<TrilinosWrappers::PreconditionBase,
3216  Preconditioner>::value,
3217  TrilinosPayload>::type
3219  Solver & solver,
3220  const Preconditioner &preconditioner) const
3221  {
3222  const auto &payload = *this;
3223 
3224  TrilinosPayload return_op(payload);
3225 
3226  // Capture by copy so the payloads are always valid
3227 
3228  return_op.inv_vmult = [payload, &solver, &preconditioner](
3229  TrilinosPayload::Domain & tril_dst,
3230  const TrilinosPayload::Range &tril_src) {
3231  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3232  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3233  Assert(&tril_src != &tril_dst,
3236  tril_src,
3237  tril_dst,
3238  !payload.UseTranspose());
3239  solver.solve(payload, tril_dst, tril_src, preconditioner);
3240  };
3241 
3242  return_op.inv_Tvmult = [payload, &solver, &preconditioner](
3243  TrilinosPayload::Range & tril_dst,
3244  const TrilinosPayload::Domain &tril_src) {
3245  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3246  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3247  Assert(&tril_src != &tril_dst,
3250  tril_src,
3251  tril_dst,
3252  payload.UseTranspose());
3253 
3254  const_cast<TrilinosPayload &>(payload).transpose();
3255  solver.solve(payload, tril_dst, tril_src, preconditioner);
3256  const_cast<TrilinosPayload &>(payload).transpose();
3257  };
3258 
3259  // If the input operator is already setup for transpose operations, then
3260  // we must do similar with its inverse.
3261  if (return_op.UseTranspose() == true)
3262  std::swap(return_op.inv_vmult, return_op.inv_Tvmult);
3263 
3264  return return_op;
3265  }
3266 
3267  template <typename Solver, typename Preconditioner>
3268  typename std::enable_if<
3269  !(std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
3270  std::is_base_of<TrilinosWrappers::PreconditionBase,
3271  Preconditioner>::value),
3272  TrilinosPayload>::type
3273  TrilinosPayload::inverse_payload(Solver &, const Preconditioner &) const
3274  {
3275  TrilinosPayload return_op(*this);
3276 
3277  return_op.inv_vmult = [](TrilinosPayload::Domain &,
3278  const TrilinosPayload::Range &) {
3279  AssertThrow(false,
3280  ExcMessage("Payload inv_vmult disabled because of "
3281  "incompatible solver/preconditioner choice."));
3282  };
3283 
3284  return_op.inv_Tvmult = [](TrilinosPayload::Range &,
3285  const TrilinosPayload::Domain &) {
3286  AssertThrow(false,
3287  ExcMessage("Payload inv_vmult disabled because of "
3288  "incompatible solver/preconditioner choice."));
3289  };
3290 
3291  return return_op;
3292  }
3293  } // namespace LinearOperatorImplementation
3294  } // namespace internal
3295 
3296  template <>
3297  void
3298  SparseMatrix::set<TrilinosScalar>(const size_type row,
3299  const size_type n_cols,
3300  const size_type * col_indices,
3301  const TrilinosScalar *values,
3302  const bool elide_zero_values);
3303 # endif // DOXYGEN
3304 
3305 } /* namespace TrilinosWrappers */
3306 
3307 
3309 
3310 
3311 # endif // DEAL_II_WITH_TRILINOS
3312 
3313 
3314 /*----------------------- trilinos_sparse_matrix.h --------------------*/
3315 
3316 #endif
3317 /*----------------------- 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:777
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
std::unique_ptr< Epetra_FECrsMatrix > matrix
const Epetra_CrsMatrix & trilinos_matrix() const
::types::global_dof_index size_type
std::enable_if<!std::is_same< typename VectorType::value_type, TrilinosScalar >::value >::type vmult(VectorType &dst, const VectorType &src) const
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< std::is_same< typename VectorType::value_type, TrilinosScalar >::value >::type vmult(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)
void vmult_add(VectorType &dst, const VectorType &src) const
const_iterator end(const size_type r) const
void compress(::VectorOperation::values operation)
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)
std::enable_if< std::is_same< typename VectorType::value_type, TrilinosScalar >::value >::type Tvmult(VectorType &dst, const VectorType &src) const
unsigned int row_length(const size_type row) const
std::enable_if<!std::is_same< typename VectorType::value_type, TrilinosScalar >::value >::type Tvmult(VectorType &dst, const VectorType &src) 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)
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
std::enable_if< std::is_base_of< TrilinosWrappers::SolverBase, Solver >::value &&std::is_base_of< TrilinosWrappers::PreconditionBase, Preconditioner >::value, TrilinosPayload >::type inverse_payload(Solver &, const Preconditioner &) const
virtual int Apply(const VectorType &X, VectorType &Y) const override
std::function< void(VectorType &, const VectorType &)> inv_vmult
std::enable_if< !(std::is_base_of< TrilinosWrappers::SolverBase, Solver >::value &&std::is_base_of< TrilinosWrappers::PreconditionBase, Preconditioner >::value), TrilinosPayload >::type inverse_payload(Solver &, const Preconditioner &) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
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:464
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:578
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:1473
#define AssertIsFinite(number)
Definition: exceptions.h:1758
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:532
static ::ExceptionBase & ExcAccessToNonLocalElement(size_type arg1, size_type arg2, size_type arg3, size_type arg4)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:487
Expression operator>(const Expression &lhs, const Expression &rhs)
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:555
static ::ExceptionBase & ExcMatrixNotCompressed()
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
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:1583
@ 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< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
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:76
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:163