Reference documentation for deal.II version Git 37fff01519 2020-04-07 14:44:07 +0200
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
trilinos_precondition.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 2019 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 
24 # include <deal.II/base/subscriptor.h>
25 
26 # include <deal.II/lac/la_parallel_vector.h>
27 # include <deal.II/lac/trilinos_vector.h>
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 
83  {
84  public:
89 
95  {};
96 
103 
108 
112  ~PreconditionBase() override = default;
113 
118  void
119  clear();
120 
124  MPI_Comm
125  get_mpi_communicator() const;
126 
136  void
137  transpose();
138 
142  virtual void
143  vmult(MPI::Vector &dst, const MPI::Vector &src) const;
144 
148  virtual void
149  Tvmult(MPI::Vector &dst, const MPI::Vector &src) const;
150 
155  virtual void
156  vmult(::Vector<double> &dst, const ::Vector<double> &src) const;
157 
162  virtual void
163  Tvmult(::Vector<double> & dst,
164  const ::Vector<double> &src) const;
165 
170  virtual void
172  const ::LinearAlgebra::distributed::Vector<double> &src) const;
173 
178  virtual void
180  const ::LinearAlgebra::distributed::Vector<double> &src) const;
181 
191  Epetra_Operator &
192  trilinos_operator() const;
194 
199 
204  IndexSet
205  locally_owned_domain_indices() const;
206 
212  IndexSet
213  locally_owned_range_indices() const;
214 
216 
225  DeclException1(ExcNonMatchingMaps,
226  std::string,
227  << "The sparse matrix the preconditioner is based on "
228  << "uses a map that is not compatible to the one in vector "
229  << arg1 << ". Check preconditioner and matrix setup.");
231 
232  friend class SolverBase;
233 
234  protected:
239  Teuchos::RCP<Epetra_Operator> preconditioner;
240 
245 # ifdef DEAL_II_WITH_MPI
246  Epetra_MpiComm communicator;
247 # else
248  Epetra_SerialComm communicator;
249 # endif
250 
255  std::shared_ptr<Epetra_Map> vector_distributor;
256  };
257 
258 
276  {
277  public:
290  {
295  AdditionalData(const double omega = 1,
296  const double min_diagonal = 0,
297  const unsigned int n_sweeps = 1);
298 
302  double omega;
303 
311  double min_diagonal;
312 
317  unsigned int n_sweeps;
318  };
319 
324  void
325  initialize(const SparseMatrix & matrix,
326  const AdditionalData &additional_data = AdditionalData());
327  };
328 
329 
330 
358  {
359  public:
374  {
381  AdditionalData(const double omega = 1,
382  const double min_diagonal = 0,
383  const unsigned int overlap = 0,
384  const unsigned int n_sweeps = 1);
385 
390  double omega;
391 
399  double min_diagonal;
400 
405  unsigned int overlap;
406 
411  unsigned int n_sweeps;
412  };
413 
419  void
420  initialize(const SparseMatrix & matrix,
421  const AdditionalData &additional_data = AdditionalData());
422  };
423 
424 
425 
453  {
454  public:
469  {
476  AdditionalData(const double omega = 1,
477  const double min_diagonal = 0,
478  const unsigned int overlap = 0,
479  const unsigned int n_sweeps = 1);
480 
485  double omega;
486 
494  double min_diagonal;
495 
500  unsigned int overlap;
501 
506  unsigned int n_sweeps;
507  };
508 
514  void
515  initialize(const SparseMatrix & matrix,
516  const AdditionalData &additional_data = AdditionalData());
517  };
518 
519 
520 
539  {
540  public:
559  {
565  AdditionalData(const unsigned int block_size = 1,
566  const std::string &block_creation_type = "linear",
567  const double omega = 1,
568  const double min_diagonal = 0,
569  const unsigned int n_sweeps = 1);
570 
574  unsigned int block_size;
575 
583  std::string block_creation_type;
584 
589  double omega;
590 
598  double min_diagonal;
599 
604  unsigned int n_sweeps;
605  };
606 
611  void
612  initialize(const SparseMatrix & matrix,
613  const AdditionalData &additional_data = AdditionalData());
614  };
615 
616 
617 
641  {
642  public:
661  {
669  AdditionalData(const unsigned int block_size = 1,
670  const std::string &block_creation_type = "linear",
671  const double omega = 1,
672  const double min_diagonal = 0,
673  const unsigned int overlap = 0,
674  const unsigned int n_sweeps = 1);
675 
679  unsigned int block_size;
680 
688  std::string block_creation_type;
689 
694  double omega;
695 
703  double min_diagonal;
704 
709  unsigned int overlap;
710 
715  unsigned int n_sweeps;
716  };
717 
723  void
724  initialize(const SparseMatrix & matrix,
725  const AdditionalData &additional_data = AdditionalData());
726  };
727 
728 
729 
753  {
754  public:
773  {
781  AdditionalData(const unsigned int block_size = 1,
782  const std::string &block_creation_type = "linear",
783  const double omega = 1,
784  const double min_diagonal = 0,
785  const unsigned int overlap = 0,
786  const unsigned int n_sweeps = 1);
787 
791  unsigned int block_size;
792 
800  std::string block_creation_type;
801 
806  double omega;
807 
815  double min_diagonal;
816 
821  unsigned int overlap;
822 
827  unsigned int n_sweeps;
828  };
829 
835  void
836  initialize(const SparseMatrix & matrix,
837  const AdditionalData &additional_data = AdditionalData());
838  };
839 
840 
841 
878  {
879  public:
897  {
907  AdditionalData(const unsigned int ic_fill = 0,
908  const double ic_atol = 0.,
909  const double ic_rtol = 1.,
910  const unsigned int overlap = 0);
911 
920  unsigned int ic_fill;
921 
927  double ic_atol;
928 
933  double ic_rtol;
934 
939  unsigned int overlap;
940  };
941 
946  void
947  initialize(const SparseMatrix & matrix,
948  const AdditionalData &additional_data = AdditionalData());
949  };
950 
951 
952 
983  {
984  public:
1023  {
1027  AdditionalData(const unsigned int ilu_fill = 0,
1028  const double ilu_atol = 0.,
1029  const double ilu_rtol = 1.,
1030  const unsigned int overlap = 0);
1031 
1035  unsigned int ilu_fill;
1036 
1041  double ilu_atol;
1042 
1047  double ilu_rtol;
1048 
1052  unsigned int overlap;
1053  };
1054 
1059  void
1060  initialize(const SparseMatrix & matrix,
1061  const AdditionalData &additional_data = AdditionalData());
1062  };
1063 
1064 
1065 
1102  {
1103  public:
1122  {
1133  AdditionalData(const double ilut_drop = 0.,
1134  const unsigned int ilut_fill = 0,
1135  const double ilut_atol = 0.,
1136  const double ilut_rtol = 1.,
1137  const unsigned int overlap = 0);
1138 
1143  double ilut_drop;
1144 
1153  unsigned int ilut_fill;
1154 
1160  double ilut_atol;
1161 
1166  double ilut_rtol;
1167 
1172  unsigned int overlap;
1173  };
1174 
1179  void
1180  initialize(const SparseMatrix & matrix,
1181  const AdditionalData &additional_data = AdditionalData());
1182  };
1183 
1184 
1185 
1205  {
1206  public:
1212  {
1216  AdditionalData(const unsigned int overlap = 0);
1217 
1218 
1223  unsigned int overlap;
1224  };
1225 
1230  void
1231  initialize(const SparseMatrix & matrix,
1232  const AdditionalData &additional_data = AdditionalData());
1233  };
1234 
1235 
1236 
1247  {
1248  public:
1254  {
1258  AdditionalData(const unsigned int degree = 1,
1259  const double max_eigenvalue = 10.,
1260  const double eigenvalue_ratio = 30.,
1261  const double min_eigenvalue = 1.,
1262  const double min_diagonal = 1e-12,
1263  const bool nonzero_starting = false);
1264 
1270  unsigned int degree;
1271 
1277 
1282 
1288 
1294 
1306  };
1307 
1312  void
1313  initialize(const SparseMatrix & matrix,
1314  const AdditionalData &additional_data = AdditionalData());
1315  };
1316 
1317 
1318 
1362  {
1363  public:
1371  {
1395  AdditionalData(const bool elliptic = true,
1396  const bool higher_order_elements = false,
1397  const unsigned int n_cycles = 1,
1398  const bool w_cyle = false,
1399  const double aggregation_threshold = 1e-4,
1400  const std::vector<std::vector<bool>> &constant_modes =
1401  std::vector<std::vector<bool>>(0),
1402  const unsigned int smoother_sweeps = 2,
1403  const unsigned int smoother_overlap = 0,
1404  const bool output_details = false,
1405  const char * smoother_type = "Chebyshev",
1406  const char * coarse_type = "Amesos-KLU");
1407 
1438  void
1439  set_parameters(
1440  Teuchos::ParameterList & parameter_list,
1441  std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1442  const Epetra_RowMatrix & matrix) const;
1443 
1451  void
1452  set_parameters(
1453  Teuchos::ParameterList & parameter_list,
1454  std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1455  const SparseMatrix & matrix) const;
1456 
1462  void
1463  set_operator_null_space(
1464  Teuchos::ParameterList & parameter_list,
1465  std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1466  const Epetra_RowMatrix & matrix) const;
1467 
1473  void
1474  set_operator_null_space(
1475  Teuchos::ParameterList & parameter_list,
1476  std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1477  const SparseMatrix & matrix) const;
1478 
1486  bool elliptic;
1487 
1493 
1498  unsigned int n_cycles;
1499 
1504  bool w_cycle;
1505 
1516 
1533  std::vector<std::vector<bool>> constant_modes;
1534 
1545  unsigned int smoother_sweeps;
1546 
1551  unsigned int smoother_overlap;
1552 
1559 
1594  const char *smoother_type;
1595 
1600  const char *coarse_type;
1601  };
1602 
1606  ~PreconditionAMG() override;
1607 
1608 
1614  void
1615  initialize(const SparseMatrix & matrix,
1616  const AdditionalData &additional_data = AdditionalData());
1617 
1636  void
1637  initialize(const Epetra_RowMatrix &matrix,
1638  const AdditionalData & additional_data = AdditionalData());
1639 
1654  void
1655  initialize(const SparseMatrix & matrix,
1656  const Teuchos::ParameterList &ml_parameters);
1657 
1665  void
1666  initialize(const Epetra_RowMatrix & matrix,
1667  const Teuchos::ParameterList &ml_parameters);
1668 
1675  template <typename number>
1676  void
1677  initialize(const ::SparseMatrix<number> &deal_ii_sparse_matrix,
1678  const AdditionalData &additional_data = AdditionalData(),
1679  const double drop_tolerance = 1e-13,
1680  const ::SparsityPattern *use_this_sparsity = nullptr);
1681 
1694  void
1695  reinit();
1696 
1701  void
1702  clear();
1703 
1707  size_type
1708  memory_consumption() const;
1709 
1710  private:
1714  std::shared_ptr<SparseMatrix> trilinos_matrix;
1715  };
1716 
1717 
1718 
1719 # if defined(DOXYGEN) || defined(DEAL_II_TRILINOS_WITH_MUELU)
1720 
1740  {
1741  public:
1749  {
1754  AdditionalData(const bool elliptic = true,
1755  const unsigned int n_cycles = 1,
1756  const bool w_cyle = false,
1757  const double aggregation_threshold = 1e-4,
1758  const std::vector<std::vector<bool>> &constant_modes =
1759  std::vector<std::vector<bool>>(0),
1760  const unsigned int smoother_sweeps = 2,
1761  const unsigned int smoother_overlap = 0,
1762  const bool output_details = false,
1763  const char * smoother_type = "Chebyshev",
1764  const char * coarse_type = "Amesos-KLU");
1765 
1773  bool elliptic;
1774 
1779  unsigned int n_cycles;
1780 
1785  bool w_cycle;
1786 
1797 
1804  std::vector<std::vector<bool>> constant_modes;
1805 
1816  unsigned int smoother_sweeps;
1817 
1822  unsigned int smoother_overlap;
1823 
1830 
1865  const char *smoother_type;
1866 
1871  const char *coarse_type;
1872  };
1873 
1878 
1882  virtual ~PreconditionAMGMueLu() override = default;
1883 
1889  void
1890  initialize(const SparseMatrix & matrix,
1891  const AdditionalData &additional_data = AdditionalData());
1892 
1899  void
1900  initialize(const Epetra_CrsMatrix &matrix,
1901  const AdditionalData & additional_data = AdditionalData());
1902 
1916  void
1917  initialize(const SparseMatrix & matrix,
1918  Teuchos::ParameterList &muelu_parameters);
1919 
1925  void
1926  initialize(const Epetra_CrsMatrix &matrix,
1927  Teuchos::ParameterList &muelu_parameters);
1928 
1935  template <typename number>
1936  void
1937  initialize(const ::SparseMatrix<number> &deal_ii_sparse_matrix,
1938  const AdditionalData &additional_data = AdditionalData(),
1939  const double drop_tolerance = 1e-13,
1940  const ::SparsityPattern *use_this_sparsity = nullptr);
1941 
1946  void
1947  clear();
1948 
1952  size_type
1953  memory_consumption() const;
1954 
1955  private:
1959  std::shared_ptr<SparseMatrix> trilinos_matrix;
1960  };
1961 # endif
1962 
1963 
1964 
1974  {
1975  public:
1981  {};
1982 
1989  void
1990  initialize(const SparseMatrix & matrix,
1991  const AdditionalData &additional_data = AdditionalData());
1992 
1996  void
1997  vmult(MPI::Vector &dst, const MPI::Vector &src) const override;
1998 
2002  void
2003  Tvmult(MPI::Vector &dst, const MPI::Vector &src) const override;
2004 
2009  void
2010  vmult(::Vector<double> & dst,
2011  const ::Vector<double> &src) const override;
2012 
2017  void
2018  Tvmult(::Vector<double> & dst,
2019  const ::Vector<double> &src) const override;
2020 
2025  void
2027  const ::LinearAlgebra::distributed::Vector<double> &src)
2028  const override;
2029 
2035  void
2037  const ::LinearAlgebra::distributed::Vector<double> &src)
2038  const override;
2039  };
2040 
2041 
2042 
2043  // ----------------------- inline and template functions --------------------
2044 
2045 
2046 # ifndef DOXYGEN
2047 
2048 
2049  inline void
2051  {
2052  // This only flips a flag that tells
2053  // Trilinos that any vmult operation
2054  // should be done with the
2055  // transpose. However, the matrix
2056  // structure is not reset.
2057  int ierr;
2058 
2059  if (!preconditioner->UseTranspose())
2060  {
2061  ierr = preconditioner->SetUseTranspose(true);
2062  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2063  }
2064  else
2065  {
2066  ierr = preconditioner->SetUseTranspose(false);
2067  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2068  }
2069  }
2070 
2071 
2072  inline void
2073  PreconditionBase::vmult(MPI::Vector &dst, const MPI::Vector &src) const
2074  {
2075  Assert(dst.trilinos_partitioner().SameAs(
2076  preconditioner->OperatorRangeMap()),
2077  ExcNonMatchingMaps("dst"));
2078  Assert(src.trilinos_partitioner().SameAs(
2079  preconditioner->OperatorDomainMap()),
2080  ExcNonMatchingMaps("src"));
2081 
2082  const int ierr = preconditioner->ApplyInverse(src.trilinos_vector(),
2083  dst.trilinos_vector());
2084  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2085  }
2086 
2087  inline void
2088  PreconditionBase::Tvmult(MPI::Vector &dst, const MPI::Vector &src) const
2089  {
2090  Assert(dst.trilinos_partitioner().SameAs(
2091  preconditioner->OperatorRangeMap()),
2092  ExcNonMatchingMaps("dst"));
2093  Assert(src.trilinos_partitioner().SameAs(
2094  preconditioner->OperatorDomainMap()),
2095  ExcNonMatchingMaps("src"));
2096 
2097  preconditioner->SetUseTranspose(true);
2098  const int ierr = preconditioner->ApplyInverse(src.trilinos_vector(),
2099  dst.trilinos_vector());
2100  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2101  preconditioner->SetUseTranspose(false);
2102  }
2103 
2104 
2105  // For the implementation of the <code>vmult</code> function with deal.II
2106  // data structures we note that invoking a call of the Trilinos
2107  // preconditioner requires us to use Epetra vectors as well. We do this by
2108  // providing a view, i.e., feed Trilinos with a pointer to the data, so we
2109  // avoid copying the content of the vectors during the iteration (this
2110  // function is only useful when used in serial anyway). In the declaration
2111  // of the right hand side, we need to cast the source vector (that is
2112  // <code>const</code> in all deal.II calls) to non-constant value, as this
2113  // is the way Trilinos wants to have them.
2114  inline void
2115  PreconditionBase::vmult(::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  const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2130  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2131  }
2132 
2133 
2134  inline void
2135  PreconditionBase::Tvmult(::Vector<double> & dst,
2136  const ::Vector<double> &src) const
2137  {
2138  AssertDimension(dst.size(),
2139  preconditioner->OperatorDomainMap().NumMyElements());
2140  AssertDimension(src.size(),
2141  preconditioner->OperatorRangeMap().NumMyElements());
2142  Epetra_Vector tril_dst(View,
2143  preconditioner->OperatorDomainMap(),
2144  dst.begin());
2145  Epetra_Vector tril_src(View,
2146  preconditioner->OperatorRangeMap(),
2147  const_cast<double *>(src.begin()));
2148 
2149  preconditioner->SetUseTranspose(true);
2150  const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2151  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2152  preconditioner->SetUseTranspose(false);
2153  }
2154 
2155 
2156 
2157  inline void
2158  PreconditionBase::vmult(
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  const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2174  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2175  }
2176 
2177  inline void
2178  PreconditionBase::Tvmult(
2181  {
2183  preconditioner->OperatorDomainMap().NumMyElements());
2185  preconditioner->OperatorRangeMap().NumMyElements());
2186  Epetra_Vector tril_dst(View,
2187  preconditioner->OperatorDomainMap(),
2188  dst.begin());
2189  Epetra_Vector tril_src(View,
2190  preconditioner->OperatorRangeMap(),
2191  const_cast<double *>(src.begin()));
2192 
2193  preconditioner->SetUseTranspose(true);
2194  const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2195  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2196  preconditioner->SetUseTranspose(false);
2197  }
2198 
2199 # endif
2200 
2201 } // namespace TrilinosWrappers
2202 
2203 
2207 DEAL_II_NAMESPACE_CLOSE
2208 
2209 # endif // DEAL_II_WITH_TRILINOS
2210 
2211 #endif // trilinos_precondition_h
2212 /*------------------------- trilinos_precondition.h -------------------------*/
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
std::shared_ptr< Epetra_Map > vector_distributor
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
constexpr SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
LinearAlgebra::distributed::Vector< Number > Vector
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:518
#define Assert(cond, exc)
Definition: exceptions.h:1419
iterator begin()
std::shared_ptr< SparseMatrix > trilinos_matrix
std::shared_ptr< SparseMatrix > trilinos_matrix
unsigned int global_dof_index
Definition: types.h:77
const Epetra_BlockMap & trilinos_partitioner() const
size_type size() const
const Epetra_MultiVector & trilinos_vector() const
Teuchos::RCP< Epetra_Operator > preconditioner