Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
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 class Ifpack_Preconditioner;
44 class Ifpack_Chebyshev;
45 namespace ML_Epetra
46 {
47  class MultiLevelPreconditioner;
48 }
49 
50 
51 DEAL_II_NAMESPACE_OPEN
52 
53 // forward declarations
54 template <typename number>
55 class SparseMatrix;
56 template <typename number>
57 class Vector;
58 class SparsityPattern;
59 
64 namespace TrilinosWrappers
65 {
66  // forward declarations
67  class SparseMatrix;
68  class BlockSparseMatrix;
69  class SolverBase;
70 
80  {
81  public:
86 
92  {};
93 
100 
105 
109  ~PreconditionBase() override = default;
110 
115  void
116  clear();
117 
121  MPI_Comm
122  get_mpi_communicator() const;
123 
133  void
134  transpose();
135 
139  virtual void
140  vmult(MPI::Vector &dst, const MPI::Vector &src) const;
141 
145  virtual void
146  Tvmult(MPI::Vector &dst, const MPI::Vector &src) const;
147 
152  virtual void
153  vmult(::Vector<double> &dst, const ::Vector<double> &src) const;
154 
159  virtual void
160  Tvmult(::Vector<double> & dst,
161  const ::Vector<double> &src) const;
162 
167  virtual void
169  const ::LinearAlgebra::distributed::Vector<double> &src) const;
170 
175  virtual void
177  const ::LinearAlgebra::distributed::Vector<double> &src) const;
178 
188  Epetra_Operator &
189  trilinos_operator() const;
191 
196 
201  IndexSet
202  locally_owned_domain_indices() const;
203 
209  IndexSet
210  locally_owned_range_indices() const;
211 
213 
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 
273  {
274  public:
287  {
292  AdditionalData(const double omega = 1,
293  const double min_diagonal = 0,
294  const unsigned int n_sweeps = 1);
295 
299  double omega;
300 
308  double min_diagonal;
309 
314  unsigned int n_sweeps;
315  };
316 
321  void
322  initialize(const SparseMatrix & matrix,
323  const AdditionalData &additional_data = AdditionalData());
324  };
325 
326 
327 
355  {
356  public:
371  {
378  AdditionalData(const double omega = 1,
379  const double min_diagonal = 0,
380  const unsigned int overlap = 0,
381  const unsigned int n_sweeps = 1);
382 
387  double omega;
388 
396  double min_diagonal;
397 
402  unsigned int overlap;
403 
408  unsigned int n_sweeps;
409  };
410 
416  void
417  initialize(const SparseMatrix & matrix,
418  const AdditionalData &additional_data = AdditionalData());
419  };
420 
421 
422 
450  {
451  public:
466  {
473  AdditionalData(const double omega = 1,
474  const double min_diagonal = 0,
475  const unsigned int overlap = 0,
476  const unsigned int n_sweeps = 1);
477 
482  double omega;
483 
491  double min_diagonal;
492 
497  unsigned int overlap;
498 
503  unsigned int n_sweeps;
504  };
505 
511  void
512  initialize(const SparseMatrix & matrix,
513  const AdditionalData &additional_data = AdditionalData());
514  };
515 
516 
517 
536  {
537  public:
556  {
562  AdditionalData(const unsigned int block_size = 1,
563  const std::string &block_creation_type = "linear",
564  const double omega = 1,
565  const double min_diagonal = 0,
566  const unsigned int n_sweeps = 1);
567 
571  unsigned int block_size;
572 
580  std::string block_creation_type;
581 
586  double omega;
587 
595  double min_diagonal;
596 
601  unsigned int n_sweeps;
602  };
603 
608  void
609  initialize(const SparseMatrix & matrix,
610  const AdditionalData &additional_data = AdditionalData());
611  };
612 
613 
614 
638  {
639  public:
658  {
666  AdditionalData(const unsigned int block_size = 1,
667  const std::string &block_creation_type = "linear",
668  const double omega = 1,
669  const double min_diagonal = 0,
670  const unsigned int overlap = 0,
671  const unsigned int n_sweeps = 1);
672 
676  unsigned int block_size;
677 
685  std::string block_creation_type;
686 
691  double omega;
692 
700  double min_diagonal;
701 
706  unsigned int overlap;
707 
712  unsigned int n_sweeps;
713  };
714 
720  void
721  initialize(const SparseMatrix & matrix,
722  const AdditionalData &additional_data = AdditionalData());
723  };
724 
725 
726 
750  {
751  public:
770  {
778  AdditionalData(const unsigned int block_size = 1,
779  const std::string &block_creation_type = "linear",
780  const double omega = 1,
781  const double min_diagonal = 0,
782  const unsigned int overlap = 0,
783  const unsigned int n_sweeps = 1);
784 
788  unsigned int block_size;
789 
797  std::string block_creation_type;
798 
803  double omega;
804 
812  double min_diagonal;
813 
818  unsigned int overlap;
819 
824  unsigned int n_sweeps;
825  };
826 
832  void
833  initialize(const SparseMatrix & matrix,
834  const AdditionalData &additional_data = AdditionalData());
835  };
836 
837 
838 
875  {
876  public:
894  {
904  AdditionalData(const unsigned int ic_fill = 0,
905  const double ic_atol = 0.,
906  const double ic_rtol = 1.,
907  const unsigned int overlap = 0);
908 
917  unsigned int ic_fill;
918 
924  double ic_atol;
925 
930  double ic_rtol;
931 
936  unsigned int overlap;
937  };
938 
943  void
944  initialize(const SparseMatrix & matrix,
945  const AdditionalData &additional_data = AdditionalData());
946  };
947 
948 
949 
980  {
981  public:
1020  {
1024  AdditionalData(const unsigned int ilu_fill = 0,
1025  const double ilu_atol = 0.,
1026  const double ilu_rtol = 1.,
1027  const unsigned int overlap = 0);
1028 
1032  unsigned int ilu_fill;
1033 
1038  double ilu_atol;
1039 
1044  double ilu_rtol;
1045 
1049  unsigned int overlap;
1050  };
1051 
1056  void
1057  initialize(const SparseMatrix & matrix,
1058  const AdditionalData &additional_data = AdditionalData());
1059  };
1060 
1061 
1062 
1099  {
1100  public:
1119  {
1130  AdditionalData(const double ilut_drop = 0.,
1131  const unsigned int ilut_fill = 0,
1132  const double ilut_atol = 0.,
1133  const double ilut_rtol = 1.,
1134  const unsigned int overlap = 0);
1135 
1140  double ilut_drop;
1141 
1150  unsigned int ilut_fill;
1151 
1157  double ilut_atol;
1158 
1163  double ilut_rtol;
1164 
1169  unsigned int overlap;
1170  };
1171 
1176  void
1177  initialize(const SparseMatrix & matrix,
1178  const AdditionalData &additional_data = AdditionalData());
1179  };
1180 
1181 
1182 
1202  {
1203  public:
1209  {
1213  AdditionalData(const unsigned int overlap = 0);
1214 
1215 
1220  unsigned int overlap;
1221  };
1222 
1227  void
1228  initialize(const SparseMatrix & matrix,
1229  const AdditionalData &additional_data = AdditionalData());
1230  };
1231 
1232 
1233 
1244  {
1245  public:
1251  {
1255  AdditionalData(const unsigned int degree = 1,
1256  const double max_eigenvalue = 10.,
1257  const double eigenvalue_ratio = 30.,
1258  const double min_eigenvalue = 1.,
1259  const double min_diagonal = 1e-12,
1260  const bool nonzero_starting = false);
1261 
1267  unsigned int degree;
1268 
1274 
1279 
1285 
1291 
1303  };
1304 
1309  void
1310  initialize(const SparseMatrix & matrix,
1311  const AdditionalData &additional_data = AdditionalData());
1312  };
1313 
1314 
1315 
1359  {
1360  public:
1368  {
1392  AdditionalData(const bool elliptic = true,
1393  const bool higher_order_elements = false,
1394  const unsigned int n_cycles = 1,
1395  const bool w_cyle = false,
1396  const double aggregation_threshold = 1e-4,
1397  const std::vector<std::vector<bool>> &constant_modes =
1398  std::vector<std::vector<bool>>(0),
1399  const unsigned int smoother_sweeps = 2,
1400  const unsigned int smoother_overlap = 0,
1401  const bool output_details = false,
1402  const char * smoother_type = "Chebyshev",
1403  const char * coarse_type = "Amesos-KLU");
1404 
1435  void
1436  set_parameters(
1437  Teuchos::ParameterList & parameter_list,
1438  std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1439  const Epetra_RowMatrix & matrix) const;
1440 
1448  void
1449  set_parameters(
1450  Teuchos::ParameterList & parameter_list,
1451  std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1452  const SparseMatrix & matrix) const;
1453 
1459  void
1460  set_operator_null_space(
1461  Teuchos::ParameterList & parameter_list,
1462  std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1463  const Epetra_RowMatrix & matrix) const;
1464 
1470  void
1471  set_operator_null_space(
1472  Teuchos::ParameterList & parameter_list,
1473  std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1474  const SparseMatrix & matrix) const;
1475 
1483  bool elliptic;
1484 
1490 
1495  unsigned int n_cycles;
1496 
1501  bool w_cycle;
1502 
1513 
1530  std::vector<std::vector<bool>> constant_modes;
1531 
1542  unsigned int smoother_sweeps;
1543 
1548  unsigned int smoother_overlap;
1549 
1556 
1591  const char *smoother_type;
1592 
1597  const char *coarse_type;
1598  };
1599 
1603  ~PreconditionAMG() override;
1604 
1605 
1611  void
1612  initialize(const SparseMatrix & matrix,
1613  const AdditionalData &additional_data = AdditionalData());
1614 
1633  void
1634  initialize(const Epetra_RowMatrix &matrix,
1635  const AdditionalData & additional_data = AdditionalData());
1636 
1651  void
1652  initialize(const SparseMatrix & matrix,
1653  const Teuchos::ParameterList &ml_parameters);
1654 
1662  void
1663  initialize(const Epetra_RowMatrix & matrix,
1664  const Teuchos::ParameterList &ml_parameters);
1665 
1672  template <typename number>
1673  void
1674  initialize(const ::SparseMatrix<number> &deal_ii_sparse_matrix,
1675  const AdditionalData &additional_data = AdditionalData(),
1676  const double drop_tolerance = 1e-13,
1677  const ::SparsityPattern *use_this_sparsity = nullptr);
1678 
1691  void
1692  reinit();
1693 
1698  void
1699  clear();
1700 
1704  size_type
1705  memory_consumption() const;
1706 
1707  private:
1711  std::shared_ptr<SparseMatrix> trilinos_matrix;
1712  };
1713 
1714 
1715 
1716 # if defined(DOXYGEN) || DEAL_II_TRILINOS_VERSION_GTE(11, 14, 0)
1717 
1735  {
1736  public:
1744  {
1749  AdditionalData(const bool elliptic = true,
1750  const unsigned int n_cycles = 1,
1751  const bool w_cyle = false,
1752  const double aggregation_threshold = 1e-4,
1753  const std::vector<std::vector<bool>> &constant_modes =
1754  std::vector<std::vector<bool>>(0),
1755  const unsigned int smoother_sweeps = 2,
1756  const unsigned int smoother_overlap = 0,
1757  const bool output_details = false,
1758  const char * smoother_type = "Chebyshev",
1759  const char * coarse_type = "Amesos-KLU");
1760 
1768  bool elliptic;
1769 
1774  unsigned int n_cycles;
1775 
1780  bool w_cycle;
1781 
1792 
1799  std::vector<std::vector<bool>> constant_modes;
1800 
1811  unsigned int smoother_sweeps;
1812 
1817  unsigned int smoother_overlap;
1818 
1825 
1860  const char *smoother_type;
1861 
1866  const char *coarse_type;
1867  };
1868 
1873 
1877  virtual ~PreconditionAMGMueLu() override = default;
1878 
1884  void
1885  initialize(const SparseMatrix & matrix,
1886  const AdditionalData &additional_data = AdditionalData());
1887 
1894  void
1895  initialize(const Epetra_CrsMatrix &matrix,
1896  const AdditionalData & additional_data = AdditionalData());
1897 
1911  void
1912  initialize(const SparseMatrix & matrix,
1913  Teuchos::ParameterList &muelu_parameters);
1914 
1920  void
1921  initialize(const Epetra_CrsMatrix &matrix,
1922  Teuchos::ParameterList &muelu_parameters);
1923 
1930  template <typename number>
1931  void
1932  initialize(const ::SparseMatrix<number> &deal_ii_sparse_matrix,
1933  const AdditionalData &additional_data = AdditionalData(),
1934  const double drop_tolerance = 1e-13,
1935  const ::SparsityPattern *use_this_sparsity = nullptr);
1936 
1941  void
1942  clear();
1943 
1947  size_type
1948  memory_consumption() const;
1949 
1950  private:
1954  std::shared_ptr<SparseMatrix> trilinos_matrix;
1955  };
1956 # endif
1957 
1958 
1959 
1969  {
1970  public:
1976  {};
1977 
1984  void
1985  initialize(const SparseMatrix & matrix,
1986  const AdditionalData &additional_data = AdditionalData());
1987 
1991  void
1992  vmult(MPI::Vector &dst, const MPI::Vector &src) const override;
1993 
1997  void
1998  Tvmult(MPI::Vector &dst, const MPI::Vector &src) const override;
1999 
2004  void
2005  vmult(::Vector<double> & dst,
2006  const ::Vector<double> &src) const override;
2007 
2012  void
2013  Tvmult(::Vector<double> & dst,
2014  const ::Vector<double> &src) const override;
2015 
2020  void
2022  const ::LinearAlgebra::distributed::Vector<double> &src)
2023  const override;
2024 
2030  void
2032  const ::LinearAlgebra::distributed::Vector<double> &src)
2033  const override;
2034  };
2035 
2036 
2037 
2038  // ----------------------- inline and template functions --------------------
2039 
2040 
2041 # ifndef DOXYGEN
2042 
2043 
2044  inline void
2046  {
2047  // This only flips a flag that tells
2048  // Trilinos that any vmult operation
2049  // should be done with the
2050  // transpose. However, the matrix
2051  // structure is not reset.
2052  int ierr;
2053 
2054  if (!preconditioner->UseTranspose())
2055  {
2056  ierr = preconditioner->SetUseTranspose(true);
2057  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2058  }
2059  else
2060  {
2061  ierr = preconditioner->SetUseTranspose(false);
2062  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2063  }
2064  }
2065 
2066 
2067  inline void
2068  PreconditionBase::vmult(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  const int ierr = preconditioner->ApplyInverse(src.trilinos_vector(),
2078  dst.trilinos_vector());
2079  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2080  }
2081 
2082  inline void
2083  PreconditionBase::Tvmult(MPI::Vector &dst, const MPI::Vector &src) const
2084  {
2085  Assert(dst.trilinos_partitioner().SameAs(
2086  preconditioner->OperatorRangeMap()),
2087  ExcNonMatchingMaps("dst"));
2088  Assert(src.trilinos_partitioner().SameAs(
2089  preconditioner->OperatorDomainMap()),
2090  ExcNonMatchingMaps("src"));
2091 
2092  preconditioner->SetUseTranspose(true);
2093  const int ierr = preconditioner->ApplyInverse(src.trilinos_vector(),
2094  dst.trilinos_vector());
2095  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2096  preconditioner->SetUseTranspose(false);
2097  }
2098 
2099 
2100  // For the implementation of the <code>vmult</code> function with deal.II
2101  // data structures we note that invoking a call of the Trilinos
2102  // preconditioner requires us to use Epetra vectors as well. We do this by
2103  // providing a view, i.e., feed Trilinos with a pointer to the data, so we
2104  // avoid copying the content of the vectors during the iteration (this
2105  // function is only useful when used in serial anyway). In the declaration
2106  // of the right hand side, we need to cast the source vector (that is
2107  // <code>const</code> in all deal.II calls) to non-constant value, as this
2108  // is the way Trilinos wants to have them.
2109  inline void
2110  PreconditionBase::vmult(::Vector<double> & dst,
2111  const ::Vector<double> &src) const
2112  {
2113  AssertDimension(dst.size(),
2114  preconditioner->OperatorDomainMap().NumMyElements());
2115  AssertDimension(src.size(),
2116  preconditioner->OperatorRangeMap().NumMyElements());
2117  Epetra_Vector tril_dst(View,
2118  preconditioner->OperatorDomainMap(),
2119  dst.begin());
2120  Epetra_Vector tril_src(View,
2121  preconditioner->OperatorRangeMap(),
2122  const_cast<double *>(src.begin()));
2123 
2124  const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2125  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2126  }
2127 
2128 
2129  inline void
2130  PreconditionBase::Tvmult(::Vector<double> & dst,
2131  const ::Vector<double> &src) const
2132  {
2133  AssertDimension(dst.size(),
2134  preconditioner->OperatorDomainMap().NumMyElements());
2135  AssertDimension(src.size(),
2136  preconditioner->OperatorRangeMap().NumMyElements());
2137  Epetra_Vector tril_dst(View,
2138  preconditioner->OperatorDomainMap(),
2139  dst.begin());
2140  Epetra_Vector tril_src(View,
2141  preconditioner->OperatorRangeMap(),
2142  const_cast<double *>(src.begin()));
2143 
2144  preconditioner->SetUseTranspose(true);
2145  const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2146  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2147  preconditioner->SetUseTranspose(false);
2148  }
2149 
2150 
2151 
2152  inline void
2153  PreconditionBase::vmult(
2156  {
2158  preconditioner->OperatorDomainMap().NumMyElements());
2160  preconditioner->OperatorRangeMap().NumMyElements());
2161  Epetra_Vector tril_dst(View,
2162  preconditioner->OperatorDomainMap(),
2163  dst.begin());
2164  Epetra_Vector tril_src(View,
2165  preconditioner->OperatorRangeMap(),
2166  const_cast<double *>(src.begin()));
2167 
2168  const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2169  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2170  }
2171 
2172  inline void
2173  PreconditionBase::Tvmult(
2176  {
2178  preconditioner->OperatorDomainMap().NumMyElements());
2180  preconditioner->OperatorRangeMap().NumMyElements());
2181  Epetra_Vector tril_dst(View,
2182  preconditioner->OperatorDomainMap(),
2183  dst.begin());
2184  Epetra_Vector tril_src(View,
2185  preconditioner->OperatorRangeMap(),
2186  const_cast<double *>(src.begin()));
2187 
2188  preconditioner->SetUseTranspose(true);
2189  const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2190  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2191  preconditioner->SetUseTranspose(false);
2192  }
2193 
2194 # endif
2195 
2196 } // namespace TrilinosWrappers
2197 
2198 
2202 DEAL_II_NAMESPACE_CLOSE
2203 
2204 # endif // DEAL_II_WITH_TRILINOS
2205 
2206 #endif // trilinos_precondition_h
2207 /*------------------------- trilinos_precondition.h -------------------------*/
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
std::shared_ptr< Epetra_Map > vector_distributor
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
LinearAlgebra::distributed::Vector< Number > Vector
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:518
#define Assert(cond, exc)
Definition: exceptions.h:1407
SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
iterator begin()
std::shared_ptr< SparseMatrix > trilinos_matrix
std::shared_ptr< SparseMatrix > trilinos_matrix
unsigned int global_dof_index
Definition: types.h:89
const Epetra_BlockMap & trilinos_partitioner() const
size_type size() const
const Epetra_MultiVector & trilinos_vector() const
Teuchos::RCP< Epetra_Operator > preconditioner