46 #ifndef MUELU_IFPACK2SMOOTHER_DEF_HPP 47 #define MUELU_IFPACK2SMOOTHER_DEF_HPP 51 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_MUELU_IFPACK2) 53 #include <Teuchos_ParameterList.hpp> 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>
97 ParameterList& paramList =
const_cast<ParameterList&
>(this->GetParameterList());
98 paramList.setParameters(list);
100 RCP<ParameterList> precList = this->RemoveFactoriesFromList(this->GetParameterList());
102 prec_->setParameters(*precList);
104 paramList.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;
160 RCP<const tRowMatrix> tA;
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) {
183 ParameterList& paramList =
const_cast<ParameterList&
>(this->GetParameterList());
185 bool isBlockedMatrix =
false;
186 RCP<Matrix> merged2Mat;
188 std::string sublistName =
"subdomain solver parameters";
189 if (paramList.isSublist(sublistName)) {
198 ParameterList& subList = paramList.sublist(sublistName);
200 std::string partName =
"partitioner: type";
201 if (subList.isParameter(partName) && subList.get<std::string>(partName) ==
"user") {
202 isBlockedMatrix =
true;
204 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
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();
212 ArrayRCP<LocalOrdinal> blockSeeds(numRows, Teuchos::OrdinalTraits<LocalOrdinal>::invalid());
214 size_t numBlocks = 0;
215 for (
size_t rowOfB = numVels; rowOfB < numVels+numPres; ++rowOfB)
216 blockSeeds[rowOfB] = numBlocks++;
218 RCP<BlockedCrsMatrix> bA2 = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
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));
242 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
244 isBlockedMatrix =
true;
245 merged2Mat = bA->Merge();
250 RCP<const tRowMatrix> tA;
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;
270 ParameterList& myparamList =
const_cast<ParameterList&
>(this->GetParameterList());
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();
295 Teuchos::ArrayRCP<LO> partitionerMap(numLocalRows, Teuchos::OrdinalTraits<LocalOrdinal>::invalid());
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;
334 typedef Teuchos::ScalarTraits<SC> STS;
335 SC negone = -STS::one();
337 SC lambdaMax = negone;
339 std::string maxEigString =
"chebyshev: max eigenvalue";
340 std::string eigRatioString =
"chebyshev: ratio eigenvalue";
342 ParameterList& paramList =
const_cast<ParameterList&
>(this->GetParameterList());
345 if (paramList.isParameter(maxEigString)) {
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;
364 if (paramList.isParameter(eigRatioString)) {
365 if (paramList.isType<
double>(eigRatioString))
366 ratio = paramList.get<
double>(eigRatioString);
368 ratio = paramList.get<
SC>(eigRatioString);
375 RCP<const Matrix> fineA = currentLevel.
GetPreviousLevel()->Get<RCP<Matrix> >(
"A");
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) {
401 A_->SetMaxEigenvalueEstimate(lambdaMax);
402 this->GetOStream(
Statistics1) <<
"chebyshev: max eigenvalue (calculated by Ifpack2)" <<
" = " << lambdaMax << std::endl;
404 TEUCHOS_TEST_FOR_EXCEPTION(lambdaMax == negone,
Exceptions::RuntimeError,
"MueLu::Ifpack2Smoother::Setup(): no maximum eigenvalue estimate");
408 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
410 typedef Tpetra::RowMatrix<SC,LO,GO,NO> tRowMatrix;
412 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
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;
424 if (!prec.is_null()) {
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>
460 Teuchos::ParameterList paramList;
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);
483 prec_->setParameters(paramList);
484 supportInitialGuess =
true;
494 if (InitialGuessIsZero || supportInitialGuess) {
497 prec_->apply(tpB, tpX);
499 typedef Teuchos::ScalarTraits<Scalar> TST;
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>
515 smoother->SetParameterList(this->GetParameterList());
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;
540 Teuchos::OSTab tab2(out);
541 out << this->GetParameterList();
542 out0 <<
"Overlap: " << overlap_ << std::endl;
546 if (prec_ != Teuchos::null) {
547 Teuchos::OSTab tab2(out);
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 toString(const T &what)
Little helper function to convert non-string types to strings.
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.
friend class Ifpack2Smoother
Constructor.
Timer to be used in factories. Similar to Monitor but with additional timers.
Print additional debugging information.
One-liner description of what is happening.
void SetParameterList(const Teuchos::ParameterList ¶mList)
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)
MatrixType::scalar_type getLambdaMaxForApply() const
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)
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)
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.
static Teuchos::ArrayRCP< const bool > DetectDirichletRows(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::magnitude(0.))
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.