43 #ifndef IFPACK2_CHEBYSHEV_DEF_HPP 44 #define IFPACK2_CHEBYSHEV_DEF_HPP 46 #include "Ifpack2_Parameters.hpp" 47 #include "Teuchos_TimeMonitor.hpp" 48 #include "Tpetra_CrsMatrix.hpp" 55 template<
class MatrixType>
56 Chebyshev<MatrixType>::
57 Chebyshev (
const Teuchos::RCP<const row_matrix_type>& A)
59 IsInitialized_ (false),
64 InitializeTime_ (0.0),
70 this->setObjectLabel (
"Ifpack2::Chebyshev");
74 template<
class MatrixType>
79 template<
class MatrixType>
82 if (A.getRawPtr () != impl_.
getMatrix ().getRawPtr ()) {
83 IsInitialized_ =
false;
90 template<
class MatrixType>
95 impl_.
setParameters (const_cast<Teuchos::ParameterList&> (List));
99 template<
class MatrixType>
100 Teuchos::RCP<const Teuchos::Comm<int> >
103 Teuchos::RCP<const row_matrix_type> A = impl_.
getMatrix ();
104 TEUCHOS_TEST_FOR_EXCEPTION(
105 A.is_null (), std::runtime_error,
"Ifpack2::Chebyshev::getComm: The input " 106 "matrix A is null. Please call setMatrix() with a nonnull input matrix " 107 "before calling this method.");
108 return A->getRowMap ()->getComm ();
112 template<
class MatrixType>
113 Teuchos::RCP<const typename Chebyshev<MatrixType>::row_matrix_type>
120 template<
class MatrixType>
121 Teuchos::RCP<
const Tpetra::CrsMatrix<
typename MatrixType::scalar_type,
122 typename MatrixType::local_ordinal_type,
123 typename MatrixType::global_ordinal_type,
124 typename MatrixType::node_type> >
129 return Teuchos::rcp_dynamic_cast<
const crs_matrix_type> (impl_.
getMatrix ());
133 template<
class MatrixType>
134 Teuchos::RCP<const typename Chebyshev<MatrixType>::map_type>
138 Teuchos::RCP<const row_matrix_type> A = impl_.
getMatrix ();
139 TEUCHOS_TEST_FOR_EXCEPTION(
140 A.is_null (), std::runtime_error,
"Ifpack2::Chebyshev::getDomainMap: The " 141 "input matrix A is null. Please call setMatrix() with a nonnull input " 142 "matrix before calling this method.");
143 return A->getDomainMap ();
147 template<
class MatrixType>
148 Teuchos::RCP<const typename Chebyshev<MatrixType>::map_type>
152 Teuchos::RCP<const row_matrix_type> A = impl_.
getMatrix ();
153 TEUCHOS_TEST_FOR_EXCEPTION(
154 A.is_null (), std::runtime_error,
"Ifpack2::Chebyshev::getRangeMap: The " 155 "input matrix A is null. Please call setMatrix() with a nonnull input " 156 "matrix before calling this method.");
157 return A->getRangeMap ();
161 template<
class MatrixType>
167 template<
class MatrixType>
169 return NumInitialize_;
173 template<
class MatrixType>
179 template<
class MatrixType>
185 template<
class MatrixType>
187 return InitializeTime_;
191 template<
class MatrixType>
197 template<
class MatrixType>
203 template<
class MatrixType>
205 return ComputeFlops_;
209 template<
class MatrixType>
215 template<
class MatrixType>
218 apply (
const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
219 Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
220 Teuchos::ETransp mode,
224 const std::string timerName (
"Ifpack2::Chebyshev::apply");
225 Teuchos::RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
226 if (timer.is_null ()) {
227 timer = Teuchos::TimeMonitor::getNewCounter (timerName);
232 Teuchos::TimeMonitor timeMon (*timer);
236 TEUCHOS_TEST_FOR_EXCEPTION(
238 "Ifpack2::Chebyshev::apply(): You must call the compute() method before " 239 "you may call apply().");
240 TEUCHOS_TEST_FOR_EXCEPTION(
241 X.getNumVectors () != Y.getNumVectors (), std::runtime_error,
242 "Ifpack2::Chebyshev::apply(): X and Y must have the same number of " 243 "columns. X.getNumVectors() = " << X.getNumVectors() <<
" != " 244 <<
"Y.getNumVectors() = " << Y.getNumVectors() <<
".");
245 applyImpl (X, Y, mode, alpha, beta);
251 ApplyTime_ = timer->totalElapsedTime ();
255 template<
class MatrixType>
258 applyMat (
const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
259 Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
260 Teuchos::ETransp mode)
const 262 TEUCHOS_TEST_FOR_EXCEPTION(
263 X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
264 "Ifpack2::Chebyshev::applyMat: X.getNumVectors() != Y.getNumVectors().");
266 Teuchos::RCP<const row_matrix_type> A = impl_.
getMatrix ();
267 TEUCHOS_TEST_FOR_EXCEPTION(
268 A.is_null (), std::runtime_error,
"Ifpack2::Chebyshev::applyMat: The input " 269 "matrix A is null. Please call setMatrix() with a nonnull input matrix " 270 "before calling this method.");
272 A->apply (X, Y, mode);
276 template<
class MatrixType>
281 const std::string timerName (
"Ifpack2::Chebyshev::initialize");
282 Teuchos::RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
283 if (timer.is_null ()) {
284 timer = Teuchos::TimeMonitor::getNewCounter (timerName);
286 IsInitialized_ =
true;
291 template<
class MatrixType>
294 const std::string timerName (
"Ifpack2::Chebyshev::compute");
295 Teuchos::RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
296 if (timer.is_null ()) {
297 timer = Teuchos::TimeMonitor::getNewCounter (timerName);
302 Teuchos::TimeMonitor timeMon (*timer);
314 ComputeTime_ = timer->totalElapsedTime ();
318 template <
class MatrixType>
320 std::ostringstream out;
325 out <<
"\"Ifpack2::Chebyshev\": {";
326 out <<
"Initialized: " << (
isInitialized () ?
"true" :
"false") <<
", " 327 <<
"Computed: " << (
isComputed () ?
"true" :
"false") <<
", ";
332 out <<
"Matrix: null";
335 out <<
"Global matrix dimensions: [" 336 << impl_.
getMatrix ()->getGlobalNumRows () <<
", " 337 << impl_.
getMatrix ()->getGlobalNumCols () <<
"]" 338 <<
", Global nnz: " << impl_.
getMatrix ()->getGlobalNumEntries();
346 template <
class MatrixType>
349 const Teuchos::EVerbosityLevel verbLevel)
const 351 const Teuchos::EVerbosityLevel vl =
352 (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
353 const int myRank = this->
getComm ()->getRank ();
355 if (vl != Teuchos::VERB_NONE && myRank == 0) {
357 Teuchos::OSTab tab0 (out);
364 using Teuchos::VERB_DEFAULT;
365 using Teuchos::VERB_NONE;
366 using Teuchos::VERB_LOW;
367 using Teuchos::VERB_MEDIUM;
368 using Teuchos::VERB_HIGH;
369 using Teuchos::VERB_EXTREME;
373 Teuchos::EVerbosityLevel vl = verbLevel;
374 if (vl == VERB_DEFAULT) {
377 RCP<const Comm<int> > comm = A_->getRowMap ()->getComm ();
379 const int myImageID = comm->getRank();
380 Teuchos::OSTab tab(out);
384 Teuchos::ArrayRCP<const scalar_type> DiagView = InvDiagonal_->get1dView();
387 for(
typename Teuchos::ArrayRCP<scalar_type>::size_type i=1; i<DiagView.size(); ++i) {
388 if (STS::magnitude(myMinVal) > STS::magnitude(DiagView[i])) myMinVal = DiagView[i];
389 if (STS::magnitude(myMaxVal) < STS::magnitude(DiagView[i])) myMaxVal = DiagView[i];
391 Teuchos::reduceAll(*comm, Teuchos::REDUCE_MIN, 1, &myMinVal, &MinVal);
392 Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, 1, &myMaxVal, &MaxVal);
400 if (vl != VERB_NONE && myImageID == 0) {
403 out <<
"===============================================================================" << std::endl;
404 out <<
"Degree of polynomial = " << PolyDegree_ << std::endl;
405 if (ZeroStartingSolution_) { out <<
"Using zero starting solution" << endl; }
406 else { out <<
"Using input starting solution" << endl; }
408 out <<
"Minimum value on stored inverse diagonal = " << MinVal << std::endl;
409 out <<
"Maximum value on stored inverse diagonal = " << MaxVal << std::endl;
412 out <<
"Phase # calls Total Time (s) Total MFlops MFlops/s " << endl;
413 out <<
"------------ ------- --------------- --------------- ---------------" << endl;
421 out <<
"===============================================================================" << std::endl;
427 template<
class MatrixType>
432 Teuchos::ETransp mode,
436 using Teuchos::ArrayRCP;
440 using Teuchos::rcp_const_cast;
441 using Teuchos::rcpFromRef;
465 Y_orig = rcp (
new MV (Y, Teuchos::Copy));
474 RCP<const MV> X_copy;
475 bool copiedInput =
false;
477 auto X_lcl_host = X.template getLocalView<Kokkos::HostSpace> ();
478 auto Y_lcl_host = Y.template getLocalView<Kokkos::HostSpace> ();
479 if (X_lcl_host.ptr_on_device () == Y_lcl_host.ptr_on_device ()) {
480 X_copy = rcp (
new MV (X, Teuchos::Copy));
483 X_copy = rcpFromRef (X);
493 RCP<MV> X_copy_nonConst = rcp_const_cast<MV> (X_copy);
495 X_copy_nonConst = rcp (
new MV (X, Teuchos::Copy));
498 X_copy_nonConst->scale (alpha);
499 X_copy = rcp_const_cast<
const MV> (X_copy_nonConst);
502 impl_.
apply (*X_copy, Y);
505 Y.update (beta, *Y_orig, one);
510 template<
class MatrixType>
512 return impl_.getLambdaMaxForApply ();
519 #define IFPACK2_CHEBYSHEV_INSTANT(S,LO,GO,N) \ 520 template class Ifpack2::Chebyshev< Tpetra::RowMatrix<S, LO, GO, N> >; 522 #endif // IFPACK2_CHEBYSHEV_DEF_HPP void setParameters(Teuchos::ParameterList &plist)
Set (or reset) parameters.
Definition: Ifpack2_Details_Chebyshev_def.hpp:316
double getApplyFlops() const
The total number of floating-point operations taken by all calls to apply().
Definition: Ifpack2_Chebyshev_def.hpp:210
MT apply(const MV &B, MV &X)
Solve Ax=b for x with Chebyshev iteration with left diagonal scaling.
Definition: Ifpack2_Details_Chebyshev_def.hpp:825
Teuchos::RCP< const Tpetra::CrsMatrix< scalar_type, local_ordinal_type, global_ordinal_type, node_type > > getCrsMatrix() const
Attempt to return the matrix A as a Tpetra::CrsMatrix.
Definition: Ifpack2_Chebyshev_def.hpp:126
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input MatrixType.
Definition: Ifpack2_Chebyshev_decl.hpp:223
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_Chebyshev_decl.hpp:217
double getApplyTime() const
The total time spent in all calls to apply().
Definition: Ifpack2_Chebyshev_def.hpp:198
void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Set the matrix.
Definition: Ifpack2_Details_Chebyshev_def.hpp:652
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Change the matrix to be preconditioned.
Definition: Ifpack2_Chebyshev_def.hpp:80
void setParameters(const Teuchos::ParameterList ¶ms)
Set (or reset) parameters.
Definition: Ifpack2_Chebyshev_def.hpp:92
Teuchos::RCP< const row_matrix_type > getMatrix() const
Get the matrix given to the constructor.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1298
double getComputeTime() const
The total time spent in all calls to compute().
Definition: Ifpack2_Chebyshev_def.hpp:192
int getNumApply() const
The total number of successful calls to apply().
Definition: Ifpack2_Chebyshev_def.hpp:180
MatrixType::node_type node_type
The Node type used by the input MatrixType.
Definition: Ifpack2_Chebyshev_decl.hpp:226
bool isInitialized() const
Definition: Ifpack2_Chebyshev_decl.hpp:432
void compute()
(Re)compute the left scaling, and (if applicable) estimate max and min eigenvalues of D_inv * A...
Definition: Ifpack2_Chebyshev_def.hpp:292
bool hasTransposeApply() const
Whether it's possible to apply the transpose of this operator.
Definition: Ifpack2_Chebyshev_def.hpp:162
double getComputeFlops() const
The total number of floating-point operations taken by all calls to compute().
Definition: Ifpack2_Chebyshev_def.hpp:204
Diagonally scaled Chebyshev iteration for Tpetra sparse matrices.
Definition: Ifpack2_Chebyshev_decl.hpp:199
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, returning the result in Y.
Definition: Ifpack2_Chebyshev_def.hpp:218
bool isComputed() const
Definition: Ifpack2_Chebyshev_decl.hpp:479
std::string description() const
A simple one-line description of this object.
Definition: Ifpack2_Chebyshev_def.hpp:319
void initialize()
Initialize the preconditioner.
Definition: Ifpack2_Chebyshev_def.hpp:277
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
The communicator over which the matrix is distributed.
Definition: Ifpack2_Chebyshev_def.hpp:101
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_Chebyshev_decl.hpp:220
int getNumInitialize() const
The total number of successful calls to initialize().
Definition: Ifpack2_Chebyshev_def.hpp:168
double getInitializeTime() const
The total time spent in all calls to initialize().
Definition: Ifpack2_Chebyshev_def.hpp:186
MatrixType::scalar_type getLambdaMaxForApply() const
The estimate of the maximum eigenvalue used in the apply().
Definition: Ifpack2_Chebyshev_def.hpp:511
std::string description() const
A single-line description of the Chebyshev solver.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1331
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object with some verbosity level to a Teuchos::FancyOStream.
Definition: Ifpack2_Chebyshev_def.hpp:348
Teuchos::RCP< const row_matrix_type > getMatrix() const
The matrix for which this is a preconditioner.
Definition: Ifpack2_Chebyshev_def.hpp:115
Teuchos::RCP< const map_type > getDomainMap() const
The Tpetra::Map representing the domain of this operator.
Definition: Ifpack2_Chebyshev_def.hpp:136
virtual ~Chebyshev()
Destructor.
Definition: Ifpack2_Chebyshev_def.hpp:75
bool hasTransposeApply() const
Whether it's possible to apply the transpose of this operator.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1305
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
Teuchos::RCP< const map_type > getRangeMap() const
The Tpetra::Map representing the range of this operator.
Definition: Ifpack2_Chebyshev_def.hpp:150
int getNumCompute() const
The total number of successful calls to compute().
Definition: Ifpack2_Chebyshev_def.hpp:174
void compute()
(Re)compute the left scaling D_inv, and estimate min and max eigenvalues of D_inv * A...
Definition: Ifpack2_Details_Chebyshev_def.hpp:665
void applyMat(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) const
Compute Y = Op(A)*X, where Op(A) is either A, , or .
Definition: Ifpack2_Chebyshev_def.hpp:258