Reference documentation for deal.II version Git 01d13113a1 2020-07-10 09:14:01 +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:
552  {
558  AdditionalData(const unsigned int block_size = 1,
559  const std::string &block_creation_type = "linear",
560  const double omega = 1,
561  const double min_diagonal = 0,
562  const unsigned int n_sweeps = 1);
563 
567  unsigned int block_size;
568 
576  std::string block_creation_type;
577 
582  double omega;
583 
591  double min_diagonal;
592 
597  unsigned int n_sweeps;
598  };
599 
604  void
605  initialize(const SparseMatrix & matrix,
606  const AdditionalData &additional_data = AdditionalData());
607  };
608 
609 
610 
633  {
634  public:
653  {
661  AdditionalData(const unsigned int block_size = 1,
662  const std::string &block_creation_type = "linear",
663  const double omega = 1,
664  const double min_diagonal = 0,
665  const unsigned int overlap = 0,
666  const unsigned int n_sweeps = 1);
667 
671  unsigned int block_size;
672 
680  std::string block_creation_type;
681 
686  double omega;
687 
695  double min_diagonal;
696 
701  unsigned int overlap;
702 
707  unsigned int n_sweeps;
708  };
709 
715  void
716  initialize(const SparseMatrix & matrix,
717  const AdditionalData &additional_data = AdditionalData());
718  };
719 
720 
721 
744  {
745  public:
764  {
772  AdditionalData(const unsigned int block_size = 1,
773  const std::string &block_creation_type = "linear",
774  const double omega = 1,
775  const double min_diagonal = 0,
776  const unsigned int overlap = 0,
777  const unsigned int n_sweeps = 1);
778 
782  unsigned int block_size;
783 
791  std::string block_creation_type;
792 
797  double omega;
798 
806  double min_diagonal;
807 
812  unsigned int overlap;
813 
818  unsigned int n_sweeps;
819  };
820 
826  void
827  initialize(const SparseMatrix & matrix,
828  const AdditionalData &additional_data = AdditionalData());
829  };
830 
831 
832 
868  {
869  public:
887  {
897  AdditionalData(const unsigned int ic_fill = 0,
898  const double ic_atol = 0.,
899  const double ic_rtol = 1.,
900  const unsigned int overlap = 0);
901 
910  unsigned int ic_fill;
911 
917  double ic_atol;
918 
923  double ic_rtol;
924 
929  unsigned int overlap;
930  };
931 
936  void
937  initialize(const SparseMatrix & matrix,
938  const AdditionalData &additional_data = AdditionalData());
939  };
940 
941 
942 
972  {
973  public:
1012  {
1016  AdditionalData(const unsigned int ilu_fill = 0,
1017  const double ilu_atol = 0.,
1018  const double ilu_rtol = 1.,
1019  const unsigned int overlap = 0);
1020 
1024  unsigned int ilu_fill;
1025 
1030  double ilu_atol;
1031 
1036  double ilu_rtol;
1037 
1041  unsigned int overlap;
1042  };
1043 
1048  void
1049  initialize(const SparseMatrix & matrix,
1050  const AdditionalData &additional_data = AdditionalData());
1051  };
1052 
1053 
1054 
1090  {
1091  public:
1110  {
1121  AdditionalData(const double ilut_drop = 0.,
1122  const unsigned int ilut_fill = 0,
1123  const double ilut_atol = 0.,
1124  const double ilut_rtol = 1.,
1125  const unsigned int overlap = 0);
1126 
1131  double ilut_drop;
1132 
1141  unsigned int ilut_fill;
1142 
1148  double ilut_atol;
1149 
1154  double ilut_rtol;
1155 
1160  unsigned int overlap;
1161  };
1162 
1167  void
1168  initialize(const SparseMatrix & matrix,
1169  const AdditionalData &additional_data = AdditionalData());
1170  };
1171 
1172 
1173 
1192  {
1193  public:
1199  {
1203  AdditionalData(const unsigned int overlap = 0);
1204 
1205 
1210  unsigned int overlap;
1211  };
1212 
1217  void
1218  initialize(const SparseMatrix & matrix,
1219  const AdditionalData &additional_data = AdditionalData());
1220  };
1221 
1222 
1223 
1233  {
1234  public:
1240  {
1244  AdditionalData(const unsigned int degree = 1,
1245  const double max_eigenvalue = 10.,
1246  const double eigenvalue_ratio = 30.,
1247  const double min_eigenvalue = 1.,
1248  const double min_diagonal = 1e-12,
1249  const bool nonzero_starting = false);
1250 
1256  unsigned int degree;
1257 
1263 
1268 
1274 
1280 
1292  };
1293 
1298  void
1299  initialize(const SparseMatrix & matrix,
1300  const AdditionalData &additional_data = AdditionalData());
1301  };
1302 
1303 
1304 
1347  {
1348  public:
1356  {
1380  AdditionalData(const bool elliptic = true,
1381  const bool higher_order_elements = false,
1382  const unsigned int n_cycles = 1,
1383  const bool w_cyle = false,
1384  const double aggregation_threshold = 1e-4,
1385  const std::vector<std::vector<bool>> &constant_modes =
1386  std::vector<std::vector<bool>>(0),
1387  const unsigned int smoother_sweeps = 2,
1388  const unsigned int smoother_overlap = 0,
1389  const bool output_details = false,
1390  const char * smoother_type = "Chebyshev",
1391  const char * coarse_type = "Amesos-KLU");
1392 
1423  void
1424  set_parameters(
1425  Teuchos::ParameterList & parameter_list,
1426  std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1427  const Epetra_RowMatrix & matrix) const;
1428 
1436  void
1437  set_parameters(
1438  Teuchos::ParameterList & parameter_list,
1439  std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1440  const SparseMatrix & matrix) const;
1441 
1447  void
1448  set_operator_null_space(
1449  Teuchos::ParameterList & parameter_list,
1450  std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1451  const Epetra_RowMatrix & matrix) const;
1452 
1458  void
1459  set_operator_null_space(
1460  Teuchos::ParameterList & parameter_list,
1461  std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1462  const SparseMatrix & matrix) const;
1463 
1471  bool elliptic;
1472 
1478 
1483  unsigned int n_cycles;
1484 
1489  bool w_cycle;
1490 
1501 
1518  std::vector<std::vector<bool>> constant_modes;
1519 
1530  unsigned int smoother_sweeps;
1531 
1536  unsigned int smoother_overlap;
1537 
1544 
1579  const char *smoother_type;
1580 
1585  const char *coarse_type;
1586  };
1587 
1591  ~PreconditionAMG() override;
1592 
1593 
1599  void
1600  initialize(const SparseMatrix & matrix,
1601  const AdditionalData &additional_data = AdditionalData());
1602 
1621  void
1622  initialize(const Epetra_RowMatrix &matrix,
1623  const AdditionalData & additional_data = AdditionalData());
1624 
1639  void
1640  initialize(const SparseMatrix & matrix,
1641  const Teuchos::ParameterList &ml_parameters);
1642 
1650  void
1651  initialize(const Epetra_RowMatrix & matrix,
1652  const Teuchos::ParameterList &ml_parameters);
1653 
1660  template <typename number>
1661  void
1662  initialize(const ::SparseMatrix<number> &deal_ii_sparse_matrix,
1663  const AdditionalData &additional_data = AdditionalData(),
1664  const double drop_tolerance = 1e-13,
1665  const ::SparsityPattern *use_this_sparsity = nullptr);
1666 
1679  void
1680  reinit();
1681 
1686  void
1687  clear();
1688 
1692  size_type
1693  memory_consumption() const;
1694 
1695  private:
1699  std::shared_ptr<SparseMatrix> trilinos_matrix;
1700  };
1701 
1702 
1703 
1704 # if defined(DOXYGEN) || defined(DEAL_II_TRILINOS_WITH_MUELU)
1705 
1724  {
1725  public:
1733  {
1738  AdditionalData(const bool elliptic = true,
1739  const unsigned int n_cycles = 1,
1740  const bool w_cyle = false,
1741  const double aggregation_threshold = 1e-4,
1742  const std::vector<std::vector<bool>> &constant_modes =
1743  std::vector<std::vector<bool>>(0),
1744  const unsigned int smoother_sweeps = 2,
1745  const unsigned int smoother_overlap = 0,
1746  const bool output_details = false,
1747  const char * smoother_type = "Chebyshev",
1748  const char * coarse_type = "Amesos-KLU");
1749 
1757  bool elliptic;
1758 
1763  unsigned int n_cycles;
1764 
1769  bool w_cycle;
1770 
1781 
1788  std::vector<std::vector<bool>> constant_modes;
1789 
1800  unsigned int smoother_sweeps;
1801 
1806  unsigned int smoother_overlap;
1807 
1814 
1849  const char *smoother_type;
1850 
1855  const char *coarse_type;
1856  };
1857 
1862 
1866  virtual ~PreconditionAMGMueLu() override = default;
1867 
1873  void
1874  initialize(const SparseMatrix & matrix,
1875  const AdditionalData &additional_data = AdditionalData());
1876 
1883  void
1884  initialize(const Epetra_CrsMatrix &matrix,
1885  const AdditionalData & additional_data = AdditionalData());
1886 
1900  void
1901  initialize(const SparseMatrix & matrix,
1902  Teuchos::ParameterList &muelu_parameters);
1903 
1909  void
1910  initialize(const Epetra_CrsMatrix &matrix,
1911  Teuchos::ParameterList &muelu_parameters);
1912 
1919  template <typename number>
1920  void
1921  initialize(const ::SparseMatrix<number> &deal_ii_sparse_matrix,
1922  const AdditionalData &additional_data = AdditionalData(),
1923  const double drop_tolerance = 1e-13,
1924  const ::SparsityPattern *use_this_sparsity = nullptr);
1925 
1930  void
1931  clear();
1932 
1936  size_type
1937  memory_consumption() const;
1938 
1939  private:
1943  std::shared_ptr<SparseMatrix> trilinos_matrix;
1944  };
1945 # endif
1946 
1947 
1948 
1956  {
1957  public:
1963  {};
1964 
1971  void
1972  initialize(const SparseMatrix & matrix,
1973  const AdditionalData &additional_data = AdditionalData());
1974 
1978  void
1979  vmult(MPI::Vector &dst, const MPI::Vector &src) const override;
1980 
1984  void
1985  Tvmult(MPI::Vector &dst, const MPI::Vector &src) const override;
1986 
1991  void
1992  vmult(::Vector<double> & dst,
1993  const ::Vector<double> &src) const override;
1994 
1999  void
2000  Tvmult(::Vector<double> & dst,
2001  const ::Vector<double> &src) const override;
2002 
2007  void
2009  const ::LinearAlgebra::distributed::Vector<double> &src)
2010  const override;
2011 
2017  void
2019  const ::LinearAlgebra::distributed::Vector<double> &src)
2020  const override;
2021  };
2022 
2023 
2024 
2025  // ----------------------- inline and template functions --------------------
2026 
2027 
2028 # ifndef DOXYGEN
2029 
2030 
2031  inline void
2033  {
2034  // This only flips a flag that tells
2035  // Trilinos that any vmult operation
2036  // should be done with the
2037  // transpose. However, the matrix
2038  // structure is not reset.
2039  int ierr;
2040 
2041  if (!preconditioner->UseTranspose())
2042  {
2043  ierr = preconditioner->SetUseTranspose(true);
2044  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2045  }
2046  else
2047  {
2048  ierr = preconditioner->SetUseTranspose(false);
2049  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2050  }
2051  }
2052 
2053 
2054  inline void
2055  PreconditionBase::vmult(MPI::Vector &dst, const MPI::Vector &src) const
2056  {
2057  Assert(dst.trilinos_partitioner().SameAs(
2058  preconditioner->OperatorRangeMap()),
2059  ExcNonMatchingMaps("dst"));
2060  Assert(src.trilinos_partitioner().SameAs(
2061  preconditioner->OperatorDomainMap()),
2062  ExcNonMatchingMaps("src"));
2063 
2064  const int ierr = preconditioner->ApplyInverse(src.trilinos_vector(),
2065  dst.trilinos_vector());
2066  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2067  }
2068 
2069  inline void
2070  PreconditionBase::Tvmult(MPI::Vector &dst, const MPI::Vector &src) const
2071  {
2072  Assert(dst.trilinos_partitioner().SameAs(
2073  preconditioner->OperatorRangeMap()),
2074  ExcNonMatchingMaps("dst"));
2075  Assert(src.trilinos_partitioner().SameAs(
2076  preconditioner->OperatorDomainMap()),
2077  ExcNonMatchingMaps("src"));
2078 
2079  preconditioner->SetUseTranspose(true);
2080  const int ierr = preconditioner->ApplyInverse(src.trilinos_vector(),
2081  dst.trilinos_vector());
2082  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2083  preconditioner->SetUseTranspose(false);
2084  }
2085 
2086 
2087  // For the implementation of the <code>vmult</code> function with deal.II
2088  // data structures we note that invoking a call of the Trilinos
2089  // preconditioner requires us to use Epetra vectors as well. We do this by
2090  // providing a view, i.e., feed Trilinos with a pointer to the data, so we
2091  // avoid copying the content of the vectors during the iteration (this
2092  // function is only useful when used in serial anyway). In the declaration
2093  // of the right hand side, we need to cast the source vector (that is
2094  // <code>const</code> in all deal.II calls) to non-constant value, as this
2095  // is the way Trilinos wants to have them.
2096  inline void
2097  PreconditionBase::vmult(::Vector<double> & dst,
2098  const ::Vector<double> &src) const
2099  {
2100  AssertDimension(dst.size(),
2101  preconditioner->OperatorDomainMap().NumMyElements());
2102  AssertDimension(src.size(),
2103  preconditioner->OperatorRangeMap().NumMyElements());
2104  Epetra_Vector tril_dst(View,
2105  preconditioner->OperatorDomainMap(),
2106  dst.begin());
2107  Epetra_Vector tril_src(View,
2108  preconditioner->OperatorRangeMap(),
2109  const_cast<double *>(src.begin()));
2110 
2111  const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2112  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2113  }
2114 
2115 
2116  inline void
2117  PreconditionBase::Tvmult(::Vector<double> & dst,
2118  const ::Vector<double> &src) const
2119  {
2120  AssertDimension(dst.size(),
2121  preconditioner->OperatorDomainMap().NumMyElements());
2122  AssertDimension(src.size(),
2123  preconditioner->OperatorRangeMap().NumMyElements());
2124  Epetra_Vector tril_dst(View,
2125  preconditioner->OperatorDomainMap(),
2126  dst.begin());
2127  Epetra_Vector tril_src(View,
2128  preconditioner->OperatorRangeMap(),
2129  const_cast<double *>(src.begin()));
2130 
2131  preconditioner->SetUseTranspose(true);
2132  const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2133  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2134  preconditioner->SetUseTranspose(false);
2135  }
2136 
2137 
2138 
2139  inline void
2140  PreconditionBase::vmult(
2143  {
2145  preconditioner->OperatorDomainMap().NumMyElements());
2147  preconditioner->OperatorRangeMap().NumMyElements());
2148  Epetra_Vector tril_dst(View,
2149  preconditioner->OperatorDomainMap(),
2150  dst.begin());
2151  Epetra_Vector tril_src(View,
2152  preconditioner->OperatorRangeMap(),
2153  const_cast<double *>(src.begin()));
2154 
2155  const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2156  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2157  }
2158 
2159  inline void
2160  PreconditionBase::Tvmult(
2163  {
2165  preconditioner->OperatorDomainMap().NumMyElements());
2167  preconditioner->OperatorRangeMap().NumMyElements());
2168  Epetra_Vector tril_dst(View,
2169  preconditioner->OperatorDomainMap(),
2170  dst.begin());
2171  Epetra_Vector tril_src(View,
2172  preconditioner->OperatorRangeMap(),
2173  const_cast<double *>(src.begin()));
2174 
2175  preconditioner->SetUseTranspose(true);
2176  const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2177  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2178  preconditioner->SetUseTranspose(false);
2179  }
2180 
2181 # endif
2182 
2183 } // namespace TrilinosWrappers
2184 
2185 
2190 
2191 # endif // DEAL_II_WITH_TRILINOS
2192 
2193 #endif // trilinos_precondition_h
2194 /*------------------------- 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:1560
typename ::internal::FEValuesViews::ViewType< dim, spacedim, Extractor >::type View
Definition: fe_values.h:1970
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:1513
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:515
#define Assert(cond, exc)
Definition: exceptions.h:1403
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
iterator begin()
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:362
size_type size() const
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