trilinos_preconditioners.cc
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//====================================================================
31 
32 namespace oomph
33 {
34 
35 
36 ///////////////////////////////////////////////////////////////////////////////
37 ///////////////////////////////////////////////////////////////////////////////
38 // functions for the TrilinosPreconditionerBase class
39 ///////////////////////////////////////////////////////////////////////////////
40 ///////////////////////////////////////////////////////////////////////////////
41 
42 
43 //=============================================================================
44 /// \short Static double that accumulates the preconditioner
45 /// solve time of all instantiations of this class. Reset
46 /// it manually, e.g. after every Newton solve.
47 //=============================================================================
49 
50 
51 //=============================================================================
52 /// \short Function to set up a preconditioner for the linear system
53 /// defined by matrix_pt. This function must be called before using
54 /// preconditioner_solve.
55 /// \b NOTE 1. matrix_pt must point to an object of class CRDoubleMatrix or
56 /// DistributedCRDoubleMatrix
57 /// This method should be called by oomph-lib solvers and preconditioners
58 //=============================================================================
60 {
61  //clean up the memory
63 
64 #ifdef PARANOID
65  // check the matrix is square
66  if ( matrix_pt()->nrow() != matrix_pt()->ncol() )
67  {
68  std::ostringstream error_message;
69  error_message << "Preconditioners require a square matrix. "
70  << "Matrix is " << matrix_pt()->nrow()
71  << " by " << matrix_pt()->ncol() << std::endl;
72  throw OomphLibError(error_message.str(),
73  OOMPH_CURRENT_FUNCTION,
74  OOMPH_EXCEPTION_LOCATION);
75  }
76 #endif
77 
78 
79  // get a pointer to the cr double matrix
80  CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt());
81 
82 #ifdef PARANOID
83  if (cr_matrix_pt == 0)
84  {
85  throw OomphLibError("Matrix must be of type CRDoubleMatrix",
86  OOMPH_CURRENT_FUNCTION,
87  OOMPH_EXCEPTION_LOCATION);
88  }
89 #endif
90 
91  // build the distribution of this preconditioner
92  this->build_distribution(cr_matrix_pt->distribution_pt());
93 
94  // create the matrices
96  (cr_matrix_pt,this->distribution_pt());
97 
98  // set up preconditioner
99  setup_trilinos_preconditioner(Epetra_matrix_pt);
100 }
101 
102 //===========================================================================
103 /// \short Function to setup a preconditioner for the linear system defined
104 /// by the oomph-lib oomph_matrix_pt and Epetra epetra_matrix_pt matrices.
105 /// This method is called by Trilinos solvers.
106 //===========================================================================
107 void TrilinosPreconditionerBase::setup(Epetra_CrsMatrix* epetra_matrix_pt)
108 {
109 // clean up old data
110  clean_up_memory();
111 
112  // first try CRDoubleMatrix
113  CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt());
114  if (cr_matrix_pt==0)
115  {
116  std::ostringstream error_message;
117  error_message << "TrilinosSolver only work with "
118  << "DistributedCRDoubleMatrix matrices" << std::endl;
119  throw OomphLibError(error_message.str(),
120  OOMPH_CURRENT_FUNCTION,
121  OOMPH_EXCEPTION_LOCATION);
122  }
123 
124  // setup the specific preconditioner
125  setup_trilinos_preconditioner(epetra_matrix_pt);
126 }
127 
128 
129 //=============================================================================
130 /// \short preconditioner_solve - applies the preconditioner to the vector r
131 /// taking distributed oomph-lib vectors (DistributedVector<double>) as
132 /// arguments.
133 //=============================================================================
136 {
137 
138  // Get ready to do cumulative timings
139 double t_start = TimingHelpers::timer();
140 
141 #ifdef PARANOID
142  // check preconditioner data exists
143  if (Epetra_preconditioner_pt == 0)
144  {
145  std::ostringstream error_message;
146  error_message << "preconditioner_solve requires that solver data has "
147  << "been set up" << std::endl;
148  throw OomphLibError(error_message.str(),
149  OOMPH_CURRENT_FUNCTION,
150  OOMPH_EXCEPTION_LOCATION);
151  }
152 
153  if (*this->distribution_pt() != r.distribution_pt())
154  {
155  std::ostringstream error_message;
156  error_message << "The rhs vector and the matrix must have the same "
157  << "distribution.\n";
158  throw OomphLibError(error_message.str(),
159  OOMPH_CURRENT_FUNCTION,
160  OOMPH_EXCEPTION_LOCATION);
161  }
162 #endif
163 
164  // create Epetra version of r
165  Epetra_Vector* epetra_r_pt =
167 
168  // create an empty Epetra vector for z
169  Epetra_Vector* epetra_z_pt =
171  (this->distribution_pt());
172 
173  // Apply preconditioner
174  Epetra_preconditioner_pt->ApplyInverse(*epetra_r_pt,*epetra_z_pt);
175 
176  // Copy result to z
177  z.build(this->distribution_pt(),0.0);
179 
180  // clean up memory
181  delete epetra_r_pt;
182  delete epetra_z_pt;
183 
184 
185  // Add to cumulative solve time
186  double t_end = TimingHelpers::timer();
187  Cumulative_preconditioner_solve_time+=(t_end-t_start);
188 
189 }
190 
191 
192 ///////////////////////////////////////////////////////////////////////////////
193 ///////////////////////////////////////////////////////////////////////////////
194 // function for the TrilinosML class
195 ///////////////////////////////////////////////////////////////////////////////
196 ///////////////////////////////////////////////////////////////////////////////
197 
198 
199  //=============================================================================
200  /// (Static) default number of V cycles (one to be consistent
201  /// with previous default). (It's an int because Trilinos wants it to be!)
202  //=============================================================================
204 
205 //=============================================================================
206 // Function to set up the ML preconditioner.
207 //=============================================================================
209 setup_trilinos_preconditioner(Epetra_CrsMatrix* epetra_matrix_pt)
210 {
211 
212 
213  // doc setup time
214  oomph_info << "Setting up TrilinosML, ";
215  double t_start = TimingHelpers::timer();
216 
217 
218  // create the preconditioner
220  new ML_Epetra::MultiLevelPreconditioner(*epetra_matrix_pt,
222  true);
223 
224  // doc setup time
225  oomph_info << "time for setup [s] : "
226  << TimingHelpers::timer()-t_start
227  << std::endl;
228 
229  // oomph_info << "In here\n";
230  // ML_Epetra::MultiLevelPreconditioner* tmp_pt=0;
231  // tmp_pt=dynamic_cast<ML_Epetra::MultiLevelPreconditioner*>(
232  // Epetra_preconditioner_pt);
233  // if (tmp_pt!=0)
234  // {
235  // oomph_info << "Doing test\n";
236  // ML_parameters.print(*(oomph_info.stream_pt()));
237  // oomph_info.stream_pt()->flush();
238  // // tmp_pt->TestSmoothers();
239  // oomph_info << "Done test\n";
240  // }
241 
242 }
243 
244 
245 ///////////////////////////////////////////////////////////////////////////////
246 ///////////////////////////////////////////////////////////////////////////////
247 // function for the TrilinosIFPACK
248 ///////////////////////////////////////////////////////////////////////////////
249 ///////////////////////////////////////////////////////////////////////////////
250 
251 
252 //=============================================================================
253 // Function to set up the IFPACK preconditioner.
254 //=============================================================================
256 setup_trilinos_preconditioner(Epetra_CrsMatrix* epetra_matrix_pt)
257 {
258  // set up a parameter list
259  Teuchos::ParameterList ifpack_parameters;
260  ifpack_parameters.set("fact: level-of-fill", ILU_fill_level);
261  ifpack_parameters.set("fact: ilut level-of-fill", ILUT_fill_level);
262  ifpack_parameters.set("fact: absolute threshold", Absolute_threshold);
263  ifpack_parameters.set("fact: relative threshold", Relative_threshold);
264 
265  // create the preconditioner
266  Ifpack ifpack_factory;
267  Ifpack_Preconditioner* tmp_pt= ifpack_factory.Create(Preconditioner_type,
268  epetra_matrix_pt,
269  Overlap);
270 
271  tmp_pt->SetParameters(ifpack_parameters);
272  tmp_pt->Initialize();
273  tmp_pt->Compute();
274 
275  // Now store the pointer to the newly created Ifpack_Preconditioner
277 }
278 
279 
280 ///////////////////////////////////////////////////////////////////////////////
281 ///////////////////////////////////////////////////////////////////////////////
282 ///////////////////////////////////////////////////////////////////////////////
283 
284 
285 }
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...
void setup()
Broken assignment operator.
double Relative_threshold
Value of relative threshold, used to pertub diagonal.
virtual DoubleMatrixBase * matrix_pt() const
Get function for matrix pointer.
unsigned nrow() const
access function to the number of global rows.
static double Cumulative_preconditioner_solve_time
Static double that accumulates the preconditioner solve time of all instantiations of this class...
Epetra_CrsMatrix * Epetra_matrix_pt
Pointer used to store the epetra matrix - only used when this preconditioner is setup using the oomph...
OomphInfo oomph_info
Epetra_CrsMatrix * create_distributed_epetra_matrix(const CRDoubleMatrix *oomph_matrix_pt, const LinearAlgebraDistribution *dist_pt)
create an Epetra_CrsMatrix from an oomph-lib CRDoubleMatrix. If oomph_matrix_pt is NOT distributed (i...
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
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 build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
applies the preconditioner
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(...)
double timer()
returns the time in seconds after some point in past
void clean_up_memory()
deletes the preconditioner, matrices and maps
string Preconditioner_type
Type of ILU preconditioner.
double ILUT_fill_level
Level of fill for "ILUT".
virtual unsigned long nrow() const =0
Return the number of rows of the matrix.
Epetra_Vector * create_distributed_epetra_vector(const DoubleVector &oomph_vec)
create an Epetra_Vector from an oomph-lib DoubleVector. If oomph_vec is NOT distributed (i...
int Overlap
Value of overlap level - used in parallel ILU.
void build_distribution(const LinearAlgebraDistribution *const dist_pt)
setup the distribution of this distributable linear algebra object
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...
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 copy_to_oomphlib_vector(const Epetra_Vector *epetra_vec_pt, DoubleVector &oomph_vec)
Helper function to copy the contents of a Trilinos vector to an oomph-lib distributed vector...
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:872
virtual unsigned long ncol() const =0
Return the number of columns of the matrix.
virtual void setup_trilinos_preconditioner(Epetra_CrsMatrix *epetra_matrix_pt)=0
Function to set up a specific Trilinos preconditioner. This is called by setup(...).