deal.II version GIT relicensing-3509-g790ffced1c 2025-06-14 07:20:00+00:00
\(\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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2008 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_trilinos_precondition_h
16#define dealii_trilinos_precondition_h
17
18
19#include <deal.II/base/config.h>
20
21#ifdef DEAL_II_WITH_TRILINOS
22
24
27
28# include <Epetra_Map.h>
29# include <Epetra_MpiComm.h>
30# include <Epetra_MultiVector.h>
31# include <Epetra_RowMatrix.h>
32# include <Epetra_Vector.h>
33# include <Teuchos_ParameterList.hpp>
34
35# include <memory>
36
37// forward declarations
38# ifndef DOXYGEN
39class Ifpack_Preconditioner;
40class Ifpack_Chebyshev;
41namespace ML_Epetra
42{
43 class MultiLevelPreconditioner;
44}
45# endif
46
48
49// forward declarations
50# ifndef DOXYGEN
51template <typename number>
52class SparseMatrix;
53template <typename number>
54class Vector;
55class SparsityPattern;
56# endif
57
63namespace TrilinosWrappers
64{
65 // forward declarations
66 class SparseMatrix;
68 class SolverBase;
69
77 {
78 public:
83
89 {};
90
97
102 void
103 clear();
104
109 get_mpi_communicator() const;
110
120 void
122
126 virtual void
127 vmult(MPI::Vector &dst, const MPI::Vector &src) const;
128
132 virtual void
133 Tvmult(MPI::Vector &dst, const MPI::Vector &src) const;
134
139 virtual void
140 vmult(::Vector<double> &dst, const ::Vector<double> &src) const;
141
146 virtual void
148 const ::Vector<double> &src) const;
149
154 virtual void
156 const ::LinearAlgebra::distributed::Vector<double> &src) const;
157
162 virtual void
164 const ::LinearAlgebra::distributed::Vector<double> &src) const;
165
176 trilinos_operator() const;
190
198
209 std::string,
210 << "The sparse matrix the preconditioner is based on "
211 << "uses a map that is not compatible to the one in vector "
212 << arg1 << ". Check preconditioner and matrix setup.");
215 friend class SolverBase;
216
217 protected:
222 Teuchos::RCP<Epetra_Operator> preconditioner;
223
228 Epetra_MpiComm communicator;
229 };
230
231
248 {
249 public:
262 {
267 AdditionalData(const double omega = 1,
268 const double min_diagonal = 0,
269 const unsigned int n_sweeps = 1);
270
274 double omega;
275
284
289 unsigned int n_sweeps;
290 };
291
296 void
297 initialize(const SparseMatrix &matrix,
298 const AdditionalData &additional_data = AdditionalData());
299 };
300
301
302
329 {
330 public:
345 {
352 AdditionalData(const double omega = 1,
353 const double min_diagonal = 0,
354 const unsigned int overlap = 0,
355 const unsigned int n_sweeps = 1);
356
361 double omega;
362
371
376 unsigned int overlap;
377
382 unsigned int n_sweeps;
383 };
384
390 void
391 initialize(const SparseMatrix &matrix,
392 const AdditionalData &additional_data = AdditionalData());
393 };
394
395
396
423 {
424 public:
439 {
446 AdditionalData(const double omega = 1,
447 const double min_diagonal = 0,
448 const unsigned int overlap = 0,
449 const unsigned int n_sweeps = 1);
450
455 double omega;
456
465
470 unsigned int overlap;
471
476 unsigned int n_sweeps;
477 };
478
484 void
485 initialize(const SparseMatrix &matrix,
486 const AdditionalData &additional_data = AdditionalData());
487 };
488
489
490
508 {
509 public:
526 {
532 AdditionalData(const unsigned int block_size = 1,
533 const std::string &block_creation_type = "linear",
534 const double omega = 1,
535 const double min_diagonal = 0,
536 const unsigned int n_sweeps = 1);
537
541 unsigned int block_size;
542
551
556 double omega;
557
566
571 unsigned int n_sweeps;
572 };
573
578 void
579 initialize(const SparseMatrix &matrix,
580 const AdditionalData &additional_data = AdditionalData());
581 };
582
583
584
607 {
608 public:
627 {
635 AdditionalData(const unsigned int block_size = 1,
636 const std::string &block_creation_type = "linear",
637 const double omega = 1,
638 const double min_diagonal = 0,
639 const unsigned int overlap = 0,
640 const unsigned int n_sweeps = 1);
641
645 unsigned int block_size;
646
655
660 double omega;
661
670
675 unsigned int overlap;
676
681 unsigned int n_sweeps;
682 };
683
689 void
690 initialize(const SparseMatrix &matrix,
691 const AdditionalData &additional_data = AdditionalData());
692 };
693
694
695
718 {
719 public:
738 {
746 AdditionalData(const unsigned int block_size = 1,
747 const std::string &block_creation_type = "linear",
748 const double omega = 1,
749 const double min_diagonal = 0,
750 const unsigned int overlap = 0,
751 const unsigned int n_sweeps = 1);
752
756 unsigned int block_size;
757
766
771 double omega;
772
781
786 unsigned int overlap;
787
792 unsigned int n_sweeps;
793 };
794
800 void
801 initialize(const SparseMatrix &matrix,
802 const AdditionalData &additional_data = AdditionalData());
803 };
804
805
806
842 {
843 public:
861 {
871 AdditionalData(const unsigned int ic_fill = 0,
872 const double ic_atol = 0.,
873 const double ic_rtol = 1.,
874 const unsigned int overlap = 0);
875
884 unsigned int ic_fill;
885
891 double ic_atol;
892
897 double ic_rtol;
898
903 unsigned int overlap;
904 };
905
910 void
911 initialize(const SparseMatrix &matrix,
912 const AdditionalData &additional_data = AdditionalData());
913 };
914
915
916
946 {
947 public:
986 {
990 AdditionalData(const unsigned int ilu_fill = 0,
991 const double ilu_atol = 0.,
992 const double ilu_rtol = 1.,
993 const unsigned int overlap = 0);
994
998 unsigned int ilu_fill;
999
1004 double ilu_atol;
1005
1010 double ilu_rtol;
1011
1015 unsigned int overlap;
1016 };
1017
1022 void
1023 initialize(const SparseMatrix &matrix,
1024 const AdditionalData &additional_data = AdditionalData());
1025 };
1026
1027
1028
1064 {
1065 public:
1084 {
1095 AdditionalData(const double ilut_drop = 0.,
1096 const unsigned int ilut_fill = 0,
1097 const double ilut_atol = 0.,
1098 const double ilut_rtol = 1.,
1099 const unsigned int overlap = 0);
1100
1106
1115 unsigned int ilut_fill;
1116
1123
1129
1134 unsigned int overlap;
1135 };
1136
1141 void
1142 initialize(const SparseMatrix &matrix,
1143 const AdditionalData &additional_data = AdditionalData());
1144 };
1145
1146
1147
1166 {
1167 public:
1173 {
1177 AdditionalData(const unsigned int overlap = 0);
1178
1179
1184 unsigned int overlap;
1185 };
1186
1191 void
1192 initialize(const SparseMatrix &matrix,
1193 const AdditionalData &additional_data = AdditionalData());
1194 };
1195
1196
1197
1207 {
1208 public:
1214 {
1218 AdditionalData(const unsigned int degree = 1,
1219 const double max_eigenvalue = 10.,
1220 const double eigenvalue_ratio = 30.,
1221 const double min_eigenvalue = 1.,
1222 const double min_diagonal = 1e-12,
1223 const bool nonzero_starting = false);
1224
1230 unsigned int degree;
1231
1237
1242
1248
1254
1266 };
1267
1272 void
1273 initialize(const SparseMatrix &matrix,
1274 const AdditionalData &additional_data = AdditionalData());
1275 };
1276
1277
1278
1321 {
1322 public:
1330 {
1354 AdditionalData(const bool elliptic = true,
1355 const bool higher_order_elements = false,
1356 const unsigned int n_cycles = 1,
1357 const bool w_cycle = false,
1358 const double aggregation_threshold = 1e-4,
1359 const std::vector<std::vector<bool>> &constant_modes =
1360 std::vector<std::vector<bool>>(0),
1361 const unsigned int smoother_sweeps = 2,
1362 const unsigned int smoother_overlap = 0,
1363 const bool output_details = false,
1364 const char *smoother_type = "Chebyshev",
1365 const char *coarse_type = "Amesos-KLU");
1366
1397 void
1399 Teuchos::ParameterList &parameter_list,
1400 std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1401 const Epetra_RowMatrix &matrix) const;
1402
1410 void
1412 Teuchos::ParameterList &parameter_list,
1413 std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1414 const SparseMatrix &matrix) const;
1415
1421 void
1423 Teuchos::ParameterList &parameter_list,
1424 std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1425 const Epetra_RowMatrix &matrix) const;
1426
1432 void
1434 Teuchos::ParameterList &parameter_list,
1435 std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
1436 const SparseMatrix &matrix) const;
1437
1446
1452
1457 unsigned int n_cycles;
1458
1464
1475
1500 std::vector<std::vector<bool>> constant_modes;
1501
1509 std::vector<std::vector<double>> constant_modes_values;
1510
1521 unsigned int smoother_sweeps;
1522
1527 unsigned int smoother_overlap;
1528
1535
1570 const char *smoother_type;
1571
1576 const char *coarse_type;
1577 };
1578
1582 ~PreconditionAMG() override;
1583
1584
1590 void
1591 initialize(const SparseMatrix &matrix,
1592 const AdditionalData &additional_data = AdditionalData());
1593
1612 void
1613 initialize(const Epetra_RowMatrix &matrix,
1614 const AdditionalData &additional_data = AdditionalData());
1615
1630 void
1631 initialize(const SparseMatrix &matrix,
1632 const Teuchos::ParameterList &ml_parameters);
1633
1641 void
1642 initialize(const Epetra_RowMatrix &matrix,
1643 const Teuchos::ParameterList &ml_parameters);
1644
1651 template <typename number>
1652 void
1653 initialize(const ::SparseMatrix<number> &deal_ii_sparse_matrix,
1654 const AdditionalData &additional_data = AdditionalData(),
1655 const double drop_tolerance = 1e-13,
1656 const ::SparsityPattern *use_this_sparsity = nullptr);
1657
1670 void
1671 reinit();
1672
1677 void
1678 clear();
1679
1683 size_type
1684 memory_consumption() const;
1685
1686 private:
1690 std::shared_ptr<SparseMatrix> trilinos_matrix;
1691 };
1692
1693
1694
1695# 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
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#else
2191
2192// Make sure the scripts that create the C++20 module input files have
2193// something to latch on if the preprocessor #ifdef above would
2194// otherwise lead to an empty content of the file.
2197
2198#endif // DEAL_II_WITH_TRILINOS
2199
2200#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
virtual void vmult(MPI::Vector &dst, const MPI::Vector &src) const
Epetra_Operator & trilinos_operator() const
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())
virtual size_type size() const override
iterator begin()
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
static ::ExceptionBase & ExcNonMatchingMaps(std::string arg1)
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define DeclException1(Exception1, type1, outsequence)
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define AssertThrow(cond, exc)
unsigned int global_dof_index
Definition types.h:94
std::vector< std::vector< double > > constant_modes_values
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