Reference documentation for deal.II version Git 24041777c5 2021-12-02 20:43:25 -0700
\(\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 - 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_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 <Epetra_Map.h>
30 # include <Epetra_MpiComm.h>
31 # include <Epetra_MultiVector.h>
32 # include <Epetra_RowMatrix.h>
33 # include <Epetra_Vector.h>
34 # include <Teuchos_ParameterList.hpp>
35 
36 # include <memory>
37 
38 // forward declarations
39 # ifndef DOXYGEN
40 class Ifpack_Preconditioner;
41 class Ifpack_Chebyshev;
42 namespace ML_Epetra
43 {
44  class MultiLevelPreconditioner;
45 }
46 # endif
47 
49 
50 // forward declarations
51 # ifndef DOXYGEN
52 template <typename number>
53 class SparseMatrix;
54 template <typename number>
55 class Vector;
56 class SparsityPattern;
57 # endif
58 
63 namespace TrilinosWrappers
64 {
65  // forward declarations
66  class SparseMatrix;
67  class BlockSparseMatrix;
68  class SolverBase;
69 
77  {
78  public:
83 
89  {};
90 
97 
102 
106  ~PreconditionBase() override = default;
107 
112  void
113  clear();
114 
118  MPI_Comm
119  get_mpi_communicator() const;
120 
130  void
131  transpose();
132 
136  virtual void
137  vmult(MPI::Vector &dst, const MPI::Vector &src) const;
138 
142  virtual void
143  Tvmult(MPI::Vector &dst, const MPI::Vector &src) const;
144 
149  virtual void
150  vmult(::Vector<double> &dst, const ::Vector<double> &src) const;
151 
156  virtual void
157  Tvmult(::Vector<double> & dst,
158  const ::Vector<double> &src) const;
159 
164  virtual void
166  const ::LinearAlgebra::distributed::Vector<double> &src) const;
167 
172  virtual void
174  const ::LinearAlgebra::distributed::Vector<double> &src) const;
175 
186  trilinos_operator() const;
188 
193 
198  IndexSet
199  locally_owned_domain_indices() const;
200 
206  IndexSet
207  locally_owned_range_indices() const;
208 
210 
218  DeclException1(ExcNonMatchingMaps,
219  std::string,
220  << "The sparse matrix the preconditioner is based on "
221  << "uses a map that is not compatible to the one in vector "
222  << arg1 << ". Check preconditioner and matrix setup.");
224 
225  friend class SolverBase;
226 
227  protected:
232  Teuchos::RCP<Epetra_Operator> preconditioner;
233 
238  Epetra_MpiComm communicator;
239 
244  std::shared_ptr<Epetra_Map> vector_distributor;
245  };
246 
247 
264  {
265  public:
278  {
283  AdditionalData(const double omega = 1,
284  const double min_diagonal = 0,
285  const unsigned int n_sweeps = 1);
286 
290  double omega;
291 
299  double min_diagonal;
300 
305  unsigned int n_sweeps;
306  };
307 
312  void
313  initialize(const SparseMatrix & matrix,
314  const AdditionalData &additional_data = AdditionalData());
315  };
316 
317 
318 
345  {
346  public:
361  {
368  AdditionalData(const double omega = 1,
369  const double min_diagonal = 0,
370  const unsigned int overlap = 0,
371  const unsigned int n_sweeps = 1);
372 
377  double omega;
378 
386  double min_diagonal;
387 
392  unsigned int overlap;
393 
398  unsigned int n_sweeps;
399  };
400 
406  void
407  initialize(const SparseMatrix & matrix,
408  const AdditionalData &additional_data = AdditionalData());
409  };
410 
411 
412 
439  {
440  public:
455  {
462  AdditionalData(const double omega = 1,
463  const double min_diagonal = 0,
464  const unsigned int overlap = 0,
465  const unsigned int n_sweeps = 1);
466 
471  double omega;
472 
480  double min_diagonal;
481 
486  unsigned int overlap;
487 
492  unsigned int n_sweeps;
493  };
494 
500  void
501  initialize(const SparseMatrix & matrix,
502  const AdditionalData &additional_data = AdditionalData());
503  };
504 
505 
506 
524  {
525  public:
542  {
548  AdditionalData(const unsigned int block_size = 1,
549  const std::string &block_creation_type = "linear",
550  const double omega = 1,
551  const double min_diagonal = 0,
552  const unsigned int n_sweeps = 1);
553 
557  unsigned int block_size;
558 
566  std::string block_creation_type;
567 
572  double omega;
573 
581  double min_diagonal;
582 
587  unsigned int n_sweeps;
588  };
589 
594  void
595  initialize(const SparseMatrix & matrix,
596  const AdditionalData &additional_data = AdditionalData());
597  };
598 
599 
600 
623  {
624  public:
643  {
651  AdditionalData(const unsigned int block_size = 1,
652  const std::string &block_creation_type = "linear",
653  const double omega = 1,
654  const double min_diagonal = 0,
655  const unsigned int overlap = 0,
656  const unsigned int n_sweeps = 1);
657 
661  unsigned int block_size;
662 
670  std::string block_creation_type;
671 
676  double omega;
677 
685  double min_diagonal;
686 
691  unsigned int overlap;
692 
697  unsigned int n_sweeps;
698  };
699 
705  void
706  initialize(const SparseMatrix & matrix,
707  const AdditionalData &additional_data = AdditionalData());
708  };
709 
710 
711 
734  {
735  public:
754  {
762  AdditionalData(const unsigned int block_size = 1,
763  const std::string &block_creation_type = "linear",
764  const double omega = 1,
765  const double min_diagonal = 0,
766  const unsigned int overlap = 0,
767  const unsigned int n_sweeps = 1);
768 
772  unsigned int block_size;
773 
781  std::string block_creation_type;
782 
787  double omega;
788 
796  double min_diagonal;
797 
802  unsigned int overlap;
803 
808  unsigned int n_sweeps;
809  };
810 
816  void
817  initialize(const SparseMatrix & matrix,
818  const AdditionalData &additional_data = AdditionalData());
819  };
820 
821 
822 
858  {
859  public:
877  {
887  AdditionalData(const unsigned int ic_fill = 0,
888  const double ic_atol = 0.,
889  const double ic_rtol = 1.,
890  const unsigned int overlap = 0);
891 
900  unsigned int ic_fill;
901 
907  double ic_atol;
908 
913  double ic_rtol;
914 
919  unsigned int overlap;
920  };
921 
926  void
927  initialize(const SparseMatrix & matrix,
928  const AdditionalData &additional_data = AdditionalData());
929  };
930 
931 
932 
962  {
963  public:
1002  {
1006  AdditionalData(const unsigned int ilu_fill = 0,
1007  const double ilu_atol = 0.,
1008  const double ilu_rtol = 1.,
1009  const unsigned int overlap = 0);
1010 
1014  unsigned int ilu_fill;
1015 
1020  double ilu_atol;
1021 
1026  double ilu_rtol;
1027 
1031  unsigned int overlap;
1032  };
1033 
1038  void
1039  initialize(const SparseMatrix & matrix,
1040  const AdditionalData &additional_data = AdditionalData());
1041  };
1042 
1043 
1044 
1080  {
1081  public:
1100  {
1111  AdditionalData(const double ilut_drop = 0.,
1112  const unsigned int ilut_fill = 0,
1113  const double ilut_atol = 0.,
1114  const double ilut_rtol = 1.,
1115  const unsigned int overlap = 0);
1116 
1121  double ilut_drop;
1122 
1131  unsigned int ilut_fill;
1132 
1138  double ilut_atol;
1139 
1144  double ilut_rtol;
1145 
1150  unsigned int overlap;
1151  };
1152 
1157  void
1158  initialize(const SparseMatrix & matrix,
1159  const AdditionalData &additional_data = AdditionalData());
1160  };
1161 
1162 
1163 
1182  {
1183  public:
1189  {
1193  AdditionalData(const unsigned int overlap = 0);
1194 
1195 
1200  unsigned int overlap;
1201  };
1202 
1207  void
1208  initialize(const SparseMatrix & matrix,
1209  const AdditionalData &additional_data = AdditionalData());
1210  };
1211 
1212 
1213 
1223  {
1224  public:
1230  {
1234  AdditionalData(const unsigned int degree = 1,
1235  const double max_eigenvalue = 10.,
1236  const double eigenvalue_ratio = 30.,
1237  const double min_eigenvalue = 1.,
1238  const double min_diagonal = 1e-12,
1239  const bool nonzero_starting = false);
1240 
1246  unsigned int degree;
1247 
1253 
1258 
1264 
1270 
1282  };
1283 
1288  void
1289  initialize(const SparseMatrix & matrix,
1290  const AdditionalData &additional_data = AdditionalData());
1291  };
1292 
1293 
1294 
1337  {
1338  public:
1346  {
1370  AdditionalData(const bool elliptic = true,
1371  const bool higher_order_elements = false,
1372  const unsigned int n_cycles = 1,
1373  const bool w_cyle = false,
1374  const double aggregation_threshold = 1e-4,
1375  const std::vector<std::vector<bool>> &constant_modes =
1376  std::vector<std::vector<bool>>(0),
1377  const unsigned int smoother_sweeps = 2,
1378  const unsigned int smoother_overlap = 0,
1379  const bool output_details = false,
1380  const char * smoother_type = "Chebyshev",
1381  const char * coarse_type = "Amesos-KLU");
1382 
1413  void
1414  set_parameters(
1415  Teuchos::ParameterList & parameter_list,
1416  std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1417  const Epetra_RowMatrix & matrix) const;
1418 
1426  void
1427  set_parameters(
1428  Teuchos::ParameterList & parameter_list,
1429  std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1430  const SparseMatrix & matrix) const;
1431 
1437  void
1438  set_operator_null_space(
1439  Teuchos::ParameterList & parameter_list,
1440  std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1441  const Epetra_RowMatrix & matrix) const;
1442 
1448  void
1449  set_operator_null_space(
1450  Teuchos::ParameterList & parameter_list,
1451  std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1452  const SparseMatrix & matrix) const;
1453 
1461  bool elliptic;
1462 
1468 
1473  unsigned int n_cycles;
1474 
1479  bool w_cycle;
1480 
1491 
1508  std::vector<std::vector<bool>> constant_modes;
1509 
1520  unsigned int smoother_sweeps;
1521 
1526  unsigned int smoother_overlap;
1527 
1534 
1569  const char *smoother_type;
1570 
1575  const char *coarse_type;
1576  };
1577 
1581  ~PreconditionAMG() override;
1582 
1583 
1589  void
1590  initialize(const SparseMatrix & matrix,
1591  const AdditionalData &additional_data = AdditionalData());
1592 
1611  void
1612  initialize(const Epetra_RowMatrix &matrix,
1613  const AdditionalData & additional_data = AdditionalData());
1614 
1629  void
1630  initialize(const SparseMatrix & matrix,
1631  const Teuchos::ParameterList &ml_parameters);
1632 
1640  void
1641  initialize(const Epetra_RowMatrix & matrix,
1642  const Teuchos::ParameterList &ml_parameters);
1643 
1650  template <typename number>
1651  void
1652  initialize(const ::SparseMatrix<number> &deal_ii_sparse_matrix,
1653  const AdditionalData &additional_data = AdditionalData(),
1654  const double drop_tolerance = 1e-13,
1655  const ::SparsityPattern *use_this_sparsity = nullptr);
1656 
1669  void
1670  reinit();
1671 
1676  void
1677  clear();
1678 
1682  size_type
1683  memory_consumption() const;
1684 
1685  private:
1689  std::shared_ptr<SparseMatrix> trilinos_matrix;
1690  };
1691 
1692 
1693 
1694 # if defined(DOXYGEN) || defined(DEAL_II_TRILINOS_WITH_MUELU)
1695 
1714  {
1715  public:
1723  {
1728  AdditionalData(const bool elliptic = true,
1729  const unsigned int n_cycles = 1,
1730  const bool w_cyle = false,
1731  const double aggregation_threshold = 1e-4,
1732  const std::vector<std::vector<bool>> &constant_modes =
1733  std::vector<std::vector<bool>>(0),
1734  const unsigned int smoother_sweeps = 2,
1735  const unsigned int smoother_overlap = 0,
1736  const bool output_details = false,
1737  const char * smoother_type = "Chebyshev",
1738  const char * coarse_type = "Amesos-KLU");
1739 
1747  bool elliptic;
1748 
1753  unsigned int n_cycles;
1754 
1759  bool w_cycle;
1760 
1771 
1778  std::vector<std::vector<bool>> constant_modes;
1779 
1790  unsigned int smoother_sweeps;
1791 
1796  unsigned int smoother_overlap;
1797 
1804 
1839  const char *smoother_type;
1840 
1845  const char *coarse_type;
1846  };
1847 
1852 
1856  virtual ~PreconditionAMGMueLu() override = default;
1857 
1863  void
1864  initialize(const SparseMatrix & matrix,
1865  const AdditionalData &additional_data = AdditionalData());
1866 
1873  void
1874  initialize(const Epetra_CrsMatrix &matrix,
1875  const AdditionalData & additional_data = AdditionalData());
1876 
1890  void
1891  initialize(const SparseMatrix & matrix,
1892  Teuchos::ParameterList &muelu_parameters);
1893 
1899  void
1900  initialize(const Epetra_CrsMatrix &matrix,
1901  Teuchos::ParameterList &muelu_parameters);
1902 
1909  template <typename number>
1910  void
1911  initialize(const ::SparseMatrix<number> &deal_ii_sparse_matrix,
1912  const AdditionalData &additional_data = AdditionalData(),
1913  const double drop_tolerance = 1e-13,
1914  const ::SparsityPattern *use_this_sparsity = nullptr);
1915 
1920  void
1921  clear();
1922 
1926  size_type
1927  memory_consumption() const;
1928 
1929  private:
1933  std::shared_ptr<SparseMatrix> trilinos_matrix;
1934  };
1935 # endif
1936 
1937 
1938 
1946  {
1947  public:
1953  {};
1954 
1961  void
1962  initialize(const SparseMatrix & matrix,
1963  const AdditionalData &additional_data = AdditionalData());
1964 
1968  void
1969  vmult(MPI::Vector &dst, const MPI::Vector &src) const override;
1970 
1974  void
1975  Tvmult(MPI::Vector &dst, const MPI::Vector &src) const override;
1976 
1981  void
1982  vmult(::Vector<double> & dst,
1983  const ::Vector<double> &src) const override;
1984 
1989  void
1990  Tvmult(::Vector<double> & dst,
1991  const ::Vector<double> &src) const override;
1992 
1997  void
1999  const ::LinearAlgebra::distributed::Vector<double> &src)
2000  const override;
2001 
2007  void
2009  const ::LinearAlgebra::distributed::Vector<double> &src)
2010  const override;
2011  };
2012 
2013 
2014 
2015  // ----------------------- inline and template functions --------------------
2016 
2017 
2018 # ifndef DOXYGEN
2019 
2020 
2021  inline void
2023  {
2024  // This only flips a flag that tells
2025  // Trilinos that any vmult operation
2026  // should be done with the
2027  // transpose. However, the matrix
2028  // structure is not reset.
2029  int ierr;
2030 
2031  if (!preconditioner->UseTranspose())
2032  {
2033  ierr = preconditioner->SetUseTranspose(true);
2034  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2035  }
2036  else
2037  {
2038  ierr = preconditioner->SetUseTranspose(false);
2039  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2040  }
2041  }
2042 
2043 
2044  inline void
2045  PreconditionBase::vmult(MPI::Vector &dst, const MPI::Vector &src) const
2046  {
2047  Assert(dst.trilinos_partitioner().SameAs(
2048  preconditioner->OperatorRangeMap()),
2049  ExcNonMatchingMaps("dst"));
2050  Assert(src.trilinos_partitioner().SameAs(
2051  preconditioner->OperatorDomainMap()),
2052  ExcNonMatchingMaps("src"));
2053 
2054  const int ierr = preconditioner->ApplyInverse(src.trilinos_vector(),
2055  dst.trilinos_vector());
2056  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2057  }
2058 
2059  inline void
2060  PreconditionBase::Tvmult(MPI::Vector &dst, const MPI::Vector &src) const
2061  {
2062  Assert(dst.trilinos_partitioner().SameAs(
2063  preconditioner->OperatorRangeMap()),
2064  ExcNonMatchingMaps("dst"));
2065  Assert(src.trilinos_partitioner().SameAs(
2066  preconditioner->OperatorDomainMap()),
2067  ExcNonMatchingMaps("src"));
2068 
2069  preconditioner->SetUseTranspose(true);
2070  const int ierr = preconditioner->ApplyInverse(src.trilinos_vector(),
2071  dst.trilinos_vector());
2072  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2073  preconditioner->SetUseTranspose(false);
2074  }
2075 
2076 
2077  // For the implementation of the <code>vmult</code> function with deal.II
2078  // data structures we note that invoking a call of the Trilinos
2079  // preconditioner requires us to use Epetra vectors as well. We do this by
2080  // providing a view, i.e., feed Trilinos with a pointer to the data, so we
2081  // avoid copying the content of the vectors during the iteration (this
2082  // function is only useful when used in serial anyway). In the declaration
2083  // of the right hand side, we need to cast the source vector (that is
2084  // <code>const</code> in all deal.II calls) to non-constant value, as this
2085  // is the way Trilinos wants to have them.
2086  inline void
2087  PreconditionBase::vmult(::Vector<double> & dst,
2088  const ::Vector<double> &src) const
2089  {
2090  AssertDimension(dst.size(),
2091  preconditioner->OperatorDomainMap().NumMyElements());
2092  AssertDimension(src.size(),
2093  preconditioner->OperatorRangeMap().NumMyElements());
2094  Epetra_Vector tril_dst(View,
2095  preconditioner->OperatorDomainMap(),
2096  dst.begin());
2097  Epetra_Vector tril_src(View,
2098  preconditioner->OperatorRangeMap(),
2099  const_cast<double *>(src.begin()));
2100 
2101  const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2102  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2103  }
2104 
2105 
2106  inline void
2107  PreconditionBase::Tvmult(::Vector<double> & dst,
2108  const ::Vector<double> &src) const
2109  {
2110  AssertDimension(dst.size(),
2111  preconditioner->OperatorDomainMap().NumMyElements());
2112  AssertDimension(src.size(),
2113  preconditioner->OperatorRangeMap().NumMyElements());
2114  Epetra_Vector tril_dst(View,
2115  preconditioner->OperatorDomainMap(),
2116  dst.begin());
2117  Epetra_Vector tril_src(View,
2118  preconditioner->OperatorRangeMap(),
2119  const_cast<double *>(src.begin()));
2120 
2121  preconditioner->SetUseTranspose(true);
2122  const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2123  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2124  preconditioner->SetUseTranspose(false);
2125  }
2126 
2127 
2128 
2129  inline void
2130  PreconditionBase::vmult(
2133  {
2135  preconditioner->OperatorDomainMap().NumMyElements());
2137  preconditioner->OperatorRangeMap().NumMyElements());
2138  Epetra_Vector tril_dst(View,
2139  preconditioner->OperatorDomainMap(),
2140  dst.begin());
2141  Epetra_Vector tril_src(View,
2142  preconditioner->OperatorRangeMap(),
2143  const_cast<double *>(src.begin()));
2144 
2145  const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2146  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2147  }
2148 
2149  inline void
2150  PreconditionBase::Tvmult(
2153  {
2155  preconditioner->OperatorDomainMap().NumMyElements());
2157  preconditioner->OperatorRangeMap().NumMyElements());
2158  Epetra_Vector tril_dst(View,
2159  preconditioner->OperatorDomainMap(),
2160  dst.begin());
2161  Epetra_Vector tril_src(View,
2162  preconditioner->OperatorRangeMap(),
2163  const_cast<double *>(src.begin()));
2164 
2165  preconditioner->SetUseTranspose(true);
2166  const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2167  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2168  preconditioner->SetUseTranspose(false);
2169  }
2170 
2171 # endif
2172 
2173 } // namespace TrilinosWrappers
2174 
2175 
2180 
2181 # endif // DEAL_II_WITH_TRILINOS
2182 
2183 #endif // trilinos_precondition_h
2184 /*------------------------- trilinos_precondition.h -------------------------*/
static ::ExceptionBase & ExcTrilinosError(int arg1)
typename ::internal::FEInterfaceViews::ViewType< dim, spacedim, Extractor >::type View
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1655
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:1571
size_type locally_owned_size() const
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
#define Assert(cond, exc)
Definition: exceptions.h:1461
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:404
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:403
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