Reference documentation for deal.II version Git 1dc1051882 2021-04-22 23:57:03 +0200
\(\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_precondition.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_precondition_h
17 # define dealii_trilinos_precondition_h
18 
19 
20 # include <deal.II/base/config.h>
21 
22 # ifdef DEAL_II_WITH_TRILINOS
23 
25 
28 
29 # include <memory>
30 
31 # ifdef DEAL_II_WITH_MPI
32 # include <Epetra_MpiComm.h>
33 # else
34 # include <Epetra_SerialComm.h>
35 # endif
36 # include <Epetra_Map.h>
37 # include <Epetra_MultiVector.h>
38 # include <Epetra_RowMatrix.h>
39 # include <Epetra_Vector.h>
40 # include <Teuchos_ParameterList.hpp>
41 
42 // forward declarations
43 # ifndef DOXYGEN
44 class Ifpack_Preconditioner;
45 class Ifpack_Chebyshev;
46 namespace ML_Epetra
47 {
48  class MultiLevelPreconditioner;
49 }
50 # endif
51 
53 
54 // forward declarations
55 # ifndef DOXYGEN
56 template <typename number>
57 class SparseMatrix;
58 template <typename number>
59 class Vector;
60 class SparsityPattern;
61 # endif
62 
67 namespace TrilinosWrappers
68 {
69  // forward declarations
70  class SparseMatrix;
71  class BlockSparseMatrix;
72  class SolverBase;
73 
81  {
82  public:
87 
93  {};
94 
101 
106 
110  ~PreconditionBase() override = default;
111 
116  void
117  clear();
118 
122  MPI_Comm
123  get_mpi_communicator() const;
124 
134  void
135  transpose();
136 
140  virtual void
141  vmult(MPI::Vector &dst, const MPI::Vector &src) const;
142 
146  virtual void
147  Tvmult(MPI::Vector &dst, const MPI::Vector &src) const;
148 
153  virtual void
154  vmult(::Vector<double> &dst, const ::Vector<double> &src) const;
155 
160  virtual void
161  Tvmult(::Vector<double> & dst,
162  const ::Vector<double> &src) const;
163 
168  virtual void
170  const ::LinearAlgebra::distributed::Vector<double> &src) const;
171 
176  virtual void
178  const ::LinearAlgebra::distributed::Vector<double> &src) const;
179 
190  trilinos_operator() const;
192 
197 
202  IndexSet
203  locally_owned_domain_indices() const;
204 
210  IndexSet
211  locally_owned_range_indices() const;
212 
214 
222  DeclException1(ExcNonMatchingMaps,
223  std::string,
224  << "The sparse matrix the preconditioner is based on "
225  << "uses a map that is not compatible to the one in vector "
226  << arg1 << ". Check preconditioner and matrix setup.");
228 
229  friend class SolverBase;
230 
231  protected:
236  Teuchos::RCP<Epetra_Operator> preconditioner;
237 
242 # ifdef DEAL_II_WITH_MPI
243  Epetra_MpiComm communicator;
244 # else
245  Epetra_SerialComm communicator;
246 # endif
247 
252  std::shared_ptr<Epetra_Map> vector_distributor;
253  };
254 
255 
272  {
273  public:
286  {
291  AdditionalData(const double omega = 1,
292  const double min_diagonal = 0,
293  const unsigned int n_sweeps = 1);
294 
298  double omega;
299 
307  double min_diagonal;
308 
313  unsigned int n_sweeps;
314  };
315 
320  void
321  initialize(const SparseMatrix & matrix,
322  const AdditionalData &additional_data = AdditionalData());
323  };
324 
325 
326 
353  {
354  public:
369  {
376  AdditionalData(const double omega = 1,
377  const double min_diagonal = 0,
378  const unsigned int overlap = 0,
379  const unsigned int n_sweeps = 1);
380 
385  double omega;
386 
394  double min_diagonal;
395 
400  unsigned int overlap;
401 
406  unsigned int n_sweeps;
407  };
408 
414  void
415  initialize(const SparseMatrix & matrix,
416  const AdditionalData &additional_data = AdditionalData());
417  };
418 
419 
420 
447  {
448  public:
463  {
470  AdditionalData(const double omega = 1,
471  const double min_diagonal = 0,
472  const unsigned int overlap = 0,
473  const unsigned int n_sweeps = 1);
474 
479  double omega;
480 
488  double min_diagonal;
489 
494  unsigned int overlap;
495 
500  unsigned int n_sweeps;
501  };
502 
508  void
509  initialize(const SparseMatrix & matrix,
510  const AdditionalData &additional_data = AdditionalData());
511  };
512 
513 
514 
532  {
533  public:
550  {
556  AdditionalData(const unsigned int block_size = 1,
557  const std::string &block_creation_type = "linear",
558  const double omega = 1,
559  const double min_diagonal = 0,
560  const unsigned int n_sweeps = 1);
561 
565  unsigned int block_size;
566 
574  std::string block_creation_type;
575 
580  double omega;
581 
589  double min_diagonal;
590 
595  unsigned int n_sweeps;
596  };
597 
602  void
603  initialize(const SparseMatrix & matrix,
604  const AdditionalData &additional_data = AdditionalData());
605  };
606 
607 
608 
631  {
632  public:
651  {
659  AdditionalData(const unsigned int block_size = 1,
660  const std::string &block_creation_type = "linear",
661  const double omega = 1,
662  const double min_diagonal = 0,
663  const unsigned int overlap = 0,
664  const unsigned int n_sweeps = 1);
665 
669  unsigned int block_size;
670 
678  std::string block_creation_type;
679 
684  double omega;
685 
693  double min_diagonal;
694 
699  unsigned int overlap;
700 
705  unsigned int n_sweeps;
706  };
707 
713  void
714  initialize(const SparseMatrix & matrix,
715  const AdditionalData &additional_data = AdditionalData());
716  };
717 
718 
719 
742  {
743  public:
762  {
770  AdditionalData(const unsigned int block_size = 1,
771  const std::string &block_creation_type = "linear",
772  const double omega = 1,
773  const double min_diagonal = 0,
774  const unsigned int overlap = 0,
775  const unsigned int n_sweeps = 1);
776 
780  unsigned int block_size;
781 
789  std::string block_creation_type;
790 
795  double omega;
796 
804  double min_diagonal;
805 
810  unsigned int overlap;
811 
816  unsigned int n_sweeps;
817  };
818 
824  void
825  initialize(const SparseMatrix & matrix,
826  const AdditionalData &additional_data = AdditionalData());
827  };
828 
829 
830 
866  {
867  public:
885  {
895  AdditionalData(const unsigned int ic_fill = 0,
896  const double ic_atol = 0.,
897  const double ic_rtol = 1.,
898  const unsigned int overlap = 0);
899 
908  unsigned int ic_fill;
909 
915  double ic_atol;
916 
921  double ic_rtol;
922 
927  unsigned int overlap;
928  };
929 
934  void
935  initialize(const SparseMatrix & matrix,
936  const AdditionalData &additional_data = AdditionalData());
937  };
938 
939 
940 
970  {
971  public:
1010  {
1014  AdditionalData(const unsigned int ilu_fill = 0,
1015  const double ilu_atol = 0.,
1016  const double ilu_rtol = 1.,
1017  const unsigned int overlap = 0);
1018 
1022  unsigned int ilu_fill;
1023 
1028  double ilu_atol;
1029 
1034  double ilu_rtol;
1035 
1039  unsigned int overlap;
1040  };
1041 
1046  void
1047  initialize(const SparseMatrix & matrix,
1048  const AdditionalData &additional_data = AdditionalData());
1049  };
1050 
1051 
1052 
1088  {
1089  public:
1108  {
1119  AdditionalData(const double ilut_drop = 0.,
1120  const unsigned int ilut_fill = 0,
1121  const double ilut_atol = 0.,
1122  const double ilut_rtol = 1.,
1123  const unsigned int overlap = 0);
1124 
1129  double ilut_drop;
1130 
1139  unsigned int ilut_fill;
1140 
1146  double ilut_atol;
1147 
1152  double ilut_rtol;
1153 
1158  unsigned int overlap;
1159  };
1160 
1165  void
1166  initialize(const SparseMatrix & matrix,
1167  const AdditionalData &additional_data = AdditionalData());
1168  };
1169 
1170 
1171 
1190  {
1191  public:
1197  {
1201  AdditionalData(const unsigned int overlap = 0);
1202 
1203 
1208  unsigned int overlap;
1209  };
1210 
1215  void
1216  initialize(const SparseMatrix & matrix,
1217  const AdditionalData &additional_data = AdditionalData());
1218  };
1219 
1220 
1221 
1231  {
1232  public:
1238  {
1242  AdditionalData(const unsigned int degree = 1,
1243  const double max_eigenvalue = 10.,
1244  const double eigenvalue_ratio = 30.,
1245  const double min_eigenvalue = 1.,
1246  const double min_diagonal = 1e-12,
1247  const bool nonzero_starting = false);
1248 
1254  unsigned int degree;
1255 
1261 
1266 
1272 
1278 
1290  };
1291 
1296  void
1297  initialize(const SparseMatrix & matrix,
1298  const AdditionalData &additional_data = AdditionalData());
1299  };
1300 
1301 
1302 
1345  {
1346  public:
1354  {
1378  AdditionalData(const bool elliptic = true,
1379  const bool higher_order_elements = false,
1380  const unsigned int n_cycles = 1,
1381  const bool w_cyle = false,
1382  const double aggregation_threshold = 1e-4,
1383  const std::vector<std::vector<bool>> &constant_modes =
1384  std::vector<std::vector<bool>>(0),
1385  const unsigned int smoother_sweeps = 2,
1386  const unsigned int smoother_overlap = 0,
1387  const bool output_details = false,
1388  const char * smoother_type = "Chebyshev",
1389  const char * coarse_type = "Amesos-KLU");
1390 
1421  void
1422  set_parameters(
1423  Teuchos::ParameterList & parameter_list,
1424  std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1425  const Epetra_RowMatrix & matrix) const;
1426 
1434  void
1435  set_parameters(
1436  Teuchos::ParameterList & parameter_list,
1437  std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1438  const SparseMatrix & matrix) const;
1439 
1445  void
1446  set_operator_null_space(
1447  Teuchos::ParameterList & parameter_list,
1448  std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1449  const Epetra_RowMatrix & matrix) const;
1450 
1456  void
1457  set_operator_null_space(
1458  Teuchos::ParameterList & parameter_list,
1459  std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1460  const SparseMatrix & matrix) const;
1461 
1469  bool elliptic;
1470 
1476 
1481  unsigned int n_cycles;
1482 
1487  bool w_cycle;
1488 
1499 
1516  std::vector<std::vector<bool>> constant_modes;
1517 
1528  unsigned int smoother_sweeps;
1529 
1534  unsigned int smoother_overlap;
1535 
1542 
1577  const char *smoother_type;
1578 
1583  const char *coarse_type;
1584  };
1585 
1589  ~PreconditionAMG() override;
1590 
1591 
1597  void
1598  initialize(const SparseMatrix & matrix,
1599  const AdditionalData &additional_data = AdditionalData());
1600 
1619  void
1620  initialize(const Epetra_RowMatrix &matrix,
1621  const AdditionalData & additional_data = AdditionalData());
1622 
1637  void
1638  initialize(const SparseMatrix & matrix,
1639  const Teuchos::ParameterList &ml_parameters);
1640 
1648  void
1649  initialize(const Epetra_RowMatrix & matrix,
1650  const Teuchos::ParameterList &ml_parameters);
1651 
1658  template <typename number>
1659  void
1660  initialize(const ::SparseMatrix<number> &deal_ii_sparse_matrix,
1661  const AdditionalData &additional_data = AdditionalData(),
1662  const double drop_tolerance = 1e-13,
1663  const ::SparsityPattern *use_this_sparsity = nullptr);
1664 
1677  void
1678  reinit();
1679 
1684  void
1685  clear();
1686 
1690  size_type
1691  memory_consumption() const;
1692 
1693  private:
1697  std::shared_ptr<SparseMatrix> trilinos_matrix;
1698  };
1699 
1700 
1701 
1702 # if defined(DOXYGEN) || defined(DEAL_II_TRILINOS_WITH_MUELU)
1703 
1722  {
1723  public:
1731  {
1736  AdditionalData(const bool elliptic = true,
1737  const unsigned int n_cycles = 1,
1738  const bool w_cyle = false,
1739  const double aggregation_threshold = 1e-4,
1740  const std::vector<std::vector<bool>> &constant_modes =
1741  std::vector<std::vector<bool>>(0),
1742  const unsigned int smoother_sweeps = 2,
1743  const unsigned int smoother_overlap = 0,
1744  const bool output_details = false,
1745  const char * smoother_type = "Chebyshev",
1746  const char * coarse_type = "Amesos-KLU");
1747 
1755  bool elliptic;
1756 
1761  unsigned int n_cycles;
1762 
1767  bool w_cycle;
1768 
1779 
1786  std::vector<std::vector<bool>> constant_modes;
1787 
1798  unsigned int smoother_sweeps;
1799 
1804  unsigned int smoother_overlap;
1805 
1812 
1847  const char *smoother_type;
1848 
1853  const char *coarse_type;
1854  };
1855 
1860 
1864  virtual ~PreconditionAMGMueLu() override = default;
1865 
1871  void
1872  initialize(const SparseMatrix & matrix,
1873  const AdditionalData &additional_data = AdditionalData());
1874 
1881  void
1882  initialize(const Epetra_CrsMatrix &matrix,
1883  const AdditionalData & additional_data = AdditionalData());
1884 
1898  void
1899  initialize(const SparseMatrix & matrix,
1900  Teuchos::ParameterList &muelu_parameters);
1901 
1907  void
1908  initialize(const Epetra_CrsMatrix &matrix,
1909  Teuchos::ParameterList &muelu_parameters);
1910 
1917  template <typename number>
1918  void
1919  initialize(const ::SparseMatrix<number> &deal_ii_sparse_matrix,
1920  const AdditionalData &additional_data = AdditionalData(),
1921  const double drop_tolerance = 1e-13,
1922  const ::SparsityPattern *use_this_sparsity = nullptr);
1923 
1928  void
1929  clear();
1930 
1934  size_type
1935  memory_consumption() const;
1936 
1937  private:
1941  std::shared_ptr<SparseMatrix> trilinos_matrix;
1942  };
1943 # endif
1944 
1945 
1946 
1954  {
1955  public:
1961  {};
1962 
1969  void
1970  initialize(const SparseMatrix & matrix,
1971  const AdditionalData &additional_data = AdditionalData());
1972 
1976  void
1977  vmult(MPI::Vector &dst, const MPI::Vector &src) const override;
1978 
1982  void
1983  Tvmult(MPI::Vector &dst, const MPI::Vector &src) const override;
1984 
1989  void
1990  vmult(::Vector<double> & dst,
1991  const ::Vector<double> &src) const override;
1992 
1997  void
1998  Tvmult(::Vector<double> & dst,
1999  const ::Vector<double> &src) const override;
2000 
2005  void
2007  const ::LinearAlgebra::distributed::Vector<double> &src)
2008  const override;
2009 
2015  void
2017  const ::LinearAlgebra::distributed::Vector<double> &src)
2018  const override;
2019  };
2020 
2021 
2022 
2023  // ----------------------- inline and template functions --------------------
2024 
2025 
2026 # ifndef DOXYGEN
2027 
2028 
2029  inline void
2031  {
2032  // This only flips a flag that tells
2033  // Trilinos that any vmult operation
2034  // should be done with the
2035  // transpose. However, the matrix
2036  // structure is not reset.
2037  int ierr;
2038 
2039  if (!preconditioner->UseTranspose())
2040  {
2041  ierr = preconditioner->SetUseTranspose(true);
2042  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2043  }
2044  else
2045  {
2046  ierr = preconditioner->SetUseTranspose(false);
2047  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2048  }
2049  }
2050 
2051 
2052  inline void
2053  PreconditionBase::vmult(MPI::Vector &dst, const MPI::Vector &src) const
2054  {
2055  Assert(dst.trilinos_partitioner().SameAs(
2056  preconditioner->OperatorRangeMap()),
2057  ExcNonMatchingMaps("dst"));
2058  Assert(src.trilinos_partitioner().SameAs(
2059  preconditioner->OperatorDomainMap()),
2060  ExcNonMatchingMaps("src"));
2061 
2062  const int ierr = preconditioner->ApplyInverse(src.trilinos_vector(),
2063  dst.trilinos_vector());
2064  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2065  }
2066 
2067  inline void
2068  PreconditionBase::Tvmult(MPI::Vector &dst, const MPI::Vector &src) const
2069  {
2070  Assert(dst.trilinos_partitioner().SameAs(
2071  preconditioner->OperatorRangeMap()),
2072  ExcNonMatchingMaps("dst"));
2073  Assert(src.trilinos_partitioner().SameAs(
2074  preconditioner->OperatorDomainMap()),
2075  ExcNonMatchingMaps("src"));
2076 
2077  preconditioner->SetUseTranspose(true);
2078  const int ierr = preconditioner->ApplyInverse(src.trilinos_vector(),
2079  dst.trilinos_vector());
2080  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2081  preconditioner->SetUseTranspose(false);
2082  }
2083 
2084 
2085  // For the implementation of the <code>vmult</code> function with deal.II
2086  // data structures we note that invoking a call of the Trilinos
2087  // preconditioner requires us to use Epetra vectors as well. We do this by
2088  // providing a view, i.e., feed Trilinos with a pointer to the data, so we
2089  // avoid copying the content of the vectors during the iteration (this
2090  // function is only useful when used in serial anyway). In the declaration
2091  // of the right hand side, we need to cast the source vector (that is
2092  // <code>const</code> in all deal.II calls) to non-constant value, as this
2093  // is the way Trilinos wants to have them.
2094  inline void
2095  PreconditionBase::vmult(::Vector<double> & dst,
2096  const ::Vector<double> &src) const
2097  {
2098  AssertDimension(dst.size(),
2099  preconditioner->OperatorDomainMap().NumMyElements());
2100  AssertDimension(src.size(),
2101  preconditioner->OperatorRangeMap().NumMyElements());
2102  Epetra_Vector tril_dst(View,
2103  preconditioner->OperatorDomainMap(),
2104  dst.begin());
2105  Epetra_Vector tril_src(View,
2106  preconditioner->OperatorRangeMap(),
2107  const_cast<double *>(src.begin()));
2108 
2109  const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2110  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2111  }
2112 
2113 
2114  inline void
2115  PreconditionBase::Tvmult(::Vector<double> & dst,
2116  const ::Vector<double> &src) const
2117  {
2118  AssertDimension(dst.size(),
2119  preconditioner->OperatorDomainMap().NumMyElements());
2120  AssertDimension(src.size(),
2121  preconditioner->OperatorRangeMap().NumMyElements());
2122  Epetra_Vector tril_dst(View,
2123  preconditioner->OperatorDomainMap(),
2124  dst.begin());
2125  Epetra_Vector tril_src(View,
2126  preconditioner->OperatorRangeMap(),
2127  const_cast<double *>(src.begin()));
2128 
2129  preconditioner->SetUseTranspose(true);
2130  const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2131  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2132  preconditioner->SetUseTranspose(false);
2133  }
2134 
2135 
2136 
2137  inline void
2138  PreconditionBase::vmult(
2141  {
2143  preconditioner->OperatorDomainMap().NumMyElements());
2145  preconditioner->OperatorRangeMap().NumMyElements());
2146  Epetra_Vector tril_dst(View,
2147  preconditioner->OperatorDomainMap(),
2148  dst.begin());
2149  Epetra_Vector tril_src(View,
2150  preconditioner->OperatorRangeMap(),
2151  const_cast<double *>(src.begin()));
2152 
2153  const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2154  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2155  }
2156 
2157  inline void
2158  PreconditionBase::Tvmult(
2161  {
2163  preconditioner->OperatorDomainMap().NumMyElements());
2165  preconditioner->OperatorRangeMap().NumMyElements());
2166  Epetra_Vector tril_dst(View,
2167  preconditioner->OperatorDomainMap(),
2168  dst.begin());
2169  Epetra_Vector tril_src(View,
2170  preconditioner->OperatorRangeMap(),
2171  const_cast<double *>(src.begin()));
2172 
2173  preconditioner->SetUseTranspose(true);
2174  const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2175  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2176  preconditioner->SetUseTranspose(false);
2177  }
2178 
2179 # endif
2180 
2181 } // namespace TrilinosWrappers
2182 
2183 
2188 
2189 # endif // DEAL_II_WITH_TRILINOS
2190 
2191 #endif // trilinos_precondition_h
2192 /*------------------------- trilinos_precondition.h -------------------------*/
static ::ExceptionBase & ExcTrilinosError(int arg1)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1630
typename ::internal::FEValuesViews::ViewType< dim, spacedim, Extractor >::type View
Definition: fe_values.h:2062
Contents is actually a matrix.
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
constexpr int block_size
Definition: cuda_size.h:29
std::shared_ptr< Epetra_Map > vector_distributor
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
size_type locally_owned_size() const
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:515
#define Assert(cond, exc)
Definition: exceptions.h:1473
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:396
std::shared_ptr< SparseMatrix > trilinos_matrix
std::shared_ptr< SparseMatrix > trilinos_matrix
unsigned int global_dof_index
Definition: types.h:76
const Epetra_BlockMap & trilinos_partitioner() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:395
const Epetra_MultiVector & trilinos_vector() const
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
Teuchos::RCP< Epetra_Operator > preconditioner