Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Ifpack_Krylov.cpp
Go to the documentation of this file.
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack: Object-Oriented Algebraic Preconditioner Package
5// Copyright (2002) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41*/
42
43#include "Ifpack_ConfigDefs.h"
44#include <iomanip>
45#include "Epetra_Operator.h"
46#include "Epetra_RowMatrix.h"
47#include "Epetra_Comm.h"
48#include "Epetra_Map.h"
49#include "Epetra_MultiVector.h"
50#include "Epetra_Vector.h"
51#include "Epetra_Time.h"
52#include "Ifpack_Krylov.h"
53#include "Ifpack_Utils.h"
54#include "Ifpack_Condest.h"
55#ifdef HAVE_IFPACK_AZTECOO
56#include "AztecOO.h"
57#endif
58
59#ifdef HAVE_IFPACK_EPETRAEXT
60#include "Epetra_CrsMatrix.h"
62#endif
63
64//==============================================================================
65// NOTE: any change to the default values should be committed to the other
66// constructor as well.
69 IsInitialized_(false),
70 IsComputed_(false),
71 NumInitialize_(0),
72 NumCompute_(0),
73 NumApplyInverse_(0),
74 InitializeTime_(0.0),
75 ComputeTime_(0.0),
76 ApplyInverseTime_(0.0),
77 ComputeFlops_(0.0),
78 ApplyInverseFlops_(0.0),
79 Iterations_(5),
80 Tolerance_(1e-6),
81 SolverType_(1),
82 PreconditionerType_(1),
83 NumSweeps_(0),
84 BlockSize_(1),
85 DampingParameter_(1.0),
86 UseTranspose_(false),
87 Condest_(-1.0),
88 /* ComputeCondest_(false), (unused; commented out to avoid build warnings) */
89 Label_(),
90 NumMyRows_(0),
91 NumMyNonzeros_(0),
92 NumGlobalRows_(0),
93 NumGlobalNonzeros_(0),
94 Operator_(Teuchos::rcp(Operator,false)),
95 IsRowMatrix_(false),
96 ZeroStartingSolution_(true)
97{
98}
99
100//==============================================================================
101// NOTE: This constructor has been introduced because SWIG does not appear
102// to appreciate dynamic_cast. An instruction of type
103// Matrix_ = dynamic_cast<const Epetra_RowMatrix*> in the
104// other construction does not work in PyTrilinos -- of course
105// it does in any C++ code (for an Epetra_Operator that is also
106// an Epetra_RowMatrix).
107//
108// FIXME: move declarations into a separate method?
111 IsInitialized_(false),
112 IsComputed_(false),
113 NumInitialize_(0),
114 NumCompute_(0),
115 NumApplyInverse_(0),
116 InitializeTime_(0.0),
117 ComputeTime_(0.0),
118 ApplyInverseTime_(0.0),
119 ComputeFlops_(0.0),
120 ApplyInverseFlops_(0.0),
121 Iterations_(5),
122 Tolerance_(1e-6),
123 SolverType_(1),
124 PreconditionerType_(1),
125 NumSweeps_(0),
126 BlockSize_(1),
127 DampingParameter_(1.0),
128 UseTranspose_(false),
129 Condest_(-1.0),
130 /* ComputeCondest_(false), (unused; commented out to avoid build warnings) */
131 Label_(),
132 NumMyRows_(0),
133 NumMyNonzeros_(0),
134 NumGlobalRows_(0),
135 NumGlobalNonzeros_(0),
136 Operator_(Teuchos::rcp(Operator,false)),
137 Matrix_(Teuchos::rcp(Operator,false)),
138 IsRowMatrix_(true),
139 ZeroStartingSolution_(true)
140{
141}
142
143//==============================================================================
144int Ifpack_Krylov::SetParameters(Teuchos::ParameterList& List)
145{
146 Iterations_ = List.get("krylov: iterations",Iterations_);
147 Tolerance_ = List.get("krylov: tolerance",Tolerance_);
148 SolverType_ = List.get("krylov: solver",SolverType_);
149 PreconditionerType_ = List.get("krylov: preconditioner",PreconditionerType_);
150 NumSweeps_ = List.get("krylov: number of sweeps",NumSweeps_);
151 BlockSize_ = List.get("krylov: block size",BlockSize_);
152 DampingParameter_ = List.get("krylov: damping parameter",DampingParameter_);
153 ZeroStartingSolution_ = List.get("krylov: zero starting solution",ZeroStartingSolution_);
154 SetLabel();
155 return(0);
156}
157
158//==============================================================================
160{
161 return(Operator_->Comm());
162}
163
164//==============================================================================
166{
167 return(Operator_->OperatorDomainMap());
168}
169
170//==============================================================================
172{
173 return(Operator_->OperatorRangeMap());
174}
175
176//==============================================================================
179{
180 if (IsComputed() == false)
181 IFPACK_CHK_ERR(-3);
182
183 if (X.NumVectors() != Y.NumVectors())
184 IFPACK_CHK_ERR(-2);
185
186 if (IsRowMatrix_)
187 {
188 IFPACK_CHK_ERR(Matrix_->Multiply(UseTranspose(),X,Y));
189 }
190 else
191 {
192 IFPACK_CHK_ERR(Operator_->Apply(X,Y));
193 }
194
195 return(0);
196}
197
198//==============================================================================
200{
201 IsInitialized_ = false;
202
203 if (Operator_ == Teuchos::null)
204 IFPACK_CHK_ERR(-2);
205
206 if (Time_ == Teuchos::null)
207 Time_ = Teuchos::rcp( new Epetra_Time(Comm()) );
208
209 if (IsRowMatrix_)
210 {
211 if (Matrix().NumGlobalRows64() != Matrix().NumGlobalCols64())
212 IFPACK_CHK_ERR(-2); // only square matrices
213
214 NumMyRows_ = Matrix_->NumMyRows();
215 NumMyNonzeros_ = Matrix_->NumMyNonzeros();
216 NumGlobalRows_ = Matrix_->NumGlobalRows64();
217 NumGlobalNonzeros_ = Matrix_->NumGlobalNonzeros64();
218 }
219 else
220 {
221 if (Operator_->OperatorDomainMap().NumGlobalElements64() !=
222 Operator_->OperatorRangeMap().NumGlobalElements64())
223 IFPACK_CHK_ERR(-2); // only square operators
224 }
225
227 InitializeTime_ += Time_->ElapsedTime();
228 IsInitialized_ = true;
229 return(0);
230}
231
232//==============================================================================
234{
235 if (!IsInitialized())
237
238 Time_->ResetStartTime();
239
240#ifdef HAVE_IFPACK_AZTECOO
241 // setup Aztec solver
242 AztecSolver_ = Teuchos::rcp( new AztecOO );
243 if(IsRowMatrix_==true) {
244 AztecSolver_ -> SetUserMatrix(&*Matrix_);
245 }
246 else {
247 AztecSolver_ -> SetUserOperator(&*Operator_);
248 }
249 if(SolverType_==0) {
250 AztecSolver_ -> SetAztecOption(AZ_solver, AZ_cg);
251 }
252 else {
253 AztecSolver_ -> SetAztecOption(AZ_solver, AZ_gmres);
254 }
255 AztecSolver_ -> SetAztecOption(AZ_output, AZ_none);
256 // setup preconditioner
257 Teuchos::ParameterList List;
258 List.set("relaxation: damping factor", DampingParameter_);
259 List.set("relaxation: sweeps",NumSweeps_);
260 if(PreconditionerType_==0) { }
261 else if(PreconditionerType_==1) { List.set("relaxation: type", "Jacobi" ); }
262 else if(PreconditionerType_==2) { List.set("relaxation: type", "Gauss-Seidel" ); }
263 else if(PreconditionerType_==3) { List.set("relaxation: type", "symmetric Gauss-Seidel"); }
264 if(BlockSize_==1) {
265 IfpackPrec_ = Teuchos::rcp( new Ifpack_PointRelaxation(&*Matrix_) );
266 }
267 else {
268 IfpackPrec_ = Teuchos::rcp( new Ifpack_BlockRelaxation< Ifpack_DenseContainer > (&*Matrix_) );
269 int NumRows;
270 if(IsRowMatrix_==true) {
271 NumRows = Matrix_->NumMyRows();
272 }
273 else {
274 long long NumRows_LL = Operator_->OperatorDomainMap().NumGlobalElements64();
275 if(NumRows_LL > std::numeric_limits<int>::max())
276 throw "Ifpack_Krylov::Compute: NumGlobalElements don't fit an int";
277 else
278 NumRows = static_cast<int>(NumRows_LL);
279 }
280 List.set("partitioner: type", "linear");
281 List.set("partitioner: local parts", NumRows/BlockSize_);
282 }
283 if(PreconditionerType_>0) {
284 IfpackPrec_ -> SetParameters(List);
286 IfpackPrec_ -> Compute();
287 AztecSolver_ -> SetPrecOperator(&*IfpackPrec_);
288 }
289#else
290 using std::cout;
291 using std::endl;
292
293 cout << "You need to configure IFPACK with support for AztecOO" << endl;
294 cout << "to use this preconditioner. This may require --enable-aztecoo" << endl;
295 cout << "in your configure script." << endl;
296 IFPACK_CHK_ERR(-1);
297#endif
298
299 // reset values
300 IsComputed_ = false;
301 Condest_ = -1.0;
302 ++NumCompute_;
303 ComputeTime_ += Time_->ElapsedTime();
304 IsComputed_ = true;
305
306 return(0);
307}
308
309//==============================================================================
310std::ostream& Ifpack_Krylov::Print(std::ostream & os) const
311{
312 using std::endl;
313
314 if (!Comm().MyPID()) {
315 os << endl;
316 os << "================================================================================" << endl;
317 os << "Ifpack_Krylov" << endl;
318 os << "Number of iterations = " << Iterations_ << endl;
319 os << "Residual Tolerance = " << Tolerance_ << endl;
320 os << "Solver type (O for CG, 1 for GMRES) = " << SolverType_ << endl;
321 os << "Preconditioner type = " << PreconditionerType_ << endl;
322 os << "(0 for none, 1 for Jacobi, 2 for GS, 3 for SGS )" << endl;
323 os << "Condition number estimate = " << Condest() << endl;
324 os << "Global number of rows = " << Operator_->OperatorRangeMap().NumGlobalElements64() << endl;
326 os << "Using zero starting solution" << endl;
327 else
328 os << "Using input starting solution" << endl;
329 os << endl;
330 os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
331 os << "----- ------- -------------- ------------ --------" << endl;
332 os << "Initialize() " << std::setw(5) << NumInitialize_
333 << " " << std::setw(15) << InitializeTime_
334 << " 0.0 0.0" << endl;
335 os << "Compute() " << std::setw(5) << NumCompute_
336 << " " << std::setw(15) << ComputeTime_
337 << " " << std::setw(15) << 1.0e-6 * ComputeFlops_;
338 if (ComputeTime_ != 0.0)
339 os << " " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
340 else
341 os << " " << std::setw(15) << 0.0 << endl;
342 os << "ApplyInverse() " << std::setw(5) << NumApplyInverse_
343 << " " << std::setw(15) << ApplyInverseTime_
344 << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
345 if (ApplyInverseTime_ != 0.0)
346 os << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
347 else
348 os << " " << std::setw(15) << 0.0 << endl;
349 os << "================================================================================" << endl;
350 os << endl;
351 }
352
353 return(os);
354}
355
356//==============================================================================
359 const int MaxIters, const double Tol,
360 Epetra_RowMatrix* Matrix_in)
361{
362 if (!IsComputed()) // cannot compute right now
363 return(-1.0);
364
365 // always computes it. Call Condest() with no parameters to get
366 // the previous estimate.
367 Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
368
369 return(Condest_);
370}
371
372//==============================================================================
374{
375 Label_ = "IFPACK (Krylov smoother), iterations=" + Ifpack_toString(Iterations_);
376}
377
378//==============================================================================
381{
382
383 if (!IsComputed())
384 IFPACK_CHK_ERR(-3);
385
386 if (Iterations_ == 0)
387 return 0;
388
389 int nVec = X.NumVectors();
390 if (nVec != Y.NumVectors())
391 IFPACK_CHK_ERR(-2);
392
393 Time_->ResetStartTime();
394
395 // AztecOO gives X and Y pointing to the same memory location,
396 // need to create an auxiliary vector, Xcopy
397 Teuchos::RCP<Epetra_MultiVector> Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
398 if(ZeroStartingSolution_==true) {
399 Y.PutScalar(0.0);
400 }
401
402#ifdef HAVE_IFPACK_AZTECOO
403 AztecSolver_ -> SetLHS(&Y);
404 AztecSolver_ -> SetRHS(&*Xcopy);
405 AztecSolver_ -> Iterate(Iterations_,Tolerance_);
406#else
407 using std::cout;
408 using std::endl;
409
410 cout << "You need to configure IFPACK with support for AztecOO" << endl;
411 cout << "to use this preconditioner. This may require --enable-aztecoo" << endl;
412 cout << "in your configure script." << endl;
413 IFPACK_CHK_ERR(-1);
414#endif
415
416 // Flops are updated in each of the following.
418 ApplyInverseTime_ += Time_->ElapsedTime();
419 return(0);
420}
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
double Ifpack_Condest(const Ifpack_Preconditioner &IFP, const Ifpack_CondestType CT, const int MaxIters, const double Tol, Epetra_RowMatrix *Matrix)
#define IFPACK_CHK_ERR(ifpack_err)
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.
int NumVectors() const
Teuchos::RefCountPtr< Epetra_Operator > Operator_
Pointers to the matrix as an Epetra_Operator.
Teuchos::RefCountPtr< Epetra_RowMatrix > Matrix_
Pointers to the matrix as an Epetra_RowMatrix.
double ComputeFlops_
Contains the number of flops for Compute().
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
Teuchos::RefCountPtr< Epetra_Time > Time_
Time object to track timing.
Ifpack_Krylov(Epetra_Operator *Matrix)
double Tolerance_
Residual Tolerance.
virtual int Compute()
Computes the preconditioners.
bool IsRowMatrix_
If true, the Operator_ is an Epetra_RowMatrix.
int SolverType_
Solver - 0 for CG, 1 for GMRES.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
double ComputeTime_
Contains the time for all successful calls to Compute().
bool IsComputed_
If true, the preconditioner has been computed successfully.
virtual double Condest() const
Returns the condition number estimate, or -1.0 if not computed.
bool ZeroStartingSolution_
If true, the starting solution is always the zero vector.
double DampingParameter_
Damping parameter for inner preconditioner.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
int NumCompute_
Contains the number of successful call to Compute().
int Iterations_
Max number of iterations.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
long long NumGlobalRows_
Number of global rows.
virtual std::ostream & Print(std::ostream &os) const
Prints object to an output stream.
int NumMyRows_
Number of local rows.
int BlockSize_
Block Size (for block relaxation)
Teuchos::RCP< Ifpack_Preconditioner > IfpackPrec_
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.
int NumInitialize_
Contains the number of successful calls to Initialize().
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized, false otherwise.
double ApplyInverseFlops_
Contain sthe number of flops for ApplyInverse().
bool IsInitialized_
If true, the preconditioner has been computed successfully.
long long NumGlobalNonzeros_
Number of global nonzeros.
virtual int Initialize()
Computes all it is necessary to initialize the preconditioner.
double ApplyInverseTime_
Contains the time for all successful calls to ApplyInverse().
virtual const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.
double InitializeTime_
Contains the time for all successful calls to Initialize().
std::string Label_
Contains the label of this object.
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
int NumApplyInverse_
Contains the number of successful call to ApplyInverse().
int NumMyNonzeros_
Number of local nonzeros.
virtual void SetLabel()
Sets the label.
double Condest_
Contains the estimated condition number.
int NumSweeps_
Number of GS or Jacobi sweeps.
int PreconditionerType_
Preconditioner - 0 for none, 1 for Jacobi, 2 for GS, 3 for SGS.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
Ifpack_PointRelaxation: a class to define point relaxation preconditioners of for Epetra_RowMatrix's.
#define true
#define false