43 #include "Ifpack_ConfigDefs.h"
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"
54 #include "Ifpack_Condest.h"
55 #ifdef HAVE_IFPACK_AZTECOO
59 #ifdef HAVE_IFPACK_EPETRAEXT
60 #include "Epetra_CrsMatrix.h"
61 #include "EpetraExt_PointToBlockDiagPermute.h"
69 IsInitialized_(false),
76 ApplyInverseTime_(0.0),
78 ApplyInverseFlops_(0.0),
82 PreconditionerType_(1),
85 DampingParameter_(1.0),
93 NumGlobalNonzeros_(0),
94 Operator_(Teuchos::rcp(Operator,false)),
96 ZeroStartingSolution_(true)
111 IsInitialized_(false),
116 InitializeTime_(0.0),
118 ApplyInverseTime_(0.0),
120 ApplyInverseFlops_(0.0),
124 PreconditionerType_(1),
127 DampingParameter_(1.0),
128 UseTranspose_(false),
135 NumGlobalNonzeros_(0),
136 Operator_(Teuchos::rcp(Operator,false)),
137 Matrix_(Teuchos::rcp(Operator,false)),
139 ZeroStartingSolution_(true)
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_);
161 return(Operator_->Comm());
167 return(Operator_->OperatorDomainMap());
173 return(Operator_->OperatorRangeMap());
183 if (X.NumVectors() != Y.NumVectors())
192 IFPACK_CHK_ERR(Operator_->Apply(X,Y));
201 IsInitialized_ =
false;
203 if (Operator_ == Teuchos::null)
206 if (Time_ == Teuchos::null)
211 if (
Matrix().NumGlobalRows64() !=
Matrix().NumGlobalCols64())
214 NumMyRows_ = Matrix_->NumMyRows();
215 NumMyNonzeros_ = Matrix_->NumMyNonzeros();
216 NumGlobalRows_ = Matrix_->NumGlobalRows64();
217 NumGlobalNonzeros_ = Matrix_->NumGlobalNonzeros64();
221 if (Operator_->OperatorDomainMap().NumGlobalElements64() !=
222 Operator_->OperatorRangeMap().NumGlobalElements64())
227 InitializeTime_ += Time_->ElapsedTime();
228 IsInitialized_ =
true;
238 Time_->ResetStartTime();
240 #ifdef HAVE_IFPACK_AZTECOO
242 AztecSolver_ = Teuchos::rcp(
new AztecOO );
243 if(IsRowMatrix_==
true) {
244 AztecSolver_ -> SetUserMatrix(&*Matrix_);
247 AztecSolver_ -> SetUserOperator(&*Operator_);
250 AztecSolver_ -> SetAztecOption(AZ_solver, AZ_cg);
253 AztecSolver_ -> SetAztecOption(AZ_solver, AZ_gmres);
255 AztecSolver_ -> SetAztecOption(AZ_output, AZ_none);
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"); }
270 if(IsRowMatrix_==
true) {
271 NumRows = Matrix_->NumMyRows();
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";
278 NumRows =
static_cast<int>(NumRows_LL);
280 List.set(
"partitioner: type",
"linear");
281 List.set(
"partitioner: local parts", NumRows/BlockSize_);
283 if(PreconditionerType_>0) {
287 AztecSolver_ -> SetPrecOperator(&*IfpackPrec_);
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;
303 ComputeTime_ += Time_->ElapsedTime();
314 if (!
Comm().MyPID()) {
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;
328 os <<
"Using input starting solution" << 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;
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;
348 os <<
" " << std::setw(15) << 0.0 << endl;
349 os <<
"================================================================================" << endl;
359 const int MaxIters,
const double Tol,
367 Condest_ = Ifpack_Condest(*
this, CT, MaxIters, Tol, Matrix_in);
373 void Ifpack_Krylov::SetLabel()
375 Label_ =
"IFPACK (Krylov smoother), iterations=" +
Ifpack_toString(Iterations_);
386 if (Iterations_ == 0)
389 int nVec = X.NumVectors();
390 if (nVec != Y.NumVectors())
393 Time_->ResetStartTime();
398 if(ZeroStartingSolution_==
true) {
402 #ifdef HAVE_IFPACK_AZTECOO
403 AztecSolver_ -> SetLHS(&Y);
404 AztecSolver_ -> SetRHS(&*Xcopy);
405 AztecSolver_ -> Iterate(Iterations_,Tolerance_);
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;
418 ApplyInverseTime_ += Time_->ElapsedTime();
Ifpack_BlockRelaxation: a class to define block relaxation preconditioners of Epetra_RowMatrix's.
virtual int Initialize()
Computes all it is necessary to initialize the preconditioner.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.
virtual double Condest() const
Returns the condition number estimate, or -1.0 if not computed.
virtual int Compute()
Computes the preconditioners.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
virtual std::ostream & Print(std::ostream &os) const
Prints object to an output stream.
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...
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized, false otherwise.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
virtual const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.