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-2023 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 
43 namespace 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.
115  using Preconditioner::setup;
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.
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)
double & value()
access function for the coefficient value
Base class for any linear algebra object that is distributable. Just contains storage for the LinearA...
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
unsigned nrow() const
access function to the number of global rows.
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...
unsigned & max_iter()
Access to max. number of iterations of the inner iteration solver.
SOLVER * Solver_pt
pointer to the underlying 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 * 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...
void operator=(const MatrixBasedLumpedPreconditioner &)=delete
Broken assignment operator.
double * inverse_lumped_vector_pt()
Access function to the inverse of the lumped vector assembled in the preconditioner setup routine.
unsigned & nrow()
Access function to number of rows for this preconditioner.
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 void preconditioner_solve(const DoubleVector &r, DoubleVector &z)=0
Apply the preconditioner. Pure virtual generic interface function. This method should apply the preco...
virtual const OomphCommunicator * comm_pt() const
Get function for comm pointer.
virtual void setup()=0
Setup the preconditioner. Pure virtual generic interface function.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...