43 #ifndef IFPACK2_KRYLOV_DEF_HPP 44 #define IFPACK2_KRYLOV_DEF_HPP 46 #include "Ifpack2_Chebyshev.hpp" 47 #include "Ifpack2_Heap.hpp" 48 #include "Ifpack2_ILUT.hpp" 49 #include "Ifpack2_Parameters.hpp" 50 #include "Ifpack2_Relaxation.hpp" 51 #include "Ifpack2_RILUK.hpp" 56 #include "Teuchos_Assert.hpp" 57 #include "Teuchos_Time.hpp" 66 template <
class MatrixType>
68 Krylov (
const Teuchos::RCP<const row_matrix_type>& A) :
71 iterationType_ (
"GMRES"),
75 ZeroStartingSolution_ (true),
76 PreconditionerType_ (1),
78 IsInitialized_ (false),
83 InitializeTime_ (0.0),
89 template <
class MatrixType>
93 template <
class MatrixType>
97 TEUCHOS_TEST_FOR_EXCEPTION(
98 ! A.is_null () && A->getComm ()->getSize () == 1 &&
99 A->getNodeNumRows () != A->getNodeNumCols (),
100 std::runtime_error,
"Ifpack2::Krylov::setMatrix: If A's communicator only " 101 "contains one process, then A must be square. Instead, you provided a " 102 "matrix A with " << A->getNodeNumRows () <<
" rows and " 103 << A->getNodeNumCols () <<
" columns.");
109 IsInitialized_ =
false;
116 template <
class MatrixType>
120 using Teuchos::ParameterList;
121 using Teuchos::Exceptions::InvalidParameterName;
122 using Teuchos::Exceptions::InvalidParameterType;
127 ParameterList params = plist;
135 int numIters = numIters_;
136 std::string iterType = iterationType_;
137 int blockSize = BlockSize_;
138 bool zeroStartingSolution = ZeroStartingSolution_;
139 int precType = PreconditionerType_;
147 bool gotIterType =
false;
149 iterType = params.get<std::string> (
"krylov: iteration type");
152 catch (InvalidParameterName) {
155 catch (InvalidParameterType) {
161 const int iterTypeInt = params.get<
int> (
"krylov: iteration type");
164 if (iterTypeInt == 1) {
166 }
else if (iterTypeInt == 2) {
169 TEUCHOS_TEST_FOR_EXCEPTION(
170 true, std::invalid_argument,
"Ifpack2::Krylov::setParameters: Invalid " 171 "\"krylov: iteration type\" value " << iterTypeInt <<
". Valid int " 172 "values are 1 (GMRES) and 2 (CG). Please prefer setting this " 173 "parameter as a string (the name of the Krylov solver to use).");
177 resTol = params.get (
"krylov: residual tolerance", resTol);
178 numIters = params.get (
"krylov: number of iterations", numIters);
179 blockSize = params.get (
"krylov: block size", blockSize);
180 zeroStartingSolution = params.get (
"krylov: zero starting solution",
181 zeroStartingSolution);
182 precType = params.get (
"krylov: preconditioner type", precType);
188 if (PreconditionerType_ == 1) {
189 precParams_.set (
"relaxation: sweeps",
190 params.get (
"relaxation: sweeps", 1));
191 precParams_.set (
"relaxation: damping factor",
192 params.get (
"relaxation: damping factor", (
scalar_type) 1.0));
193 precParams_.set (
"relaxation: min diagonal value",
194 params.get (
"relaxation: min diagonal value", STS::one ()));
195 precParams_.set (
"relaxation: zero starting solution",
196 params.get (
"relaxation: zero starting solution",
true));
197 precParams_.set (
"relaxation: backward mode",
198 params.get (
"relaxation: backward mode",
false));
203 if (PreconditionerType_ == 2 || PreconditionerType_ == 3) {
207 precParams_.set (
"fact: ilut level-of-fill",
208 params.get (
"fact: ilut level-of-fill", (
double) 1.0));
209 precParams_.set (
"fact: iluk level-of-fill",
210 params.get (
"fact: iluk level-of-fill", (
double) 1.0));
213 precParams_.set (
"fact: absolute threshold",
214 params.get (
"fact: absolute threshold", (
double) 0.0));
217 precParams_.set (
"fact: relative threshold",
218 params.get(
"fact: relative threshold", (
double) 1.0));
221 precParams_.set (
"fact: relax value",
222 params.get (
"fact: relax value", (
double) 0.0));
226 iterationType_ = iterType;
227 numIters_ = numIters;
229 BlockSize_ = blockSize;
230 ZeroStartingSolution_ = zeroStartingSolution;
231 PreconditionerType_ = precType;
235 template <
class MatrixType>
236 Teuchos::RCP<const Teuchos::Comm<int> >
238 TEUCHOS_TEST_FOR_EXCEPTION(
239 A_.is_null (), std::runtime_error,
"Ifpack2::Krylov::getComm: " 240 "The input matrix A is null. Please call setMatrix() with a nonnull " 241 "input matrix before calling this method.");
242 return A_->getComm ();
246 template <
class MatrixType>
247 Teuchos::RCP<const Tpetra::RowMatrix<typename MatrixType::scalar_type,typename MatrixType::local_ordinal_type,typename MatrixType::global_ordinal_type,typename MatrixType::node_type> >
253 template <
class MatrixType>
254 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,typename MatrixType::global_ordinal_type,typename MatrixType::node_type> >
257 TEUCHOS_TEST_FOR_EXCEPTION(
258 A_.is_null (), std::runtime_error,
"Ifpack2::Krylov::getDomainMap: " 259 "The input matrix A is null. Please call setMatrix() with a nonnull " 260 "input matrix before calling this method.");
261 return A_->getDomainMap ();
265 template <
class MatrixType>
266 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,typename MatrixType::global_ordinal_type,typename MatrixType::node_type> >
269 TEUCHOS_TEST_FOR_EXCEPTION(
270 A_.is_null (), std::runtime_error,
"Ifpack2::Krylov::getRangeMap: " 271 "The input matrix A is null. Please call setMatrix() with a nonnull " 272 "input matrix before calling this method.");
273 return A_->getRangeMap ();
277 template <
class MatrixType>
285 template <
class MatrixType>
287 return NumInitialize_;
291 template <
class MatrixType>
297 template <
class MatrixType>
303 template <
class MatrixType>
305 return InitializeTime_;
309 template <
class MatrixType>
315 template <
class MatrixType>
321 template <
class MatrixType>
324 using Teuchos::ParameterList;
332 TEUCHOS_TEST_FOR_EXCEPTION(
333 A_.is_null (), std::runtime_error,
"Ifpack2::Krylov::initialize: " 334 "The input matrix A is null. Please call setMatrix() with a nonnull " 335 "input matrix before calling this method.");
338 IsInitialized_ =
false;
341 Teuchos::Time timer (
"initialize");
343 Teuchos::TimeMonitor timeMon (timer);
346 RCP<ParameterList> belosList = rcp (
new ParameterList (
"GMRES"));
347 belosList->set (
"Maximum Iterations", numIters_);
348 belosList->set (
"Convergence Tolerance", resTol_);
355 if (PreconditionerType_ == 0) {
358 else if (PreconditionerType_==1) {
361 else if (PreconditionerType_==2) {
364 else if (PreconditionerType_==3) {
367 else if (PreconditionerType_==4) {
370 if (PreconditionerType_>0) {
371 ifpack2_prec_->initialize();
372 ifpack2_prec_->setParameters(precParams_);
375 belosProblem_->setOperator (A_);
377 if (iterationType_ ==
"GMRES") {
387 IsInitialized_ =
true;
389 InitializeTime_ += timer.totalElapsedTime ();
393 template <
class MatrixType>
396 TEUCHOS_TEST_FOR_EXCEPTION(
397 A_.is_null (), std::runtime_error,
"Ifpack2::Krylov::compute: " 398 "The input matrix A is null. Please call setMatrix() with a nonnull " 399 "input matrix before calling this method.");
406 Teuchos::Time timer (
"compute");
408 Teuchos::TimeMonitor timeMon (timer);
409 if (PreconditionerType_ > 0) {
410 ifpack2_prec_->compute ();
411 belosProblem_->setLeftPrec (ifpack2_prec_);
416 ComputeTime_ += timer.totalElapsedTime ();
420 template <
class MatrixType>
422 apply (
const Tpetra::MultiVector<
typename MatrixType::scalar_type,
423 typename MatrixType::local_ordinal_type,
424 typename MatrixType::global_ordinal_type,
425 typename MatrixType::node_type>& X,
426 Tpetra::MultiVector<
typename MatrixType::scalar_type,
427 typename MatrixType::local_ordinal_type,
428 typename MatrixType::global_ordinal_type,
429 typename MatrixType::node_type>& Y,
430 Teuchos::ETransp mode,
431 typename MatrixType::scalar_type alpha,
432 typename MatrixType::scalar_type beta)
const 436 using Teuchos::rcpFromRef;
439 TEUCHOS_TEST_FOR_EXCEPTION(
441 "Ifpack2::Krylov::apply: You must call compute() before you may call apply().");
442 TEUCHOS_TEST_FOR_EXCEPTION(
443 X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
444 "Ifpack2::Krylov::apply: The MultiVector inputs X and Y do not have the " 445 "same number of columns. X.getNumVectors() = " << X.getNumVectors ()
446 <<
" != Y.getNumVectors() = " << Y.getNumVectors () <<
".");
449 TEUCHOS_TEST_FOR_EXCEPTION(
450 alpha != STS::one (), std::logic_error,
451 "Ifpack2::Krylov::apply: alpha != 1 has not been implemented.");
452 TEUCHOS_TEST_FOR_EXCEPTION(
453 beta != STS::zero (), std::logic_error,
454 "Ifpack2::Krylov::apply: zero != 0 has not been implemented.");
455 TEUCHOS_TEST_FOR_EXCEPTION(
456 mode != Teuchos::NO_TRANS, std::logic_error,
457 "Ifpack2::Krylov::apply: mode != Teuchos::NO_TRANS has not been implemented.");
459 Teuchos::Time timer (
"apply");
461 Teuchos::TimeMonitor timeMon (timer);
467 auto X_lcl_host = X.template getLocalView<Kokkos::HostSpace> ();
468 auto Y_lcl_host = Y.template getLocalView<Kokkos::HostSpace> ();
469 if (X_lcl_host.ptr_on_device () == Y_lcl_host.ptr_on_device ()) {
470 Xcopy = rcp (
new MV (X, Teuchos::Copy));
472 Xcopy = rcpFromRef (X);
476 RCP<MV> Ycopy = rcpFromRef (Y);
477 if (ZeroStartingSolution_) {
478 Ycopy->putScalar (STS::zero ());
482 belosProblem_->setProblem (Ycopy, Xcopy);
483 belosSolver_->solve ();
486 ApplyTime_ += timer.totalElapsedTime ();
490 template <
class MatrixType>
493 std::ostringstream os;
498 os <<
"\"Ifpack2::Krylov\": {";
499 if (this->getObjectLabel () !=
"") {
500 os <<
"Label: \"" << this->getObjectLabel () <<
"\", ";
502 os <<
"Initialized: " << (
isInitialized () ?
"true" :
"false") <<
", " 503 <<
"Computed: " << (
isComputed () ?
"true" :
"false") <<
", ";
506 os <<
"Matrix: null";
509 os <<
"Matrix: not null" 510 <<
", Global matrix dimensions: [" 511 << A_->getGlobalNumRows () <<
", " << A_->getGlobalNumCols () <<
"]";
519 template <
class MatrixType>
522 const Teuchos::EVerbosityLevel verbLevel)
const 526 using Teuchos::VERB_DEFAULT;
527 using Teuchos::VERB_NONE;
528 using Teuchos::VERB_LOW;
529 using Teuchos::VERB_MEDIUM;
530 using Teuchos::VERB_HIGH;
531 using Teuchos::VERB_EXTREME;
533 const Teuchos::EVerbosityLevel vl =
534 (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
536 if (vl != VERB_NONE) {
538 Teuchos::OSTab tab0 (out);
539 out <<
"\"Ifpack2::Krylov\":";
541 Teuchos::OSTab tab1 (out);
542 if (this->getObjectLabel () !=
"") {
543 out <<
"Label: " << this->getObjectLabel () << endl;
545 out <<
"Initialized: " << (
isInitialized () ?
"true" :
"false") << endl
546 <<
"Computed: " << (
isComputed () ?
"true" :
"false") << endl
547 <<
"Global number of rows: " << A_->getGlobalNumRows () << endl
548 <<
"Global number of columns: " << A_->getGlobalNumCols () << endl
551 out <<
" null" << endl;
553 A_->describe (out, vl);
569 #define IFPACK2_KRYLOV_INSTANT(S,LO,GO,N) \ 570 template class Ifpack2::Krylov< Tpetra::RowMatrix<S, LO, GO, N> >; double getApplyTime() const
Returns the time spent in apply().
Definition: Ifpack2_Krylov_def.hpp:316
void initialize()
Do any initialization that depends on the input matrix's structure.
Definition: Ifpack2_Krylov_def.hpp:322
int getNumApply() const
Returns the number of calls to apply().
Definition: Ifpack2_Krylov_def.hpp:298
double getInitializeTime() const
Returns the time spent in Initialize().
Definition: Ifpack2_Krylov_def.hpp:304
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_Krylov_decl.hpp:110
Teuchos::RCP< const Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > > getDomainMap() const
Tpetra::Map representing the domain of this operator.
Definition: Ifpack2_Krylov_def.hpp:255
void apply(const Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &X, Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, scalar_type alpha=Teuchos::ScalarTraits< scalar_type >::one(), scalar_type beta=Teuchos::ScalarTraits< scalar_type >::zero()) const
Apply the preconditioner to X, putting the result in Y.
Definition: Ifpack2_Krylov_def.hpp:422
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Change the matrix to be preconditioned.
Definition: Ifpack2_Krylov_def.hpp:94
int getNumCompute() const
Returns the number of calls to Compute().
Definition: Ifpack2_Krylov_def.hpp:292
Teuchos::ScalarTraits< scalar_type >::magnitudeType magnitude_type
The type of the magnitude (absolute value) of a matrix entry.
Definition: Ifpack2_Krylov_decl.hpp:125
Krylov(const Teuchos::RCP< const row_matrix_type > &A)
Constructor that takes a Tpetra::RowMatrix.
Definition: Ifpack2_Krylov_def.hpp:68
std::string description() const
Return a simple one-line description of this object.
Definition: Ifpack2_Krylov_def.hpp:491
ILU(k) factorization of a given Tpetra::RowMatrix.
Definition: Ifpack2_RILUK_decl.hpp:243
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input MatrixType.
Definition: Ifpack2_Krylov_decl.hpp:119
MatrixType::node_type node_type
The Node type used by the input MatrixType.
Definition: Ifpack2_Krylov_decl.hpp:122
ILUT (incomplete LU factorization with threshold) of a Tpetra sparse matrix.
Definition: Ifpack2_ILUT_decl.hpp:91
virtual ~Krylov()
Destructor.
Definition: Ifpack2_Krylov_def.hpp:90
Teuchos::RCP< const Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > > getRangeMap() const
Tpetra::Map representing the range of this operator.
Definition: Ifpack2_Krylov_def.hpp:267
Teuchos::RCP< const Tpetra::RowMatrix< scalar_type, local_ordinal_type, global_ordinal_type, node_type > > getMatrix() const
Returns a reference to the matrix to be preconditioned.
Definition: Ifpack2_Krylov_def.hpp:248
void compute()
Do any initialization that depends on the input matrix's values.
Definition: Ifpack2_Krylov_def.hpp:394
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_Krylov_decl.hpp:116
int getNumInitialize() const
Returns the number of calls to Initialize().
Definition: Ifpack2_Krylov_def.hpp:286
Diagonally scaled Chebyshev iteration for Tpetra sparse matrices.
Definition: Ifpack2_Chebyshev_decl.hpp:199
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Return the operator's communicator.
Definition: Ifpack2_Krylov_def.hpp:237
bool hasTransposeApply() const
Whether this object's apply() method can apply the transpose (or conjugate transpose, if applicable).
Definition: Ifpack2_Krylov_def.hpp:278
bool isComputed() const
Return true if compute() completed successfully, else false.
Definition: Ifpack2_Krylov_decl.hpp:183
Relaxation preconditioners for Tpetra::RowMatrix and Tpetra::CrsMatrix sparse matrices.
Definition: Ifpack2_Relaxation_decl.hpp:222
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
void setParameters(const Teuchos::ParameterList ¶ms)
Set the preconditioner's parameters.
Definition: Ifpack2_Krylov_def.hpp:117
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object with some verbosity level to an FancyOStream object.
Definition: Ifpack2_Krylov_def.hpp:521
bool isInitialized() const
Return true if initialize() completed successfully, else false.
Definition: Ifpack2_Krylov_decl.hpp:175
double getComputeTime() const
Returns the time spent in Compute().
Definition: Ifpack2_Krylov_def.hpp:310