IFPACK Development
Loading...
Searching...
No Matches
Ifpack_Krylov.cpp
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"
61#include "EpetraExt_PointToBlockDiagPermute.h"
62#endif
63
64//==============================================================================
65// NOTE: any change to the default values should be committed to the other
66// constructor as well.
67Ifpack_Krylov::
68Ifpack_Krylov(Epetra_Operator* Operator) :
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?
109Ifpack_Krylov::
110Ifpack_Krylov(Epetra_RowMatrix* Operator) :
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
226 ++NumInitialize_;
227 InitializeTime_ += Time_->ElapsedTime();
228 IsInitialized_ = true;
229 return(0);
230}
231
232//==============================================================================
234{
235 if (!IsInitialized())
236 IFPACK_CHK_ERR(Initialize());
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);
285 IfpackPrec_ -> Initialize();
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;
325 if (ZeroStartingSolution_)
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//==============================================================================
358Condest(const Ifpack_CondestType CT,
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//==============================================================================
373void Ifpack_Krylov::SetLabel()
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.
417 ++NumApplyInverse_;
418 ApplyInverseTime_ += Time_->ElapsedTime();
419 return(0);
420}
int NumVectors() const
int PutScalar(double ScalarConstant)
Ifpack_BlockRelaxation: a class to define block relaxation preconditioners of Epetra_RowMatrix's.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
virtual int Compute()
Computes the preconditioners.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
virtual double Condest() const
Returns the condition number estimate, or -1.0 if not computed.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
virtual std::ostream & Print(std::ostream &os) const
Prints object to an output stream.
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.
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized, false otherwise.
virtual int Initialize()
Computes all it is necessary to initialize the preconditioner.
virtual const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
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.