Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
trilinos_sparse_matrix.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_trilinos_sparse_matrix_h
17 # define dealii_trilinos_sparse_matrix_h
18 
19 
20 # include <deal.II/base/config.h>
21 
22 # ifdef DEAL_II_WITH_TRILINOS
23 
24 # include <deal.II/base/index_set.h>
26 
27 # include <deal.II/lac/exceptions.h>
28 # include <deal.II/lac/full_matrix.h>
35 
36 # include <Epetra_Comm.h>
37 # include <Epetra_CrsGraph.h>
38 # include <Epetra_Export.h>
39 # include <Epetra_FECrsMatrix.h>
40 # include <Epetra_Map.h>
41 # include <Epetra_MultiVector.h>
42 # include <Epetra_Operator.h>
43 
44 # include <cmath>
45 # include <memory>
46 # include <type_traits>
47 # include <vector>
48 # ifdef DEAL_II_WITH_MPI
49 # include <Epetra_MpiComm.h>
50 # include <mpi.h>
51 # else
52 # include <Epetra_SerialComm.h>
53 # endif
54 
56 
57 // forward declarations
58 # ifndef DOXYGEN
59 template <typename MatrixType>
60 class BlockMatrixBase;
61 
62 template <typename number>
63 class SparseMatrix;
64 class SparsityPattern;
66 
67 namespace TrilinosWrappers
68 {
69  class SparseMatrix;
70  class SparsityPattern;
71 
72  namespace SparseMatrixIterators
73  {
74  template <bool Constness>
75  class Iterator;
76  }
77 } // namespace TrilinosWrappers
78 # endif
79 
80 namespace TrilinosWrappers
81 {
86  {
91 
96  std::size_t,
97  std::size_t,
98  std::size_t,
99  << "You tried to access row " << arg1
100  << " of a distributed sparsity pattern, "
101  << " but only rows " << arg2 << " through " << arg3
102  << " are stored locally and can be accessed.");
103 
115  {
116  public:
121 
126  const size_type row,
127  const size_type index);
128 
132  size_type
133  row() const;
134 
138  size_type
139  index() const;
140 
144  size_type
145  column() const;
146 
147  protected:
159 
164 
170  void
171  visit_present_row();
172 
185  std::shared_ptr<std::vector<size_type>> colnum_cache;
186 
190  std::shared_ptr<std::vector<TrilinosScalar>> value_cache;
191  };
192 
203  template <bool Constess>
204  class Accessor : public AccessorBase
205  {
210  value() const;
211 
216  value();
217  };
218 
222  template <>
223  class Accessor<true> : public AccessorBase
224  {
225  public:
230  using MatrixType = const SparseMatrix;
231 
236  Accessor(MatrixType *matrix, const size_type row, const size_type index);
237 
242  template <bool Other>
243  Accessor(const Accessor<Other> &a);
244 
249  value() const;
250 
251  private:
252  // Make iterator class a friend.
253  template <bool>
254  friend class Iterator;
255  };
256 
260  template <>
261  class Accessor<false> : public AccessorBase
262  {
263  class Reference
264  {
265  public:
269  Reference(const Accessor<false> &accessor);
270 
274  operator TrilinosScalar() const;
275 
279  const Reference &
280  operator=(const TrilinosScalar n) const;
281 
285  const Reference &
286  operator+=(const TrilinosScalar n) const;
287 
291  const Reference &
292  operator-=(const TrilinosScalar n) const;
293 
297  const Reference &
298  operator*=(const TrilinosScalar n) const;
299 
303  const Reference &
304  operator/=(const TrilinosScalar n) const;
305 
306  private:
312  };
313 
314  public:
320 
325  Accessor(MatrixType *matrix, const size_type row, const size_type index);
326 
330  Reference
331  value() const;
332 
333  private:
334  // Make iterator class a friend.
335  template <bool>
336  friend class Iterator;
337 
338  // Make Reference object a friend.
339  friend class Reference;
340  };
341 
356  template <bool Constness>
357  class Iterator
358  {
359  public:
364 
370 
375  Iterator(MatrixType *matrix, const size_type row, const size_type index);
376 
380  template <bool Other>
381  Iterator(const Iterator<Other> &other);
382 
387  operator++();
388 
393  operator++(int);
394 
398  const Accessor<Constness> &operator*() const;
399 
403  const Accessor<Constness> *operator->() const;
404 
409  bool
410  operator==(const Iterator<Constness> &) const;
411 
415  bool
416  operator!=(const Iterator<Constness> &) const;
417 
423  bool
424  operator<(const Iterator<Constness> &) const;
425 
429  bool
430  operator>(const Iterator<Constness> &) const;
431 
435  DeclException2(ExcInvalidIndexWithinRow,
436  size_type,
437  size_type,
438  << "Attempt to access element " << arg2 << " of row "
439  << arg1 << " which doesn't have that many elements.");
440 
441  private:
446 
447  template <bool Other>
448  friend class Iterator;
449  };
450 
451  } // namespace SparseMatrixIterators
452 
453 
515  class SparseMatrix : public Subscriptor
516  {
517  public:
522 
527  std::size_t,
528  << "You tried to access row " << arg1
529  << " of a non-contiguous locally owned row set."
530  << " The row " << arg1
531  << " is not stored locally and can't be accessed.");
532 
540  struct Traits
541  {
546  static const bool zero_addition_can_be_elided = true;
547  };
548 
553 
558 
563 
571  SparseMatrix();
572 
580  SparseMatrix(const size_type m,
581  const size_type n,
582  const unsigned int n_max_entries_per_row);
583 
591  SparseMatrix(const size_type m,
592  const size_type n,
593  const std::vector<unsigned int> &n_entries_per_row);
594 
598  SparseMatrix(const SparsityPattern &InputSparsityPattern);
599 
604  SparseMatrix(SparseMatrix &&other) noexcept;
605 
609  SparseMatrix(const SparseMatrix &) = delete;
610 
614  SparseMatrix &
615  operator=(const SparseMatrix &) = delete;
616 
620  virtual ~SparseMatrix() override = default;
621 
637  template <typename SparsityPatternType>
638  void
639  reinit(const SparsityPatternType &sparsity_pattern);
640 
653  void
654  reinit(const SparsityPattern &sparsity_pattern);
655 
664  void
665  reinit(const SparseMatrix &sparse_matrix);
666 
687  template <typename number>
688  void
689  reinit(const ::SparseMatrix<number> &dealii_sparse_matrix,
690  const double drop_tolerance = 1e-13,
691  const bool copy_values = true,
692  const ::SparsityPattern * use_this_sparsity = nullptr);
693 
699  void
700  reinit(const Epetra_CrsMatrix &input_matrix, const bool copy_values = true);
702 
719  SparseMatrix(const IndexSet & parallel_partitioning,
720  const MPI_Comm & communicator = MPI_COMM_WORLD,
721  const unsigned int n_max_entries_per_row = 0);
722 
730  SparseMatrix(const IndexSet & parallel_partitioning,
731  const MPI_Comm & communicator,
732  const std::vector<unsigned int> &n_entries_per_row);
733 
748  SparseMatrix(const IndexSet &row_parallel_partitioning,
749  const IndexSet &col_parallel_partitioning,
750  const MPI_Comm &communicator = MPI_COMM_WORLD,
751  const size_type n_max_entries_per_row = 0);
752 
767  SparseMatrix(const IndexSet & row_parallel_partitioning,
768  const IndexSet & col_parallel_partitioning,
769  const MPI_Comm & communicator,
770  const std::vector<unsigned int> &n_entries_per_row);
771 
792  template <typename SparsityPatternType>
793  void
794  reinit(const IndexSet & parallel_partitioning,
795  const SparsityPatternType &sparsity_pattern,
796  const MPI_Comm & communicator = MPI_COMM_WORLD,
797  const bool exchange_data = false);
798 
811  template <typename SparsityPatternType>
812  typename std::enable_if<
813  !std::is_same<SparsityPatternType,
815  reinit(const IndexSet & row_parallel_partitioning,
816  const IndexSet & col_parallel_partitioning,
817  const SparsityPatternType &sparsity_pattern,
818  const MPI_Comm & communicator = MPI_COMM_WORLD,
819  const bool exchange_data = false);
820 
837  template <typename number>
838  void
839  reinit(const IndexSet & parallel_partitioning,
840  const ::SparseMatrix<number> &dealii_sparse_matrix,
841  const MPI_Comm & communicator = MPI_COMM_WORLD,
842  const double drop_tolerance = 1e-13,
843  const bool copy_values = true,
844  const ::SparsityPattern * use_this_sparsity = nullptr);
845 
859  template <typename number>
860  void
861  reinit(const IndexSet & row_parallel_partitioning,
862  const IndexSet & col_parallel_partitioning,
863  const ::SparseMatrix<number> &dealii_sparse_matrix,
864  const MPI_Comm & communicator = MPI_COMM_WORLD,
865  const double drop_tolerance = 1e-13,
866  const bool copy_values = true,
867  const ::SparsityPattern * use_this_sparsity = nullptr);
869 
873 
877  size_type
878  m() const;
879 
883  size_type
884  n() const;
885 
894  unsigned int
895  local_size() const;
896 
905  std::pair<size_type, size_type>
906  local_range() const;
907 
912  bool
913  in_local_range(const size_type index) const;
914 
919  size_type
920  n_nonzero_elements() const;
921 
925  unsigned int
926  row_length(const size_type row) const;
927 
934  bool
935  is_compressed() const;
936 
942  size_type
943  memory_consumption() const;
944 
948  MPI_Comm
949  get_mpi_communicator() const;
950 
952 
956 
966  SparseMatrix &
967  operator=(const double d);
968 
976  void
977  clear();
978 
1006  void
1007  compress(::VectorOperation::values operation);
1008 
1030  void
1031  set(const size_type i, const size_type j, const TrilinosScalar value);
1032 
1065  void
1066  set(const std::vector<size_type> & indices,
1067  const FullMatrix<TrilinosScalar> &full_matrix,
1068  const bool elide_zero_values = false);
1069 
1075  void
1076  set(const std::vector<size_type> & row_indices,
1077  const std::vector<size_type> & col_indices,
1078  const FullMatrix<TrilinosScalar> &full_matrix,
1079  const bool elide_zero_values = false);
1080 
1108  void
1109  set(const size_type row,
1110  const std::vector<size_type> & col_indices,
1111  const std::vector<TrilinosScalar> &values,
1112  const bool elide_zero_values = false);
1113 
1141  template <typename Number>
1142  void
1143  set(const size_type row,
1144  const size_type n_cols,
1145  const size_type *col_indices,
1146  const Number * values,
1147  const bool elide_zero_values = false);
1148 
1158  void
1159  add(const size_type i, const size_type j, const TrilinosScalar value);
1160 
1179  void
1180  add(const std::vector<size_type> & indices,
1181  const FullMatrix<TrilinosScalar> &full_matrix,
1182  const bool elide_zero_values = true);
1183 
1189  void
1190  add(const std::vector<size_type> & row_indices,
1191  const std::vector<size_type> & col_indices,
1192  const FullMatrix<TrilinosScalar> &full_matrix,
1193  const bool elide_zero_values = true);
1194 
1208  void
1209  add(const size_type row,
1210  const std::vector<size_type> & col_indices,
1211  const std::vector<TrilinosScalar> &values,
1212  const bool elide_zero_values = true);
1213 
1227  void
1228  add(const size_type row,
1229  const size_type n_cols,
1230  const size_type * col_indices,
1231  const TrilinosScalar *values,
1232  const bool elide_zero_values = true,
1233  const bool col_indices_are_sorted = false);
1234 
1238  SparseMatrix &
1239  operator*=(const TrilinosScalar factor);
1240 
1244  SparseMatrix &
1245  operator/=(const TrilinosScalar factor);
1246 
1250  void
1251  copy_from(const SparseMatrix &source);
1252 
1260  void
1261  add(const TrilinosScalar factor, const SparseMatrix &matrix);
1262 
1289  void
1290  clear_row(const size_type row, const TrilinosScalar new_diag_value = 0);
1291 
1312  void
1313  clear_rows(const std::vector<size_type> &rows,
1314  const TrilinosScalar new_diag_value = 0);
1315 
1325  void
1326  transpose();
1327 
1329 
1333 
1343  operator()(const size_type i, const size_type j) const;
1344 
1362  el(const size_type i, const size_type j) const;
1363 
1371  diag_element(const size_type i) const;
1372 
1374 
1378 
1409  template <typename VectorType>
1410  typename std::enable_if<std::is_same<typename VectorType::value_type,
1411  TrilinosScalar>::value>::type
1412  vmult(VectorType &dst, const VectorType &src) const;
1413 
1420  template <typename VectorType>
1421  typename std::enable_if<!std::is_same<typename VectorType::value_type,
1422  TrilinosScalar>::value>::type
1423  vmult(VectorType &dst, const VectorType &src) const;
1424 
1439  template <typename VectorType>
1440  typename std::enable_if<std::is_same<typename VectorType::value_type,
1441  TrilinosScalar>::value>::type
1442  Tvmult(VectorType &dst, const VectorType &src) const;
1443 
1450  template <typename VectorType>
1451  typename std::enable_if<!std::is_same<typename VectorType::value_type,
1452  TrilinosScalar>::value>::type
1453  Tvmult(VectorType &dst, const VectorType &src) const;
1454 
1464  template <typename VectorType>
1465  void
1466  vmult_add(VectorType &dst, const VectorType &src) const;
1467 
1478  template <typename VectorType>
1479  void
1480  Tvmult_add(VectorType &dst, const VectorType &src) const;
1481 
1504  matrix_norm_square(const MPI::Vector &v) const;
1505 
1526  matrix_scalar_product(const MPI::Vector &u, const MPI::Vector &v) const;
1527 
1545  residual(MPI::Vector & dst,
1546  const MPI::Vector &x,
1547  const MPI::Vector &b) const;
1548 
1563  void
1564  mmult(SparseMatrix & C,
1565  const SparseMatrix &B,
1566  const MPI::Vector & V = MPI::Vector()) const;
1567 
1568 
1585  void
1586  Tmmult(SparseMatrix & C,
1587  const SparseMatrix &B,
1588  const MPI::Vector & V = MPI::Vector()) const;
1589 
1591 
1595 
1604  l1_norm() const;
1605 
1615  linfty_norm() const;
1616 
1622  frobenius_norm() const;
1623 
1625 
1629 
1634  const Epetra_CrsMatrix &
1635  trilinos_matrix() const;
1636 
1641  const Epetra_CrsGraph &
1642  trilinos_sparsity_pattern() const;
1643 
1645 
1650 
1655  IndexSet
1656  locally_owned_domain_indices() const;
1657 
1663  IndexSet
1664  locally_owned_range_indices() const;
1665 
1667 
1672 
1692  begin() const;
1693 
1697  iterator
1698  begin();
1699 
1705  end() const;
1706 
1710  iterator
1711  end();
1712 
1742  begin(const size_type r) const;
1743 
1747  iterator
1748  begin(const size_type r);
1749 
1760  end(const size_type r) const;
1761 
1765  iterator
1766  end(const size_type r);
1767 
1769 
1773 
1779  void
1780  write_ascii();
1781 
1789  void
1790  print(std::ostream &out,
1791  const bool write_extended_trilinos_info = false) const;
1792 
1794 
1803  int,
1804  << "An error with error number " << arg1
1805  << " occurred while calling a Trilinos function");
1806 
1810  DeclException2(ExcInvalidIndex,
1811  size_type,
1812  size_type,
1813  << "The entry with index <" << arg1 << ',' << arg2
1814  << "> does not exist.");
1815 
1819  DeclExceptionMsg(ExcSourceEqualsDestination,
1820  "You are attempting an operation on two matrices that "
1821  "are the same object, but the operation requires that the "
1822  "two objects are in fact different.");
1823 
1827  DeclException0(ExcMatrixNotCompressed);
1828 
1832  DeclException4(ExcAccessToNonLocalElement,
1833  size_type,
1834  size_type,
1835  size_type,
1836  size_type,
1837  << "You tried to access element (" << arg1 << "/" << arg2
1838  << ")"
1839  << " of a distributed matrix, but only rows " << arg3
1840  << " through " << arg4
1841  << " are stored locally and can be accessed.");
1842 
1846  DeclException2(ExcAccessToNonPresentElement,
1847  size_type,
1848  size_type,
1849  << "You tried to access element (" << arg1 << "/" << arg2
1850  << ")"
1851  << " of a sparse matrix, but it appears to not"
1852  << " exist in the Trilinos sparsity pattern.");
1854 
1855 
1856 
1857  protected:
1868  void
1869  prepare_add();
1870 
1876  void
1877  prepare_set();
1878 
1879 
1880 
1881  private:
1886  std::unique_ptr<Epetra_Map> column_space_map;
1887 
1893  std::unique_ptr<Epetra_FECrsMatrix> matrix;
1894 
1900  std::unique_ptr<Epetra_CrsMatrix> nonlocal_matrix;
1901 
1905  std::unique_ptr<Epetra_Export> nonlocal_matrix_exporter;
1906 
1918  Epetra_CombineMode last_action;
1919 
1925 
1926  // To allow calling protected prepare_add() and prepare_set().
1928  };
1929 
1930 
1931 
1932  // forwards declarations
1933  class SolverBase;
1934  class PreconditionBase;
1935 
1936  namespace internal
1937  {
1938  inline void
1939  check_vector_map_equality(const Epetra_CrsMatrix & mtrx,
1940  const Epetra_MultiVector &src,
1941  const Epetra_MultiVector &dst,
1942  const bool transpose)
1943  {
1944  if (transpose == false)
1945  {
1946  Assert(src.Map().SameAs(mtrx.DomainMap()) == true,
1947  ExcMessage(
1948  "Column map of matrix does not fit with vector map!"));
1949  Assert(dst.Map().SameAs(mtrx.RangeMap()) == true,
1950  ExcMessage("Row map of matrix does not fit with vector map!"));
1951  }
1952  else
1953  {
1954  Assert(src.Map().SameAs(mtrx.RangeMap()) == true,
1955  ExcMessage(
1956  "Column map of matrix does not fit with vector map!"));
1957  Assert(dst.Map().SameAs(mtrx.DomainMap()) == true,
1958  ExcMessage("Row map of matrix does not fit with vector map!"));
1959  }
1960  (void)mtrx; // removes -Wunused-variable in optimized mode
1961  (void)src;
1962  (void)dst;
1963  }
1964 
1965  inline void
1967  const Epetra_MultiVector &src,
1968  const Epetra_MultiVector &dst,
1969  const bool transpose)
1970  {
1971  if (transpose == false)
1972  {
1973  Assert(src.Map().SameAs(op.OperatorDomainMap()) == true,
1974  ExcMessage(
1975  "Column map of operator does not fit with vector map!"));
1976  Assert(dst.Map().SameAs(op.OperatorRangeMap()) == true,
1977  ExcMessage(
1978  "Row map of operator does not fit with vector map!"));
1979  }
1980  else
1981  {
1982  Assert(src.Map().SameAs(op.OperatorRangeMap()) == true,
1983  ExcMessage(
1984  "Column map of operator does not fit with vector map!"));
1985  Assert(dst.Map().SameAs(op.OperatorDomainMap()) == true,
1986  ExcMessage(
1987  "Row map of operator does not fit with vector map!"));
1988  }
1989  (void)op; // removes -Wunused-variable in optimized mode
1990  (void)src;
1991  (void)dst;
1992  }
1993 
1994  namespace LinearOperatorImplementation
1995  {
2017  {
2018  public:
2022  using VectorType = Epetra_MultiVector;
2023 
2028 
2033 
2038 
2046  TrilinosPayload();
2047 
2051  TrilinosPayload(const TrilinosWrappers::SparseMatrix &matrix_exemplar,
2053 
2058  const TrilinosWrappers::SparseMatrix & matrix_exemplar,
2059  const TrilinosWrappers::PreconditionBase &preconditioner);
2060 
2065  const TrilinosWrappers::PreconditionBase &preconditioner_exemplar,
2066  const TrilinosWrappers::PreconditionBase &preconditioner);
2067 
2071  TrilinosPayload(const TrilinosPayload &payload);
2072 
2080  TrilinosPayload(const TrilinosPayload &first_op,
2081  const TrilinosPayload &second_op);
2082 
2086  virtual ~TrilinosPayload() override = default;
2087 
2092  identity_payload() const;
2093 
2098  null_payload() const;
2099 
2104  transpose_payload() const;
2105 
2122  template <typename Solver, typename Preconditioner>
2123  typename std::enable_if<
2125  std::is_base_of<TrilinosWrappers::PreconditionBase,
2126  Preconditioner>::value,
2127  TrilinosPayload>::type
2128  inverse_payload(Solver &, const Preconditioner &) const;
2129 
2147  template <typename Solver, typename Preconditioner>
2148  typename std::enable_if<
2150  std::is_base_of<TrilinosWrappers::PreconditionBase,
2151  Preconditioner>::value),
2152  TrilinosPayload>::type
2153  inverse_payload(Solver &, const Preconditioner &) const;
2154 
2156 
2161 
2167  IndexSet
2168  locally_owned_domain_indices() const;
2169 
2175  IndexSet
2176  locally_owned_range_indices() const;
2177 
2181  MPI_Comm
2182  get_mpi_communicator() const;
2183 
2190  void
2191  transpose();
2192 
2200  std::function<void(VectorType &, const VectorType &)> vmult;
2201 
2209  std::function<void(VectorType &, const VectorType &)> Tvmult;
2210 
2219  std::function<void(VectorType &, const VectorType &)> inv_vmult;
2220 
2229  std::function<void(VectorType &, const VectorType &)> inv_Tvmult;
2230 
2232 
2237 
2244  virtual bool
2245  UseTranspose() const override;
2246 
2262  virtual int
2263  SetUseTranspose(bool UseTranspose) override;
2264 
2276  virtual int
2277  Apply(const VectorType &X, VectorType &Y) const override;
2278 
2297  virtual int
2298  ApplyInverse(const VectorType &Y, VectorType &X) const override;
2300 
2305 
2312  virtual const char *
2313  Label() const override;
2314 
2322  virtual const Epetra_Comm &
2323  Comm() const override;
2324 
2332  virtual const Epetra_Map &
2333  OperatorDomainMap() const override;
2334 
2343  virtual const Epetra_Map &
2344  OperatorRangeMap() const override;
2346 
2347  private:
2353 
2358 # ifdef DEAL_II_WITH_MPI
2359  Epetra_MpiComm communicator;
2360 # else
2361  Epetra_SerialComm communicator;
2362 # endif
2363 
2368  Epetra_Map domain_map;
2369 
2374  Epetra_Map range_map;
2375 
2384  virtual bool
2385  HasNormInf() const override;
2386 
2394  virtual double
2395  NormInf() const override;
2396  };
2397 
2403  operator+(const TrilinosPayload &first_op,
2404  const TrilinosPayload &second_op);
2405 
2410  TrilinosPayload operator*(const TrilinosPayload &first_op,
2411  const TrilinosPayload &second_op);
2412 
2413  } // namespace LinearOperatorImplementation
2414  } /* namespace internal */
2415 
2416 
2417 
2418  // ----------------------- inline and template functions --------------------
2419 
2420 # ifndef DOXYGEN
2421 
2422  namespace SparseMatrixIterators
2423  {
2424  inline AccessorBase::AccessorBase(SparseMatrix *matrix,
2425  size_type row,
2426  size_type index)
2427  : matrix(matrix)
2428  , a_row(row)
2429  , a_index(index)
2430  {
2431  visit_present_row();
2432  }
2433 
2434 
2436  AccessorBase::row() const
2437  {
2438  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2439  return a_row;
2440  }
2441 
2442 
2444  AccessorBase::column() const
2445  {
2446  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2447  return (*colnum_cache)[a_index];
2448  }
2449 
2450 
2452  AccessorBase::index() const
2453  {
2454  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2455  return a_index;
2456  }
2457 
2458 
2459  inline Accessor<true>::Accessor(MatrixType * matrix,
2460  const size_type row,
2461  const size_type index)
2462  : AccessorBase(const_cast<SparseMatrix *>(matrix), row, index)
2463  {}
2464 
2465 
2466  template <bool Other>
2467  inline Accessor<true>::Accessor(const Accessor<Other> &other)
2468  : AccessorBase(other)
2469  {}
2470 
2471 
2472  inline TrilinosScalar
2473  Accessor<true>::value() const
2474  {
2475  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2476  return (*value_cache)[a_index];
2477  }
2478 
2479 
2480  inline Accessor<false>::Reference::Reference(const Accessor<false> &acc)
2481  : accessor(const_cast<Accessor<false> &>(acc))
2482  {}
2483 
2484 
2485  inline Accessor<false>::Reference::operator TrilinosScalar() const
2486  {
2487  return (*accessor.value_cache)[accessor.a_index];
2488  }
2489 
2490  inline const Accessor<false>::Reference &
2491  Accessor<false>::Reference::operator=(const TrilinosScalar n) const
2492  {
2493  (*accessor.value_cache)[accessor.a_index] = n;
2494  accessor.matrix->set(accessor.row(),
2495  accessor.column(),
2496  static_cast<TrilinosScalar>(*this));
2497  return *this;
2498  }
2499 
2500 
2501  inline const Accessor<false>::Reference &
2502  Accessor<false>::Reference::operator+=(const TrilinosScalar n) const
2503  {
2504  (*accessor.value_cache)[accessor.a_index] += n;
2505  accessor.matrix->set(accessor.row(),
2506  accessor.column(),
2507  static_cast<TrilinosScalar>(*this));
2508  return *this;
2509  }
2510 
2511 
2512  inline const Accessor<false>::Reference &
2513  Accessor<false>::Reference::operator-=(const TrilinosScalar n) const
2514  {
2515  (*accessor.value_cache)[accessor.a_index] -= n;
2516  accessor.matrix->set(accessor.row(),
2517  accessor.column(),
2518  static_cast<TrilinosScalar>(*this));
2519  return *this;
2520  }
2521 
2522 
2523  inline const Accessor<false>::Reference &
2524  Accessor<false>::Reference::operator*=(const TrilinosScalar n) const
2525  {
2526  (*accessor.value_cache)[accessor.a_index] *= n;
2527  accessor.matrix->set(accessor.row(),
2528  accessor.column(),
2529  static_cast<TrilinosScalar>(*this));
2530  return *this;
2531  }
2532 
2533 
2534  inline const Accessor<false>::Reference &
2535  Accessor<false>::Reference::operator/=(const TrilinosScalar n) const
2536  {
2537  (*accessor.value_cache)[accessor.a_index] /= n;
2538  accessor.matrix->set(accessor.row(),
2539  accessor.column(),
2540  static_cast<TrilinosScalar>(*this));
2541  return *this;
2542  }
2543 
2544 
2545  inline Accessor<false>::Accessor(MatrixType * matrix,
2546  const size_type row,
2547  const size_type index)
2548  : AccessorBase(matrix, row, index)
2549  {}
2550 
2551 
2552  inline Accessor<false>::Reference
2553  Accessor<false>::value() const
2554  {
2555  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2556  return {*this};
2557  }
2558 
2559 
2560 
2561  template <bool Constness>
2562  inline Iterator<Constness>::Iterator(MatrixType * matrix,
2563  const size_type row,
2564  const size_type index)
2565  : accessor(matrix, row, index)
2566  {}
2567 
2568 
2569  template <bool Constness>
2570  template <bool Other>
2571  inline Iterator<Constness>::Iterator(const Iterator<Other> &other)
2572  : accessor(other.accessor)
2573  {}
2574 
2575 
2576  template <bool Constness>
2577  inline Iterator<Constness> &
2579  {
2580  Assert(accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
2581 
2582  ++accessor.a_index;
2583 
2584  // If at end of line: do one
2585  // step, then cycle until we
2586  // find a row with a nonzero
2587  // number of entries.
2588  if (accessor.a_index >= accessor.colnum_cache->size())
2589  {
2590  accessor.a_index = 0;
2591  ++accessor.a_row;
2592 
2593  while ((accessor.a_row < accessor.matrix->m()) &&
2594  ((accessor.matrix->in_local_range(accessor.a_row) == false) ||
2595  (accessor.matrix->row_length(accessor.a_row) == 0)))
2596  ++accessor.a_row;
2597 
2598  accessor.visit_present_row();
2599  }
2600  return *this;
2601  }
2602 
2603 
2604  template <bool Constness>
2605  inline Iterator<Constness>
2607  {
2608  const Iterator<Constness> old_state = *this;
2609  ++(*this);
2610  return old_state;
2611  }
2612 
2613 
2614 
2615  template <bool Constness>
2616  inline const Accessor<Constness> &Iterator<Constness>::operator*() const
2617  {
2618  return accessor;
2619  }
2620 
2621 
2622 
2623  template <bool Constness>
2624  inline const Accessor<Constness> *Iterator<Constness>::operator->() const
2625  {
2626  return &accessor;
2627  }
2628 
2629 
2630 
2631  template <bool Constness>
2632  inline bool
2633  Iterator<Constness>::operator==(const Iterator<Constness> &other) const
2634  {
2635  return (accessor.a_row == other.accessor.a_row &&
2636  accessor.a_index == other.accessor.a_index);
2637  }
2638 
2639 
2640 
2641  template <bool Constness>
2642  inline bool
2643  Iterator<Constness>::operator!=(const Iterator<Constness> &other) const
2644  {
2645  return !(*this == other);
2646  }
2647 
2648 
2649 
2650  template <bool Constness>
2651  inline bool
2652  Iterator<Constness>::operator<(const Iterator<Constness> &other) const
2653  {
2654  return (accessor.row() < other.accessor.row() ||
2655  (accessor.row() == other.accessor.row() &&
2656  accessor.index() < other.accessor.index()));
2657  }
2658 
2659 
2660  template <bool Constness>
2661  inline bool
2662  Iterator<Constness>::operator>(const Iterator<Constness> &other) const
2663  {
2664  return (other < *this);
2665  }
2666 
2667  } // namespace SparseMatrixIterators
2668 
2669 
2670 
2672  SparseMatrix::begin() const
2673  {
2674  return begin(0);
2675  }
2676 
2677 
2678 
2680  SparseMatrix::end() const
2681  {
2682  return const_iterator(this, m(), 0);
2683  }
2684 
2685 
2686 
2688  SparseMatrix::begin(const size_type r) const
2689  {
2690  AssertIndexRange(r, m());
2691  if (in_local_range(r) && (row_length(r) > 0))
2692  return const_iterator(this, r, 0);
2693  else
2694  return end(r);
2695  }
2696 
2697 
2698 
2700  SparseMatrix::end(const size_type r) const
2701  {
2702  AssertIndexRange(r, m());
2703 
2704  // place the iterator on the first entry
2705  // past this line, or at the end of the
2706  // matrix
2707  for (size_type i = r + 1; i < m(); ++i)
2708  if (in_local_range(i) && (row_length(i) > 0))
2709  return const_iterator(this, i, 0);
2710 
2711  // if there is no such line, then take the
2712  // end iterator of the matrix
2713  return end();
2714  }
2715 
2716 
2717 
2718  inline SparseMatrix::iterator
2720  {
2721  return begin(0);
2722  }
2723 
2724 
2725 
2726  inline SparseMatrix::iterator
2728  {
2729  return iterator(this, m(), 0);
2730  }
2731 
2732 
2733 
2734  inline SparseMatrix::iterator
2736  {
2737  AssertIndexRange(r, m());
2738  if (in_local_range(r) && (row_length(r) > 0))
2739  return iterator(this, r, 0);
2740  else
2741  return end(r);
2742  }
2743 
2744 
2745 
2746  inline SparseMatrix::iterator
2747  SparseMatrix::end(const size_type r)
2748  {
2749  AssertIndexRange(r, m());
2750 
2751  // place the iterator on the first entry
2752  // past this line, or at the end of the
2753  // matrix
2754  for (size_type i = r + 1; i < m(); ++i)
2755  if (in_local_range(i) && (row_length(i) > 0))
2756  return iterator(this, i, 0);
2757 
2758  // if there is no such line, then take the
2759  // end iterator of the matrix
2760  return end();
2761  }
2762 
2763 
2764 
2765  inline bool
2766  SparseMatrix::in_local_range(const size_type index) const
2767  {
2769 # ifndef DEAL_II_WITH_64BIT_INDICES
2770  begin = matrix->RowMap().MinMyGID();
2771  end = matrix->RowMap().MaxMyGID() + 1;
2772 # else
2773  begin = matrix->RowMap().MinMyGID64();
2774  end = matrix->RowMap().MaxMyGID64() + 1;
2775 # endif
2776 
2777  return ((index >= static_cast<size_type>(begin)) &&
2778  (index < static_cast<size_type>(end)));
2779  }
2780 
2781 
2782 
2783  inline bool
2785  {
2786  return compressed;
2787  }
2788 
2789 
2790 
2791  // Inline the set() and add() functions, since they will be called
2792  // frequently, and the compiler can optimize away some unnecessary loops
2793  // when the sizes are given at compile time.
2794  template <>
2795  void
2796  SparseMatrix::set<TrilinosScalar>(const size_type row,
2797  const size_type n_cols,
2798  const size_type * col_indices,
2799  const TrilinosScalar *values,
2800  const bool elide_zero_values);
2801 
2802 
2803 
2804  template <typename Number>
2805  void
2806  SparseMatrix::set(const size_type row,
2807  const size_type n_cols,
2808  const size_type *col_indices,
2809  const Number * values,
2810  const bool elide_zero_values)
2811  {
2812  std::vector<TrilinosScalar> trilinos_values(n_cols);
2813  std::copy(values, values + n_cols, trilinos_values.begin());
2814  this->set(
2815  row, n_cols, col_indices, trilinos_values.data(), elide_zero_values);
2816  }
2817 
2818 
2819 
2820  inline void
2821  SparseMatrix::set(const size_type i,
2822  const size_type j,
2823  const TrilinosScalar value)
2824  {
2825  AssertIsFinite(value);
2826 
2827  set(i, 1, &j, &value, false);
2828  }
2829 
2830 
2831 
2832  inline void
2833  SparseMatrix::set(const std::vector<size_type> & indices,
2834  const FullMatrix<TrilinosScalar> &values,
2835  const bool elide_zero_values)
2836  {
2837  Assert(indices.size() == values.m(),
2838  ExcDimensionMismatch(indices.size(), values.m()));
2839  Assert(values.m() == values.n(), ExcNotQuadratic());
2840 
2841  for (size_type i = 0; i < indices.size(); ++i)
2842  set(indices[i],
2843  indices.size(),
2844  indices.data(),
2845  &values(i, 0),
2846  elide_zero_values);
2847  }
2848 
2849 
2850 
2851  inline void
2852  SparseMatrix::add(const size_type i,
2853  const size_type j,
2854  const TrilinosScalar value)
2855  {
2856  AssertIsFinite(value);
2857 
2858  if (value == 0)
2859  {
2860  // we have to check after Insert/Add in any case to be consistent
2861  // with the MPI communication model, but we can save some
2862  // work if the addend is zero. However, these actions are done in case
2863  // we pass on to the other function.
2864 
2865  // TODO: fix this (do not run compress here, but fail)
2866  if (last_action == Insert)
2867  {
2868  int ierr;
2869  ierr = matrix->GlobalAssemble(*column_space_map,
2870  matrix->RowMap(),
2871  false);
2872 
2873  Assert(ierr == 0, ExcTrilinosError(ierr));
2874  (void)ierr; // removes -Wunused-but-set-variable in optimized mode
2875  }
2876 
2877  last_action = Add;
2878 
2879  return;
2880  }
2881  else
2882  add(i, 1, &j, &value, false);
2883  }
2884 
2885 
2886 
2887  // inline "simple" functions that are called frequently and do only involve
2888  // a call to some Trilinos function.
2890  SparseMatrix::m() const
2891  {
2892 # ifndef DEAL_II_WITH_64BIT_INDICES
2893  return matrix->NumGlobalRows();
2894 # else
2895  return matrix->NumGlobalRows64();
2896 # endif
2897  }
2898 
2899 
2900 
2902  SparseMatrix::n() const
2903  {
2904  // If the matrix structure has not been fixed (i.e., we did not have a
2905  // sparsity pattern), it does not know about the number of columns so we
2906  // must always take this from the additional column space map
2907  Assert(column_space_map.get() != nullptr, ExcInternalError());
2908  return n_global_elements(*column_space_map);
2909  }
2910 
2911 
2912 
2913  inline unsigned int
2914  SparseMatrix::local_size() const
2915  {
2916  return matrix->NumMyRows();
2917  }
2918 
2919 
2920 
2921  inline std::pair<SparseMatrix::size_type, SparseMatrix::size_type>
2923  {
2924  size_type begin, end;
2925 # ifndef DEAL_II_WITH_64BIT_INDICES
2926  begin = matrix->RowMap().MinMyGID();
2927  end = matrix->RowMap().MaxMyGID() + 1;
2928 # else
2929  begin = matrix->RowMap().MinMyGID64();
2930  end = matrix->RowMap().MaxMyGID64() + 1;
2931 # endif
2932 
2933  return std::make_pair(begin, end);
2934  }
2935 
2936 
2937 
2940  {
2941 # ifndef DEAL_II_WITH_64BIT_INDICES
2942  return matrix->NumGlobalNonzeros();
2943 # else
2944  return matrix->NumGlobalNonzeros64();
2945 # endif
2946  }
2947 
2948 
2949 
2950  template <typename SparsityPatternType>
2951  inline void
2952  SparseMatrix::reinit(const IndexSet & parallel_partitioning,
2953  const SparsityPatternType &sparsity_pattern,
2954  const MPI_Comm & communicator,
2955  const bool exchange_data)
2956  {
2957  reinit(parallel_partitioning,
2958  parallel_partitioning,
2959  sparsity_pattern,
2960  communicator,
2961  exchange_data);
2962  }
2963 
2964 
2965 
2966  template <typename number>
2967  inline void
2968  SparseMatrix::reinit(const IndexSet &parallel_partitioning,
2969  const ::SparseMatrix<number> &sparse_matrix,
2970  const MPI_Comm & communicator,
2971  const double drop_tolerance,
2972  const bool copy_values,
2973  const ::SparsityPattern * use_this_sparsity)
2974  {
2975  Epetra_Map map =
2976  parallel_partitioning.make_trilinos_map(communicator, false);
2977  reinit(parallel_partitioning,
2978  parallel_partitioning,
2979  sparse_matrix,
2980  drop_tolerance,
2981  copy_values,
2982  use_this_sparsity);
2983  }
2984 
2985 
2986 
2987  inline const Epetra_CrsMatrix &
2989  {
2990  return static_cast<const Epetra_CrsMatrix &>(*matrix);
2991  }
2992 
2993 
2994 
2995  inline const Epetra_CrsGraph &
2997  {
2998  return matrix->Graph();
2999  }
3000 
3001 
3002 
3003  inline IndexSet
3005  {
3006  return IndexSet(matrix->DomainMap());
3007  }
3008 
3009 
3010 
3011  inline IndexSet
3013  {
3014  return IndexSet(matrix->RangeMap());
3015  }
3016 
3017 
3018 
3019  inline void
3021  {
3022  // nothing to do here
3023  }
3024 
3025 
3026 
3027  inline void
3029  {
3030  // nothing to do here
3031  }
3032 
3033 
3034  namespace internal
3035  {
3036  namespace LinearOperatorImplementation
3037  {
3038  template <typename Solver, typename Preconditioner>
3039  typename std::enable_if<
3041  std::is_base_of<TrilinosWrappers::PreconditionBase,
3042  Preconditioner>::value,
3043  TrilinosPayload>::type
3045  Solver & solver,
3046  const Preconditioner &preconditioner) const
3047  {
3048  const auto &payload = *this;
3049 
3050  TrilinosPayload return_op(payload);
3051 
3052  // Capture by copy so the payloads are always valid
3053 
3054  return_op.inv_vmult = [payload, &solver, &preconditioner](
3055  TrilinosPayload::Domain & tril_dst,
3056  const TrilinosPayload::Range &tril_src) {
3057  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3058  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3059  Assert(&tril_src != &tril_dst,
3062  tril_src,
3063  tril_dst,
3064  !payload.UseTranspose());
3065  solver.solve(payload, tril_dst, tril_src, preconditioner);
3066  };
3067 
3068  return_op.inv_Tvmult = [payload, &solver, &preconditioner](
3069  TrilinosPayload::Range & tril_dst,
3070  const TrilinosPayload::Domain &tril_src) {
3071  // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3072  // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3073  Assert(&tril_src != &tril_dst,
3076  tril_src,
3077  tril_dst,
3078  payload.UseTranspose());
3079 
3080  const_cast<TrilinosPayload &>(payload).transpose();
3081  solver.solve(payload, tril_dst, tril_src, preconditioner);
3082  const_cast<TrilinosPayload &>(payload).transpose();
3083  };
3084 
3085  // If the input operator is already setup for transpose operations, then
3086  // we must do similar with its inverse.
3087  if (return_op.UseTranspose() == true)
3088  std::swap(return_op.inv_vmult, return_op.inv_Tvmult);
3089 
3090  return return_op;
3091  }
3092 
3093  template <typename Solver, typename Preconditioner>
3094  typename std::enable_if<
3096  std::is_base_of<TrilinosWrappers::PreconditionBase,
3097  Preconditioner>::value),
3098  TrilinosPayload>::type
3099  TrilinosPayload::inverse_payload(Solver &, const Preconditioner &) const
3100  {
3101  TrilinosPayload return_op(*this);
3102 
3103  return_op.inv_vmult = [](TrilinosPayload::Domain &,
3104  const TrilinosPayload::Range &) {
3105  AssertThrow(false,
3106  ExcMessage("Payload inv_vmult disabled because of "
3107  "incompatible solver/preconditioner choice."));
3108  };
3109 
3110  return_op.inv_Tvmult = [](TrilinosPayload::Range &,
3111  const TrilinosPayload::Domain &) {
3112  AssertThrow(false,
3113  ExcMessage("Payload inv_vmult disabled because of "
3114  "incompatible solver/preconditioner choice."));
3115  };
3116 
3117  return return_op;
3118  }
3119  } // namespace LinearOperatorImplementation
3120  } // namespace internal
3121 
3122  template <>
3123  void
3124  SparseMatrix::set<TrilinosScalar>(const size_type row,
3125  const size_type n_cols,
3126  const size_type * col_indices,
3127  const TrilinosScalar *values,
3128  const bool elide_zero_values);
3129 # endif // DOXYGEN
3130 
3131 } /* namespace TrilinosWrappers */
3132 
3133 
3135 
3136 
3137 # endif // DEAL_II_WITH_TRILINOS
3138 
3139 
3140 /*----------------------- trilinos_sparse_matrix.h --------------------*/
3141 
3142 #endif
3143 /*----------------------- 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:621
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:541
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:1649
static const char V
TrilinosPayload operator+(const TrilinosPayload &first_op, const TrilinosPayload &second_op)
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)
static ::ExceptionBase & ExcAccessToNonlocalRow(std::size_t arg1, std::size_t arg2, std::size_t arg3)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:596
std::shared_ptr< std::vector< size_type > > colnum_cache
std::unique_ptr< Epetra_FECrsMatrix > matrix
::types::global_dof_index size_type
bool in_local_range(const size_type index) const
static ::ExceptionBase & ExcMessage(std::string arg1)
std::string compress(const std::string &input)
Definition: utilities.cc:392
size_type n() const
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:518
void add(const size_type i, const size_type j, const TrilinosScalar value)
#define Assert(cond, exc)
Definition: exceptions.h:1419
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:496
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2732
#define DeclException0(Exception0)
Definition: exceptions.h:473
typename Accessor< Constness >::MatrixType MatrixType
MatrixTableIterators::AccessorBase< TransposeTable< T >, Constness, MatrixTableIterators::Storage::column_major > AccessorBase
Definition: table.h:1901
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
std::function< void(VectorType &, const VectorType &)> inv_vmult
VectorType::value_type * end(VectorType &V)
double TrilinosScalar
Definition: types.h:158
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
MatrixTableIterators::Accessor< TransposeTable< T >, Constness, MatrixTableIterators::Storage::column_major > Accessor
Definition: table.h:1907
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
MatrixTableIterators::Iterator< TransposeTable< T >, Constness, MatrixTableIterators::Storage::column_major > Iterator
Definition: table.h:1913
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
Definition: memory_space.h:103
void reinit(const SparsityPatternType &sparsity_pattern)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
VectorType::value_type * begin(VectorType &V)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:587
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2705
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:564
std::unique_ptr< Epetra_Map > column_space_map
std::unique_ptr< Epetra_CrsMatrix > nonlocal_matrix
static const bool value
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:1681
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