Reference documentation for deal.II version GIT relicensing-422-gb369f187d8 2024-04-17 18:10:02+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
1512 unsigned int smoother_sweeps;
1513
1518 unsigned int smoother_overlap;
1519
1526
1561 const char *smoother_type;
1562
1567 const char *coarse_type;
1568 };
1569
1573 ~PreconditionAMG() override;
1574
1575
1581 void
1582 initialize(const SparseMatrix &matrix,
1583 const AdditionalData &additional_data = AdditionalData());
1584
1603 void
1604 initialize(const Epetra_RowMatrix &matrix,
1605 const AdditionalData &additional_data = AdditionalData());
1606
1621 void
1622 initialize(const SparseMatrix &matrix,
1623 const Teuchos::ParameterList &ml_parameters);
1624
1632 void
1633 initialize(const Epetra_RowMatrix &matrix,
1634 const Teuchos::ParameterList &ml_parameters);
1635
1642 template <typename number>
1643 void
1644 initialize(const ::SparseMatrix<number> &deal_ii_sparse_matrix,
1645 const AdditionalData &additional_data = AdditionalData(),
1646 const double drop_tolerance = 1e-13,
1647 const ::SparsityPattern *use_this_sparsity = nullptr);
1648
1661 void
1662 reinit();
1663
1668 void
1669 clear();
1670
1674 size_type
1675 memory_consumption() const;
1676
1677 private:
1681 std::shared_ptr<SparseMatrix> trilinos_matrix;
1682 };
1683
1684
1685
1686# if defined(DOXYGEN) || defined(DEAL_II_TRILINOS_WITH_MUELU)
1706 {
1707 public:
1715 {
1720 AdditionalData(const bool elliptic = true,
1721 const unsigned int n_cycles = 1,
1722 const bool w_cycle = false,
1723 const double aggregation_threshold = 1e-4,
1724 const std::vector<std::vector<bool>> &constant_modes =
1725 std::vector<std::vector<bool>>(0),
1726 const unsigned int smoother_sweeps = 2,
1727 const unsigned int smoother_overlap = 0,
1728 const bool output_details = false,
1729 const char *smoother_type = "Chebyshev",
1730 const char *coarse_type = "Amesos-KLU");
1731
1740
1745 unsigned int n_cycles;
1746
1752
1763
1770 std::vector<std::vector<bool>> constant_modes;
1771
1782 unsigned int smoother_sweeps;
1783
1788 unsigned int smoother_overlap;
1789
1796
1831 const char *smoother_type;
1832
1837 const char *coarse_type;
1838 };
1839
1844
1848 virtual ~PreconditionAMGMueLu() override = default;
1849
1855 void
1856 initialize(const SparseMatrix &matrix,
1857 const AdditionalData &additional_data = AdditionalData());
1858
1865 void
1866 initialize(const Epetra_CrsMatrix &matrix,
1867 const AdditionalData &additional_data = AdditionalData());
1868
1882 void
1883 initialize(const SparseMatrix &matrix,
1884 Teuchos::ParameterList &muelu_parameters);
1885
1891 void
1892 initialize(const Epetra_CrsMatrix &matrix,
1893 Teuchos::ParameterList &muelu_parameters);
1894
1901 template <typename number>
1902 void
1903 initialize(const ::SparseMatrix<number> &deal_ii_sparse_matrix,
1904 const AdditionalData &additional_data = AdditionalData(),
1905 const double drop_tolerance = 1e-13,
1906 const ::SparsityPattern *use_this_sparsity = nullptr);
1907
1912 void
1913 clear();
1914
1918 size_type
1919 memory_consumption() const;
1920
1921 private:
1925 std::shared_ptr<SparseMatrix> trilinos_matrix;
1926 };
1927# endif
1928
1929
1930
1938 {
1939 public:
1945 {};
1946
1953 void
1954 initialize(const SparseMatrix &matrix,
1955 const AdditionalData &additional_data = AdditionalData());
1956
1960 void
1961 vmult(MPI::Vector &dst, const MPI::Vector &src) const override;
1962
1966 void
1967 Tvmult(MPI::Vector &dst, const MPI::Vector &src) const override;
1968
1973 void
1975 const ::Vector<double> &src) const override;
1976
1981 void
1983 const ::Vector<double> &src) const override;
1984
1989 void
1991 const ::LinearAlgebra::distributed::Vector<double> &src)
1992 const override;
1993
1999 void
2001 const ::LinearAlgebra::distributed::Vector<double> &src)
2002 const override;
2003 };
2004
2005
2006
2007 // ----------------------- inline and template functions --------------------
2008
2009
2010# ifndef DOXYGEN
2011
2012
2013 inline void
2015 {
2016 // This only flips a flag that tells
2017 // Trilinos that any vmult operation
2018 // should be done with the
2019 // transpose. However, the matrix
2020 // structure is not reset.
2021 int ierr;
2022
2023 if (!preconditioner->UseTranspose())
2024 {
2025 ierr = preconditioner->SetUseTranspose(true);
2026 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2027 }
2028 else
2029 {
2030 ierr = preconditioner->SetUseTranspose(false);
2031 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2032 }
2033 }
2034
2035
2036 inline void
2037 PreconditionBase::vmult(MPI::Vector &dst, const MPI::Vector &src) const
2038 {
2039 Assert(dst.trilinos_partitioner().SameAs(
2040 preconditioner->OperatorRangeMap()),
2041 ExcNonMatchingMaps("dst"));
2042 Assert(src.trilinos_partitioner().SameAs(
2043 preconditioner->OperatorDomainMap()),
2044 ExcNonMatchingMaps("src"));
2045
2046 const int ierr = preconditioner->ApplyInverse(src.trilinos_vector(),
2047 dst.trilinos_vector());
2048 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2049 }
2050
2051 inline void
2052 PreconditionBase::Tvmult(MPI::Vector &dst, const MPI::Vector &src) const
2053 {
2054 Assert(dst.trilinos_partitioner().SameAs(
2055 preconditioner->OperatorRangeMap()),
2056 ExcNonMatchingMaps("dst"));
2057 Assert(src.trilinos_partitioner().SameAs(
2058 preconditioner->OperatorDomainMap()),
2059 ExcNonMatchingMaps("src"));
2060
2061 preconditioner->SetUseTranspose(true);
2062 const int ierr = preconditioner->ApplyInverse(src.trilinos_vector(),
2063 dst.trilinos_vector());
2064 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2065 preconditioner->SetUseTranspose(false);
2066 }
2067
2068
2069 // For the implementation of the <code>vmult</code> function with deal.II
2070 // data structures we note that invoking a call of the Trilinos
2071 // preconditioner requires us to use Epetra vectors as well. We do this by
2072 // providing a view, i.e., feed Trilinos with a pointer to the data, so we
2073 // avoid copying the content of the vectors during the iteration (this
2074 // function is only useful when used in serial anyway). In the declaration
2075 // of the right hand side, we need to cast the source vector (that is
2076 // <code>const</code> in all deal.II calls) to non-constant value, as this
2077 // is the way Trilinos wants to have them.
2078 inline void
2080 const ::Vector<double> &src) const
2081 {
2082 AssertDimension(dst.size(),
2083 preconditioner->OperatorDomainMap().NumMyElements());
2084 AssertDimension(src.size(),
2085 preconditioner->OperatorRangeMap().NumMyElements());
2086 Epetra_Vector tril_dst(View,
2087 preconditioner->OperatorDomainMap(),
2088 dst.begin());
2089 Epetra_Vector tril_src(View,
2090 preconditioner->OperatorRangeMap(),
2091 const_cast<double *>(src.begin()));
2092
2093 const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2094 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2095 }
2096
2097
2098 inline void
2100 const ::Vector<double> &src) const
2101 {
2102 AssertDimension(dst.size(),
2103 preconditioner->OperatorDomainMap().NumMyElements());
2104 AssertDimension(src.size(),
2105 preconditioner->OperatorRangeMap().NumMyElements());
2106 Epetra_Vector tril_dst(View,
2107 preconditioner->OperatorDomainMap(),
2108 dst.begin());
2109 Epetra_Vector tril_src(View,
2110 preconditioner->OperatorRangeMap(),
2111 const_cast<double *>(src.begin()));
2112
2113 preconditioner->SetUseTranspose(true);
2114 const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2115 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2116 preconditioner->SetUseTranspose(false);
2117 }
2118
2119
2120
2121 inline void
2125 {
2127 preconditioner->OperatorDomainMap().NumMyElements());
2129 preconditioner->OperatorRangeMap().NumMyElements());
2130 Epetra_Vector tril_dst(View,
2131 preconditioner->OperatorDomainMap(),
2132 dst.begin());
2133 Epetra_Vector tril_src(View,
2134 preconditioner->OperatorRangeMap(),
2135 const_cast<double *>(src.begin()));
2136
2137 const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2138 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2139 }
2140
2141 inline void
2145 {
2147 preconditioner->OperatorDomainMap().NumMyElements());
2149 preconditioner->OperatorRangeMap().NumMyElements());
2150 Epetra_Vector tril_dst(View,
2151 preconditioner->OperatorDomainMap(),
2152 dst.begin());
2153 Epetra_Vector tril_src(View,
2154 preconditioner->OperatorRangeMap(),
2155 const_cast<double *>(src.begin()));
2156
2157 preconditioner->SetUseTranspose(true);
2158 const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2159 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2160 preconditioner->SetUseTranspose(false);
2161 }
2162
2163# endif
2164
2165} // namespace TrilinosWrappers
2166
2167
2172
2173#endif // DEAL_II_WITH_TRILINOS
2174
2175#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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcNonMatchingMaps(std::string arg1)
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:516
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define AssertThrow(cond, exc)
unsigned int global_dof_index
Definition types.h:81
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