Reference documentation for deal.II version Git 8f0f58a98a 2021-10-21 17:29:01 -0600
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
trilinos_sparse_matrix.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 2021 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_trilinos_sparse_matrix_h
17 # define dealii_trilinos_sparse_matrix_h
18 
19 
20 # include <deal.II/base/config.h>
21 
22 # ifdef DEAL_II_WITH_TRILINOS
23 
24 # include <deal.II/base/index_set.h>
26 
27 # include <deal.II/lac/exceptions.h>
28 # include <deal.II/lac/full_matrix.h>
35 
36 # include <Epetra_Comm.h>
37 # include <Epetra_CrsGraph.h>
38 # include <Epetra_Export.h>
39 # include <Epetra_FECrsMatrix.h>
40 # include <Epetra_Map.h>
41 # include <Epetra_MpiComm.h>
42 # include <Epetra_MultiVector.h>
43 # include <Epetra_Operator.h>
44 # include <mpi.h>
45 
46 # include <cmath>
47 # include <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
165  visit_present_row();
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 
232  Accessor(MatrixType *matrix, const size_type row, const size_type index);
233 
238  template <bool Other>
239  Accessor(const Accessor<Other> &a);
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 
321  Accessor(MatrixType *matrix, const size_type row, const size_type index);
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 
394  operator++();
395 
400  operator++(int);
401 
405  const Accessor<Constness> &
406  operator*() const;
407 
411  const Accessor<Constness> *
412  operator->() const;
413 
418  bool
419  operator==(const Iterator<Constness> &) const;
420 
424  bool
425  operator!=(const Iterator<Constness> &) const;
426 
432  bool
433  operator<(const Iterator<Constness> &) const;
434 
438  bool
439  operator>(const Iterator<Constness> &) const;
440 
444  DeclException2(ExcInvalidIndexWithinRow,
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;
475  using difference_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  size_type
951  n_nonzero_elements() const;
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 &
1666  trilinos_matrix() const;
1667 
1672  const Epetra_CrsGraph &
1673  trilinos_sparsity_pattern() const;
1674 
1676 
1681 
1686  IndexSet
1687  locally_owned_domain_indices() const;
1688 
1694  IndexSet
1695  locally_owned_range_indices() const;
1696 
1698 
1703 
1723  begin() const;
1724 
1728  iterator
1729  begin();
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 
1833  int,
1834  << "An error with error number " << arg1
1835  << " occurred while calling a Trilinos function");
1836 
1840  DeclException2(ExcInvalidIndex,
1841  size_type,
1842  size_type,
1843  << "The entry with index <" << arg1 << ',' << arg2
1844  << "> does not exist.");
1845 
1849  DeclExceptionMsg(ExcSourceEqualsDestination,
1850  "You are attempting an operation on two matrices that "
1851  "are the same object, but the operation requires that the "
1852  "two objects are in fact different.");
1853 
1857  DeclException0(ExcMatrixNotCompressed);
1858 
1862  DeclException4(ExcAccessToNonLocalElement,
1863  size_type,
1864  size_type,
1865  size_type,
1866  size_type,
1867  << "You tried to access element (" << arg1 << "/" << arg2
1868  << ")"
1869  << " of a distributed matrix, but only rows in range ["
1870  << arg3 << "," << arg4
1871  << "] are stored locally and can be accessed.");
1872 
1876  DeclException2(ExcAccessToNonPresentElement,
1877  size_type,
1878  size_type,
1879  << "You tried to access element (" << arg1 << "/" << arg2
1880  << ")"
1881  << " of a sparse matrix, but it appears to not"
1882  << " exist in the Trilinos sparsity pattern.");
1884 
1885 
1886 
1887  protected:
1898  void
1899  prepare_add();
1900 
1906  void
1907  prepare_set();
1908 
1909 
1910 
1911  private:
1916  std::unique_ptr<Epetra_Map> column_space_map;
1917 
1923  std::unique_ptr<Epetra_FECrsMatrix> matrix;
1924 
1930  std::unique_ptr<Epetra_CrsMatrix> nonlocal_matrix;
1931 
1935  std::unique_ptr<Epetra_Export> nonlocal_matrix_exporter;
1936 
1948  Epetra_CombineMode last_action;
1949 
1955 
1956  // To allow calling protected prepare_add() and prepare_set().
1958  };
1959 
1960 
1961 
1962  // forwards declarations
1963  class SolverBase;
1964  class PreconditionBase;
1965 
1966  namespace internal
1967  {
1968  inline void
1969  check_vector_map_equality(const Epetra_CrsMatrix & mtrx,
1970  const Epetra_MultiVector &src,
1971  const Epetra_MultiVector &dst,
1972  const bool transpose)
1973  {
1974  if (transpose == false)
1975  {
1976  Assert(src.Map().SameAs(mtrx.DomainMap()) == true,
1977  ExcMessage(
1978  "Column map of matrix does not fit with vector map!"));
1979  Assert(dst.Map().SameAs(mtrx.RangeMap()) == true,
1980  ExcMessage("Row map of matrix does not fit with vector map!"));
1981  }
1982  else
1983  {
1984  Assert(src.Map().SameAs(mtrx.RangeMap()) == true,
1985  ExcMessage(
1986  "Column map of matrix does not fit with vector map!"));
1987  Assert(dst.Map().SameAs(mtrx.DomainMap()) == true,
1988  ExcMessage("Row map of matrix does not fit with vector map!"));
1989  }
1990  (void)mtrx; // removes -Wunused-variable in optimized mode
1991  (void)src;
1992  (void)dst;
1993  }
1994 
1995  inline void
1997  const Epetra_MultiVector &src,
1998  const Epetra_MultiVector &dst,
1999  const bool transpose)
2000  {
2001  if (transpose == false)
2002  {
2003  Assert(src.Map().SameAs(op.OperatorDomainMap()) == true,
2004  ExcMessage(
2005  "Column map of operator does not fit with vector map!"));
2006  Assert(dst.Map().SameAs(op.OperatorRangeMap()) == true,
2007  ExcMessage(
2008  "Row map of operator does not fit with vector map!"));
2009  }
2010  else
2011  {
2012  Assert(src.Map().SameAs(op.OperatorRangeMap()) == true,
2013  ExcMessage(
2014  "Column map of operator does not fit with vector map!"));
2015  Assert(dst.Map().SameAs(op.OperatorDomainMap()) == true,
2016  ExcMessage(
2017  "Row map of operator does not fit with vector map!"));
2018  }
2019  (void)op; // removes -Wunused-variable in optimized mode
2020  (void)src;
2021  (void)dst;
2022  }
2023 
2024  namespace LinearOperatorImplementation
2025  {
2046  {
2047  public:
2051  using VectorType = Epetra_MultiVector;
2052 
2057 
2062 
2067 
2075  TrilinosPayload();
2076 
2080  TrilinosPayload(const TrilinosWrappers::SparseMatrix &matrix_exemplar,
2082 
2086  TrilinosPayload(const TrilinosPayload & payload_exemplar,
2088 
2093  const TrilinosWrappers::SparseMatrix & matrix_exemplar,
2094  const TrilinosWrappers::PreconditionBase &preconditioner);
2095 
2100  const TrilinosWrappers::PreconditionBase &preconditioner_exemplar,
2101  const TrilinosWrappers::PreconditionBase &preconditioner);
2102 
2107  const TrilinosPayload & payload_exemplar,
2108  const TrilinosWrappers::PreconditionBase &preconditioner);
2109 
2113  TrilinosPayload(const TrilinosPayload &payload);
2114 
2122  TrilinosPayload(const TrilinosPayload &first_op,
2123  const TrilinosPayload &second_op);
2124 
2128  virtual ~TrilinosPayload() override = default;
2129 
2134  identity_payload() const;
2135 
2140  null_payload() const;
2141 
2146  transpose_payload() const;
2147 
2164  template <typename Solver, typename Preconditioner>
2165  typename std::enable_if<
2166  std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
2167  std::is_base_of<TrilinosWrappers::PreconditionBase,
2168  Preconditioner>::value,
2169  TrilinosPayload>::type
2170  inverse_payload(Solver &, const Preconditioner &) const;
2171 
2189  template <typename Solver, typename Preconditioner>
2190  typename std::enable_if<
2191  !(std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
2192  std::is_base_of<TrilinosWrappers::PreconditionBase,
2193  Preconditioner>::value),
2194  TrilinosPayload>::type
2195  inverse_payload(Solver &, const Preconditioner &) const;
2196 
2198 
2203 
2209  IndexSet
2210  locally_owned_domain_indices() const;
2211 
2217  IndexSet
2218  locally_owned_range_indices() const;
2219 
2223  MPI_Comm
2224  get_mpi_communicator() const;
2225 
2232  void
2233  transpose();
2234 
2242  std::function<void(VectorType &, const VectorType &)> vmult;
2243 
2251  std::function<void(VectorType &, const VectorType &)> Tvmult;
2252 
2261  std::function<void(VectorType &, const VectorType &)> inv_vmult;
2262 
2271  std::function<void(VectorType &, const VectorType &)> inv_Tvmult;
2272 
2274 
2279 
2286  virtual bool
2287  UseTranspose() const override;
2288 
2304  virtual int
2305  SetUseTranspose(bool UseTranspose) override;
2306 
2318  virtual int
2319  Apply(const VectorType &X, VectorType &Y) const override;
2320 
2339  virtual int
2340  ApplyInverse(const VectorType &Y, VectorType &X) const override;
2342 
2347 
2354  virtual const char *
2355  Label() const override;
2356 
2364  virtual const Epetra_Comm &
2365  Comm() const override;
2366 
2374  virtual const Epetra_Map &
2375  OperatorDomainMap() const override;
2376 
2385  virtual const Epetra_Map &
2386  OperatorRangeMap() const override;
2388 
2389  private:
2399  template <typename EpetraOpType>
2400  TrilinosPayload(EpetraOpType & op,
2401  const bool supports_inverse_operations,
2402  const bool use_transpose,
2403  const MPI_Comm &mpi_communicator,
2404  const IndexSet &locally_owned_domain_indices,
2405  const IndexSet &locally_owned_range_indices);
2406 
2412 
2417  Epetra_MpiComm communicator;
2418 
2423  Epetra_Map domain_map;
2424 
2429  Epetra_Map range_map;
2430 
2439  virtual bool
2440  HasNormInf() const override;
2441 
2449  virtual double
2450  NormInf() const override;
2451  };
2452 
2458  operator+(const TrilinosPayload &first_op,
2459  const TrilinosPayload &second_op);
2460 
2466  operator*(const TrilinosPayload &first_op,
2467  const TrilinosPayload &second_op);
2468 
2469  } // namespace LinearOperatorImplementation
2470  } /* namespace internal */
2471 
2472 
2473 
2474  // ----------------------- inline and template functions --------------------
2475 
2476 # ifndef DOXYGEN
2477 
2478  namespace SparseMatrixIterators
2479  {
2480  inline AccessorBase::AccessorBase(SparseMatrix *matrix,
2481  size_type row,
2482  size_type index)
2483  : matrix(matrix)
2484  , a_row(row)
2485  , a_index(index)
2486  {
2487  visit_present_row();
2488  }
2489 
2490 
2492  AccessorBase::row() const
2493  {
2494  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2495  return a_row;
2496  }
2497 
2498 
2500  AccessorBase::column() const
2501  {
2502  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2503  return (*colnum_cache)[a_index];
2504  }
2505 
2506 
2508  AccessorBase::index() const
2509  {
2510  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2511  return a_index;
2512  }
2513 
2514 
2515  inline Accessor<true>::Accessor(MatrixType * matrix,
2516  const size_type row,
2517  const size_type index)
2518  : AccessorBase(const_cast<SparseMatrix *>(matrix), row, index)
2519  {}
2520 
2521 
2522  template <bool Other>
2523  inline Accessor<true>::Accessor(const Accessor<Other> &other)
2524  : AccessorBase(other)
2525  {}
2526 
2527 
2528  inline TrilinosScalar
2529  Accessor<true>::value() const
2530  {
2531  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2532  return (*value_cache)[a_index];
2533  }
2534 
2535 
2536  inline Accessor<false>::Reference::Reference(const Accessor<false> &acc)
2537  : accessor(const_cast<Accessor<false> &>(acc))
2538  {}
2539 
2540 
2541  inline Accessor<false>::Reference::operator TrilinosScalar() const
2542  {
2543  return (*accessor.value_cache)[accessor.a_index];
2544  }
2545 
2546  inline const Accessor<false>::Reference &
2547  Accessor<false>::Reference::operator=(const TrilinosScalar n) const
2548  {
2549  (*accessor.value_cache)[accessor.a_index] = n;
2550  accessor.matrix->set(accessor.row(),
2551  accessor.column(),
2552  static_cast<TrilinosScalar>(*this));
2553  return *this;
2554  }
2555 
2556 
2557  inline const Accessor<false>::Reference &
2558  Accessor<false>::Reference::operator+=(const TrilinosScalar n) const
2559  {
2560  (*accessor.value_cache)[accessor.a_index] += n;
2561  accessor.matrix->set(accessor.row(),
2562  accessor.column(),
2563  static_cast<TrilinosScalar>(*this));
2564  return *this;
2565  }
2566 
2567 
2568  inline const Accessor<false>::Reference &
2569  Accessor<false>::Reference::operator-=(const TrilinosScalar n) const
2570  {
2571  (*accessor.value_cache)[accessor.a_index] -= n;
2572  accessor.matrix->set(accessor.row(),
2573  accessor.column(),
2574  static_cast<TrilinosScalar>(*this));
2575  return *this;
2576  }
2577 
2578 
2579  inline const Accessor<false>::Reference &
2580  Accessor<false>::Reference::operator*=(const TrilinosScalar n) const
2581  {
2582  (*accessor.value_cache)[accessor.a_index] *= n;
2583  accessor.matrix->set(accessor.row(),
2584  accessor.column(),
2585  static_cast<TrilinosScalar>(*this));
2586  return *this;
2587  }
2588 
2589 
2590  inline const Accessor<false>::Reference &
2591  Accessor<false>::Reference::operator/=(const TrilinosScalar n) const
2592  {
2593  (*accessor.value_cache)[accessor.a_index] /= n;
2594  accessor.matrix->set(accessor.row(),
2595  accessor.column(),
2596  static_cast<TrilinosScalar>(*this));
2597  return *this;
2598  }
2599 
2600 
2601  inline Accessor<false>::Accessor(MatrixType * matrix,
2602  const size_type row,
2603  const size_type index)
2604  : AccessorBase(matrix, row, index)
2605  {}
2606 
2607 
2608  inline Accessor<false>::Reference
2609  Accessor<false>::value() const
2610  {
2611  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2612  return {*this};
2613  }
2614 
2615 
2616 
2617  template <bool Constness>
2618  inline Iterator<Constness>::Iterator(MatrixType * matrix,
2619  const size_type row,
2620  const size_type index)
2621  : accessor(matrix, row, index)
2622  {}
2623 
2624 
2625  template <bool Constness>
2626  template <bool Other>
2627  inline Iterator<Constness>::Iterator(const Iterator<Other> &other)
2628  : accessor(other.accessor)
2629  {}
2630 
2631 
2632  template <bool Constness>
2633  inline Iterator<Constness> &
2635  {
2636  Assert(accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
2637 
2638  ++accessor.a_index;
2639 
2640  // If at end of line: do one
2641  // step, then cycle until we
2642  // find a row with a nonzero
2643  // number of entries.
2644  if (accessor.a_index >= accessor.colnum_cache->size())
2645  {
2646  accessor.a_index = 0;
2647  ++accessor.a_row;
2648 
2649  while ((accessor.a_row < accessor.matrix->m()) &&
2650  ((accessor.matrix->in_local_range(accessor.a_row) == false) ||
2651  (accessor.matrix->row_length(accessor.a_row) == 0)))
2652  ++accessor.a_row;
2653 
2654  accessor.visit_present_row();
2655  }
2656  return *this;
2657  }
2658 
2659 
2660  template <bool Constness>
2661  inline Iterator<Constness>
2663  {
2664  const Iterator<Constness> old_state = *this;
2665  ++(*this);
2666  return old_state;
2667  }
2668 
2669 
2670 
2671  template <bool Constness>
2672  inline const Accessor<Constness> &
2674  {
2675  return accessor;
2676  }
2677 
2678 
2679 
2680  template <bool Constness>
2681  inline const Accessor<Constness> *
2682  Iterator<Constness>::operator->() const
2683  {
2684  return &accessor;
2685  }
2686 
2687 
2688 
2689  template <bool Constness>
2690  inline bool
2691  Iterator<Constness>::operator==(const Iterator<Constness> &other) const
2692  {
2693  return (accessor.a_row == other.accessor.a_row &&
2694  accessor.a_index == other.accessor.a_index);
2695  }
2696 
2697 
2698 
2699  template <bool Constness>
2700  inline bool
2701  Iterator<Constness>::operator!=(const Iterator<Constness> &other) const
2702  {
2703  return !(*this == other);
2704  }
2705 
2706 
2707 
2708  template <bool Constness>
2709  inline bool
2710  Iterator<Constness>::operator<(const Iterator<Constness> &other) const
2711  {
2712  return (accessor.row() < other.accessor.row() ||
2713  (accessor.row() == other.accessor.row() &&
2714  accessor.index() < other.accessor.index()));
2715  }
2716 
2717 
2718  template <bool Constness>
2719  inline bool
2720  Iterator<Constness>::operator>(const Iterator<Constness> &other) const
2721  {
2722  return (other < *this);
2723  }
2724 
2725  } // namespace SparseMatrixIterators
2726 
2727 
2728 
2730  SparseMatrix::begin() const
2731  {
2732  return begin(0);
2733  }
2734 
2735 
2736 
2738  SparseMatrix::end() const
2739  {
2740  return const_iterator(this, m(), 0);
2741  }
2742 
2743 
2744 
2746  SparseMatrix::begin(const size_type r) const
2747  {
2748  AssertIndexRange(r, m());
2749  if (in_local_range(r) && (row_length(r) > 0))
2750  return const_iterator(this, r, 0);
2751  else
2752  return end(r);
2753  }
2754 
2755 
2756 
2758  SparseMatrix::end(const size_type r) const
2759  {
2760  AssertIndexRange(r, m());
2761 
2762  // place the iterator on the first entry
2763  // past this line, or at the end of the
2764  // matrix
2765  for (size_type i = r + 1; i < m(); ++i)
2766  if (in_local_range(i) && (row_length(i) > 0))
2767  return const_iterator(this, i, 0);
2768 
2769  // if there is no such line, then take the
2770  // end iterator of the matrix
2771  return end();
2772  }
2773 
2774 
2775 
2776  inline SparseMatrix::iterator
2778  {
2779  return begin(0);
2780  }
2781 
2782 
2783 
2784  inline SparseMatrix::iterator
2786  {
2787  return iterator(this, m(), 0);
2788  }
2789 
2790 
2791 
2792  inline SparseMatrix::iterator
2794  {
2795  AssertIndexRange(r, m());
2796  if (in_local_range(r) && (row_length(r) > 0))
2797  return iterator(this, r, 0);
2798  else
2799  return end(r);
2800  }
2801 
2802 
2803 
2804  inline SparseMatrix::iterator
2805  SparseMatrix::end(const size_type r)
2806  {
2807  AssertIndexRange(r, m());
2808 
2809  // place the iterator on the first entry
2810  // past this line, or at the end of the
2811  // matrix
2812  for (size_type i = r + 1; i < m(); ++i)
2813  if (in_local_range(i) && (row_length(i) > 0))
2814  return iterator(this, i, 0);
2815 
2816  // if there is no such line, then take the
2817  // end iterator of the matrix
2818  return end();
2819  }
2820 
2821 
2822 
2823  inline bool
2824  SparseMatrix::in_local_range(const size_type index) const
2825  {
2827 # ifndef DEAL_II_WITH_64BIT_INDICES
2828  begin = matrix->RowMap().MinMyGID();
2829  end = matrix->RowMap().MaxMyGID() + 1;
2830 # else
2831  begin = matrix->RowMap().MinMyGID64();
2832  end = matrix->RowMap().MaxMyGID64() + 1;
2833 # endif
2834 
2835  return ((index >= static_cast<size_type>(begin)) &&
2836  (index < static_cast<size_type>(end)));
2837  }
2838 
2839 
2840 
2841  inline bool
2843  {
2844  return compressed;
2845  }
2846 
2847 
2848 
2849  // Inline the set() and add() functions, since they will be called
2850  // frequently, and the compiler can optimize away some unnecessary loops
2851  // when the sizes are given at compile time.
2852  template <>
2853  void
2854  SparseMatrix::set<TrilinosScalar>(const size_type row,
2855  const size_type n_cols,
2856  const size_type * col_indices,
2857  const TrilinosScalar *values,
2858  const bool elide_zero_values);
2859 
2860 
2861 
2862  template <typename Number>
2863  void
2864  SparseMatrix::set(const size_type row,
2865  const size_type n_cols,
2866  const size_type *col_indices,
2867  const Number * values,
2868  const bool elide_zero_values)
2869  {
2870  std::vector<TrilinosScalar> trilinos_values(n_cols);
2871  std::copy(values, values + n_cols, trilinos_values.begin());
2872  this->set(
2873  row, n_cols, col_indices, trilinos_values.data(), elide_zero_values);
2874  }
2875 
2876 
2877 
2878  inline void
2879  SparseMatrix::set(const size_type i,
2880  const size_type j,
2881  const TrilinosScalar value)
2882  {
2883  AssertIsFinite(value);
2884 
2885  set(i, 1, &j, &value, false);
2886  }
2887 
2888 
2889 
2890  inline void
2891  SparseMatrix::set(const std::vector<size_type> & indices,
2892  const FullMatrix<TrilinosScalar> &values,
2893  const bool elide_zero_values)
2894  {
2895  Assert(indices.size() == values.m(),
2896  ExcDimensionMismatch(indices.size(), values.m()));
2897  Assert(values.m() == values.n(), ExcNotQuadratic());
2898 
2899  for (size_type i = 0; i < indices.size(); ++i)
2900  set(indices[i],
2901  indices.size(),
2902  indices.data(),
2903  &values(i, 0),
2904  elide_zero_values);
2905  }
2906 
2907 
2908 
2909  inline void
2910  SparseMatrix::add(const size_type i,
2911  const size_type j,
2912  const TrilinosScalar value)
2913  {
2914  AssertIsFinite(value);
2915 
2916  if (value == 0)
2917  {
2918  // we have to check after Insert/Add in any case to be consistent
2919  // with the MPI communication model, but we can save some
2920  // work if the addend is zero. However, these actions are done in case
2921  // we pass on to the other function.
2922 
2923  // TODO: fix this (do not run compress here, but fail)
2924  if (last_action == Insert)
2925  {
2926 # ifdef DEBUG
2927  int ierr;
2928  ierr =
2929 # endif
2930  matrix->GlobalAssemble(*column_space_map,
2931  matrix->RowMap(),
2932  false);
2933 
2934  Assert(ierr == 0, ExcTrilinosError(ierr));
2935  }
2936 
2937  last_action = Add;
2938 
2939  return;
2940  }
2941  else
2942  add(i, 1, &j, &value, false);
2943  }
2944 
2945 
2946 
2947  // inline "simple" functions that are called frequently and do only involve
2948  // a call to some Trilinos function.
2950  SparseMatrix::m() const
2951  {
2952 # ifndef DEAL_II_WITH_64BIT_INDICES
2953  return matrix->NumGlobalRows();
2954 # else
2955  return matrix->NumGlobalRows64();
2956 # endif
2957  }
2958 
2959 
2960 
2962  SparseMatrix::n() const
2963  {
2964  // If the matrix structure has not been fixed (i.e., we did not have a
2965  // sparsity pattern), it does not know about the number of columns so we
2966  // must always take this from the additional column space map
2967  Assert(column_space_map.get() != nullptr, ExcInternalError());
2968  return n_global_elements(*column_space_map);
2969  }
2970 
2971 
2972 
2973  inline unsigned int
2974  SparseMatrix::local_size() const
2975  {
2976  return matrix->NumMyRows();
2977  }
2978 
2979 
2980 
2981  inline std::pair<SparseMatrix::size_type, SparseMatrix::size_type>
2983  {
2984  size_type begin, end;
2985 # ifndef DEAL_II_WITH_64BIT_INDICES
2986  begin = matrix->RowMap().MinMyGID();
2987  end = matrix->RowMap().MaxMyGID() + 1;
2988 # else
2989  begin = matrix->RowMap().MinMyGID64();
2990  end = matrix->RowMap().MaxMyGID64() + 1;
2991 # endif
2992 
2993  return std::make_pair(begin, end);
2994  }
2995 
2996 
2997 
3000  {
3001 # ifndef DEAL_II_WITH_64BIT_INDICES
3002  return matrix->NumGlobalNonzeros();
3003 # else
3004  return matrix->NumGlobalNonzeros64();
3005 # endif
3006  }
3007 
3008 
3009 
3010  template <typename SparsityPatternType>
3011  inline void
3012  SparseMatrix::reinit(const IndexSet & parallel_partitioning,
3013  const SparsityPatternType &sparsity_pattern,
3014  const MPI_Comm & communicator,
3015  const bool exchange_data)
3016  {
3017  reinit(parallel_partitioning,
3018  parallel_partitioning,
3019  sparsity_pattern,
3020  communicator,
3021  exchange_data);
3022  }
3023 
3024 
3025 
3026  template <typename number>
3027  inline void
3028  SparseMatrix::reinit(const IndexSet &parallel_partitioning,
3029  const ::SparseMatrix<number> &sparse_matrix,
3030  const MPI_Comm & communicator,
3031  const double drop_tolerance,
3032  const bool copy_values,
3033  const ::SparsityPattern * use_this_sparsity)
3034  {
3035  Epetra_Map map =
3036  parallel_partitioning.make_trilinos_map(communicator, false);
3037  reinit(parallel_partitioning,
3038  parallel_partitioning,
3039  sparse_matrix,
3040  drop_tolerance,
3041  copy_values,
3042  use_this_sparsity);
3043  }
3044 
3045 
3046 
3047  inline const Epetra_CrsMatrix &
3049  {
3050  return static_cast<const Epetra_CrsMatrix &>(*matrix);
3051  }
3052 
3053 
3054 
3055  inline const Epetra_CrsGraph &
3057  {
3058  return matrix->Graph();
3059  }
3060 
3061 
3062 
3063  inline IndexSet
3065  {
3066  return IndexSet(matrix->DomainMap());
3067  }
3068 
3069 
3070 
3071  inline IndexSet
3073  {
3074  return IndexSet(matrix->RangeMap());
3075  }
3076 
3077 
3078 
3079  inline void
3081  {
3082  // nothing to do here
3083  }
3084 
3085 
3086 
3087  inline void
3089  {
3090  // nothing to do here
3091  }
3092 
3093 
3094  namespace internal
3095  {
3096  namespace LinearOperatorImplementation
3097  {
3098  template <typename EpetraOpType>
3100  EpetraOpType & op,
3101  const bool supports_inverse_operations,
3102  const bool use_transpose,
3103  const MPI_Comm &mpi_communicator,
3106  : use_transpose(use_transpose)
3107  , communicator(mpi_communicator)
3108  , domain_map(
3109  locally_owned_domain_indices.make_trilinos_map(communicator.Comm()))
3110  , range_map(
3111  locally_owned_range_indices.make_trilinos_map(communicator.Comm()))
3112  {
3113  vmult = [&op](Range &tril_dst, const Domain &tril_src) {
3114  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3115  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3116  Assert(&tril_src != &tril_dst,
3119  tril_src,
3120  tril_dst,
3121  op.UseTranspose());
3122 
3123  const int ierr = op.Apply(tril_src, tril_dst);
3124  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3125  };
3126 
3127  Tvmult = [&op](Domain &tril_dst, const Range &tril_src) {
3128  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3129  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3130  Assert(&tril_src != &tril_dst,
3133  tril_src,
3134  tril_dst,
3135  !op.UseTranspose());
3136 
3137  op.SetUseTranspose(!op.UseTranspose());
3138  const int ierr = op.Apply(tril_src, tril_dst);
3139  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3140  op.SetUseTranspose(!op.UseTranspose());
3141  };
3142 
3143  if (supports_inverse_operations)
3144  {
3145  inv_vmult = [&op](Domain &tril_dst, const Range &tril_src) {
3146  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3147  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3148  Assert(
3149  &tril_src != &tril_dst,
3152  tril_src,
3153  tril_dst,
3154  !op.UseTranspose());
3155 
3156  const int ierr = op.ApplyInverse(tril_src, tril_dst);
3157  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3158  };
3159 
3160  inv_Tvmult = [&op](Range &tril_dst, const Domain &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  op.SetUseTranspose(!op.UseTranspose());
3172  const int ierr = op.ApplyInverse(tril_src, tril_dst);
3173  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3174  op.SetUseTranspose(!op.UseTranspose());
3175  };
3176  }
3177  else
3178  {
3179  inv_vmult = [](Domain &, const Range &) {
3180  Assert(false,
3181  ExcMessage(
3182  "Uninitialized TrilinosPayload::inv_vmult called. "
3183  "The operator does not support inverse operations."));
3184  };
3185 
3186  inv_Tvmult = [](Range &, const Domain &) {
3187  Assert(false,
3188  ExcMessage(
3189  "Uninitialized TrilinosPayload::inv_Tvmult called. "
3190  "The operator does not support inverse operations."));
3191  };
3192  }
3193  }
3194 
3195 
3196  template <typename Solver, typename Preconditioner>
3197  typename std::enable_if<
3198  std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
3199  std::is_base_of<TrilinosWrappers::PreconditionBase,
3200  Preconditioner>::value,
3201  TrilinosPayload>::type
3203  Solver & solver,
3204  const Preconditioner &preconditioner) const
3205  {
3206  const auto &payload = *this;
3207 
3208  TrilinosPayload return_op(payload);
3209 
3210  // Capture by copy so the payloads are always valid
3211 
3212  return_op.inv_vmult = [payload, &solver, &preconditioner](
3213  TrilinosPayload::Domain & tril_dst,
3214  const TrilinosPayload::Range &tril_src) {
3215  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3216  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3217  Assert(&tril_src != &tril_dst,
3220  tril_src,
3221  tril_dst,
3222  !payload.UseTranspose());
3223  solver.solve(payload, tril_dst, tril_src, preconditioner);
3224  };
3225 
3226  return_op.inv_Tvmult = [payload, &solver, &preconditioner](
3227  TrilinosPayload::Range & tril_dst,
3228  const TrilinosPayload::Domain &tril_src) {
3229  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3230  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3231  Assert(&tril_src != &tril_dst,
3234  tril_src,
3235  tril_dst,
3236  payload.UseTranspose());
3237 
3238  const_cast<TrilinosPayload &>(payload).transpose();
3239  solver.solve(payload, tril_dst, tril_src, preconditioner);
3240  const_cast<TrilinosPayload &>(payload).transpose();
3241  };
3242 
3243  // If the input operator is already setup for transpose operations, then
3244  // we must do similar with its inverse.
3245  if (return_op.UseTranspose() == true)
3246  std::swap(return_op.inv_vmult, return_op.inv_Tvmult);
3247 
3248  return return_op;
3249  }
3250 
3251  template <typename Solver, typename Preconditioner>
3252  typename std::enable_if<
3253  !(std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
3254  std::is_base_of<TrilinosWrappers::PreconditionBase,
3255  Preconditioner>::value),
3256  TrilinosPayload>::type
3257  TrilinosPayload::inverse_payload(Solver &, const Preconditioner &) const
3258  {
3259  TrilinosPayload return_op(*this);
3260 
3261  return_op.inv_vmult = [](TrilinosPayload::Domain &,
3262  const TrilinosPayload::Range &) {
3263  AssertThrow(false,
3264  ExcMessage("Payload inv_vmult disabled because of "
3265  "incompatible solver/preconditioner choice."));
3266  };
3267 
3268  return_op.inv_Tvmult = [](TrilinosPayload::Range &,
3269  const TrilinosPayload::Domain &) {
3270  AssertThrow(false,
3271  ExcMessage("Payload inv_vmult disabled because of "
3272  "incompatible solver/preconditioner choice."));
3273  };
3274 
3275  return return_op;
3276  }
3277  } // namespace LinearOperatorImplementation
3278  } // namespace internal
3279 
3280  template <>
3281  void
3282  SparseMatrix::set<TrilinosScalar>(const size_type row,
3283  const size_type n_cols,
3284  const size_type * col_indices,
3285  const TrilinosScalar *values,
3286  const bool elide_zero_values);
3287 # endif // DOXYGEN
3288 
3289 } /* namespace TrilinosWrappers */
3290 
3291 
3293 
3294 
3295 # endif // DEAL_II_WITH_TRILINOS
3296 
3297 
3298 /*----------------------- trilinos_sparse_matrix.h --------------------*/
3299 
3300 #endif
3301 /*----------------------- trilinos_sparse_matrix.h --------------------*/
static ::ExceptionBase & ExcTrilinosError(int arg1)
size_type m() const
const_iterator end() const
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:532
Contents is actually a matrix.
static ::ExceptionBase & ExcSourceEqualsDestination()
TrilinosPayload operator*(const TrilinosPayload &first_op, const TrilinosPayload &second_op)
types::global_dof_index size_type
Definition: cuda_kernels.h:45
static ::ExceptionBase & ExcBeyondEndOfMatrix()
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
std::function< void(VectorType &, const VectorType &)> inv_Tvmult
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1720
static const char V
TrilinosPayload operator+(const TrilinosPayload &first_op, const TrilinosPayload &second_op)
STL namespace.
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)
static ::ExceptionBase & ExcAccessToNonlocalRow(std::size_t arg1, std::size_t arg2, std::size_t arg3)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1571
typename ::TrilinosWrappers::SparseMatrixIterators::Iterator< Constness >::value_type value_type
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:604
std::shared_ptr< std::vector< size_type > > colnum_cache
std::unique_ptr< Epetra_FECrsMatrix > matrix
::types::global_dof_index size_type
bool in_local_range(const size_type index) const
static ::ExceptionBase & ExcMessage(std::string arg1)
std::string compress(const std::string &input)
Definition: utilities.cc:392
size_type n() const
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
void add(const size_type i, const size_type j, const TrilinosScalar value)
#define Assert(cond, exc)
Definition: exceptions.h:1461
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
std::unique_ptr< Epetra_Export > nonlocal_matrix_exporter
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:487
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:3025
#define DeclException0(Exception0)
Definition: exceptions.h:464
typename Accessor< Constness >::MatrixType MatrixType
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
std::function< void(VectorType &, const VectorType &)> inv_vmult
VectorType::value_type * end(VectorType &V)
double TrilinosScalar
Definition: types.h:163
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
std::enable_if< std::is_base_of< TrilinosWrappers::SolverBase, Solver >::value &&std::is_base_of< TrilinosWrappers::PreconditionBase, Preconditioner >::value, TrilinosPayload >::type inverse_payload(Solver &, const Preconditioner &) const
Expression operator>(const Expression &lhs, const Expression &rhs)
const Epetra_CrsMatrix & trilinos_matrix() const
AccessorBase(SparseMatrix *matrix, const size_type row, const size_type index)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SparseMatrix< double > SparseMatrix
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
IndexSet locally_owned_range_indices() const
size_type n_nonzero_elements() const
static ::ExceptionBase & ExcIteratorPastEnd()
const_iterator begin() const
static ::ExceptionBase & ExcNotQuadratic()
TrilinosWrappers::types::int_type n_global_elements(const Epetra_BlockMap &map)
unsigned int global_dof_index
Definition: types.h:76
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
void reinit(const SparsityPatternType &sparsity_pattern)
typename ::TrilinosWrappers::SparseMatrixIterators::Iterator< Constness >::difference_type difference_type
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
VectorType::value_type * begin(VectorType &V)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:578
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2999
std::pair< size_type, size_type > local_range() const
IndexSet locally_owned_domain_indices() const
const Epetra_CrsGraph & trilinos_sparsity_pattern() const
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:555
std::unique_ptr< Epetra_Map > column_space_map
std::unique_ptr< Epetra_CrsMatrix > nonlocal_matrix
void set(const size_type i, const size_type j, const TrilinosScalar value)
void check_vector_map_equality(const Epetra_Operator &op, const Epetra_MultiVector &src, const Epetra_MultiVector &dst, const bool transpose)
void copy(const T *begin, const T *end, U *dest)
#define AssertIsFinite(number)
Definition: exceptions.h:1746
unsigned int local_size() const
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
void check_vector_map_equality(const Epetra_CrsMatrix &mtrx, const Epetra_MultiVector &src, const Epetra_MultiVector &dst, const bool transpose)
static ::ExceptionBase & ExcInternalError()
std::shared_ptr< std::vector< TrilinosScalar > > value_cache