43 #ifndef IFPACK2_HIPTMAIR_DEF_HPP 44 #define IFPACK2_HIPTMAIR_DEF_HPP 46 #include "Ifpack2_Details_OneLevelFactory.hpp" 47 #include "Ifpack2_Parameters.hpp" 48 #include "Teuchos_TimeMonitor.hpp" 49 #include "Tpetra_MultiVector.hpp" 56 template <
class MatrixType>
58 Hiptmair (
const Teuchos::RCP<const row_matrix_type>& A,
59 const Teuchos::RCP<const row_matrix_type>& PtAP,
60 const Teuchos::RCP<const row_matrix_type>& P) :
65 precType1_ (
"CHEBYSHEV"),
66 precType2_ (
"CHEBYSHEV"),
68 ZeroStartingSolution_ (true),
70 IsInitialized_ (false),
75 InitializeTime_ (0.0),
81 template <
class MatrixType>
84 template <
class MatrixType>
88 using Teuchos::ParameterList;
89 using Teuchos::Exceptions::InvalidParameterName;
90 using Teuchos::Exceptions::InvalidParameterType;
92 ParameterList params = plist;
99 std::string precType1 = precType1_;
100 std::string precType2 = precType2_;
101 std::string preOrPost = preOrPost_;
102 Teuchos::ParameterList precList1 = precList1_;
103 Teuchos::ParameterList precList2 = precList2_;
104 bool zeroStartingSolution = ZeroStartingSolution_;
106 precType1 = params.get(
"hiptmair: smoother type 1", precType1);
107 precType2 = params.get(
"hiptmair: smoother type 2", precType2);
108 precList1 = params.get(
"hiptmair: smoother list 1", precList1);
109 precList2 = params.get(
"hiptmair: smoother list 2", precList2);
110 preOrPost = params.get(
"hiptmair: pre or post", preOrPost);
111 zeroStartingSolution = params.get(
"hiptmair: zero starting solution",
112 zeroStartingSolution);
115 precType1_ = precType1;
116 precType2_ = precType2;
117 precList1_ = precList1;
118 precList2_ = precList2;
119 preOrPost_ = preOrPost;
120 ZeroStartingSolution_ = zeroStartingSolution;
124 template <
class MatrixType>
125 Teuchos::RCP<const Teuchos::Comm<int> >
127 TEUCHOS_TEST_FOR_EXCEPTION(
128 A_.is_null (), std::runtime_error,
"Ifpack2::Hiptmair::getComm: " 129 "The input matrix A is null. Please call setMatrix() with a nonnull " 130 "input matrix before calling this method.");
131 return A_->getComm ();
135 template <
class MatrixType>
136 Teuchos::RCP<const Tpetra::RowMatrix<typename MatrixType::scalar_type,typename MatrixType::local_ordinal_type,typename MatrixType::global_ordinal_type,typename MatrixType::node_type> >
142 template <
class MatrixType>
143 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,typename MatrixType::global_ordinal_type,typename MatrixType::node_type> >
146 TEUCHOS_TEST_FOR_EXCEPTION(
147 A_.is_null (), std::runtime_error,
"Ifpack2::Hiptmair::getDomainMap: " 148 "The input matrix A is null. Please call setMatrix() with a nonnull " 149 "input matrix before calling this method.");
150 return A_->getDomainMap ();
154 template <
class MatrixType>
155 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,typename MatrixType::global_ordinal_type,typename MatrixType::node_type> >
158 TEUCHOS_TEST_FOR_EXCEPTION(
159 A_.is_null (), std::runtime_error,
"Ifpack2::Hiptmair::getRangeMap: " 160 "The input matrix A is null. Please call setMatrix() with a nonnull " 161 "input matrix before calling this method.");
162 return A_->getRangeMap ();
166 template <
class MatrixType>
174 template <
class MatrixType>
176 return NumInitialize_;
180 template <
class MatrixType>
186 template <
class MatrixType>
192 template <
class MatrixType>
194 return InitializeTime_;
198 template <
class MatrixType>
204 template <
class MatrixType>
210 template <
class MatrixType>
213 using Teuchos::ParameterList;
217 TEUCHOS_TEST_FOR_EXCEPTION(
218 A_.is_null (), std::runtime_error,
"Ifpack2::Hiptmair::initialize: " 219 "The input matrix A is null. Please call setMatrix() with a nonnull " 220 "input matrix before calling this method.");
223 IsInitialized_ =
false;
226 Teuchos::Time timer (
"initialize");
228 Teuchos::TimeMonitor timeMon (timer);
232 ifpack2_prec1_=factory.
create(precType1_,A_);
233 ifpack2_prec1_->initialize();
234 ifpack2_prec1_->setParameters(precList1_);
236 ifpack2_prec2_=factory.
create(precType2_,PtAP_);
237 ifpack2_prec2_->initialize();
238 ifpack2_prec2_->setParameters(precList2_);
241 IsInitialized_ =
true;
243 InitializeTime_ += timer.totalElapsedTime ();
247 template <
class MatrixType>
250 TEUCHOS_TEST_FOR_EXCEPTION(
251 A_.is_null (), std::runtime_error,
"Ifpack2::Hiptmair::compute: " 252 "The input matrix A is null. Please call setMatrix() with a nonnull " 253 "input matrix before calling this method.");
256 if (! isInitialized ()) {
260 Teuchos::Time timer (
"compute");
262 Teuchos::TimeMonitor timeMon (timer);
263 ifpack2_prec1_->compute();
264 ifpack2_prec2_->compute();
268 ComputeTime_ += timer.totalElapsedTime ();
272 template <
class MatrixType>
274 apply (
const Tpetra::MultiVector<
typename MatrixType::scalar_type,
275 typename MatrixType::local_ordinal_type,
276 typename MatrixType::global_ordinal_type,
277 typename MatrixType::node_type>& X,
278 Tpetra::MultiVector<
typename MatrixType::scalar_type,
279 typename MatrixType::local_ordinal_type,
280 typename MatrixType::global_ordinal_type,
281 typename MatrixType::node_type>& Y,
282 Teuchos::ETransp mode,
283 typename MatrixType::scalar_type alpha,
284 typename MatrixType::scalar_type beta)
const 288 using Teuchos::rcpFromRef;
291 TEUCHOS_TEST_FOR_EXCEPTION(
292 ! isComputed (), std::runtime_error,
293 "Ifpack2::Hiptmair::apply: You must call compute() before you may call apply().");
294 TEUCHOS_TEST_FOR_EXCEPTION(
295 X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
296 "Ifpack2::Hiptmair::apply: The MultiVector inputs X and Y do not have the " 297 "same number of columns. X.getNumVectors() = " << X.getNumVectors ()
298 <<
" != Y.getNumVectors() = " << Y.getNumVectors () <<
".");
301 TEUCHOS_TEST_FOR_EXCEPTION(
302 alpha != STS::one (), std::logic_error,
303 "Ifpack2::Hiptmair::apply: alpha != 1 has not been implemented.");
304 TEUCHOS_TEST_FOR_EXCEPTION(
305 beta != STS::zero (), std::logic_error,
306 "Ifpack2::Hiptmair::apply: zero != 0 has not been implemented.");
307 TEUCHOS_TEST_FOR_EXCEPTION(
308 mode != Teuchos::NO_TRANS, std::logic_error,
309 "Ifpack2::Hiptmair::apply: mode != Teuchos::NO_TRANS has not been implemented.");
311 Teuchos::Time timer (
"apply");
313 Teuchos::TimeMonitor timeMon (timer);
319 auto X_lcl_host = X.getLocalViewHost ();
320 auto Y_lcl_host = Y.getLocalViewHost ();
322 if (X_lcl_host.data () == Y_lcl_host.data ()) {
323 Xcopy = rcp (
new MV (X, Teuchos::Copy));
325 Xcopy = rcpFromRef (X);
329 RCP<MV> Ycopy = rcpFromRef (Y);
330 if (ZeroStartingSolution_) {
331 Ycopy->putScalar (STS::zero ());
335 applyHiptmairSmoother (*Xcopy, *Ycopy);
339 ApplyTime_ += timer.totalElapsedTime ();
342 template <
class MatrixType>
345 typename MatrixType::local_ordinal_type,
346 typename MatrixType::global_ordinal_type,
347 typename MatrixType::node_type>& X,
348 Tpetra::MultiVector<
typename MatrixType::scalar_type,
349 typename MatrixType::local_ordinal_type,
350 typename MatrixType::global_ordinal_type,
351 typename MatrixType::node_type>& Y)
const 355 using Teuchos::rcpFromRef;
356 typedef Tpetra::MultiVector<scalar_type, local_ordinal_type,
357 global_ordinal_type, node_type> MV;
358 const scalar_type ZERO = STS::zero ();
359 const scalar_type ONE = STS::one ();
361 RCP<MV> res1 = rcp (
new MV (A_->getRowMap (), X.getNumVectors ()));
362 RCP<MV> vec1 = rcp (
new MV (A_->getRowMap (), X.getNumVectors ()));
363 RCP<MV> res2 = rcp (
new MV (PtAP_->getRowMap (), X.getNumVectors ()));
364 RCP<MV> vec2 = rcp (
new MV (PtAP_->getRowMap (), X.getNumVectors ()));
366 if (preOrPost_ ==
"pre" || preOrPost_ ==
"both") {
368 A_->apply (Y, *res1);
369 res1->update (ONE, X, -ONE);
370 vec1->putScalar (ZERO);
371 ifpack2_prec1_->apply (*res1, *vec1);
372 Y.update (ONE, *vec1, ONE);
376 A_->apply (Y, *res1);
377 res1->update (ONE, X, -ONE);
378 P_->apply (*res1, *res2, Teuchos::TRANS);
379 vec2->putScalar (ZERO);
380 ifpack2_prec2_->apply (*res2, *vec2);
381 P_->apply (*vec2, *vec1, Teuchos::NO_TRANS);
382 Y.update (ONE,*vec1,ONE);
384 if (preOrPost_ ==
"post" || preOrPost_ ==
"both") {
386 A_->apply (Y, *res1);
387 res1->update (ONE, X, -ONE);
388 vec1->putScalar (ZERO);
389 ifpack2_prec1_->apply (*res1, *vec1);
390 Y.update (ONE, *vec1, ONE);
394 template <
class MatrixType>
397 std::ostringstream os;
402 os <<
"\"Ifpack2::Hiptmair\": {";
403 if (this->getObjectLabel () !=
"") {
404 os <<
"Label: \"" << this->getObjectLabel () <<
"\", ";
406 os <<
"Initialized: " << (isInitialized () ?
"true" :
"false") <<
", " 407 <<
"Computed: " << (isComputed () ?
"true" :
"false") <<
", ";
410 os <<
"Matrix: null";
413 os <<
"Matrix: not null" 414 <<
", Global matrix dimensions: [" 415 << A_->getGlobalNumRows () <<
", " << A_->getGlobalNumCols () <<
"]";
423 template <
class MatrixType>
426 const Teuchos::EVerbosityLevel verbLevel)
const 430 using Teuchos::VERB_DEFAULT;
431 using Teuchos::VERB_NONE;
432 using Teuchos::VERB_LOW;
433 using Teuchos::VERB_MEDIUM;
434 using Teuchos::VERB_HIGH;
435 using Teuchos::VERB_EXTREME;
437 const Teuchos::EVerbosityLevel vl =
438 (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
440 if (vl != VERB_NONE) {
442 Teuchos::OSTab tab0 (out);
443 out <<
"\"Ifpack2::Hiptmair\":";
445 Teuchos::OSTab tab1 (out);
446 if (this->getObjectLabel () !=
"") {
447 out <<
"Label: " << this->getObjectLabel () << endl;
449 out <<
"Initialized: " << (isInitialized () ?
"true" :
"false") << endl
450 <<
"Computed: " << (isComputed () ?
"true" :
"false") << endl
451 <<
"Global number of rows: " << A_->getGlobalNumRows () << endl
452 <<
"Global number of columns: " << A_->getGlobalNumCols () << endl
455 out <<
" null" << endl;
457 A_->describe (out, vl);
464 #define IFPACK2_HIPTMAIR_INSTANT(S,LO,GO,N) \ 465 template class Ifpack2::Hiptmair< Tpetra::RowMatrix<S, LO, GO, N> >; void initialize()
Do any initialization that depends on the input matrix's structure.
Definition: Ifpack2_Hiptmair_def.hpp:211
int getNumInitialize() const
Returns the number of calls to Initialize().
Definition: Ifpack2_Hiptmair_def.hpp:175
double getComputeTime() const
Returns the time spent in Compute().
Definition: Ifpack2_Hiptmair_def.hpp:199
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_Hiptmair_decl.hpp:85
double getApplyTime() const
Returns the time spent in apply().
Definition: Ifpack2_Hiptmair_def.hpp:205
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_Hiptmair_decl.hpp:82
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_Hiptmair_def.hpp:156
int getNumCompute() const
Returns the number of calls to Compute().
Definition: Ifpack2_Hiptmair_def.hpp:181
void setParameters(const Teuchos::ParameterList ¶ms)
Set the preconditioner's parameters.
Definition: Ifpack2_Hiptmair_def.hpp:85
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_Hiptmair_def.hpp:144
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_Hiptmair_def.hpp:137
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input MatrixType.
Definition: Ifpack2_Hiptmair_decl.hpp:88
Wrapper for Hiptmair smoothers.
Definition: Ifpack2_Hiptmair_decl.hpp:71
Hiptmair(const Teuchos::RCP< const row_matrix_type > &A, const Teuchos::RCP< const row_matrix_type > &PtAP, const Teuchos::RCP< const row_matrix_type > &P)
Constructor that takes 3 Tpetra matrices.
Definition: Ifpack2_Hiptmair_def.hpp:58
MatrixType::node_type node_type
The Node type used by the input MatrixType.
Definition: Ifpack2_Hiptmair_decl.hpp:91
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_Hiptmair_def.hpp:425
void compute()
Do any initialization that depends on the input matrix's values.
Definition: Ifpack2_Hiptmair_def.hpp:248
virtual ~Hiptmair()
Destructor.
Definition: Ifpack2_Hiptmair_def.hpp:82
double getInitializeTime() const
Returns the time spent in Initialize().
Definition: Ifpack2_Hiptmair_def.hpp:193
std::string description() const
Return a simple one-line description of this object.
Definition: Ifpack2_Hiptmair_def.hpp:395
bool hasTransposeApply() const
Whether this object's apply() method can apply the transpose (or conjugate transpose, if applicable).
Definition: Ifpack2_Hiptmair_def.hpp:167
int getNumApply() const
Returns the number of calls to apply().
Definition: Ifpack2_Hiptmair_def.hpp:187
"Factory" for creating single-level preconditioners.
Definition: Ifpack2_Details_OneLevelFactory_decl.hpp:123
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:73
Teuchos::RCP< prec_type > create(const std::string &precType, const Teuchos::RCP< const row_matrix_type > &matrix) const
Create an instance of Preconditioner given the string name of the preconditioner type.
Definition: Ifpack2_Details_OneLevelFactory_def.hpp:78
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Returns the operator's communicator.
Definition: Ifpack2_Hiptmair_def.hpp:126
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_Hiptmair_def.hpp:274