preconditioner.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_PRECONDITION_HEADER
28#define OOMPH_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 "matrices.h"
37
38
39namespace oomph
40{
41 // Forward declaration of block preconditioners (for turn into subs.)
42 template<typename MATRIX>
43 class BlockPreconditioner;
44
45 // Forward declaration of problem class (for obsolete setup function)
46 class Problem;
47
48 //=============================================================================
49 /// Preconditioner base class. Gives an interface to call all other
50 /// preconditioners through and stores the matrix and communicator
51 /// pointers. All preconditioners should be derived from this class.
52 //=============================================================================
54 {
55 public:
56 /// Constructor
59 Stream_pt(0),
60 Matrix_pt(0),
61 Comm_pt(0),
62 Setup_time(0){};
63
64 /// Broken copy constructor
66
67 /// Broken assignment operator
68 void operator=(const Preconditioner&) = delete;
69
70 /// Destructor (empty)
71 virtual ~Preconditioner() {}
72
73 /// Apply the preconditioner. Pure virtual generic interface
74 /// function. This method should apply the preconditioner operator to the
75 /// vector r and return the vector z.
76 virtual void preconditioner_solve(const DoubleVector& r,
77 DoubleVector& z) = 0;
78
79 /// Apply the preconditioner. Pure virtual generic interface
80 /// function. This method should apply the preconditioner operator to the
81 /// vector r and return the vector z. (broken virtual)
83 DoubleVector& z)
84 {
85 // Throw an error
86 throw OomphLibError("This function hasn't been implemented yet!",
87 OOMPH_CURRENT_FUNCTION,
88 OOMPH_EXCEPTION_LOCATION);
89 }
90
91 /// Setup the preconditioner: store the matrix pointer and the
92 /// communicator pointer then call preconditioner specific setup()
93 /// function.
95 {
96 // Store matrix pointer
98
99 // Extract and store communicator pointer
102 if (dist_obj_pt != 0)
103 {
105 }
106 else
107 {
108 set_comm_pt(0);
109 }
110
111 double setup_time_start = TimingHelpers::timer();
112 setup();
113 double setup_time_finish = TimingHelpers::timer();
114 Setup_time = setup_time_finish - setup_time_start;
115 }
116
117 /// Compatability layer for old preconditioners where problem
118 /// pointers were needed. The problem pointer is only used to get a
119 /// communicator pointer.
120 void setup(const Problem* problem_pt, DoubleMatrixBase* matrix_pt)
121 {
124 }
125
126 /// Set up the block preconditioner quietly!
128 {
129 // Set the appropriate (silencing) boolean to true
131 } // End of enable_silent_preconditioner_setup
132
133
134 /// Be verbose in the block preconditioner setup
136 {
137 // Set the appropriate (silencing) boolean to false
139 } // End of disable_silent_preconditioner_setup
140
141
142 /// Setup the preconditioner. Pure virtual generic interface
143 /// function.
144 virtual void setup() = 0;
145
146 /// Clean up memory (empty). Generic interface function.
147 virtual void clean_up_memory() {}
148
149 /// Get function for matrix pointer.
151 {
152#ifdef PARANOID
153 if (Matrix_pt == 0)
154 {
155 std::ostringstream error_msg;
156 error_msg << "Matrix pointer is null.";
157 throw OomphLibError(
158 error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
159 }
160#endif
161 return Matrix_pt;
162 }
163
164 /// Set the matrix pointer.
166 {
168 }
169
170 /// Get function for comm pointer.
171 virtual const OomphCommunicator* comm_pt() const
172 {
173 // If we are using MPI then the comm pointer should not be null.
174#ifdef OOMPH_HAS_MPI
175#ifdef PARANOID
176 if (Comm_pt == 0)
177 {
178 std::ostringstream error_msg;
179 error_msg << "Tried to access a null communicator pointer. This might "
180 "mean you are\n"
181 << "trying to use it in a non-parallel case. Or it might "
182 "mean you haven't\n"
183 << "set it properly.";
184 throw OomphLibError(
185 error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
186 }
187#endif
188#endif
189 return Comm_pt;
190 }
191
192 /// Set the communicator pointer
193 virtual void set_comm_pt(const OomphCommunicator* const comm_pt)
194 {
196 }
197
198 /// Returns the time to setup the preconditioner.
199 double setup_time() const
200 {
201 return Setup_time;
202 }
203
204 /// Virtual interface function for making a preconditioner a subsidiary
205 /// of a block preconditioner. By default nothing is needed, but if this
206 /// preconditioner is also a block preconditioner then things need to
207 /// happen. There's an assumption here that the block preconditioner will
208 /// be in CR form but since that assumption is hard coded all over
209 /// BlockPreconditioner we're safe.
211 BlockPreconditioner<CRDoubleMatrix>* master_block_prec_pt,
212 const Vector<unsigned>& doftype_in_master_preconditioner_coarse)
213 {
214 }
215
216 /// Virtual interface function for making a preconditioner a subsidiary
217 /// of a block preconditioner. By default nothing is needed, but if this
218 /// preconditioner is also a block preconditioner then things need to
219 /// happen. Version for coarsening dof-types.
221 BlockPreconditioner<CRDoubleMatrix>* master_block_prec_pt,
222 const Vector<unsigned>& doftype_in_master_preconditioner_coarse,
223 const Vector<Vector<unsigned>>& doftype_coarsen_map_coarse)
224 {
225 }
226
227 protected:
228 /// Boolean to indicate whether or not the build should be done silently
230
231 /// Pointer to the output stream -- defaults to std::cout
232 std::ostream* Stream_pt;
233
234 private:
235 /// Storage for a pointer to the matrix.
237
238 /// Storage for a pointer to the communicator. Null
239 /// if the preconditioner should not be distributed.
241
242 /// The time it takes to set up this preconditioner.
244
245 }; // Preconditioner
246
247
248 //=============================================================================
249 /// The Identity Preconditioner
250 //=============================================================================
252 {
253 public:
254 // Constructor - nothing to construct
256
257 /// Broken copy constructor
259
260 /// Broken assignment operator
261 void operator=(const IdentityPreconditioner&) = delete;
262
263 /// Destructor (empty)
265
266 /// setup method - just sets the distribution
267 virtual void setup()
268 {
269 // first attempt to cast to DistributableLinearAlgebraObject
270 DistributableLinearAlgebraObject* dist_matrix_pt =
272
273 // if it is a distributable matrix
274 if (dist_matrix_pt != 0)
275 {
276 // store the distribution
277 this->build_distribution(dist_matrix_pt->distribution_pt());
278 }
279 // else it is not a distributable matrix
280 else
281 {
282 // # of rows in the matrix
283 unsigned n_row = matrix_pt()->nrow();
284
285 LinearAlgebraDistribution dist(comm_pt(), n_row, false);
286 this->build_distribution(dist);
287 }
288 }
289
290 /// Apply the preconditioner. This method should apply the
291 /// preconditioner operator to the vector r and return the vector z.
293 {
294#ifdef PARANOID
295 if (*r.distribution_pt() != *this->distribution_pt())
296 {
297 std::ostringstream error_message_stream;
298 error_message_stream
299 << "The r vector must have the same distribution as the "
300 "preconditioner. "
301 << "(this is the same as the matrix passed to setup())";
302 throw OomphLibError(error_message_stream.str(),
303 OOMPH_CURRENT_FUNCTION,
304 OOMPH_EXCEPTION_LOCATION);
305 }
306 if (z.built())
307 {
308 if (*z.distribution_pt() != *this->distribution_pt())
309 {
310 std::ostringstream error_message_stream;
311 error_message_stream
312 << "The z vector distribution has been setup; it must have the "
313 << "same distribution as the r vector (and preconditioner).";
314 throw OomphLibError(error_message_stream.str(),
315 OOMPH_CURRENT_FUNCTION,
316 OOMPH_EXCEPTION_LOCATION);
317 }
318 }
319#endif
320
321 // apply
322 z = r;
323 }
324
325
326 /// Apply the preconditioner. This method should apply the
327 /// preconditioner operator to the vector r and return the vector z.
329 {
330 // Applying the preconditioner to the transposed system is exactly the
331 // same as applying the preconditioner to the original system
333 }
334 };
335} // namespace oomph
336
337#endif
Base class for any linear algebra object that is distributable. Just contains storage for the LinearA...
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
Abstract base class for matrices of doubles – adds abstract interfaces for solving,...
Definition: matrices.h:261
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
bool built() const
The Identity Preconditioner.
void preconditioner_solve_transpose(const DoubleVector &r, DoubleVector &z)
Apply the preconditioner. This method should apply the preconditioner operator to the vector r and re...
IdentityPreconditioner(const IdentityPreconditioner &)=delete
Broken copy constructor.
virtual ~IdentityPreconditioner()
Destructor (empty)
void operator=(const IdentityPreconditioner &)=delete
Broken assignment operator.
virtual void setup()
setup method - just sets the distribution
virtual void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply the preconditioner. This method should apply the preconditioner operator to the vector r and re...
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
OomphCommunicator * communicator_pt() const
const access to the communicator pointer
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
Definition: communicator.h:54
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...
double Setup_time
The time it takes to set up this preconditioner.
DoubleMatrixBase * Matrix_pt
Storage for a pointer to the matrix.
virtual DoubleMatrixBase * matrix_pt() const
Get function for matrix pointer.
void setup(DoubleMatrixBase *matrix_pt)
Setup the preconditioner: store the matrix pointer and the communicator pointer then call preconditio...
virtual void turn_into_subsidiary_block_preconditioner(BlockPreconditioner< CRDoubleMatrix > *master_block_prec_pt, const Vector< unsigned > &doftype_in_master_preconditioner_coarse)
Virtual interface function for making a preconditioner a subsidiary of a block preconditioner....
virtual void preconditioner_solve_transpose(const DoubleVector &r, DoubleVector &z)
Apply the preconditioner. Pure virtual generic interface function. This method should apply the preco...
void disable_silent_preconditioner_setup()
Be verbose in the block preconditioner setup.
std::ostream * Stream_pt
Pointer to the output stream – defaults to std::cout.
virtual void set_comm_pt(const OomphCommunicator *const comm_pt)
Set the communicator pointer.
virtual void clean_up_memory()
Clean up memory (empty). Generic interface function.
virtual void turn_into_subsidiary_block_preconditioner(BlockPreconditioner< CRDoubleMatrix > *master_block_prec_pt, const Vector< unsigned > &doftype_in_master_preconditioner_coarse, const Vector< Vector< unsigned > > &doftype_coarsen_map_coarse)
Virtual interface function for making a preconditioner a subsidiary of a block preconditioner....
void operator=(const Preconditioner &)=delete
Broken assignment operator.
void setup(const Problem *problem_pt, DoubleMatrixBase *matrix_pt)
Compatability layer for old preconditioners where problem pointers were needed. The problem pointer i...
Preconditioner()
Constructor.
virtual const OomphCommunicator * comm_pt() const
Get function for comm pointer.
double setup_time() const
Returns the time to setup the preconditioner.
const OomphCommunicator * Comm_pt
Storage for a pointer to the communicator. Null if the preconditioner should not be distributed.
void enable_silent_preconditioner_setup()
Set up the block preconditioner quietly!
virtual ~Preconditioner()
Destructor (empty)
Preconditioner(const Preconditioner &)=delete
Broken copy constructor.
virtual void preconditioner_solve(const DoubleVector &r, DoubleVector &z)=0
Apply the preconditioner. Pure virtual generic interface function. This method should apply the preco...
bool Silent_preconditioner_setup
Boolean to indicate whether or not the build should be done silently.
virtual void set_matrix_pt(DoubleMatrixBase *matrix_pt)
Set the matrix pointer.
virtual void setup()=0
Setup the preconditioner. Pure virtual generic interface function.
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: problem.h:151
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
void obsolete()
Output warning message.
double timer()
returns the time in seconds after some point in past
//////////////////////////////////////////////////////////////////// ////////////////////////////////...