trilinos_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: 1182 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-05-20 16:50:20 +0100 (Fri, 20 May 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 #ifndef OOMPH_TRILINOS_OPERATORS_HEADER
31 #define OOMPH_TRILINOS_OPERATORS_HEADER
32 
33 // Config header generated by autoconfig
34 #ifdef HAVE_CONFIG_H
35  #include <oomph-lib-config.h>
36 #endif
37 
38 #ifdef OOMPH_HAS_MPI
39 #include "mpi.h"
40 #endif
41 
42 // Needed in trilinos (as of gcc 4.6.* or so...)
43 #include<cstddef>
44 
45 // trilinos headers
46 #include "ml_include.h"
47 #include "ml_MultiLevelPreconditioner.h"
48 #include "Ifpack.h"
49 
50 //oomph-lib headers
51 #include "trilinos_helpers.h"
52 #include "preconditioner.h"
53 
54 namespace oomph
55 {
56 
57 //=============================================================================
58 /// Base class for Trilinos preconditioners as oomph-lib preconditioner.
59 //=============================================================================
61  {
62  public:
63 
64  /// Constructor.
66  {
67  // Initialise pointers
70  }
71 
72  /// \short Destructor.
74  {
76  }
77 
78  /// \short Static double that accumulates the preconditioner
79  /// solve time of all instantiations of this class. Reset
80  /// it manually, e.g. after every Newton solve.
82 
83  /// deletes the preconditioner, matrices and maps
85  {
86  // delete the Epetra preconditioner
89 
90  // delete the epetra matrix
91  delete Epetra_matrix_pt;
92  Epetra_matrix_pt = 0;
93  }
94 
95  /// Broken copy constructor.
97  {
98  BrokenCopy::broken_copy("TrilinosPreconditionerBase");
99  }
100 
101 
102  /// Broken assignment operator.
103 //Commented out broken assignment operator because this can lead to a conflict warning
104 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
105 //realise that two separate implementations of the broken function are the same and so,
106 //quite rightly, it shouts.
107  /*void operator=(const TrilinosPreconditionerBase&)
108  {
109  BrokenCopy::broken_assign("TrilinosPreconditionerBase");
110  }*/
111 
112  /// \short Function to set up a preconditioner for the linear system
113  /// defined by matrix_pt. This function must be called before using
114  /// preconditioner_solve.
115  /// \b NOTE 1. matrix_pt must point to an object of class CRDoubleMatrix or
116  /// DistributedCRDoubleMatrix
117  /// This method should be called by oomph-lib solvers and preconditioners
118  void setup();
119 
120  /// \short Function to setup a preconditioner for the linear system defined
121  /// by the oomph-lib oomph_matrix_pt and Epetra epetra_matrix_pt matrices.
122  /// This method is called by Trilinos solvers.
123  void setup(Epetra_CrsMatrix* epetra_matrix_pt);
124 
125  /// \short applies the preconditioner
126  void preconditioner_solve(const DoubleVector& r,
127  DoubleVector &z);
128 
129  /// Access function to Epetra_preconditioner_pt.
130  /// For use with \c TrilinosAztecOOSolver
131  Epetra_Operator*& epetra_operator_pt()
132  {
134  }
135 
136  /// Access function to Epetra_preconditioner_pt (const version)
137  /// For use with \c TrilinosAztecOOSolver
138  Epetra_Operator* epetra_operator_pt() const
139  {
141  }
142 
143  protected:
144 
145  /// \short Function to set up a specific Trilinos preconditioner.
146  /// This is called by setup(...).
147  virtual void setup_trilinos_preconditioner
148  (Epetra_CrsMatrix* epetra_matrix_pt)=0;
149 
150  /// \short The preconditioner which will be set up using function
151  /// setup_trilinos_preconditioner(...)
152  Epetra_Operator* Epetra_preconditioner_pt;
153 
154  /// \short Pointer used to store the epetra matrix - only used when this
155  /// preconditioner is setup using the oomph-lib interface
156  Epetra_CrsMatrix* Epetra_matrix_pt;
157  };
158 
159 
160  //============================================================================
161  /// \short An interface to the Trilinos ML class - provides a function
162  /// to construct a serial ML object, and functions to modify some
163  /// of the ML paramaters.
164  //============================================================================
166  {
167  public:
168 
169  /// \short Constructor. Build with Smooth Aggretation (SA) default
170  /// settings, but our own default number of V cycles (initialised
171  /// to 1 to replicate TrilinosML's own behaviour).
173  {
174  // set default values
175  ML_Epetra::SetDefaults("SA", ML_parameters);
176 
177  // Set number of MG cycles performed in preconditioner
178  ML_parameters.set("cycle applications", Default_n_cycles);
179  }
180 
181  /// Destructor empty -- clean up is done in base class
183 
184  /// Broken copy constructor.
186  {
187  BrokenCopy::broken_copy("TrilinosMLPreconditioner");
188  }
189 
190  /// Broken assignment operator.
191  /*void operator=(const TrilinosMLPreconditioner&)
192  {
193  BrokenCopy::broken_assign("TrilinosMLPreconditioner");
194  }*/
195 
196  /// \short Set control flags to values for Petrov-Galerkin
197  /// preconditioning - for non symmetric systems
199  {
200  ML_Epetra::SetDefaults("NSSA", ML_parameters);
201  }
202 
203 
204  /// \short Set control flags to values for classical smoothed aggregation-
205  /// based 2-level domain decomposition
207  {
208  ML_Epetra::SetDefaults("DD", ML_parameters);
209  }
210 
211 
212  /// \short Set control flags to values 3-level algebraic domain decomposition
214  {
215  ML_Epetra::SetDefaults("DD-ML", ML_parameters);
216  }
217 
218  /// \short Set control flags to values for classical smoothed
219  /// aggregation preconditioning
221  {
222  ML_Epetra::SetDefaults("SA", ML_parameters);
223  }
224 
225  /// Function to set maximum number of levels
226  void set_max_levels(int max_levels)
227  {
228  ML_parameters.set("max levels", max_levels);
229  }
230 
231  /// Function to set the number of cycles used
232  void set_n_cycles(int n_cycles)
233  {
234  ML_parameters.set("cycle applications", n_cycles);
235  }
236 
237  /// Function to set Smoother_damping
238  void set_smoother_damping(double smoother_damping)
239  {
240  ML_parameters.set("smoother: damping factor", smoother_damping);
241  }
242 
243  /// Function to set Smoother_sweeps
244  void set_smoother_sweeps(int smoother_sweeps)
245  {
246  ML_parameters.set("smoother: sweeps", smoother_sweeps);
247  }
248 
249  /// Function to set smoother type to "Jacobi"
251  {
252  ML_parameters.set("smoother: type", "Jacobi");
253  }
254 
255  /// Function to set smoother type to "symmetric Gauss-Seidel"
257  {
258  ML_parameters.set("smoother: type", "symmetric Gauss-Seidel");
259  }
260 
261  /// Function to set output - controls level of information output by ML
262  void set_output(int output)
263  {
264  ML_parameters.set("output", output);
265  }
266 
267 
268  /// \short Default number of V cycles (one to be consistent with
269  /// previous default) (It's an int because Trilinos wants it to be!)
270  static int Default_n_cycles;
271 
272  protected:
273 
274  /// \short Function to set up the ML preconditioner. It is assumed
275  /// Trilinos_matrix_pt points to a suitable matrix.
276  void setup_trilinos_preconditioner(Epetra_CrsMatrix* epetra_matrix_pt);
277 
278  // Parameter list of control flags for the preconditioner
279  Teuchos::ParameterList ML_parameters;
280  };
281 
282 
283  //============================================================================
284  /// \short An interface to the Trilinos IFPACK class- provides a function
285  /// to construct an IFPACK object, and functions to modify some
286  /// of the IFPACK paramaters.
287  //============================================================================
289  {
290  public:
291 
292  /// Constructor.
294  {
295  // set default values
296  Preconditioner_type = "ILU";
297  ILU_fill_level = 0;
298  ILUT_fill_level = 1.0;
299  Overlap = 0;
300  Absolute_threshold = 0.0;
301  Relative_threshold = 1.0;
302  }
303 
304  /// Destructor -- empty, cleanup is done in base class
306 
307  /// Broken copy constructor.
309  {
310  BrokenCopy::broken_copy("TrilinosIFPACKPreconditioner");
311  }
312 
313 
314  /// Broken assignment operator.
315  /*void operator=(const TrilinosIFPACKPreconditioner&)
316  {
317  BrokenCopy::broken_assign("TrilinosIFPACKPreconditioner");
318  }*/
319 
320  /// Function to set Preconditioner_type to "ILU"
322  {
323  Preconditioner_type="ILU";
324  }
325 
326  /// Function to set Preconditioner_type to "ILUT"
328  {
329  Preconditioner_type="ILUT";
330  }
331 
332  /// Access function for ILU_fill_level
334 
335  /// Access function for ILUT_fill_level
336  double& ilut_fill_level() {return ILUT_fill_level;}
337 
338  /// Access function for the absolute threshold
340 
341  /// Access function for the relative threshold
343 
344  protected:
345 
346  /// \short Function to set up an IFPACK preconditioner. It is assumed
347  /// Trilinos_matrix_pt points to a suitable matrix.
348  void setup_trilinos_preconditioner(Epetra_CrsMatrix* epetra_matrix_pt);
349 
350  /// Type of ILU preconditioner
352 
353  /// Level of fill for "ILU"
355 
356  /// Level of fill for "ILUT"
358 
359  /// Value of overlap level - used in parallel ILU
360  int Overlap;
361 
362  /// Value of absolute threshold, used to peturb diagonal
364 
365  /// Value of relative threshold, used to pertub diagonal
367 
368  };
369 }
370 #endif
An interface to the Trilinos ML class - provides a function to construct a serial ML object...
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
void set_max_levels(int max_levels)
Function to set maximum number of levels.
void set_DD_default_values()
Set control flags to values for classical smoothed aggregation- based 2-level domain decomposition...
void set_smoother_damping(double smoother_damping)
Function to set Smoother_damping.
TrilinosMLPreconditioner(const TrilinosMLPreconditioner &)
Broken copy constructor.
void setup_trilinos_preconditioner(Epetra_CrsMatrix *epetra_matrix_pt)
Function to set up the ML preconditioner. It is assumed Trilinos_matrix_pt points to a suitable matri...
virtual ~TrilinosMLPreconditioner()
Destructor empty – clean up is done in base class.
void setup()
Broken assignment operator.
double Relative_threshold
Value of relative threshold, used to pertub diagonal.
virtual ~TrilinosPreconditionerBase()
Destructor.
void set_DDML_default_values()
Set control flags to values 3-level algebraic domain decomposition.
void set_preconditioner_ILUT()
Function to set Preconditioner_type to "ILUT".
Epetra_Operator * epetra_operator_pt() const
void set_smoother_jacobi()
Function to set smoother type to "Jacobi".
static double Cumulative_preconditioner_solve_time
Static double that accumulates the preconditioner solve time of all instantiations of this class...
void set_smoother_sweeps(int smoother_sweeps)
Function to set Smoother_sweeps.
Epetra_CrsMatrix * Epetra_matrix_pt
Pointer used to store the epetra matrix - only used when this preconditioner is setup using the oomph...
int & ilu_fill_level()
Access function for ILU_fill_level.
Base class for Trilinos preconditioners as oomph-lib preconditioner.
TrilinosMLPreconditioner()
Constructor. Build with Smooth Aggretation (SA) default settings, but our own default number of V cyc...
double Absolute_threshold
Value of absolute threshold, used to peturb diagonal.
static int Default_n_cycles
Default number of V cycles (one to be consistent with previous default) (It's an int because Trilinos...
void set_n_cycles(int n_cycles)
Function to set the number of cycles used.
TrilinosIFPACKPreconditioner(const TrilinosIFPACKPreconditioner &)
Broken copy constructor.
void set_smoother_gauss_seidel()
Function to set smoother type to "symmetric Gauss-Seidel".
double & absolute_threshold()
Access function for the absolute threshold.
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
applies the preconditioner
double & relative_threshold()
Access function for the relative threshold.
int ILU_fill_level
Level of fill for "ILU".
Epetra_Operator * Epetra_preconditioner_pt
The preconditioner which will be set up using function setup_trilinos_preconditioner(...)
void set_SA_default_values()
Set control flags to values for classical smoothed aggregation preconditioning.
void set_NSSA_default_values()
Broken assignment operator.
void output(std::ostream &outfile)
Output with default number of plot points.
void clean_up_memory()
deletes the preconditioner, matrices and maps
string Preconditioner_type
Type of ILU preconditioner.
Preconditioner base class. Gives an interface to call all other preconditioners through and stores th...
virtual ~TrilinosIFPACKPreconditioner()
Destructor – empty, cleanup is done in base class.
void set_preconditioner_ILU()
Broken assignment operator.
double ILUT_fill_level
Level of fill for "ILUT".
double & ilut_fill_level()
Access function for ILUT_fill_level.
int Overlap
Value of overlap level - used in parallel ILU.
void setup_trilinos_preconditioner(Epetra_CrsMatrix *epetra_matrix_pt)
Function to set up an IFPACK preconditioner. It is assumed Trilinos_matrix_pt points to a suitable ma...
TrilinosPreconditionerBase(const TrilinosPreconditionerBase &)
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
An interface to the Trilinos IFPACK class- provides a function to construct an IFPACK object...
void set_output(int output)
Function to set output - controls level of information output by ML.
virtual void setup_trilinos_preconditioner(Epetra_CrsMatrix *epetra_matrix_pt)=0
Function to set up a specific Trilinos preconditioner. This is called by setup(...).