Reference documentation for deal.II version Git 7baca85e57 2019-08-23 17:29:05 -0400
\(\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 # ifndef DOXYGEN
44 class Ifpack_Preconditioner;
45 class Ifpack_Chebyshev;
46 namespace ML_Epetra
47 {
48  class MultiLevelPreconditioner;
49 }
50 # endif
51 
52 DEAL_II_NAMESPACE_OPEN
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) || DEAL_II_TRILINOS_VERSION_GTE(11, 14, 0)
1720 
1738  {
1739  public:
1747  {
1752  AdditionalData(const bool elliptic = true,
1753  const unsigned int n_cycles = 1,
1754  const bool w_cyle = false,
1755  const double aggregation_threshold = 1e-4,
1756  const std::vector<std::vector<bool>> &constant_modes =
1757  std::vector<std::vector<bool>>(0),
1758  const unsigned int smoother_sweeps = 2,
1759  const unsigned int smoother_overlap = 0,
1760  const bool output_details = false,
1761  const char * smoother_type = "Chebyshev",
1762  const char * coarse_type = "Amesos-KLU");
1763 
1771  bool elliptic;
1772 
1777  unsigned int n_cycles;
1778 
1783  bool w_cycle;
1784 
1795 
1802  std::vector<std::vector<bool>> constant_modes;
1803 
1814  unsigned int smoother_sweeps;
1815 
1820  unsigned int smoother_overlap;
1821 
1828 
1863  const char *smoother_type;
1864 
1869  const char *coarse_type;
1870  };
1871 
1876 
1880  virtual ~PreconditionAMGMueLu() override = default;
1881 
1887  void
1888  initialize(const SparseMatrix & matrix,
1889  const AdditionalData &additional_data = AdditionalData());
1890 
1897  void
1898  initialize(const Epetra_CrsMatrix &matrix,
1899  const AdditionalData & additional_data = AdditionalData());
1900 
1914  void
1915  initialize(const SparseMatrix & matrix,
1916  Teuchos::ParameterList &muelu_parameters);
1917 
1923  void
1924  initialize(const Epetra_CrsMatrix &matrix,
1925  Teuchos::ParameterList &muelu_parameters);
1926 
1933  template <typename number>
1934  void
1935  initialize(const ::SparseMatrix<number> &deal_ii_sparse_matrix,
1936  const AdditionalData &additional_data = AdditionalData(),
1937  const double drop_tolerance = 1e-13,
1938  const ::SparsityPattern *use_this_sparsity = nullptr);
1939 
1944  void
1945  clear();
1946 
1950  size_type
1951  memory_consumption() const;
1952 
1953  private:
1957  std::shared_ptr<SparseMatrix> trilinos_matrix;
1958  };
1959 # endif
1960 
1961 
1962 
1972  {
1973  public:
1979  {};
1980 
1987  void
1988  initialize(const SparseMatrix & matrix,
1989  const AdditionalData &additional_data = AdditionalData());
1990 
1994  void
1995  vmult(MPI::Vector &dst, const MPI::Vector &src) const override;
1996 
2000  void
2001  Tvmult(MPI::Vector &dst, const MPI::Vector &src) const override;
2002 
2007  void
2008  vmult(::Vector<double> & dst,
2009  const ::Vector<double> &src) const override;
2010 
2015  void
2016  Tvmult(::Vector<double> & dst,
2017  const ::Vector<double> &src) const override;
2018 
2023  void
2025  const ::LinearAlgebra::distributed::Vector<double> &src)
2026  const override;
2027 
2033  void
2035  const ::LinearAlgebra::distributed::Vector<double> &src)
2036  const override;
2037  };
2038 
2039 
2040 
2041  // ----------------------- inline and template functions --------------------
2042 
2043 
2044 # ifndef DOXYGEN
2045 
2046 
2047  inline void
2049  {
2050  // This only flips a flag that tells
2051  // Trilinos that any vmult operation
2052  // should be done with the
2053  // transpose. However, the matrix
2054  // structure is not reset.
2055  int ierr;
2056 
2057  if (!preconditioner->UseTranspose())
2058  {
2059  ierr = preconditioner->SetUseTranspose(true);
2060  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2061  }
2062  else
2063  {
2064  ierr = preconditioner->SetUseTranspose(false);
2065  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2066  }
2067  }
2068 
2069 
2070  inline void
2071  PreconditionBase::vmult(MPI::Vector &dst, const MPI::Vector &src) const
2072  {
2073  Assert(dst.trilinos_partitioner().SameAs(
2074  preconditioner->OperatorRangeMap()),
2075  ExcNonMatchingMaps("dst"));
2076  Assert(src.trilinos_partitioner().SameAs(
2077  preconditioner->OperatorDomainMap()),
2078  ExcNonMatchingMaps("src"));
2079 
2080  const int ierr = preconditioner->ApplyInverse(src.trilinos_vector(),
2081  dst.trilinos_vector());
2082  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2083  }
2084 
2085  inline void
2086  PreconditionBase::Tvmult(MPI::Vector &dst, const MPI::Vector &src) const
2087  {
2088  Assert(dst.trilinos_partitioner().SameAs(
2089  preconditioner->OperatorRangeMap()),
2090  ExcNonMatchingMaps("dst"));
2091  Assert(src.trilinos_partitioner().SameAs(
2092  preconditioner->OperatorDomainMap()),
2093  ExcNonMatchingMaps("src"));
2094 
2095  preconditioner->SetUseTranspose(true);
2096  const int ierr = preconditioner->ApplyInverse(src.trilinos_vector(),
2097  dst.trilinos_vector());
2098  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2099  preconditioner->SetUseTranspose(false);
2100  }
2101 
2102 
2103  // For the implementation of the <code>vmult</code> function with deal.II
2104  // data structures we note that invoking a call of the Trilinos
2105  // preconditioner requires us to use Epetra vectors as well. We do this by
2106  // providing a view, i.e., feed Trilinos with a pointer to the data, so we
2107  // avoid copying the content of the vectors during the iteration (this
2108  // function is only useful when used in serial anyway). In the declaration
2109  // of the right hand side, we need to cast the source vector (that is
2110  // <code>const</code> in all deal.II calls) to non-constant value, as this
2111  // is the way Trilinos wants to have them.
2112  inline void
2113  PreconditionBase::vmult(::Vector<double> & dst,
2114  const ::Vector<double> &src) const
2115  {
2116  AssertDimension(dst.size(),
2117  preconditioner->OperatorDomainMap().NumMyElements());
2118  AssertDimension(src.size(),
2119  preconditioner->OperatorRangeMap().NumMyElements());
2120  Epetra_Vector tril_dst(View,
2121  preconditioner->OperatorDomainMap(),
2122  dst.begin());
2123  Epetra_Vector tril_src(View,
2124  preconditioner->OperatorRangeMap(),
2125  const_cast<double *>(src.begin()));
2126 
2127  const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2128  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2129  }
2130 
2131 
2132  inline void
2133  PreconditionBase::Tvmult(::Vector<double> & dst,
2134  const ::Vector<double> &src) const
2135  {
2136  AssertDimension(dst.size(),
2137  preconditioner->OperatorDomainMap().NumMyElements());
2138  AssertDimension(src.size(),
2139  preconditioner->OperatorRangeMap().NumMyElements());
2140  Epetra_Vector tril_dst(View,
2141  preconditioner->OperatorDomainMap(),
2142  dst.begin());
2143  Epetra_Vector tril_src(View,
2144  preconditioner->OperatorRangeMap(),
2145  const_cast<double *>(src.begin()));
2146 
2147  preconditioner->SetUseTranspose(true);
2148  const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2149  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2150  preconditioner->SetUseTranspose(false);
2151  }
2152 
2153 
2154 
2155  inline void
2156  PreconditionBase::vmult(
2159  {
2161  preconditioner->OperatorDomainMap().NumMyElements());
2163  preconditioner->OperatorRangeMap().NumMyElements());
2164  Epetra_Vector tril_dst(View,
2165  preconditioner->OperatorDomainMap(),
2166  dst.begin());
2167  Epetra_Vector tril_src(View,
2168  preconditioner->OperatorRangeMap(),
2169  const_cast<double *>(src.begin()));
2170 
2171  const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2172  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2173  }
2174 
2175  inline void
2176  PreconditionBase::Tvmult(
2179  {
2181  preconditioner->OperatorDomainMap().NumMyElements());
2183  preconditioner->OperatorRangeMap().NumMyElements());
2184  Epetra_Vector tril_dst(View,
2185  preconditioner->OperatorDomainMap(),
2186  dst.begin());
2187  Epetra_Vector tril_src(View,
2188  preconditioner->OperatorRangeMap(),
2189  const_cast<double *>(src.begin()));
2190 
2191  preconditioner->SetUseTranspose(true);
2192  const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2193  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2194  preconditioner->SetUseTranspose(false);
2195  }
2196 
2197 # endif
2198 
2199 } // namespace TrilinosWrappers
2200 
2201 
2205 DEAL_II_NAMESPACE_CLOSE
2206 
2207 # endif // DEAL_II_WITH_TRILINOS
2208 
2209 #endif // trilinos_precondition_h
2210 /*------------------------- 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
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:1407
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