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// Version 1.0; svn revision $LastChangedRevision: 1132 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-02-23 05:43:26 +0000 (Tue, 23 Feb 2016) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //Include guards
31 #ifndef OOMPH_GENERAL_PRECONDITION_HEADER
32 #define OOMPH_GENERAL_PRECONDITION_HEADER
33 
34 
35 // Config header generated by autoconfig
36 #ifdef HAVE_CONFIG_H
37  #include <oomph-lib-config.h>
38 #endif
39 
40 #include "preconditioner.h"
42 #include "matrices.h"
43 #include "problem.h"
44 #include <algorithm>
45 
46 
47 
48 namespace oomph
49 {
50 
51 
52 
53 //=====================================================================
54 /// Matrix-based diagonal preconditioner
55 //=====================================================================
57 {
58  public:
59 
60  ///Constructor (empty)
62 
63  ///Destructor (empty)
65 
66  /// Broken copy constructor
68  {
69  BrokenCopy::broken_copy("MatrixBasedDiagPreconditioner");
70  }
71 
72  /// Broken assignment operator
74  {
75  BrokenCopy::broken_assign("MatrixBasedDiagPreconditioner");
76  }
77 
78  /// Apply preconditioner to z, i.e. z=D^-1
80 
81  /// \short Setup the preconditioner (store diagonal) from the fully
82  /// assembled matrix.
83  void setup();
84 
85  private:
86 
87  /// Vector of inverse diagonal entries
89 };
90 
91 //=============================================================================
92 /// Matrix-based lumped preconditioner
93 //=============================================================================
94 template<typename MATRIX>
96 {
97  public:
98 
99  ///Constructor
101  {
102  // default the positive matrix boolean to false
103  Positive_matrix = false;
104 
105  // set the pointers to the lumped vector to 0
106  Inv_lumped_diag_pt = 0;
107  };
108 
109  ///Destructor
111  {
112  this->clean_up_memory();
113  }
114 
115  /// Broken copy constructor
117  {
118  BrokenCopy::broken_copy("MatrixBasedDiagPreconditioner");
119  }
120 
121  /// Broken assignment operator
123  {
124  BrokenCopy::broken_assign("MatrixBasedDiagPreconditioner");
125  }
126 
127  /// Apply preconditioner to z, i.e. z=D^-1
129 
130  /// \short Setup the preconditioner (store diagonal) from the fully
131  /// assembled matrix. Problem pointer is ignored.
132  void setup();
133 
134  /// \short For some reason we need to remind the compiler that there is
135  /// also a function named setup in the base class.
136  using Preconditioner::setup;
137 
138  /// \short Access function to the Positive_matrix which indicates whether
139  /// lumped matrix was positive
140  bool positive_matrix() const
141  {
142  /// paranoid check that preconditioner has been setup
143 #ifdef PARANOID
144  if (Inv_lumped_diag_pt == 0)
145  {
146  throw OomphLibError(
147  "The preconditioner has not been setup.",
148  OOMPH_CURRENT_FUNCTION,
149  OOMPH_EXCEPTION_LOCATION);
150  }
151 #endif
152  return Positive_matrix;
153  }
154 
155 
156  /// \short Access function to the inverse of the lumped vector assembled in
157  /// the preconditioner setup routine
159  {
160  /// paranoid check that vector has been created
161 #ifdef PARANOID
162  if (Inv_lumped_diag_pt == 0)
163  {
164  throw OomphLibError(
165  "The inverse lumped vector has not been created. Created in setup(...)",
166  OOMPH_CURRENT_FUNCTION,
167  OOMPH_EXCEPTION_LOCATION);
168  }
169 #endif
170  return Inv_lumped_diag_pt;
171  }
172 
173 
174  /// \short Access function to number of rows for this preconditioner
175  unsigned& nrow() { return Nrow; }
176 
177  /// clean up memory - just delete the inverse lumped vector
179  {
180  delete[] Inv_lumped_diag_pt;
181  }
182 
183  private:
184 
185 
186  /// Vector of inverse diagonal entries
188 
189  // boolean to indicate whether the lumped matrix was positive
191 
192  // number of rows in preconditioner
193  unsigned Nrow;
194 };
195 
196 
197 
198 //=============================================================================
199 /// \short Class for a compressed-matrix coefficent (for either CC or CR
200 /// matrices). Contains the (row or column) index and value of a
201 /// coefficient in a compressed row or column.
202 /// Currently only used in ILU(0) for CCDoubleMatrices to allow the
203 /// coefficients in each compressed column [row] to be sorted by
204 /// their row [column] index.
205 //=============================================================================
207 {
208 
209  public:
210 
211  /// Constructor (no arguments)
213 
214  /// Constructor (takes the index and value as arguments)
215  CompressedMatrixCoefficient(const unsigned& index, const double& value)
216  {
217  // store the index and value
218  Index = index;
219  Value = value;
220  }
221 
222  /// Destructor (does nothing)
224 
225  /// Copy Constructor. Not Broken. Required for STL sort function.
227  {
228  Index = a.index();
229  Value = a.value();
230  }
231 
232  /// Assignment Operator. Not Broken. Required for STL sort function.
234  {
235  Index = a.index();
236  Value = a.value();
237  }
238 
239  /// Less Than Operator (for the STL sort function)
241  {
242  return Index < a.index();
243  }
244 
245  /// access function for the coefficient's (row or column) index
246  unsigned& index() { return Index; }
247 
248  /// access function for the coefficient value
249  double& value() { return Value; }
250 
251  /// \short Access function for the coefficient's (row or column_ index
252  /// (const version)
253  unsigned index() const { return Index; }
254 
255  /// access function for the coefficient's value (const version)
256  double value() const { return Value; }
257 
258  private:
259 
260  /// the row or column index of the compressed-matrix coefficient
261  unsigned Index;
262 
263  /// the value of the compressed-matrix coefficient
264  double Value;
265 
266 };
267 
268 
269 
270 
271 
272 //=============================================================================
273 /// ILU(0) Preconditioner
274 //=============================================================================
275 template <typename MATRIX>
277 {
278 };
279 
280 
281 //=============================================================================
282 /// ILU(0) Preconditioner for matrices of CCDoubleMatrix Format
283 //=============================================================================
284 template <>
286 {
287 
288  public :
289 
290  ///Constructor (empty)
292 
293  ///Destructor (empty)
295 
296 
297  /// Broken copy constructor
299  {
300  BrokenCopy::broken_copy("ILUZeroPreconditioner");
301  }
302 
303  /// Broken assignment operator
305  {
306  BrokenCopy::broken_assign("ILUZeroPreconditioner");
307  }
308 
309 
310  /// Apply preconditioner to r
312 
313  /// \short Setup the preconditioner (store diagonal) from the fully
314  /// assembled matrix. Problem pointer is ignored.
315  void setup();
316 
317  private:
318 
319  /// Column start for upper triangular matrix
321 
322  /// \short Row entry for the upper triangular matrix (each element of the
323  /// vector contains the row index and coefficient)
325 
326  /// Column start for lower triangular matrix
328 
329  /// \short Row entry for the lower triangular matrix (each element of the
330  /// vector contains the row index and coefficient)
332 
333 };
334 
335 
336 //=============================================================================
337 /// ILU(0) Preconditioner for matrices of CRDoubleMatrix Format
338 //=============================================================================
339 template <>
341 {
342 
343  public :
344 
345  ///Constructor (empty)
347 
348 
349  /// Broken copy constructor
351  {
352  BrokenCopy::broken_copy("ILUZeroPreconditioner");
353  }
354 
355  /// Broken assignment operator
357  {
358  BrokenCopy::broken_assign("ILUZeroPreconditioner");
359  }
360 
361  ///Destructor (empty)
363 
364  /// Apply preconditioner to r
366 
367  /// \short Setup the preconditioner (store diagonal) from the fully
368  /// assembled matrix. Problem pointer is ignored.
369  void setup();
370 
371 
372  private:
373 
374 
375  /// Row start for upper triangular matrix
377 
378  /// \short column entry for the upper triangular matrix (each element of the
379  /// vector contains the column index and coefficient)
381 
382  /// Row start for lower triangular matrix
384 
385  /// \short column entry for the lower triangular matrix (each element of the
386  /// vector contains the column index and coefficient)
388 };
389 
390 //=============================================================================
391 /// \short A preconditioner for performing inner iteration preconditioner
392 /// solves. The template argument SOLVER specifies the inner iteration
393 /// solver (which must be derived from IterativeLinearSolver) and the
394 /// template argument PRECONDITIONER specifies the preconditioner for the
395 /// inner iteration iterative solver.
396 /// Note: For no preconditioning use the IdentityPreconditioner.
397 //=============================================================================
398 template <class SOLVER, class PRECONDITIONER>
400 {
401  public:
402 
403  /// Constructor
405  {
406  // create the solver
407  Solver_pt = new SOLVER;
408 
409  // create the preconditioner
410  Preconditioner_pt = new PRECONDITIONER;
411 
412 #ifdef PARANOID
413  // paranoid check that the solver is an iterative solver and
414  // the preconditioner is a preconditioner
415  if (dynamic_cast<IterativeLinearSolver* >(Solver_pt) == 0)
416  {
417  throw OomphLibError(
418  "The template argument SOLVER must be of type IterativeLinearSolver",
419  OOMPH_CURRENT_FUNCTION,
420  OOMPH_EXCEPTION_LOCATION);
421  }
422  if (dynamic_cast<Preconditioner*>(Preconditioner_pt) == 0)
423  {
424  throw OomphLibError(
425  "The template argument PRECONDITIONER must be of type Preconditioner",
426  OOMPH_CURRENT_FUNCTION,
427  OOMPH_EXCEPTION_LOCATION);
428  }
429 #endif
430 
431  // ensure the solver does not re-setup the preconditioner
432  Solver_pt->disable_setup_preconditioner_before_solve();
433 
434  // pass the preconditioner to the solver
435  Solver_pt->preconditioner_pt() = Preconditioner_pt;
436  }
437 
438  // destructor
440  {
441  delete Solver_pt;
442  delete Preconditioner_pt;
443  }
444 
445  // clean the memory
447  {
448  Preconditioner_pt->clean_up_memory();
449  Solver_pt->clean_up_memory();
450  }
451 
452  /// \short Preconditioner setup method. Setup the preconditioner for the inner
453  /// iteration solver.
454  void setup()
455  {
456 
457  // set the distribution
459  dynamic_cast<DistributableLinearAlgebraObject*>
460  (matrix_pt());
461  if (dist_pt != 0)
462  {
463  this->build_distribution(dist_pt->distribution_pt());
464  }
465  else
466  {
468  matrix_pt()->nrow(),false);
469  this->build_distribution(dist);
470  }
471 
472  // Setup the inner iteration preconditioner (For some reason we need to
473  // remind the compiler that there is also a function named setup in the
474  // base class.)
475  Preconditioner_pt->Preconditioner::setup(matrix_pt());
476 
477  // setup the solverready for resolve
478  unsigned max_iter = Solver_pt->max_iter();
479  Solver_pt->max_iter() = 1;
480  DoubleVector x(this->distribution_pt(),0.0);
481  DoubleVector y(x);
482  Solver_pt->enable_resolve();
483  Solver_pt->solve(matrix_pt(),x,y);
484  Solver_pt->max_iter() = max_iter;
485  }
486 
487  /// \short Preconditioner solve method. Performs the specified number
488  /// of Krylov iterations preconditioned with the specified preconditioner
490  {
491  Solver_pt->resolve(r,z);
492  }
493 
494  /// Access to convergence tolerance of the inner iteration solver
495  double& tolerance() {return Solver_pt->tolerance();}
496 
497  /// Access to max. number of iterations of the inner iteration solver
498  unsigned& max_iter() {return Solver_pt->max_iter();}
499 
500  ///
501  SOLVER* solver_pt() { return Solver_pt; }
502 
503  ///
504  PRECONDITIONER* preconditioner_pt() { return Preconditioner_pt; }
505 
506  private:
507 
508  /// pointer to the underlying solver
509  SOLVER* Solver_pt;
510 
511  /// pointer to the underlying preconditioner
512  PRECONDITIONER* Preconditioner_pt;
513 };
514 }
515 #endif
MatrixBasedLumpedPreconditioner(const MatrixBasedDiagPreconditioner &)
Broken copy constructor.
Vector< unsigned > L_column_start
Column start for lower triangular matrix.
ILUZeroPreconditioner(const ILUZeroPreconditioner &)
Broken copy constructor.
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
Vector< unsigned > U_row_start
Row start for upper triangular matrix.
virtual void preconditioner_solve(const DoubleVector &r, DoubleVector &z)=0
Apply the preconditioner. Pure virtual generic interface function. This method should apply the preco...
unsigned & index()
access function for the coefficient's (row or column) index
Vector< CompressedMatrixCoefficient > L_row_entry
column entry for the lower triangular matrix (each element of the vector contains the column index an...
bool positive_matrix() const
Access function to the Positive_matrix which indicates whether lumped matrix was positive.
unsigned Index
the row or column index of the compressed-matrix coefficient
double value() const
access function for the coefficient's value (const version)
Vector< CompressedMatrixCoefficient > L_row_entry
Row entry for the lower triangular matrix (each element of the vector contains the row index and coef...
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply preconditioner to z, i.e. z=D^-1.
double & tolerance()
Access to convergence tolerance of the inner iteration solver.
virtual DoubleMatrixBase * matrix_pt() const
Get function for matrix pointer.
Vector< CompressedMatrixCoefficient > U_row_entry
column entry for the upper triangular matrix (each element of the vector contains the column index an...
unsigned nrow() const
access function to the number of global rows.
void setup()
Setup the preconditioner (store diagonal) from the fully assembled matrix.
CompressedMatrixCoefficient(const CompressedMatrixCoefficient &a)
Copy Constructor. Not Broken. Required for STL sort function.
CompressedMatrixCoefficient(const unsigned &index, const double &value)
Constructor (takes the index and value as arguments)
Vector< unsigned > U_column_start
Column start for upper triangular matrix.
double & value()
access function for the coefficient value
Vector< double > Inv_diag
Vector of inverse diagonal entries.
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 &)
Broken assignment operator.
Vector< unsigned > L_row_start
Row start for lower triangular matrix.
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
ILUZeroPreconditioner(const ILUZeroPreconditioner &)
Broken copy constructor.
Matrix-based diagonal preconditioner.
double Value
the value of the compressed-matrix coefficient
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
Base class for any linear algebra object that is distributable. Just contains storage for the LinearA...
void setup()
Preconditioner setup method. Setup the preconditioner for the inner iteration solver.
unsigned & nrow()
Access function to number of rows for this preconditioner.
CompressedMatrixCoefficient()
Constructor (no arguments)
A preconditioner for performing inner iteration preconditioner solves. The template argument SOLVER s...
unsigned index() const
Access function for the coefficient's (row or column_ index (const version)
void setup()
Setup the preconditioner (store diagonal) from the fully assembled matrix. Problem pointer is ignored...
bool operator<(const CompressedMatrixCoefficient &a) const
Less Than Operator (for the STL sort function)
void operator=(const CompressedMatrixCoefficient &a)
Assignment Operator. Not Broken. Required for STL sort function.
SOLVER * Solver_pt
pointer to the underlying solver
double * Inv_lumped_diag_pt
Vector of inverse diagonal entries.
virtual const OomphCommunicator * comm_pt() const
Get function for comm pointer.
PRECONDITIONER * Preconditioner_pt
pointer to the underlying preconditioner
A class for compressed column matrices that store doubles.
Definition: matrices.h:2573
Preconditioner base class. Gives an interface to call all other preconditioners through and stores th...
unsigned & max_iter()
Access to max. number of iterations of the inner iteration solver.
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
void operator=(const ILUZeroPreconditioner &)
Broken assignment operator.
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Preconditioner solve method. Performs the specified number of Krylov iterations preconditioned with t...
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply preconditioner to z, i.e. z=D^-1.
void operator=(const MatrixBasedLumpedPreconditioner &)
Broken assignment operator.
void build_distribution(const LinearAlgebraDistribution *const dist_pt)
setup the distribution of this distributable linear algebra object
virtual void setup()=0
Setup the preconditioner. Pure virtual generic interface function.
MatrixBasedDiagPreconditioner(const MatrixBasedDiagPreconditioner &)
Broken copy constructor.
A vector in the mathematical sense, initially developed for linear algebra type applications. If MPI then this vector can be distributed - its distribution is described by the LinearAlgebraDistribution object at Distribution_pt. Data is stored in a C-style pointer vector (double*)
Definition: double_vector.h:61
Class for a compressed-matrix coefficent (for either CC or CR matrices). Contains the (row or column)...
~CompressedMatrixCoefficient()
Destructor (does nothing)
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:872
void operator=(const MatrixBasedDiagPreconditioner &)
Broken assignment operator.
void clean_up_memory()
Clean up memory (empty). Generic interface function.
void clean_up_memory()
clean up memory - just delete the inverse lumped vector
double * inverse_lumped_vector_pt()
Access function to the inverse of the lumped vector assembled in the preconditioner setup routine...