general_purpose_preconditioners.h
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//====================================================================
26// Include guards
27#ifndef OOMPH_GENERAL_PRECONDITION_HEADER
28#define OOMPH_GENERAL_PRECONDITION_HEADER
29
30
31// Config header generated by autoconfig
32#ifdef HAVE_CONFIG_H
33#include <oomph-lib-config.h>
34#endif
35
36#include "preconditioner.h"
38#include "matrices.h"
39#include "problem.h"
40#include <algorithm>
41
42
43namespace oomph
44{
45 //=====================================================================
46 /// Matrix-based diagonal preconditioner
47 //=====================================================================
49 {
50 public:
51 /// Constructor (empty)
53
54 /// Destructor (empty)
56
57 /// Broken copy constructor
59 delete;
60
61 /// Broken assignment operator
63
64 /// Apply preconditioner to z, i.e. z=D^-1
66
67 /// Setup the preconditioner (store diagonal) from the fully
68 /// assembled matrix.
69 void setup();
70
71 private:
72 /// Vector of inverse diagonal entries
74 };
75
76 //=============================================================================
77 /// Matrix-based lumped preconditioner
78 //=============================================================================
79 template<typename MATRIX>
81 {
82 public:
83 /// Constructor
85 {
86 // default the positive matrix boolean to false
87 Positive_matrix = false;
88
89 // set the pointers to the lumped vector to 0
91 };
92
93 /// Destructor
95 {
96 this->clean_up_memory();
97 }
98
99 /// Broken copy constructor
101 delete;
102
103 /// Broken assignment operator
105
106 /// Apply preconditioner to z, i.e. z=D^-1
108
109 /// Setup the preconditioner (store diagonal) from the fully
110 /// assembled matrix. Problem pointer is ignored.
111 void setup();
112
113 /// For some reason we need to remind the compiler that there is
114 /// also a function named setup in the base class.
116
117 /// Access function to the Positive_matrix which indicates whether
118 /// lumped matrix was positive
119 bool positive_matrix() const
120 {
121 /// paranoid check that preconditioner has been setup
122#ifdef PARANOID
123 if (Inv_lumped_diag_pt == 0)
124 {
125 throw OomphLibError("The preconditioner has not been setup.",
126 OOMPH_CURRENT_FUNCTION,
127 OOMPH_EXCEPTION_LOCATION);
128 }
129#endif
130 return Positive_matrix;
131 }
132
133
134 /// Access function to the inverse of the lumped vector assembled in
135 /// the preconditioner setup routine
137 {
138 /// paranoid check that vector has been created
139#ifdef PARANOID
140 if (Inv_lumped_diag_pt == 0)
141 {
142 throw OomphLibError("The inverse lumped vector has not been created. "
143 "Created in setup(...)",
144 OOMPH_CURRENT_FUNCTION,
145 OOMPH_EXCEPTION_LOCATION);
146 }
147#endif
148 return Inv_lumped_diag_pt;
149 }
150
151
152 /// Access function to number of rows for this preconditioner
153 unsigned& nrow()
154 {
155 return Nrow;
156 }
157
158 /// clean up memory - just delete the inverse lumped vector
160 {
161 delete[] Inv_lumped_diag_pt;
162 }
163
164 private:
165 /// Vector of inverse diagonal entries
167
168 // boolean to indicate whether the lumped matrix was positive
170
171 // number of rows in preconditioner
172 unsigned Nrow;
173 };
174
175
176 //=============================================================================
177 /// Class for a compressed-matrix coefficent (for either CC or CR
178 /// matrices). Contains the (row or column) index and value of a
179 /// coefficient in a compressed row or column.
180 /// Currently only used in ILU(0) for CCDoubleMatrices to allow the
181 /// coefficients in each compressed column [row] to be sorted by
182 /// their row [column] index.
183 //=============================================================================
185 {
186 public:
187 /// Constructor (no arguments)
189
190 /// Constructor (takes the index and value as arguments)
191 CompressedMatrixCoefficient(const unsigned& index, const double& value)
192 {
193 // store the index and value
194 Index = index;
195 Value = value;
196 }
197
198 /// Destructor (does nothing)
200
201 /// Copy Constructor. Not Broken. Required for STL sort function.
203 {
204 Index = a.index();
205 Value = a.value();
206 }
207
208 /// Assignment Operator. Not Broken. Required for STL sort function.
210 {
211 Index = a.index();
212 Value = a.value();
213 }
214
215 /// Less Than Operator (for the STL sort function)
217 {
218 return Index < a.index();
219 }
220
221 /// access function for the coefficient's (row or column) index
222 unsigned& index()
223 {
224 return Index;
225 }
226
227 /// access function for the coefficient value
228 double& value()
229 {
230 return Value;
231 }
232
233 /// Access function for the coefficient's (row or column_ index
234 /// (const version)
235 unsigned index() const
236 {
237 return Index;
238 }
239
240 /// access function for the coefficient's value (const version)
241 double value() const
242 {
243 return Value;
244 }
245
246 private:
247 /// the row or column index of the compressed-matrix coefficient
248 unsigned Index;
249
250 /// the value of the compressed-matrix coefficient
251 double Value;
252 };
253
254
255 //=============================================================================
256 /// ILU(0) Preconditioner
257 //=============================================================================
258 template<typename MATRIX>
260 {
261 };
262
263
264 //=============================================================================
265 /// ILU(0) Preconditioner for matrices of CCDoubleMatrix Format
266 //=============================================================================
267 template<>
269 {
270 public:
271 /// Constructor (empty)
273
274 /// Destructor (empty)
276
277
278 /// Broken copy constructor
280
281 /// Broken assignment operator
282 void operator=(const ILUZeroPreconditioner&) = delete;
283
284 /// Apply preconditioner to r
286
287 /// Setup the preconditioner (store diagonal) from the fully
288 /// assembled matrix. Problem pointer is ignored.
289 void setup();
290
291 private:
292 /// Column start for upper triangular matrix
294
295 /// Row entry for the upper triangular matrix (each element of the
296 /// vector contains the row index and coefficient)
298
299 /// Column start for lower triangular matrix
301
302 /// Row entry for the lower triangular matrix (each element of the
303 /// vector contains the row index and coefficient)
305 };
306
307
308 //=============================================================================
309 /// ILU(0) Preconditioner for matrices of CRDoubleMatrix Format
310 //=============================================================================
311 template<>
313 {
314 public:
315 /// Constructor (empty)
317
318
319 /// Broken copy constructor
321
322 /// Broken assignment operator
323 void operator=(const ILUZeroPreconditioner&) = delete;
324
325 /// Destructor (empty)
327
328 /// Apply preconditioner to r
330
331 /// Setup the preconditioner (store diagonal) from the fully
332 /// assembled matrix. Problem pointer is ignored.
333 void setup();
334
335
336 private:
337 /// Row start for upper triangular matrix
339
340 /// column entry for the upper triangular matrix (each element of the
341 /// vector contains the column index and coefficient)
343
344 /// Row start for lower triangular matrix
346
347 /// column entry for the lower triangular matrix (each element of the
348 /// vector contains the column index and coefficient)
350 };
351
352 //=============================================================================
353 /// A preconditioner for performing inner iteration preconditioner
354 /// solves. The template argument SOLVER specifies the inner iteration
355 /// solver (which must be derived from IterativeLinearSolver) and the
356 /// template argument PRECONDITIONER specifies the preconditioner for the
357 /// inner iteration iterative solver.
358 /// Note: For no preconditioning use the IdentityPreconditioner.
359 //=============================================================================
360 template<class SOLVER, class PRECONDITIONER>
362 {
363 public:
364 /// Constructor
366 {
367 // create the solver
368 Solver_pt = new SOLVER;
369
370 // create the preconditioner
371 Preconditioner_pt = new PRECONDITIONER;
372
373#ifdef PARANOID
374 // paranoid check that the solver is an iterative solver and
375 // the preconditioner is a preconditioner
376 if (dynamic_cast<IterativeLinearSolver*>(Solver_pt) == 0)
377 {
378 throw OomphLibError(
379 "The template argument SOLVER must be of type IterativeLinearSolver",
380 OOMPH_CURRENT_FUNCTION,
381 OOMPH_EXCEPTION_LOCATION);
382 }
383 if (dynamic_cast<Preconditioner*>(Preconditioner_pt) == 0)
384 {
385 throw OomphLibError(
386 "The template argument PRECONDITIONER must be of type Preconditioner",
387 OOMPH_CURRENT_FUNCTION,
388 OOMPH_EXCEPTION_LOCATION);
389 }
390#endif
391
392 // ensure the solver does not re-setup the preconditioner
393 Solver_pt->disable_setup_preconditioner_before_solve();
394
395 // pass the preconditioner to the solver
396 Solver_pt->preconditioner_pt() = Preconditioner_pt;
397 }
398
399 // destructor
401 {
402 delete Solver_pt;
403 delete Preconditioner_pt;
404 }
405
406 // clean the memory
408 {
409 Preconditioner_pt->clean_up_memory();
410 Solver_pt->clean_up_memory();
411 }
412
413 /// Preconditioner setup method. Setup the preconditioner for the
414 /// inner iteration solver.
415 void setup()
416 {
417 // set the distribution
420 if (dist_pt != 0)
421 {
422 this->build_distribution(dist_pt->distribution_pt());
423 }
424 else
425 {
426 LinearAlgebraDistribution dist(comm_pt(), matrix_pt()->nrow(), false);
427 this->build_distribution(dist);
428 }
429
430 // Setup the inner iteration preconditioner (For some reason we need to
431 // remind the compiler that there is also a function named setup in the
432 // base class.)
433 Preconditioner_pt->Preconditioner::setup(matrix_pt());
434
435 // setup the solverready for resolve
436 unsigned max_iter = Solver_pt->max_iter();
437 Solver_pt->max_iter() = 1;
438 DoubleVector x(this->distribution_pt(), 0.0);
439 DoubleVector y(x);
440 Solver_pt->enable_resolve();
441 Solver_pt->solve(matrix_pt(), x, y);
442 Solver_pt->max_iter() = max_iter;
443 }
444
445 /// Preconditioner solve method. Performs the specified number
446 /// of Krylov iterations preconditioned with the specified preconditioner
448 {
449 Solver_pt->resolve(r, z);
450 }
451
452 /// Access to convergence tolerance of the inner iteration solver
453 double& tolerance()
454 {
455 return Solver_pt->tolerance();
456 }
457
458 /// Access to max. number of iterations of the inner iteration solver
459 unsigned& max_iter()
460 {
461 return Solver_pt->max_iter();
462 }
463
464 ///
465 SOLVER* solver_pt()
466 {
467 return Solver_pt;
468 }
469
470 ///
471 PRECONDITIONER* preconditioner_pt()
472 {
473 return Preconditioner_pt;
474 }
475
476 private:
477 /// pointer to the underlying solver
478 SOLVER* Solver_pt;
479
480 /// pointer to the underlying preconditioner
481 PRECONDITIONER* Preconditioner_pt;
482 };
483} // namespace oomph
484#endif
//////////////////////////////////////////////////////////////// ////////////////////////////////////...
Definition: matrices.h:2791
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:888
Class for a compressed-matrix coefficent (for either CC or CR matrices). Contains the (row or column)...
unsigned & index()
access function for the coefficient's (row or column) index
bool operator<(const CompressedMatrixCoefficient &a) const
Less Than Operator (for the STL sort function)
CompressedMatrixCoefficient()
Constructor (no arguments)
CompressedMatrixCoefficient(const unsigned &index, const double &value)
Constructor (takes the index and value as arguments)
double value() const
access function for the coefficient's value (const version)
double Value
the value of the compressed-matrix coefficient
CompressedMatrixCoefficient(const CompressedMatrixCoefficient &a)
Copy Constructor. Not Broken. Required for STL sort function.
double & value()
access function for the coefficient value
void operator=(const CompressedMatrixCoefficient &a)
Assignment Operator. Not Broken. Required for STL sort function.
unsigned Index
the row or column index of the compressed-matrix coefficient
unsigned index() const
Access function for the coefficient's (row or column_ index (const version)
Base class for any linear algebra object that is distributable. Just contains storage for the LinearA...
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
A vector in the mathematical sense, initially developed for linear algebra type applications....
Definition: double_vector.h:58
Vector< unsigned > L_column_start
Column start for lower triangular matrix.
Vector< CompressedMatrixCoefficient > U_row_entry
Row entry for the upper triangular matrix (each element of the vector contains the row index and coef...
void operator=(const ILUZeroPreconditioner &)=delete
Broken assignment operator.
Vector< CompressedMatrixCoefficient > L_row_entry
Row entry for the lower triangular matrix (each element of the vector contains the row index and coef...
Vector< unsigned > U_column_start
Column start for upper triangular matrix.
ILUZeroPreconditioner(const ILUZeroPreconditioner &)=delete
Broken copy constructor.
ILUZeroPreconditioner(const ILUZeroPreconditioner &)=delete
Broken copy constructor.
Vector< CompressedMatrixCoefficient > U_row_entry
column entry for the upper triangular matrix (each element of the vector contains the column index an...
Vector< unsigned > U_row_start
Row start for upper triangular matrix.
Vector< unsigned > L_row_start
Row start for lower triangular matrix.
void operator=(const ILUZeroPreconditioner &)=delete
Broken assignment operator.
Vector< CompressedMatrixCoefficient > L_row_entry
column entry for the lower triangular matrix (each element of the vector contains the column index an...
A preconditioner for performing inner iteration preconditioner solves. The template argument SOLVER s...
void setup()
Preconditioner setup method. Setup the preconditioner for the inner iteration solver.
void clean_up_memory()
Clean up memory (empty). Generic interface function.
double & tolerance()
Access to convergence tolerance of the inner iteration solver.
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Preconditioner solve method. Performs the specified number of Krylov iterations preconditioned with t...
SOLVER * Solver_pt
pointer to the underlying solver
unsigned & max_iter()
Access to max. number of iterations of the inner iteration solver.
PRECONDITIONER * Preconditioner_pt
pointer to the underlying preconditioner
Base class for all linear iterative solvers. This merely defines standard interfaces for linear itera...
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
Matrix-based diagonal preconditioner.
void operator=(const MatrixBasedDiagPreconditioner &)=delete
Broken assignment operator.
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply preconditioner to z, i.e. z=D^-1.
Vector< double > Inv_diag
Vector of inverse diagonal entries.
MatrixBasedDiagPreconditioner(const MatrixBasedDiagPreconditioner &)=delete
Broken copy constructor.
void setup()
Setup the preconditioner (store diagonal) from the fully assembled matrix.
void clean_up_memory()
clean up memory - just delete the inverse lumped vector
MatrixBasedLumpedPreconditioner(const MatrixBasedDiagPreconditioner &)=delete
Broken copy constructor.
bool positive_matrix() const
Access function to the Positive_matrix which indicates whether lumped matrix was positive.
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply preconditioner to z, i.e. z=D^-1.
double * inverse_lumped_vector_pt()
Access function to the inverse of the lumped vector assembled in the preconditioner setup routine.
double * Inv_lumped_diag_pt
Vector of inverse diagonal entries.
void setup()
Setup the preconditioner (store diagonal) from the fully assembled matrix. Problem pointer is ignored...
unsigned & nrow()
Access function to number of rows for this preconditioner.
void operator=(const MatrixBasedLumpedPreconditioner &)=delete
Broken assignment operator.
An OomphLibError object which should be thrown when an run-time error is encountered....
Preconditioner base class. Gives an interface to call all other preconditioners through and stores th...
virtual DoubleMatrixBase * matrix_pt() const
Get function for matrix pointer.
virtual const OomphCommunicator * comm_pt() const
Get function for comm pointer.
virtual void preconditioner_solve(const DoubleVector &r, DoubleVector &z)=0
Apply the preconditioner. Pure virtual generic interface function. This method should apply the preco...
virtual void setup()=0
Setup the preconditioner. Pure virtual generic interface function.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...