trilinos_preconditioners.cc
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
27
28namespace oomph
29{
30 /// ////////////////////////////////////////////////////////////////////////////
31 /// ////////////////////////////////////////////////////////////////////////////
32 // functions for the TrilinosPreconditionerBase class
33 /// ////////////////////////////////////////////////////////////////////////////
34 /// ////////////////////////////////////////////////////////////////////////////
35
36
37 //=============================================================================
38 /// Static double that accumulates the preconditioner
39 /// solve time of all instantiations of this class. Reset
40 /// it manually, e.g. after every Newton solve.
41 //=============================================================================
43
44
45 //=============================================================================
46 /// Function to set up a preconditioner for the linear system
47 /// defined by matrix_pt. This function must be called before using
48 /// preconditioner_solve.
49 /// \b NOTE 1. matrix_pt must point to an object of class CRDoubleMatrix or
50 /// DistributedCRDoubleMatrix
51 /// This method should be called by oomph-lib solvers and preconditioners
52 //=============================================================================
54 {
55 // clean up the memory
57
58#ifdef PARANOID
59 // check the matrix is square
60 if (matrix_pt()->nrow() != matrix_pt()->ncol())
61 {
62 std::ostringstream error_message;
63 error_message << "Preconditioners require a square matrix. "
64 << "Matrix is " << matrix_pt()->nrow() << " by "
65 << matrix_pt()->ncol() << std::endl;
66 throw OomphLibError(
67 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
68 }
69#endif
70
71
72 // get a pointer to the cr double matrix
73 CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt());
74
75#ifdef PARANOID
76 if (cr_matrix_pt == 0)
77 {
78 throw OomphLibError("Matrix must be of type CRDoubleMatrix",
79 OOMPH_CURRENT_FUNCTION,
80 OOMPH_EXCEPTION_LOCATION);
81 }
82#endif
83
84 // build the distribution of this preconditioner
85 this->build_distribution(cr_matrix_pt->distribution_pt());
86
87 // create the matrices
89 cr_matrix_pt, this->distribution_pt());
90
91 // set up preconditioner
93 }
94
95 //===========================================================================
96 /// Function to setup a preconditioner for the linear system defined
97 /// by the oomph-lib oomph_matrix_pt and Epetra epetra_matrix_pt matrices.
98 /// This method is called by Trilinos solvers.
99 //===========================================================================
100 void TrilinosPreconditionerBase::setup(Epetra_CrsMatrix* epetra_matrix_pt)
101 {
102 // clean up old data
104
105 // first try CRDoubleMatrix
106 CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt());
107 if (cr_matrix_pt == 0)
108 {
109 std::ostringstream error_message;
110 error_message << "TrilinosSolver only work with "
111 << "DistributedCRDoubleMatrix matrices" << std::endl;
112 throw OomphLibError(
113 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
114 }
115
116 // setup the specific preconditioner
117 setup_trilinos_preconditioner(epetra_matrix_pt);
118 }
119
120
121 //=============================================================================
122 /// preconditioner_solve - applies the preconditioner to the vector r
123 /// taking distributed oomph-lib vectors (DistributedVector<double>) as
124 /// arguments.
125 //=============================================================================
127 DoubleVector& z)
128 {
129 // Get ready to do cumulative timings
130 double t_start = TimingHelpers::timer();
131
132#ifdef PARANOID
133 // check preconditioner data exists
135 {
136 std::ostringstream error_message;
137 error_message << "preconditioner_solve requires that solver data has "
138 << "been set up" << std::endl;
139 throw OomphLibError(
140 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
141 }
142
143 if (*this->distribution_pt() != r.distribution_pt())
144 {
145 std::ostringstream error_message;
146 error_message << "The rhs vector and the matrix must have the same "
147 << "distribution.\n";
148 throw OomphLibError(
149 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
150 }
151#endif
152
153 // create Epetra version of r
154 Epetra_Vector* epetra_r_pt =
156
157 // create an empty Epetra vector for z
158 Epetra_Vector* epetra_z_pt =
160 this->distribution_pt());
161
162 // Apply preconditioner
163 Epetra_preconditioner_pt->ApplyInverse(*epetra_r_pt, *epetra_z_pt);
164
165 // Copy result to z
166 z.build(this->distribution_pt(), 0.0);
168
169 // clean up memory
170 delete epetra_r_pt;
171 delete epetra_z_pt;
172
173
174 // Add to cumulative solve time
175 double t_end = TimingHelpers::timer();
176 Cumulative_preconditioner_solve_time += (t_end - t_start);
177 }
178
179
180 /// ////////////////////////////////////////////////////////////////////////////
181 /// ////////////////////////////////////////////////////////////////////////////
182 // function for the TrilinosML class
183 /// ////////////////////////////////////////////////////////////////////////////
184 /// ////////////////////////////////////////////////////////////////////////////
185
186
187 //=============================================================================
188 /// (Static) default number of V cycles (one to be consistent
189 /// with previous default). (It's an int because Trilinos wants it to be!)
190 //=============================================================================
192
193 //=============================================================================
194 // Function to set up the ML preconditioner.
195 //=============================================================================
197 Epetra_CrsMatrix* epetra_matrix_pt)
198 {
199 // doc setup time
200 oomph_info << "Setting up TrilinosML, ";
201 double t_start = TimingHelpers::timer();
202
203
204 // create the preconditioner
205 Epetra_preconditioner_pt = new ML_Epetra::MultiLevelPreconditioner(
206 *epetra_matrix_pt, ML_parameters, true);
207
208 // doc setup time
209 oomph_info << "time for setup [s] : " << TimingHelpers::timer() - t_start
210 << std::endl;
211
212 // oomph_info << "In here\n";
213 // ML_Epetra::MultiLevelPreconditioner* tmp_pt=0;
214 // tmp_pt=dynamic_cast<ML_Epetra::MultiLevelPreconditioner*>(
215 // Epetra_preconditioner_pt);
216 // if (tmp_pt!=0)
217 // {
218 // oomph_info << "Doing test\n";
219 // ML_parameters.print(*(oomph_info.stream_pt()));
220 // oomph_info.stream_pt()->flush();
221 // // tmp_pt->TestSmoothers();
222 // oomph_info << "Done test\n";
223 // }
224 }
225
226
227 /// ////////////////////////////////////////////////////////////////////////////
228 /// ////////////////////////////////////////////////////////////////////////////
229 // function for the TrilinosIFPACK
230 /// ////////////////////////////////////////////////////////////////////////////
231 /// ////////////////////////////////////////////////////////////////////////////
232
233
234 //=============================================================================
235 // Function to set up the IFPACK preconditioner.
236 //=============================================================================
238 Epetra_CrsMatrix* epetra_matrix_pt)
239 {
240 // set up a parameter list
241 Teuchos::ParameterList ifpack_parameters;
242 ifpack_parameters.set("fact: level-of-fill", ILU_fill_level);
243 ifpack_parameters.set("fact: ilut level-of-fill", ILUT_fill_level);
244 ifpack_parameters.set("fact: absolute threshold", Absolute_threshold);
245 ifpack_parameters.set("fact: relative threshold", Relative_threshold);
246
247 // create the preconditioner
248 Ifpack ifpack_factory;
249 Ifpack_Preconditioner* tmp_pt =
250 ifpack_factory.Create(Preconditioner_type, epetra_matrix_pt, Overlap);
251
252 tmp_pt->SetParameters(ifpack_parameters);
253 tmp_pt->Initialize();
254 tmp_pt->Compute();
255
256 // Now store the pointer to the newly created Ifpack_Preconditioner
258 }
259
260
261 /// ////////////////////////////////////////////////////////////////////////////
262 /// ////////////////////////////////////////////////////////////////////////////
263 /// ////////////////////////////////////////////////////////////////////////////
264
265
266} // namespace oomph
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:888
unsigned nrow() const
access function to the number of global rows.
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
void build_distribution(const LinearAlgebraDistribution *const dist_pt)
setup the distribution of this distributable linear algebra object
virtual unsigned long ncol() const =0
Return the number of columns of the matrix.
virtual unsigned long nrow() const =0
Return the number of rows of the matrix.
A vector in the mathematical sense, initially developed for linear algebra type applications....
Definition: double_vector.h:58
void build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
An OomphLibError object which should be thrown when an run-time error is encountered....
virtual DoubleMatrixBase * matrix_pt() const
Get function for matrix pointer.
double Absolute_threshold
Value of absolute threshold, used to peturb diagonal.
string Preconditioner_type
Type of ILU preconditioner.
int ILU_fill_level
Level of fill for "ILU".
void setup_trilinos_preconditioner(Epetra_CrsMatrix *epetra_matrix_pt)
Function to set up an IFPACK preconditioner. It is assumed Trilinos_matrix_pt points to a suitable ma...
double ILUT_fill_level
Level of fill for "ILUT".
double Relative_threshold
Value of relative threshold, used to pertub diagonal.
int Overlap
Value of overlap level - used in parallel ILU.
static int Default_n_cycles
Default number of V cycles (one to be consistent with previous default) (It's an int because Trilinos...
void setup_trilinos_preconditioner(Epetra_CrsMatrix *epetra_matrix_pt)
Function to set up the ML preconditioner. It is assumed Trilinos_matrix_pt points to a suitable matri...
Epetra_Operator * Epetra_preconditioner_pt
The preconditioner which will be set up using function setup_trilinos_preconditioner(....
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
applies the preconditioner
Epetra_CrsMatrix * Epetra_matrix_pt
Pointer used to store the epetra matrix - only used when this preconditioner is setup using the oomph...
void clean_up_memory()
deletes the preconditioner, matrices and maps
static double Cumulative_preconditioner_solve_time
Static double that accumulates the preconditioner solve time of all instantiations of this class....
void setup()
Broken assignment operator.
virtual void setup_trilinos_preconditioner(Epetra_CrsMatrix *epetra_matrix_pt)=0
Function to set up a specific Trilinos preconditioner. This is called by setup(......
double timer()
returns the time in seconds after some point in past
Epetra_CrsMatrix * create_distributed_epetra_matrix(const CRDoubleMatrix *oomph_matrix_pt, const LinearAlgebraDistribution *dist_pt)
create an Epetra_CrsMatrix from an oomph-lib CRDoubleMatrix. If oomph_matrix_pt is NOT distributed (i...
void copy_to_oomphlib_vector(const Epetra_Vector *epetra_vec_pt, DoubleVector &oomph_vec)
Helper function to copy the contents of a Trilinos vector to an oomph-lib distributed vector....
Epetra_Vector * create_distributed_epetra_vector(const DoubleVector &oomph_vec)
create an Epetra_Vector from an oomph-lib DoubleVector. If oomph_vec is NOT distributed (i....
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...