46 #ifndef MUELU_IFPACK2SMOOTHER_DEF_HPP 47 #define MUELU_IFPACK2SMOOTHER_DEF_HPP 51 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_MUELU_IFPACK2) 55 #include <Tpetra_RowMatrix.hpp> 57 #include <Ifpack2_Chebyshev.hpp> 58 #include <Ifpack2_Factory.hpp> 59 #include <Ifpack2_Parameters.hpp> 70 #include "MueLu_Utilities.hpp" 77 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
79 : type_(type), overlap_(overlap)
84 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
95 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
100 RCP<ParameterList> precList = this->RemoveFactoriesFromList(this->GetParameterList());
102 prec_->setParameters(*precList);
107 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
109 this->Input(currentLevel,
"A");
111 if (type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
112 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
113 type_ ==
"LINESMOOTHING_BANDEDRELAXATION" ||
114 type_ ==
"LINESMOOTHING_BLOCK_RELAXATION" ||
115 type_ ==
"LINESMOOTHING_BLOCK RELAXATION" ||
116 type_ ==
"LINESMOOTHING_BLOCKRELAXATION") {
117 this->Input(currentLevel,
"CoarseNumZLayers");
118 this->Input(currentLevel,
"LineDetection_VertLineIds");
122 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
126 A_ = Factory::Get< RCP<Matrix> >(currentLevel,
"A");
128 if (type_ ==
"SCHWARZ")
129 SetupSchwarz(currentLevel);
131 else if (type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
132 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
133 type_ ==
"LINESMOOTHING_BANDEDRELAXATION" ||
134 type_ ==
"LINESMOOTHING_BLOCK_RELAXATION" ||
135 type_ ==
"LINESMOOTHING_BLOCK RELAXATION" ||
136 type_ ==
"LINESMOOTHING_BLOCKRELAXATION")
137 SetupLineSmoothing(currentLevel);
139 else if (type_ ==
"CHEBYSHEV")
140 SetupChebyshev(currentLevel);
143 SetupGeneric(currentLevel);
147 this->GetOStream(
Statistics1) << description() << std::endl;
150 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
152 typedef Tpetra::RowMatrix<SC,LO,GO,NO> tRowMatrix;
154 bool reusePreconditioner =
false;
155 if (this->IsSetup() ==
true) {
157 this->GetOStream(
Runtime1) <<
"MueLu::Ifpack2Smoother::SetupSchwarz(): Setup() has already been called, assuming reuse" << std::endl;
159 bool isTRowMatrix =
true;
164 isTRowMatrix =
false;
168 if (!prec.
is_null() && isTRowMatrix) {
169 #ifdef IFPACK2_HAS_PROPER_REUSE 170 prec->resetMatrix(tA);
171 reusePreconditioner =
true;
173 this->GetOStream(
Errors) <<
"Ifpack2 does not have proper reuse yet." << std::endl;
177 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupSchwarz(): reuse of this type is not available " 178 "(either failed cast to CanChangeMatrix, or to Tpetra Row Matrix), reverting to full construction" << std::endl;
182 if (!reusePreconditioner) {
185 bool isBlockedMatrix =
false;
188 std::string sublistName =
"subdomain solver parameters";
200 std::string partName =
"partitioner: type";
201 if (subList.
isParameter(partName) && subList.
get<std::string>(partName) ==
"user") {
202 isBlockedMatrix =
true;
206 "Matrix A must be of type BlockedCrsMatrix.");
208 size_t numVels = bA->getMatrix(0,0)->getNodeNumRows();
209 size_t numPres = bA->getMatrix(1,0)->getNodeNumRows();
210 size_t numRows = A_->getNodeNumRows();
214 size_t numBlocks = 0;
215 for (
size_t rowOfB = numVels; rowOfB < numVels+numPres; ++rowOfB)
216 blockSeeds[rowOfB] = numBlocks++;
220 "Matrix A must be of type BlockedCrsMatrix.");
222 merged2Mat = bA2->Merge();
226 bool haveBoundary =
false;
227 for (
LO i = 0; i < boundaryNodes.
size(); i++)
228 if (boundaryNodes[i]) {
232 blockSeeds[i] = numBlocks;
238 subList.
set(
"partitioner: map", blockSeeds);
239 subList.
set(
"partitioner: local parts", as<int>(numBlocks));
244 isBlockedMatrix =
true;
245 merged2Mat = bA->Merge();
263 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
265 if (this->IsSetup() ==
true) {
266 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupLineSmoothing(): Setup() has already been called" << std::endl;
267 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupLineSmoothing(): reuse of this type is not available, reverting to full construction" << std::endl;
272 LO CoarseNumZLayers = Factory::Get<LO>(currentLevel,
"CoarseNumZLayers");
273 if (CoarseNumZLayers > 0) {
274 Teuchos::ArrayRCP<LO> TVertLineIdSmoo = Factory::Get< Teuchos::ArrayRCP<LO> >(currentLevel,
"LineDetection_VertLineIds");
278 for(
size_t k = 0; k < Teuchos::as<size_t>(TVertLineIdSmoo.
size()); k++) {
279 if(maxPart < TVertLineIdSmoo[k]) maxPart = TVertLineIdSmoo[k];
282 size_t numLocalRows = A_->getNodeNumRows();
284 "MueLu::Ifpack2Smoother::Setup(): the number of local nodes is incompatible with the TVertLineIdsSmoo.");
286 if (numLocalRows == Teuchos::as<size_t>(TVertLineIdSmoo.
size())) {
287 myparamList.
set(
"partitioner: type",
"user");
288 myparamList.
set(
"partitioner: map",TVertLineIdSmoo);
289 myparamList.
set(
"partitioner: local parts",maxPart+1);
292 size_t numDofsPerNode = numLocalRows / TVertLineIdSmoo.
size();
296 for (
size_t blockRow = 0; blockRow < Teuchos::as<size_t>(TVertLineIdSmoo.
size()); ++blockRow)
297 for (
size_t dof = 0; dof < numDofsPerNode; dof++)
298 partitionerMap[blockRow * numDofsPerNode + dof] = TVertLineIdSmoo[blockRow];
299 myparamList.
set(
"partitioner: type",
"user");
300 myparamList.
set(
"partitioner: map",partitionerMap);
301 myparamList.
set(
"partitioner: local parts",maxPart + 1);
304 if (type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
305 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
306 type_ ==
"LINESMOOTHING_BANDEDRELAXATION")
307 type_ =
"BANDEDRELAXATION";
309 type_ =
"BLOCKRELAXATION";
312 this->GetOStream(
Runtime0) <<
"Line detection failed: fall back to point-wise relaxation" << std::endl;
313 myparamList.
remove(
"partitioner: type",
false);
314 myparamList.
remove(
"partitioner: map",
false);
315 myparamList.
remove(
"partitioner: local parts",
false);
316 type_ =
"RELAXATION";
327 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
329 if (this->IsSetup() ==
true) {
330 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupChebyshev(): SetupChebyshev() has already been called" << std::endl;
331 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupChebyshev(): reuse of this type is not available, reverting to full construction" << std::endl;
335 SC negone = -STS::one();
337 SC lambdaMax = negone;
339 std::string maxEigString =
"chebyshev: max eigenvalue";
340 std::string eigRatioString =
"chebyshev: ratio eigenvalue";
346 if (paramList.
isType<
double>(maxEigString))
347 lambdaMax = paramList.
get<
double>(maxEigString);
349 lambdaMax = paramList.
get<
SC>(maxEigString);
350 this->GetOStream(
Statistics1) << maxEigString <<
" (cached with smoother parameter list) = " << lambdaMax << std::endl;
353 lambdaMax = A_->GetMaxEigenvalueEstimate();
354 if (lambdaMax != negone) {
355 this->GetOStream(
Statistics1) << maxEigString <<
" (cached with matrix) = " << lambdaMax << std::endl;
356 paramList.
set(maxEigString, lambdaMax);
361 const SC defaultEigRatio = 20;
363 SC ratio = defaultEigRatio;
365 if (paramList.
isType<
double>(eigRatioString))
366 ratio = paramList.
get<
double>(eigRatioString);
368 ratio = paramList.
get<
SC>(eigRatioString);
376 size_t nRowsFine = fineA->getGlobalNumRows();
377 size_t nRowsCoarse = A_->getGlobalNumRows();
379 SC levelRatio = as<SC>(as<float>(nRowsFine)/nRowsCoarse);
380 if (STS::magnitude(levelRatio) > STS::magnitude(ratio))
384 this->GetOStream(
Statistics1) << eigRatioString <<
" (computed) = " << ratio << std::endl;
385 paramList.
set(eigRatioString, ratio);
395 if (lambdaMax == negone) {
396 typedef Tpetra::RowMatrix<SC, LO, GO, NO> MatrixType;
399 if (chebyPrec != Teuchos::null) {
400 lambdaMax = chebyPrec->getLambdaMaxForApply();
401 A_->SetMaxEigenvalueEstimate(lambdaMax);
402 this->GetOStream(
Statistics1) <<
"chebyshev: max eigenvalue (calculated by Ifpack2)" <<
" = " << lambdaMax << std::endl;
408 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
410 typedef Tpetra::RowMatrix<SC,LO,GO,NO> tRowMatrix;
418 bool reusePreconditioner =
false;
419 if (this->IsSetup() ==
true) {
421 this->GetOStream(
Runtime1) <<
"MueLu::Ifpack2Smoother::SetupGeneric(): Setup() has already been called, assuming reuse" << std::endl;
425 #ifdef IFPACK2_HAS_PROPER_REUSE 426 prec->resetMatrix(tA);
427 reusePreconditioner =
true;
429 this->GetOStream(
Errors) <<
"Ifpack2 does not have proper reuse yet." << std::endl;
433 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupSchwarz(): reuse of this type is not available (failed cast to CanChangeMatrix), " 434 "reverting to full construction" << std::endl;
438 if (!reusePreconditioner) {
447 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
461 bool supportInitialGuess =
false;
462 if (type_ ==
"CHEBYSHEV") {
463 paramList.
set(
"chebyshev: zero starting solution", InitialGuessIsZero);
464 SetPrecParameters(paramList);
465 supportInitialGuess =
true;
467 }
else if (type_ ==
"RELAXATION") {
468 paramList.
set(
"relaxation: zero starting solution", InitialGuessIsZero);
469 SetPrecParameters(paramList);
470 supportInitialGuess =
true;
472 }
else if (type_ ==
"KRYLOV") {
473 paramList.
set(
"krylov: zero starting solution", InitialGuessIsZero);
474 SetPrecParameters(paramList);
475 supportInitialGuess =
true;
477 }
else if (type_ ==
"SCHWARZ") {
478 paramList.
set(
"schwarz: zero starting solution", InitialGuessIsZero);
484 supportInitialGuess =
true;
494 if (InitialGuessIsZero || supportInitialGuess) {
497 prec_->apply(tpB, tpX);
501 RCP<MultiVector> Correction = MultiVectorFactory::Build(A_->getDomainMap(), X.getNumVectors());
506 prec_->apply(tpB, tpX);
508 X.update(TST::one(), *Correction, TST::one());
512 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
519 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
521 std::ostringstream out;
523 out << prec_->description();
526 out <<
"{type = " << type_ <<
"}";
531 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
536 out0 <<
"Prec. type: " << type_ << std::endl;
539 out0 <<
"Parameter list: " << std::endl;
541 out << this->GetParameterList();
542 out0 <<
"Overlap: " << overlap_ << std::endl;
546 if (prec_ != Teuchos::null) {
548 out << *prec_ << std::endl;
551 if (verbLevel &
Debug) {
554 <<
"RCP<prec_>: " << prec_ << std::endl;
560 #endif // HAVE_MUELU_TPETRA && HAVE_MUELU_IFPACK2 561 #endif // MUELU_IFPACK2SMOOTHER_DEF_HPP Important warning messages (one line)
void SetupGeneric(Level ¤tLevel)
RCP< SmootherPrototype > Copy() const
Exception indicating invalid cast attempted.
RCP< Level > & GetPreviousLevel()
Previous level.
void SetPrecParameters(const Teuchos::ParameterList &list=Teuchos::ParameterList()) const
static RCP< const Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2TpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > const vec)
Helper utility to pull out the underlying Tpetra objects from an Xpetra object.
std::string description() const
Return a simple one-line description of this object.
Print external lib objects.
bool IsSetup() const
Get the state of a smoother prototype.
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
friend class Ifpack2Smoother
Constructor.
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
bool isSublist(const std::string &name) const
Print additional debugging information.
One-liner description of what is happening.
void SetParameterList(const Teuchos::ParameterList ¶mList)
Set parameters from a parameter list and return with default values.
Namespace for MueLu classes and methods.
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
int GetLevelID() const
Return level number.
void SetupSchwarz(Level ¤tLevel)
void SetupLineSmoothing(Level ¤tLevel)
bool remove(std::string const &name, bool throwIfNotExists=true)
static RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &RHS)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual void SetParameterList(const ParameterList ¶mList)
Set parameters from a parameter list and return with default values.
Class that holds all level-specific information.
void Setup(Level ¤tLevel)
Set up the smoother.
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the preconditioner.
static Teuchos::RCP< Preconditioner< typename MatrixType::scalar_type, typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type > > create(const std::string &precType, const Teuchos::RCP< const MatrixType > &matrix)
bool isType(const std::string &name) const
void DeclareInput(Level ¤tLevel) const
Input.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
ParameterList & setParameters(const ParameterList &source)
bool isParameter(const std::string &name) const
static Teuchos::ArrayRCP< const bool > DetectDirichletRows(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::magnitude(0.))
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Print class parameters (more parameters, more verbose)
Exception throws to report errors in the internal logical of the program.
void SetupChebyshev(Level ¤tLevel)
Description of what is happening (more verbose)
Class that encapsulates Ifpack2 smoothers.
static RCP< Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraRow(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
virtual std::string description() const
Return a simple one-line description of this object.
std::string toString(const T &t)