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// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
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_PRECONDITION_HEADER
32 #define OOMPH_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 "matrices.h"
41 
42 
43 
44 namespace oomph
45 {
46 
47  // Forward declaration of block preconditioners (for turn into subs.)
48  template<typename MATRIX> class BlockPreconditioner;
49 
50  // Forward declaration of problem class (for obsolete setup function)
51  class Problem;
52 
53  //=============================================================================
54  /// \short Preconditioner base class. Gives an interface to call all other
55  /// preconditioners through and stores the matrix and communicator
56  /// pointers. All preconditioners should be derived from this class.
57  //=============================================================================
59  {
60 
61  public:
62 
63  /// Constructor
65 
66  /// Broken copy constructor
68  {
69  BrokenCopy::broken_copy("Preconditioner");
70  }
71 
72  /// Broken assignment operator
74  {
75  BrokenCopy::broken_assign("Preconditioner");
76  }
77 
78  /// Destructor (empty)
79  virtual ~Preconditioner(){}
80 
81  /// \short Apply the preconditioner. Pure virtual generic interface
82  /// function. This method should apply the preconditioner operator to the
83  /// vector r and return the vector z.
84  virtual void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
85  = 0;
86 
87  /// \short Setup the preconditioner: store the matrix pointer and the
88  /// communicator pointer then call preconditioner specific setup()
89  /// function.
91  {
92  // Store matrix pointer
93  set_matrix_pt(matrix_pt);
94 
95  // Extract and store communicator pointer
98  if (dist_obj_pt!=0)
99  {
100  set_comm_pt(dist_obj_pt->distribution_pt()->communicator_pt());
101  }
102  else
103  {
104  set_comm_pt(0);
105  }
106 
107  double setup_time_start = TimingHelpers::timer();
108  setup();
109  double setup_time_finish = TimingHelpers::timer();
110  Setup_time = setup_time_finish - setup_time_start;
111  }
112 
113  /// \short Compatability layer for old preconditioners where problem
114  /// pointers were needed. The problem pointer is only used to get a
115  /// communicator pointer.
116  void setup(const Problem* problem_pt, DoubleMatrixBase* matrix_pt)
117  {
119  setup(matrix_pt);
120  }
121 
122  /// \short Setup the preconditioner. Pure virtual generic interface
123  /// function.
124  virtual void setup() = 0;
125 
126  /// \short Clean up memory (empty). Generic interface function.
127  virtual void clean_up_memory(){}
128 
129  /// Get function for matrix pointer.
130  virtual DoubleMatrixBase* matrix_pt() const
131  {
132 #ifdef PARANOID
133  if(Matrix_pt == 0)
134  {
135  std::ostringstream error_msg;
136  error_msg << "Matrix pointer is null.";
137  throw OomphLibError(error_msg.str(),
138  OOMPH_CURRENT_FUNCTION,
139  OOMPH_EXCEPTION_LOCATION);
140  }
141 #endif
142  return Matrix_pt;
143  }
144 
145  /// \short Set the matrix pointer.
147 
148  /// Get function for comm pointer.
149  virtual const OomphCommunicator* comm_pt() const
150  {
151  // If we are using MPI then the comm pointer should not be null.
152 #ifdef OOMPH_HAS_MPI
153 #ifdef PARANOID
154  if(Comm_pt == 0)
155  {
156  std::ostringstream error_msg;
157  error_msg
158  << "Tried to access a null communicator pointer. This might mean you are\n"
159  << "trying to use it in a non-parallel case. Or it might mean you haven't\n"
160  << "set it properly.";
161  throw OomphLibError(error_msg.str(),
162  OOMPH_CURRENT_FUNCTION,
163  OOMPH_EXCEPTION_LOCATION);
164  }
165 #endif
166 #endif
167  return Comm_pt;
168  }
169 
170  /// \short Set the communicator pointer
171  virtual void set_comm_pt(const OomphCommunicator* const comm_pt)
172  {Comm_pt = comm_pt;}
173 
174  /// \short Returns the time to setup the preconditioner.
175  double setup_time() const
176  {return Setup_time;}
177 
178  /// Virtual interface function for making a preconditioner a subsidiary
179  /// of a block preconditioner. By default nothing is needed, but if this
180  /// preconditioner is also a block preconditioner then things need to
181  /// happen. There's an assumption here that the block preconditioner will
182  /// be in CR form but since that assumption is hard coded all over
183  /// BlockPreconditioner we're safe.
185  (BlockPreconditioner<CRDoubleMatrix>* master_block_prec_pt,
186  const Vector<unsigned>& doftype_in_master_preconditioner_coarse) {}
187 
188  /// Virtual interface function for making a preconditioner a subsidiary
189  /// of a block preconditioner. By default nothing is needed, but if this
190  /// preconditioner is also a block preconditioner then things need to
191  /// happen. Version for coarsening dof-types.
193  (BlockPreconditioner<CRDoubleMatrix>* master_block_prec_pt,
194  const Vector<unsigned>& doftype_in_master_preconditioner_coarse,
195  const Vector<Vector<unsigned> > & doftype_coarsen_map_coarse) {}
196 
197  private:
198 
199  /// Storage for a pointer to the matrix.
201 
202  /// \short Storage for a pointer to the communicator. Null
203  /// if the preconditioner should not be distributed.
205 
206  /// The time it takes to set up this preconditioner.
207  double Setup_time;
208 
209  };//Preconditioner
210 
211 
212  //=============================================================================
213  /// The Identity Preconditioner
214  //=============================================================================
216  {
217 
218  public:
219 
220  // Constructor - nothing to construct
222 
223  /// Broken copy constructor
225  {
226  BrokenCopy::broken_copy("IdentityPreconditioner");
227  }
228 
229  /// Broken assignment operator
231  {
232  BrokenCopy::broken_assign("IdentityPreconditioner");
233  }
234 
235  /// Destructor (empty)
237 
238  /// setup method - just sets the distribution
239  virtual void setup()
240  {
241  // first attempt to cast to DistributableLinearAlgebraObject
242  DistributableLinearAlgebraObject* dist_matrix_pt =
244 
245  // if it is a distributable matrix
246  if (dist_matrix_pt != 0)
247  {
248  // store the distribution
249  this->build_distribution(dist_matrix_pt->distribution_pt());
250  }
251  // else it is not a distributable matrix
252  else
253  {
254  // # of rows in the matrix
255  unsigned n_row=matrix_pt()->nrow();
256 
257  LinearAlgebraDistribution dist(comm_pt(),n_row,false);
258  this->build_distribution(dist);
259  }
260  }
261 
262  /// \short Apply the preconditioner. This method should apply the
263  /// preconditioner operator to the vector r and return the vector z.
265  {
266 #ifdef PARANOID
267  if (*r.distribution_pt() != *this->distribution_pt())
268  {
269  std::ostringstream error_message_stream;
270  error_message_stream
271  << "The r vector must have the same distribution as the preconditioner. "
272  << "(this is the same as the matrix passed to setup())";
273  throw OomphLibError(error_message_stream.str(), OOMPH_CURRENT_FUNCTION,
274  OOMPH_EXCEPTION_LOCATION);
275  }
276  if (z.built())
277  {
278  if (*z.distribution_pt() != *this->distribution_pt())
279  {
280  std::ostringstream error_message_stream;
281  error_message_stream
282  << "The z vector distribution has been setup; it must have the "
283  << "same distribution as the r vector (and preconditioner).";
284  throw OomphLibError(error_message_stream.str(), OOMPH_CURRENT_FUNCTION,
285  OOMPH_EXCEPTION_LOCATION);
286  }
287  }
288 #endif
289 
290  // apply
291  z=r;
292  }
293  };
294 }
295 
296 #endif
virtual void set_matrix_pt(DoubleMatrixBase *matrix_pt)
Set the matrix pointer.
void setup(const Problem *problem_pt, DoubleMatrixBase *matrix_pt)
Compatability layer for old preconditioners where problem pointers were needed. The problem pointer i...
bool built() const
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
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...
The Identity Preconditioner.
virtual void preconditioner_solve(const DoubleVector &r, DoubleVector &z)=0
Apply the preconditioner. Pure virtual generic interface function. This method should apply the preco...
const OomphCommunicator * Comm_pt
Storage for a pointer to the communicator. Null if the preconditioner should not be distributed...
DoubleMatrixBase * Matrix_pt
Storage for a pointer to the matrix.
virtual DoubleMatrixBase * matrix_pt() const
Get function for matrix pointer.
void operator=(const Preconditioner &)
Broken assignment operator.
IdentityPreconditioner(const IdentityPreconditioner &)
Broken copy constructor.
void obsolete()
Output warning message.
The Problem class.
Definition: problem.h:152
virtual void turn_into_subsidiary_block_preconditioner(BlockPreconditioner< CRDoubleMatrix > *master_block_prec_pt, const Vector< unsigned > &doftype_in_master_preconditioner_coarse)
virtual ~IdentityPreconditioner()
Destructor (empty)
void operator=(const IdentityPreconditioner &)
Broken assignment operator.
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
double Setup_time
The time it takes to set up this preconditioner.
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...
virtual ~Preconditioner()
Destructor (empty)
virtual void setup()
setup method - just sets the distribution
Preconditioner()
Constructor.
double timer()
returns the time in seconds after some point in past
void setup(DoubleMatrixBase *matrix_pt)
Setup the preconditioner: store the matrix pointer and the communicator pointer then call preconditio...
virtual const OomphCommunicator * comm_pt() const
Get function for comm pointer.
virtual void set_comm_pt(const OomphCommunicator *const comm_pt)
Set the communicator pointer.
Preconditioner base class. Gives an interface to call all other preconditioners through and stores th...
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
virtual unsigned long nrow() const =0
Return the number of rows of the matrix.
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.
OomphCommunicator * communicator_pt() const
const access to the communicator pointer
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
virtual void clean_up_memory()
Clean up memory (empty). Generic interface function.
Abstract base class for matrices of doubles – adds abstract interfaces for solving, LU decomposition and multiplication by vectors.
Definition: matrices.h:275
double setup_time() const
Returns the time to setup the preconditioner.
Preconditioner(const Preconditioner &)
Broken copy constructor.
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
Definition: communicator.h:57