Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
trilinos_precondition.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2008 - 2023 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
40class Ifpack_Preconditioner;
41class Ifpack_Chebyshev;
42namespace ML_Epetra
43{
44 class MultiLevelPreconditioner;
45}
46# endif
47
49
50// forward declarations
51# ifndef DOXYGEN
52template <typename number>
53class SparseMatrix;
54template <typename number>
55class Vector;
56class SparsityPattern;
57# endif
58
64namespace TrilinosWrappers
65{
66 // forward declarations
67 class SparseMatrix;
69 class SolverBase;
70
78 {
79 public:
84
90 {};
91
98
103
107 ~PreconditionBase() override = default;
108
113 void
114 clear();
115
120 get_mpi_communicator() const;
121
131 void
133
137 virtual void
138 vmult(MPI::Vector &dst, const MPI::Vector &src) const;
139
143 virtual void
144 Tvmult(MPI::Vector &dst, const MPI::Vector &src) const;
145
150 virtual void
151 vmult(::Vector<double> &dst, const ::Vector<double> &src) const;
152
157 virtual void
159 const ::Vector<double> &src) const;
160
165 virtual void
167 const ::LinearAlgebra::distributed::Vector<double> &src) const;
168
173 virtual void
175 const ::LinearAlgebra::distributed::Vector<double> &src) const;
176
187 trilinos_operator() const;
201
209
220 std::string,
221 << "The sparse matrix the preconditioner is based on "
222 << "uses a map that is not compatible to the one in vector "
223 << arg1 << ". Check preconditioner and matrix setup.");
226 friend class SolverBase;
227
228 protected:
233 Teuchos::RCP<Epetra_Operator> preconditioner;
234
239 Epetra_MpiComm communicator;
240
245 std::shared_ptr<Epetra_Map> vector_distributor;
246 };
247
248
265 {
266 public:
279 {
284 AdditionalData(const double omega = 1,
285 const double min_diagonal = 0,
286 const unsigned int n_sweeps = 1);
287
291 double omega;
292
301
306 unsigned int n_sweeps;
307 };
308
313 void
314 initialize(const SparseMatrix & matrix,
315 const AdditionalData &additional_data = AdditionalData());
316 };
317
318
319
346 {
347 public:
362 {
369 AdditionalData(const double omega = 1,
370 const double min_diagonal = 0,
371 const unsigned int overlap = 0,
372 const unsigned int n_sweeps = 1);
373
378 double omega;
379
388
393 unsigned int overlap;
394
399 unsigned int n_sweeps;
400 };
401
407 void
408 initialize(const SparseMatrix & matrix,
409 const AdditionalData &additional_data = AdditionalData());
410 };
411
412
413
440 {
441 public:
456 {
463 AdditionalData(const double omega = 1,
464 const double min_diagonal = 0,
465 const unsigned int overlap = 0,
466 const unsigned int n_sweeps = 1);
467
472 double omega;
473
482
487 unsigned int overlap;
488
493 unsigned int n_sweeps;
494 };
495
501 void
502 initialize(const SparseMatrix & matrix,
503 const AdditionalData &additional_data = AdditionalData());
504 };
505
506
507
525 {
526 public:
543 {
549 AdditionalData(const unsigned int block_size = 1,
550 const std::string &block_creation_type = "linear",
551 const double omega = 1,
552 const double min_diagonal = 0,
553 const unsigned int n_sweeps = 1);
554
558 unsigned int block_size;
559
568
573 double omega;
574
583
588 unsigned int n_sweeps;
589 };
590
595 void
596 initialize(const SparseMatrix & matrix,
597 const AdditionalData &additional_data = AdditionalData());
598 };
599
600
601
624 {
625 public:
644 {
652 AdditionalData(const unsigned int block_size = 1,
653 const std::string &block_creation_type = "linear",
654 const double omega = 1,
655 const double min_diagonal = 0,
656 const unsigned int overlap = 0,
657 const unsigned int n_sweeps = 1);
658
662 unsigned int block_size;
663
672
677 double omega;
678
687
692 unsigned int overlap;
693
698 unsigned int n_sweeps;
699 };
700
706 void
707 initialize(const SparseMatrix & matrix,
708 const AdditionalData &additional_data = AdditionalData());
709 };
710
711
712
735 {
736 public:
755 {
763 AdditionalData(const unsigned int block_size = 1,
764 const std::string &block_creation_type = "linear",
765 const double omega = 1,
766 const double min_diagonal = 0,
767 const unsigned int overlap = 0,
768 const unsigned int n_sweeps = 1);
769
773 unsigned int block_size;
774
783
788 double omega;
789
798
803 unsigned int overlap;
804
809 unsigned int n_sweeps;
810 };
811
817 void
818 initialize(const SparseMatrix & matrix,
819 const AdditionalData &additional_data = AdditionalData());
820 };
821
822
823
859 {
860 public:
878 {
888 AdditionalData(const unsigned int ic_fill = 0,
889 const double ic_atol = 0.,
890 const double ic_rtol = 1.,
891 const unsigned int overlap = 0);
892
901 unsigned int ic_fill;
902
908 double ic_atol;
909
914 double ic_rtol;
915
920 unsigned int overlap;
921 };
922
927 void
928 initialize(const SparseMatrix & matrix,
929 const AdditionalData &additional_data = AdditionalData());
930 };
931
932
933
963 {
964 public:
1003 {
1007 AdditionalData(const unsigned int ilu_fill = 0,
1008 const double ilu_atol = 0.,
1009 const double ilu_rtol = 1.,
1010 const unsigned int overlap = 0);
1011
1015 unsigned int ilu_fill;
1016
1021 double ilu_atol;
1022
1027 double ilu_rtol;
1028
1032 unsigned int overlap;
1033 };
1034
1039 void
1040 initialize(const SparseMatrix & matrix,
1041 const AdditionalData &additional_data = AdditionalData());
1042 };
1043
1044
1045
1081 {
1082 public:
1101 {
1112 AdditionalData(const double ilut_drop = 0.,
1113 const unsigned int ilut_fill = 0,
1114 const double ilut_atol = 0.,
1115 const double ilut_rtol = 1.,
1116 const unsigned int overlap = 0);
1117
1123
1132 unsigned int ilut_fill;
1133
1140
1146
1151 unsigned int overlap;
1152 };
1153
1158 void
1159 initialize(const SparseMatrix & matrix,
1160 const AdditionalData &additional_data = AdditionalData());
1161 };
1162
1163
1164
1183 {
1184 public:
1190 {
1194 AdditionalData(const unsigned int overlap = 0);
1195
1196
1201 unsigned int overlap;
1202 };
1203
1208 void
1209 initialize(const SparseMatrix & matrix,
1210 const AdditionalData &additional_data = AdditionalData());
1211 };
1212
1213
1214
1224 {
1225 public:
1231 {
1235 AdditionalData(const unsigned int degree = 1,
1236 const double max_eigenvalue = 10.,
1237 const double eigenvalue_ratio = 30.,
1238 const double min_eigenvalue = 1.,
1239 const double min_diagonal = 1e-12,
1240 const bool nonzero_starting = false);
1241
1247 unsigned int degree;
1248
1254
1259
1265
1271
1283 };
1284
1289 void
1290 initialize(const SparseMatrix & matrix,
1291 const AdditionalData &additional_data = AdditionalData());
1292 };
1293
1294
1295
1338 {
1339 public:
1347 {
1371 AdditionalData(const bool elliptic = true,
1372 const bool higher_order_elements = false,
1373 const unsigned int n_cycles = 1,
1374 const bool w_cycle = false,
1375 const double aggregation_threshold = 1e-4,
1376 const std::vector<std::vector<bool>> &constant_modes =
1377 std::vector<std::vector<bool>>(0),
1378 const unsigned int smoother_sweeps = 2,
1379 const unsigned int smoother_overlap = 0,
1380 const bool output_details = false,
1381 const char * smoother_type = "Chebyshev",
1382 const char * coarse_type = "Amesos-KLU");
1383
1414 void
1416 Teuchos::ParameterList & parameter_list,
1417 std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1418 const Epetra_RowMatrix & matrix) const;
1419
1427 void
1429 Teuchos::ParameterList & parameter_list,
1430 std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1431 const SparseMatrix & matrix) const;
1432
1438 void
1440 Teuchos::ParameterList & parameter_list,
1441 std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1442 const Epetra_RowMatrix & matrix) const;
1443
1449 void
1451 Teuchos::ParameterList & parameter_list,
1452 std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1453 const SparseMatrix & matrix) const;
1454
1463
1469
1474 unsigned int n_cycles;
1475
1481
1492
1517 std::vector<std::vector<bool>> constant_modes;
1518
1529 unsigned int smoother_sweeps;
1530
1535 unsigned int smoother_overlap;
1536
1543
1578 const char *smoother_type;
1579
1584 const char *coarse_type;
1585 };
1586
1590 ~PreconditionAMG() override;
1591
1592
1598 void
1599 initialize(const SparseMatrix & matrix,
1600 const AdditionalData &additional_data = AdditionalData());
1601
1620 void
1621 initialize(const Epetra_RowMatrix &matrix,
1622 const AdditionalData & additional_data = AdditionalData());
1623
1638 void
1639 initialize(const SparseMatrix & matrix,
1640 const Teuchos::ParameterList &ml_parameters);
1641
1649 void
1650 initialize(const Epetra_RowMatrix & matrix,
1651 const Teuchos::ParameterList &ml_parameters);
1652
1659 template <typename number>
1660 void
1661 initialize(const ::SparseMatrix<number> &deal_ii_sparse_matrix,
1662 const AdditionalData &additional_data = AdditionalData(),
1663 const double drop_tolerance = 1e-13,
1664 const ::SparsityPattern *use_this_sparsity = nullptr);
1665
1678 void
1679 reinit();
1680
1685 void
1686 clear();
1687
1691 size_type
1692 memory_consumption() const;
1693
1694 private:
1698 std::shared_ptr<SparseMatrix> trilinos_matrix;
1699 };
1700
1701
1702
1703# if defined(DOXYGEN) || defined(DEAL_II_TRILINOS_WITH_MUELU)
1723 {
1724 public:
1732 {
1737 AdditionalData(const bool elliptic = true,
1738 const unsigned int n_cycles = 1,
1739 const bool w_cycle = false,
1740 const double aggregation_threshold = 1e-4,
1741 const std::vector<std::vector<bool>> &constant_modes =
1742 std::vector<std::vector<bool>>(0),
1743 const unsigned int smoother_sweeps = 2,
1744 const unsigned int smoother_overlap = 0,
1745 const bool output_details = false,
1746 const char * smoother_type = "Chebyshev",
1747 const char * coarse_type = "Amesos-KLU");
1748
1757
1762 unsigned int n_cycles;
1763
1769
1780
1787 std::vector<std::vector<bool>> constant_modes;
1788
1799 unsigned int smoother_sweeps;
1800
1805 unsigned int smoother_overlap;
1806
1813
1848 const char *smoother_type;
1849
1854 const char *coarse_type;
1855 };
1856
1861
1865 virtual ~PreconditionAMGMueLu() override = default;
1866
1872 void
1873 initialize(const SparseMatrix & matrix,
1874 const AdditionalData &additional_data = AdditionalData());
1875
1882 void
1883 initialize(const Epetra_CrsMatrix &matrix,
1884 const AdditionalData & additional_data = AdditionalData());
1885
1899 void
1900 initialize(const SparseMatrix & matrix,
1901 Teuchos::ParameterList &muelu_parameters);
1902
1908 void
1909 initialize(const Epetra_CrsMatrix &matrix,
1910 Teuchos::ParameterList &muelu_parameters);
1911
1918 template <typename number>
1919 void
1920 initialize(const ::SparseMatrix<number> &deal_ii_sparse_matrix,
1921 const AdditionalData &additional_data = AdditionalData(),
1922 const double drop_tolerance = 1e-13,
1923 const ::SparsityPattern *use_this_sparsity = nullptr);
1924
1929 void
1930 clear();
1931
1935 size_type
1936 memory_consumption() const;
1937
1938 private:
1942 std::shared_ptr<SparseMatrix> trilinos_matrix;
1943 };
1944# endif
1945
1946
1947
1955 {
1956 public:
1962 {};
1963
1970 void
1971 initialize(const SparseMatrix & matrix,
1972 const AdditionalData &additional_data = AdditionalData());
1973
1977 void
1978 vmult(MPI::Vector &dst, const MPI::Vector &src) const override;
1979
1983 void
1984 Tvmult(MPI::Vector &dst, const MPI::Vector &src) const override;
1985
1990 void
1991 vmult(::Vector<double> & dst,
1992 const ::Vector<double> &src) const override;
1993
1998 void
2000 const ::Vector<double> &src) const override;
2001
2006 void
2008 const ::LinearAlgebra::distributed::Vector<double> &src)
2009 const override;
2010
2016 void
2018 const ::LinearAlgebra::distributed::Vector<double> &src)
2019 const override;
2020 };
2021
2022
2023
2024 // ----------------------- inline and template functions --------------------
2025
2026
2027# ifndef DOXYGEN
2028
2029
2030 inline void
2032 {
2033 // This only flips a flag that tells
2034 // Trilinos that any vmult operation
2035 // should be done with the
2036 // transpose. However, the matrix
2037 // structure is not reset.
2038 int ierr;
2039
2040 if (!preconditioner->UseTranspose())
2041 {
2042 ierr = preconditioner->SetUseTranspose(true);
2043 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2044 }
2045 else
2046 {
2047 ierr = preconditioner->SetUseTranspose(false);
2048 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2049 }
2050 }
2051
2052
2053 inline void
2054 PreconditionBase::vmult(MPI::Vector &dst, const MPI::Vector &src) const
2055 {
2056 Assert(dst.trilinos_partitioner().SameAs(
2057 preconditioner->OperatorRangeMap()),
2058 ExcNonMatchingMaps("dst"));
2059 Assert(src.trilinos_partitioner().SameAs(
2060 preconditioner->OperatorDomainMap()),
2061 ExcNonMatchingMaps("src"));
2062
2063 const int ierr = preconditioner->ApplyInverse(src.trilinos_vector(),
2064 dst.trilinos_vector());
2065 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2066 }
2067
2068 inline void
2069 PreconditionBase::Tvmult(MPI::Vector &dst, const MPI::Vector &src) const
2070 {
2071 Assert(dst.trilinos_partitioner().SameAs(
2072 preconditioner->OperatorRangeMap()),
2073 ExcNonMatchingMaps("dst"));
2074 Assert(src.trilinos_partitioner().SameAs(
2075 preconditioner->OperatorDomainMap()),
2076 ExcNonMatchingMaps("src"));
2077
2078 preconditioner->SetUseTranspose(true);
2079 const int ierr = preconditioner->ApplyInverse(src.trilinos_vector(),
2080 dst.trilinos_vector());
2081 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2082 preconditioner->SetUseTranspose(false);
2083 }
2084
2085
2086 // For the implementation of the <code>vmult</code> function with deal.II
2087 // data structures we note that invoking a call of the Trilinos
2088 // preconditioner requires us to use Epetra vectors as well. We do this by
2089 // providing a view, i.e., feed Trilinos with a pointer to the data, so we
2090 // avoid copying the content of the vectors during the iteration (this
2091 // function is only useful when used in serial anyway). In the declaration
2092 // of the right hand side, we need to cast the source vector (that is
2093 // <code>const</code> in all deal.II calls) to non-constant value, as this
2094 // is the way Trilinos wants to have them.
2095 inline void
2097 const ::Vector<double> &src) const
2098 {
2099 AssertDimension(dst.size(),
2100 preconditioner->OperatorDomainMap().NumMyElements());
2101 AssertDimension(src.size(),
2102 preconditioner->OperatorRangeMap().NumMyElements());
2103 Epetra_Vector tril_dst(View,
2104 preconditioner->OperatorDomainMap(),
2105 dst.begin());
2106 Epetra_Vector tril_src(View,
2107 preconditioner->OperatorRangeMap(),
2108 const_cast<double *>(src.begin()));
2109
2110 const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2111 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2112 }
2113
2114
2115 inline void
2117 const ::Vector<double> &src) const
2118 {
2119 AssertDimension(dst.size(),
2120 preconditioner->OperatorDomainMap().NumMyElements());
2121 AssertDimension(src.size(),
2122 preconditioner->OperatorRangeMap().NumMyElements());
2123 Epetra_Vector tril_dst(View,
2124 preconditioner->OperatorDomainMap(),
2125 dst.begin());
2126 Epetra_Vector tril_src(View,
2127 preconditioner->OperatorRangeMap(),
2128 const_cast<double *>(src.begin()));
2129
2130 preconditioner->SetUseTranspose(true);
2131 const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2132 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2133 preconditioner->SetUseTranspose(false);
2134 }
2135
2136
2137
2138 inline void
2142 {
2144 preconditioner->OperatorDomainMap().NumMyElements());
2146 preconditioner->OperatorRangeMap().NumMyElements());
2147 Epetra_Vector tril_dst(View,
2148 preconditioner->OperatorDomainMap(),
2149 dst.begin());
2150 Epetra_Vector tril_src(View,
2151 preconditioner->OperatorRangeMap(),
2152 const_cast<double *>(src.begin()));
2153
2154 const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2155 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2156 }
2157
2158 inline void
2162 {
2164 preconditioner->OperatorDomainMap().NumMyElements());
2166 preconditioner->OperatorRangeMap().NumMyElements());
2167 Epetra_Vector tril_dst(View,
2168 preconditioner->OperatorDomainMap(),
2169 dst.begin());
2170 Epetra_Vector tril_src(View,
2171 preconditioner->OperatorRangeMap(),
2172 const_cast<double *>(src.begin()));
2173
2174 preconditioner->SetUseTranspose(true);
2175 const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2176 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2177 preconditioner->SetUseTranspose(false);
2178 }
2179
2180# endif
2181
2182} // namespace TrilinosWrappers
2183
2184
2189
2190#endif // DEAL_II_WITH_TRILINOS
2191
2192#endif
size_type locally_owned_size() const
std::shared_ptr< SparseMatrix > trilinos_matrix
virtual ~PreconditionAMGMueLu() override=default
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
std::shared_ptr< SparseMatrix > trilinos_matrix
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
virtual void Tvmult(::LinearAlgebra::distributed::Vector< double > &dst, const ::LinearAlgebra::distributed::Vector< double > &src) const
virtual void Tvmult(MPI::Vector &dst, const MPI::Vector &src) const
~PreconditionBase() override=default
virtual void vmult(MPI::Vector &dst, const MPI::Vector &src) const
Epetra_Operator & trilinos_operator() const
std::shared_ptr< Epetra_Map > vector_distributor
Teuchos::RCP< Epetra_Operator > preconditioner
virtual void Tvmult(::Vector< double > &dst, const ::Vector< double > &src) const
virtual void vmult(::Vector< double > &dst, const ::Vector< double > &src) const
virtual void vmult(::LinearAlgebra::distributed::Vector< double > &dst, const ::LinearAlgebra::distributed::Vector< double > &src) const
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
void vmult(LinearAlgebra::distributed::Vector< double > &dst, const ::LinearAlgebra::distributed::Vector< double > &src) const override
void vmult(MPI::Vector &dst, const MPI::Vector &src) const override
void Tvmult(MPI::Vector &dst, const MPI::Vector &src) const override
void Tvmult(LinearAlgebra::distributed::Vector< double > &dst, const ::LinearAlgebra::distributed::Vector< double > &src) const override
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
size_type size() const
iterator begin()
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNonMatchingMaps(std::string arg1)
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:510
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define AssertThrow(cond, exc)
unsigned int global_dof_index
Definition types.h:82
void set_operator_null_space(Teuchos::ParameterList &parameter_list, std::unique_ptr< Epetra_MultiVector > &distributed_constant_modes, const Epetra_RowMatrix &matrix) const
void set_parameters(Teuchos::ParameterList &parameter_list, std::unique_ptr< Epetra_MultiVector > &distributed_constant_modes, const Epetra_RowMatrix &matrix) const