deal.II version GIT relicensing-1939-gfdca551f47 2024-10-07 12:00: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)
1715 {
1716 public:
1724 {
1729 AdditionalData(const bool elliptic = true,
1730 const unsigned int n_cycles = 1,
1731 const bool w_cycle = false,
1732 const double aggregation_threshold = 1e-4,
1733 const std::vector<std::vector<bool>> &constant_modes =
1734 std::vector<std::vector<bool>>(0),
1735 const unsigned int smoother_sweeps = 2,
1736 const unsigned int smoother_overlap = 0,
1737 const bool output_details = false,
1738 const char *smoother_type = "Chebyshev",
1739 const char *coarse_type = "Amesos-KLU");
1740
1749
1754 unsigned int n_cycles;
1755
1761
1772
1779 std::vector<std::vector<bool>> constant_modes;
1780
1791 unsigned int smoother_sweeps;
1792
1797 unsigned int smoother_overlap;
1798
1805
1840 const char *smoother_type;
1841
1846 const char *coarse_type;
1847 };
1848
1853
1857 virtual ~PreconditionAMGMueLu() override = default;
1858
1864 void
1865 initialize(const SparseMatrix &matrix,
1866 const AdditionalData &additional_data = AdditionalData());
1867
1874 void
1875 initialize(const Epetra_CrsMatrix &matrix,
1876 const AdditionalData &additional_data = AdditionalData());
1877
1891 void
1892 initialize(const SparseMatrix &matrix,
1893 Teuchos::ParameterList &muelu_parameters);
1894
1900 void
1901 initialize(const Epetra_CrsMatrix &matrix,
1902 Teuchos::ParameterList &muelu_parameters);
1903
1910 template <typename number>
1911 void
1912 initialize(const ::SparseMatrix<number> &deal_ii_sparse_matrix,
1913 const AdditionalData &additional_data = AdditionalData(),
1914 const double drop_tolerance = 1e-13,
1915 const ::SparsityPattern *use_this_sparsity = nullptr);
1916
1921 void
1922 clear();
1923
1927 size_type
1928 memory_consumption() const;
1929
1930 private:
1934 std::shared_ptr<SparseMatrix> trilinos_matrix;
1935 };
1936# endif
1937
1938
1939
1947 {
1948 public:
1954 {};
1955
1962 void
1963 initialize(const SparseMatrix &matrix,
1964 const AdditionalData &additional_data = AdditionalData());
1965
1969 void
1970 vmult(MPI::Vector &dst, const MPI::Vector &src) const override;
1971
1975 void
1976 Tvmult(MPI::Vector &dst, const MPI::Vector &src) const override;
1977
1982 void
1984 const ::Vector<double> &src) const override;
1985
1990 void
1992 const ::Vector<double> &src) const override;
1993
1998 void
2000 const ::LinearAlgebra::distributed::Vector<double> &src)
2001 const override;
2002
2008 void
2010 const ::LinearAlgebra::distributed::Vector<double> &src)
2011 const override;
2012 };
2013
2014
2015
2016 // ----------------------- inline and template functions --------------------
2017
2018
2019# ifndef DOXYGEN
2020
2021
2022 inline void
2024 {
2025 // This only flips a flag that tells
2026 // Trilinos that any vmult operation
2027 // should be done with the
2028 // transpose. However, the matrix
2029 // structure is not reset.
2030 int ierr;
2031
2032 if (!preconditioner->UseTranspose())
2033 {
2034 ierr = preconditioner->SetUseTranspose(true);
2035 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2036 }
2037 else
2038 {
2039 ierr = preconditioner->SetUseTranspose(false);
2040 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2041 }
2042 }
2043
2044
2045 inline void
2046 PreconditionBase::vmult(MPI::Vector &dst, const MPI::Vector &src) const
2047 {
2048 Assert(dst.trilinos_partitioner().SameAs(
2049 preconditioner->OperatorRangeMap()),
2050 ExcNonMatchingMaps("dst"));
2051 Assert(src.trilinos_partitioner().SameAs(
2052 preconditioner->OperatorDomainMap()),
2053 ExcNonMatchingMaps("src"));
2054
2055 const int ierr = preconditioner->ApplyInverse(src.trilinos_vector(),
2056 dst.trilinos_vector());
2057 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2058 }
2059
2060 inline void
2061 PreconditionBase::Tvmult(MPI::Vector &dst, const MPI::Vector &src) const
2062 {
2063 Assert(dst.trilinos_partitioner().SameAs(
2064 preconditioner->OperatorRangeMap()),
2065 ExcNonMatchingMaps("dst"));
2066 Assert(src.trilinos_partitioner().SameAs(
2067 preconditioner->OperatorDomainMap()),
2068 ExcNonMatchingMaps("src"));
2069
2070 preconditioner->SetUseTranspose(true);
2071 const int ierr = preconditioner->ApplyInverse(src.trilinos_vector(),
2072 dst.trilinos_vector());
2073 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2074 preconditioner->SetUseTranspose(false);
2075 }
2076
2077
2078 // For the implementation of the <code>vmult</code> function with deal.II
2079 // data structures we note that invoking a call of the Trilinos
2080 // preconditioner requires us to use Epetra vectors as well. We do this by
2081 // providing a view, i.e., feed Trilinos with a pointer to the data, so we
2082 // avoid copying the content of the vectors during the iteration (this
2083 // function is only useful when used in serial anyway). In the declaration
2084 // of the right hand side, we need to cast the source vector (that is
2085 // <code>const</code> in all deal.II calls) to non-constant value, as this
2086 // is the way Trilinos wants to have them.
2087 inline void
2089 const ::Vector<double> &src) const
2090 {
2091 AssertDimension(dst.size(),
2092 preconditioner->OperatorDomainMap().NumMyElements());
2093 AssertDimension(src.size(),
2094 preconditioner->OperatorRangeMap().NumMyElements());
2095 Epetra_Vector tril_dst(View,
2096 preconditioner->OperatorDomainMap(),
2097 dst.begin());
2098 Epetra_Vector tril_src(View,
2099 preconditioner->OperatorRangeMap(),
2100 const_cast<double *>(src.begin()));
2101
2102 const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2103 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2104 }
2105
2106
2107 inline void
2109 const ::Vector<double> &src) const
2110 {
2111 AssertDimension(dst.size(),
2112 preconditioner->OperatorDomainMap().NumMyElements());
2113 AssertDimension(src.size(),
2114 preconditioner->OperatorRangeMap().NumMyElements());
2115 Epetra_Vector tril_dst(View,
2116 preconditioner->OperatorDomainMap(),
2117 dst.begin());
2118 Epetra_Vector tril_src(View,
2119 preconditioner->OperatorRangeMap(),
2120 const_cast<double *>(src.begin()));
2121
2122 preconditioner->SetUseTranspose(true);
2123 const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2124 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2125 preconditioner->SetUseTranspose(false);
2126 }
2127
2128
2129
2130 inline void
2134 {
2136 preconditioner->OperatorDomainMap().NumMyElements());
2138 preconditioner->OperatorRangeMap().NumMyElements());
2139 Epetra_Vector tril_dst(View,
2140 preconditioner->OperatorDomainMap(),
2141 dst.begin());
2142 Epetra_Vector tril_src(View,
2143 preconditioner->OperatorRangeMap(),
2144 const_cast<double *>(src.begin()));
2145
2146 const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2147 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2148 }
2149
2150 inline void
2154 {
2156 preconditioner->OperatorDomainMap().NumMyElements());
2158 preconditioner->OperatorRangeMap().NumMyElements());
2159 Epetra_Vector tril_dst(View,
2160 preconditioner->OperatorDomainMap(),
2161 dst.begin());
2162 Epetra_Vector tril_src(View,
2163 preconditioner->OperatorRangeMap(),
2164 const_cast<double *>(src.begin()));
2165
2166 preconditioner->SetUseTranspose(true);
2167 const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
2168 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2169 preconditioner->SetUseTranspose(false);
2170 }
2171
2172# endif
2173
2174} // namespace TrilinosWrappers
2175
2176
2181
2182#endif // DEAL_II_WITH_TRILINOS
2183
2184#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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
static ::ExceptionBase & ExcNonMatchingMaps(std::string arg1)
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:511
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define AssertThrow(cond, exc)
unsigned int global_dof_index
Definition types.h:81
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