mumps_solver.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 //This is the header file for the C++ wrapper functions for the
31 //mumps solver
32 
33 //Include guards to prevent multiple inclusions of the header
34 #ifndef NEW_MUMPS_SOLVER_HEADER
35 #define NEW_MUMPS_SOLVER_HEADER
36 
37 
38 // Config header generated by autoconfig
39 #ifdef HAVE_CONFIG_H
40  #include <oomph-lib-config.h>
41 #endif
42 
43 #ifdef OOMPH_HAS_MPI
44 #include "mpi.h"
45 #endif
46 
47 #include "linear_solver.h"
48 #include "preconditioner.h"
49 
50 #include <mumps_c_types.h>
51 #include <dmumps_c.h>
52 
53 namespace oomph
54 {
55 
56 
57 //====================================================================
58 /// \short Wrapper to Mumps solver
59 //====================================================================
60 class MumpsSolver : public LinearSolver
61 {
62 
63  public:
64 
65  /// \short Static flag that determines whether the warning about
66  /// incorrect distribution of RHSs will be printed or not
68 
69  /// \short Constructor: Call setup
70  MumpsSolver();
71 
72  /// Broken copy constructor
73  MumpsSolver(const MumpsSolver& dummy)
74  {
75  BrokenCopy::broken_copy("MumpsSolver");
76  }
77 
78  /// Broken assignment operator
79  void operator=(const MumpsSolver&)
80  {
81  BrokenCopy::broken_assign("MumpsSolver");
82  }
83 
84  ///Destructor: Cleanup
85  ~MumpsSolver();
86 
87  /// Overload disable resolve so that it cleans up memory too
89  {
92  }
93 
94  /// \short Set boolean to suppress warning about communicator
95  /// not equal to MPI_COMM_WORLD
98 
99  /// Don't suppress warning about communicator not equal to MPI_COMM_WORLD
102 
103  /// \short Set boolean to suppress info being printed to screen
104  /// during MUMPS solve
107 
108  /// Don't suppress info being printed to screen during MUMPS solve
111 
112  /// \short Solver: Takes pointer to problem and returns the results Vector
113  /// which contains the solution of the linear system defined by
114  /// the problem's fully assembled Jacobian and residual Vector.
115  void solve(Problem* const &problem_pt, DoubleVector &result);
116 
117  /// \short Linear-algebra-type solver: Takes pointer to a matrix and rhs
118  /// vector and returns the solution of the linear system.
119  /// The function returns the global result Vector.
120  /// Note: if Delete_matrix_data is true the function
121  /// matrix_pt->clean_up_memory() will be used to wipe the matrix data.
122  void solve(DoubleMatrixBase* const &matrix_pt,
123  const DoubleVector &rhs,
124  DoubleVector &result);
125 
126  /// \short Resolve the system defined by the last assembled Jacobian
127  /// and the specified rhs vector if resolve has been enabled.
128  /// Note: returns the global result Vector.
129  void resolve(const DoubleVector &rhs, DoubleVector &result);
130 
131  /// Enable documentation of statistics
132  void enable_doc_stats() {Doc_stats = true;}
133 
134  /// Disable documentation of statistics
135  void disable_doc_stats() {Doc_stats = false;}
136 
137  /// \short Returns the time taken to assemble the Jacobian matrix and
138  /// residual vector
140  {
141  return Jacobian_setup_time;
142  }
143 
144  /// \short Return the time taken to solve the linear system (needs to be
145  /// overloaded for each linear solver)
147  {
148  return Solution_time;
149  }
150 
151  /// \short Set the flag to avoid solution of the system, so
152  /// just assemble the Jacobian and the rhs.
153  /// (Used only for timing runs, obviously)
155 
156  /// \short Unset the flag so that the system is actually solved again
157  /// This is the default (obviously)
159 
160  /// \short Set Delete_matrix_data flag. MumpsSolver needs its own copy
161  /// of the input matrix, therefore a copy must be made if any matrix
162  /// used with this solver is to be preserved. If the input matrix can be
163  /// deleted the flag can be set to true to reduce the amount of memory
164  /// required, and the matrix data will be wiped using its clean_up_memory()
165  /// function. Default value is false.
167 
168  /// \short Unset Delete_matrix_data flag. MumpsSolver needs its own copy
169  /// of the input matrix, therefore a copy must be made if any matrix
170  /// used with this solver is to be preserved. If the input matrix can be
171  /// deleted the flag can be set to true to reduce the amount of memory
172  /// required, and the matrix data will be wiped using its clean_up_memory()
173  /// function. Default value is false.
175 
176  /// \short Do the factorisation stage
177  /// Note: if Delete_matrix_data is true the function
178  /// matrix_pt->clean_up_memory() will be used to wipe the matrix data.
179  void factorise(DoubleMatrixBase* const &matrix_pt);
180 
181  /// \short Do the backsubstitution for mumps solver
182  /// Note: returns the global result Vector.
183  void backsub(const DoubleVector &rhs,
184  DoubleVector &result);
185 
186  /// Clean up the memory allocated by the mumps solver
187  void clean_up_memory();
188 
189  /// \short Default factor for workspace -- static so it can be overwritten
190  /// globally.
192 
193  private:
194 
195  /// Initialise instance of mumps data structure
196  void initialise_mumps();
197 
198  /// Shutdown mumps
199  void shutdown_mumps();
200 
201  /// Jacobian setup time
203 
204  /// Solution time
206 
207  /// Suppress solve?
209 
210  /// Set to true for MumpsSolver to output statistics (false by default).
211  bool Doc_stats;
212 
213  /// Boolean to suppress warning about communicator not equal to MPI_COMM_WORLD
215 
216  /// Boolean to suppress info being printed to screen during MUMPS solve
218 
219  /// Has mumps been initialised?
221 
222  // Work space scaling factor
224 
225  /// \short Delete_matrix_data flag. MumpsSolver needs its own copy
226  /// of the input matrix, therefore a copy must be made if any matrix
227  /// used with this solver is to be preserved. If the input matrix can be
228  /// deleted the flag can be set to true to reduce the amount of memory
229  /// required, and the matrix data will be wiped using its clean_up_memory()
230  /// function. Default value is false.
232 
233  /// Vector for row numbers
235 
236  // Vector for column numbers
238 
239  // Vector for entries
241 
242  /// Pointer to MUMPS struct that contains the solver data
243  DMUMPS_STRUC_C* Mumps_struc_pt;
244 };
245 
246 
247 
248 
249 ///////////////////////////////////////////////////////////////////////
250 ///////////////////////////////////////////////////////////////////////
251 ///////////////////////////////////////////////////////////////////////
252 
253 
254 
255 //====================================================================
256 /// An interface to allow Mumps to be used as an (exact) Preconditioner
257 //====================================================================
259 {
260  public:
261 
262  /// Constructor.
264  {}
265 
266  /// Destructor.
268  {}
269 
270  /// Broken copy constructor.
272  {
273  BrokenCopy::broken_copy("NewMumpsPreconditioner");
274  }
275 
276 
277  /// Broken assignment operator.
279  {
280  BrokenCopy::broken_assign("NewMumpsPreconditioner");
281  }
282 
283  /// \short Function to set up a preconditioner for the linear
284  /// system defined by matrix_pt. This function must be called
285  /// before using preconditioner_solve.
286  /// Note: matrix_pt must point to an object of class
287  /// CRDoubleMatrix or CCDoubleMatrix
288  void setup()
289  {
290  oomph_info << "Setting up Mumps (exact) preconditioner"
291  << std::endl;
292 
293  DistributableLinearAlgebraObject* dist_matrix_pt=
295  if (dist_matrix_pt!=0)
296  {
297  LinearAlgebraDistribution dist(dist_matrix_pt->distribution_pt());
298  this->build_distribution(dist);
300  }
301  else
302  {
303  std::ostringstream error_message_stream;
304  error_message_stream
305  << "NewMumpsPreconditioner can only be applied to matrices derived \n"
306  << "DistributableLinearAlgebraObject.\n";
307  throw OomphLibError(error_message_stream.str(),
308  OOMPH_CURRENT_FUNCTION,
309  OOMPH_EXCEPTION_LOCATION);
310  }
311  }
312 
313  /// \short Function applies Mumps to vector r for (exact)
314  /// preconditioning, this requires a call to setup(...) first.
316  {
317  Solver.resolve(r, z);
318  }
319 
320 
321  /// \short Clean up memory -- forward the call to the version in
322  /// Mumps in its LinearSolver incarnation.
324  {
326  }
327 
328 
329  /// Enable documentation of timings
331 
332  /// Disable the documentation of timings
334 
335  private:
336 
337  /// \short the Mumps solver emplyed by this preconditioner
339 
340 };
341 
342 }
343 
344 #endif
345 
void setup()
Function to set up a preconditioner for the linear system defined by matrix_pt. This function must be...
Definition: mumps_solver.h:288
NewMumpsPreconditioner(const NewMumpsPreconditioner &)
Broken copy constructor.
Definition: mumps_solver.h:271
double Solution_time
Solution time.
Definition: mumps_solver.h:205
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
virtual double linear_solver_solution_time()
Return the time taken to solve the linear system (needs to be overloaded for each linear solver) ...
Definition: mumps_solver.h:146
unsigned Workspace_scaling_factor
Definition: mumps_solver.h:223
void disable_suppress_mumps_info_during_solve()
Don't suppress info being printed to screen during MUMPS solve.
Definition: mumps_solver.h:109
bool Doc_stats
Set to true for MumpsSolver to output statistics (false by default).
Definition: mumps_solver.h:211
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Function applies Mumps to vector r for (exact) preconditioning, this requires a call to setup(...
Definition: mumps_solver.h:315
void enable_doc_stats()
Enable documentation of statistics.
Definition: mumps_solver.h:132
static bool Suppress_incorrect_rhs_distribution_warning_in_resolve
Static flag that determines whether the warning about incorrect distribution of RHSs will be printed ...
Definition: mumps_solver.h:67
static int Default_workspace_scaling_factor
Default factor for workspace – static so it can be overwritten globally.
Definition: mumps_solver.h:191
virtual DoubleMatrixBase * matrix_pt() const
Get function for matrix pointer.
Wrapper to Mumps solver.
Definition: mumps_solver.h:60
void enable_doc_time()
Enable documentation of timings.
Definition: mumps_solver.h:330
void disable_resolve()
Overload disable resolve so that it cleans up memory too.
Definition: mumps_solver.h:88
The Problem class.
Definition: problem.h:152
void initialise_mumps()
Initialise instance of mumps data structure.
Definition: mumps_solver.cc:93
MumpsSolver Solver
the Mumps solver emplyed by this preconditioner
Definition: mumps_solver.h:338
Vector< double > A_loc
Definition: mumps_solver.h:240
void disable_suppress_solve()
Unset the flag so that the system is actually solved again This is the default (obviously) ...
Definition: mumps_solver.h:158
~NewMumpsPreconditioner()
Destructor.
Definition: mumps_solver.h:267
MumpsSolver(const MumpsSolver &dummy)
Broken copy constructor.
Definition: mumps_solver.h:73
void clean_up_memory()
Clean up memory – forward the call to the version in Mumps in its LinearSolver incarnation.
Definition: mumps_solver.h:323
Vector< int > Jcn_loc
Definition: mumps_solver.h:237
MumpsSolver()
Constructor: Call setup.
Definition: mumps_solver.cc:77
OomphInfo oomph_info
double Jacobian_setup_time
Jacobian setup time.
Definition: mumps_solver.h:202
void enable_delete_matrix_data()
Set Delete_matrix_data flag. MumpsSolver needs its own copy of the input matrix, therefore a copy mus...
Definition: mumps_solver.h:166
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
bool Delete_matrix_data
Delete_matrix_data flag. MumpsSolver needs its own copy of the input matrix, therefore a copy must be...
Definition: mumps_solver.h:231
void backsub(const DoubleVector &rhs, DoubleVector &result)
Do the backsubstitution for mumps solver Note: returns the global result Vector.
DMUMPS_STRUC_C * Mumps_struc_pt
Pointer to MUMPS struct that contains the solver data.
Definition: mumps_solver.h:243
Vector< int > Irn_loc
Vector for row numbers.
Definition: mumps_solver.h:234
bool Mumps_is_initialised
Has mumps been initialised?
Definition: mumps_solver.h:220
void resolve(const DoubleVector &rhs, DoubleVector &result)
Resolve the system defined by the last assembled Jacobian and the specified rhs vector if resolve has...
An interface to allow Mumps to be used as an (exact) Preconditioner.
Definition: mumps_solver.h:258
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 disable_delete_matrix_data()
Unset Delete_matrix_data flag. MumpsSolver needs its own copy of the input matrix, therefore a copy must be made if any matrix used with this solver is to be preserved. If the input matrix can be deleted the flag can be set to true to reduce the amount of memory required, and the matrix data will be wiped using its clean_up_memory() function. Default value is false.
Definition: mumps_solver.h:174
void factorise(DoubleMatrixBase *const &matrix_pt)
Do the factorisation stage Note: if Delete_matrix_data is true the function matrix_pt->clean_up_memor...
bool Suppress_mumps_info_during_solve
Boolean to suppress info being printed to screen during MUMPS solve.
Definition: mumps_solver.h:217
void disable_doc_time()
Disable documentation of solve times.
bool Suppress_warning_about_MPI_COMM_WORLD
Boolean to suppress warning about communicator not equal to MPI_COMM_WORLD.
Definition: mumps_solver.h:214
NewMumpsPreconditioner()
Constructor.
Definition: mumps_solver.h:263
void clean_up_memory()
Clean up the memory allocated by the mumps solver.
virtual void disable_resolve()
Disable resolve (i.e. store matrix and/or LU decomposition, say) This function simply resets an inter...
~MumpsSolver()
Destructor: Cleanup.
double jacobian_setup_time()
Returns the time taken to assemble the Jacobian matrix and residual vector.
Definition: mumps_solver.h:139
Preconditioner base class. Gives an interface to call all other preconditioners through and stores th...
void solve(Problem *const &problem_pt, DoubleVector &result)
Solver: Takes pointer to problem and returns the results Vector which contains the solution of the li...
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
void enable_suppress_solve()
Set the flag to avoid solution of the system, so just assemble the Jacobian and the rhs...
Definition: mumps_solver.h:154
void disable_suppress_warning_about_MPI_COMM_WORLD()
Don't suppress warning about communicator not equal to MPI_COMM_WORLD.
Definition: mumps_solver.h:100
void operator=(const MumpsSolver &)
Broken assignment operator.
Definition: mumps_solver.h:79
void disable_doc_stats()
Disable documentation of statistics.
Definition: mumps_solver.h:135
void build_distribution(const LinearAlgebraDistribution *const dist_pt)
setup the distribution of this distributable linear algebra object
void enable_suppress_warning_about_MPI_COMM_WORLD()
Set boolean to suppress warning about communicator not equal to MPI_COMM_WORLD.
Definition: mumps_solver.h:96
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
void enable_doc_time()
Enable documentation of solve times.
Abstract base class for matrices of doubles – adds abstract interfaces for solving, LU decomposition and multiplication by vectors.
Definition: matrices.h:275
void enable_suppress_mumps_info_during_solve()
Set boolean to suppress info being printed to screen during MUMPS solve.
Definition: mumps_solver.h:105
bool Suppress_solve
Suppress solve?
Definition: mumps_solver.h:208
void operator=(const NewMumpsPreconditioner &)
Broken assignment operator.
Definition: mumps_solver.h:278
void disable_doc_time()
Disable the documentation of timings.
Definition: mumps_solver.h:333
void shutdown_mumps()
Shutdown mumps.