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> 61 #include <Xpetra_BlockedCrsMatrix.hpp> 62 #include <Xpetra_CrsMatrix.hpp> 63 #include <Xpetra_CrsMatrixWrap.hpp> 64 #include <Xpetra_Matrix.hpp> 65 #include <Xpetra_MultiVectorFactory.hpp> 70 #include "MueLu_Utilities.hpp" 75 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
77 : type_(type), overlap_(overlap)
82 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
93 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
95 ParameterList& paramList =
const_cast<ParameterList&
>(this->
GetParameterList());
96 paramList.setParameters(list);
100 prec_->setParameters(*precList);
102 paramList.setParameters(*precList);
105 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
107 this->
Input(currentLevel,
"A");
109 if (
type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
110 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
111 type_ ==
"LINESMOOTHING_BANDEDRELAXATION" ||
112 type_ ==
"LINESMOOTHING_BLOCK_RELAXATION" ||
113 type_ ==
"LINESMOOTHING_BLOCK RELAXATION" ||
114 type_ ==
"LINESMOOTHING_BLOCKRELAXATION") {
115 this->
Input(currentLevel,
"CoarseNumZLayers");
116 this->
Input(currentLevel,
"LineDetection_VertLineIds");
120 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
125 this->
GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::Setup(): Setup() has already been called";
127 A_ = Factory::Get< RCP<Matrix> >(currentLevel,
"A");
129 typedef Teuchos::ScalarTraits<SC> STS;
130 SC negone = -STS::one();
132 SC lambdaMax = negone;
143 bool isBlockedMatrix =
false;
144 RCP<Matrix> merged2Mat;
145 if (
type_ ==
"SCHWARZ") {
146 ParameterList& paramList =
const_cast<ParameterList&
>(this->
GetParameterList());
148 std::string sublistName =
"subdomain solver parameters";
149 if (paramList.isSublist(sublistName)) {
150 ParameterList& subList = paramList.sublist(sublistName);
152 std::string partName =
"partitioner: type";
153 if (subList.isParameter(partName) && subList.get<std::string>(partName) ==
"user") {
154 isBlockedMatrix =
true;
156 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(
A_);
158 "Matrix A must be of type BlockedCrsMatrix.");
160 size_t numVels = bA->getMatrix(0,0)->getNodeNumRows();
161 size_t numPres = bA->getMatrix(1,0)->getNodeNumRows();
162 size_t numRows =
A_->getNodeNumRows();
164 ArrayRCP<LocalOrdinal> blockSeeds(numRows, Teuchos::OrdinalTraits<LocalOrdinal>::invalid());
166 size_t numBlocks = 0;
167 for (
size_t rowOfB = numVels; rowOfB < numVels+numPres; ++rowOfB)
168 blockSeeds[rowOfB] = numBlocks++;
170 RCP<BlockedCrsMatrix> bA2 = rcp_dynamic_cast<BlockedCrsMatrix>(
A_);
172 "Matrix A must be of type BlockedCrsMatrix.");
174 RCP<CrsMatrix> mergedMat = bA2->Merge();
175 merged2Mat = rcp(
new CrsMatrixWrap(mergedMat));
178 ArrayRCP<const bool> boundaryNodes;
180 bool haveBoundary =
false;
181 for (LO i = 0; i < boundaryNodes.size(); i++)
182 if (boundaryNodes[i]) {
186 blockSeeds[i] = numBlocks;
192 subList.set(
"partitioner: map", blockSeeds);
193 subList.set(
"partitioner: local parts", as<int>(numBlocks));
198 if (
type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
199 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
200 type_ ==
"LINESMOOTHING_BANDEDRELAXATION" ||
201 type_ ==
"LINESMOOTHING_BLOCK_RELAXATION" ||
202 type_ ==
"LINESMOOTHING_BLOCK RELAXATION" ||
203 type_ ==
"LINESMOOTHING_BLOCKRELAXATION" ) {
204 ParameterList& myparamList =
const_cast<ParameterList&
>(this->
GetParameterList());
206 LO CoarseNumZLayers = Factory::Get<LO>(currentLevel,
"CoarseNumZLayers");
207 if (CoarseNumZLayers > 0) {
208 Teuchos::ArrayRCP<LO> TVertLineIdSmoo = Factory::Get< Teuchos::ArrayRCP<LO> >(currentLevel,
"LineDetection_VertLineIds");
212 for(
size_t k = 0; k < Teuchos::as<size_t>(TVertLineIdSmoo.size()); k++) {
213 if(maxPart < TVertLineIdSmoo[k]) maxPart = TVertLineIdSmoo[k];
216 size_t numLocalRows =
A_->getNodeNumRows();
217 TEUCHOS_TEST_FOR_EXCEPTION(numLocalRows % TVertLineIdSmoo.size() != 0,
Exceptions::RuntimeError,
"MueLu::Ifpack2Smoother::Setup(): the number of local nodes is incompatible with the TVertLineIdsSmoo.");
219 if (numLocalRows == Teuchos::as<size_t>(TVertLineIdSmoo.size())) {
220 myparamList.set(
"partitioner: type",
"user");
221 myparamList.set(
"partitioner: map",TVertLineIdSmoo);
222 myparamList.set(
"partitioner: local parts",maxPart+1);
225 size_t numDofsPerNode = numLocalRows / TVertLineIdSmoo.size();
228 Teuchos::ArrayRCP<LO> partitionerMap(numLocalRows, Teuchos::OrdinalTraits<LocalOrdinal>::invalid());
229 for (
size_t blockRow = 0; blockRow < Teuchos::as<size_t>(TVertLineIdSmoo.size()); ++blockRow)
230 for (
size_t dof = 0; dof < numDofsPerNode; dof++)
231 partitionerMap[blockRow * numDofsPerNode + dof] = TVertLineIdSmoo[blockRow];
232 myparamList.set(
"partitioner: type",
"user");
233 myparamList.set(
"partitioner: map",partitionerMap);
234 myparamList.set(
"partitioner: local parts",maxPart + 1);
237 if (
type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
238 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
239 type_ ==
"LINESMOOTHING_BANDEDRELAXATION")
240 type_ =
"BANDEDRELAXATION";
242 type_ =
"BLOCKRELAXATION";
245 this->
GetOStream(
Runtime0) <<
"Line detection failed: fall back to point-wise relaxation" << std::endl;
246 myparamList.remove(
"partitioner: type",
false);
247 myparamList.remove(
"partitioner: map",
false);
248 myparamList.remove(
"partitioner: local parts",
false);
249 type_ =
"RELAXATION";
254 if (
type_ ==
"CHEBYSHEV") {
255 std::string maxEigString =
"chebyshev: max eigenvalue";
256 std::string eigRatioString =
"chebyshev: ratio eigenvalue";
258 ParameterList& paramList =
const_cast<ParameterList&
>(this->
GetParameterList());
261 if (paramList.isParameter(maxEigString)) {
262 if (paramList.isType<
double>(maxEigString))
263 lambdaMax = paramList.get<
double>(maxEigString);
265 lambdaMax = paramList.get<SC>(maxEigString);
266 this->
GetOStream(
Statistics1) << maxEigString <<
" (cached with smoother parameter list) = " << lambdaMax << std::endl;
269 lambdaMax =
A_->GetMaxEigenvalueEstimate();
270 if (lambdaMax != negone) {
271 this->
GetOStream(
Statistics1) << maxEigString <<
" (cached with matrix) = " << lambdaMax << std::endl;
272 paramList.set(maxEigString, lambdaMax);
277 const SC defaultEigRatio = 20;
279 SC ratio = defaultEigRatio;
280 if (paramList.isParameter(eigRatioString)) {
281 if (paramList.isType<
double>(eigRatioString))
282 ratio = paramList.get<
double>(eigRatioString);
284 ratio = paramList.get<SC>(eigRatioString);
286 if (currentLevel.GetLevelID()) {
291 RCP<const Matrix> fineA = currentLevel.GetPreviousLevel()->Get<RCP<Matrix> >(
"A");
292 size_t nRowsFine = fineA->getGlobalNumRows();
293 size_t nRowsCoarse =
A_->getGlobalNumRows();
295 SC levelRatio = as<SC>(as<float>(nRowsFine)/nRowsCoarse);
296 if (STS::magnitude(levelRatio) > STS::magnitude(ratio))
301 paramList.set(eigRatioString, ratio);
304 RCP<const Tpetra::RowMatrix<SC, LO, GO, NO> > tpA;
314 if (
type_ ==
"CHEBYSHEV" && lambdaMax == negone) {
315 typedef Tpetra::RowMatrix<SC, LO, GO, NO> MatrixType;
318 if (chebyPrec != Teuchos::null) {
320 A_->SetMaxEigenvalueEstimate(lambdaMax);
321 this->
GetOStream(
Statistics1) <<
"chebyshev: max eigenvalue (calculated by Ifpack2)" <<
" = " << lambdaMax << std::endl;
323 TEUCHOS_TEST_FOR_EXCEPTION(lambdaMax == negone,
Exceptions::RuntimeError,
"MueLu::IfpackSmoother::Setup(): no maximum eigenvalue estimate");
329 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
342 Teuchos::ParameterList paramList;
343 bool supportInitialGuess =
false;
344 if (
type_ ==
"CHEBYSHEV") {
345 paramList.set(
"chebyshev: zero starting solution", InitialGuessIsZero);
347 supportInitialGuess =
true;
349 }
else if (
type_ ==
"RELAXATION") {
350 paramList.set(
"relaxation: zero starting solution", InitialGuessIsZero);
352 supportInitialGuess =
true;
354 }
else if (
type_ ==
"KRYLOV") {
355 paramList.set(
"krylov: zero starting solution", InitialGuessIsZero);
357 supportInitialGuess =
true;
359 }
else if (
type_ ==
"SCHWARZ") {
360 paramList.set(
"schwarz: zero starting solution", InitialGuessIsZero);
365 prec_->setParameters(paramList);
366 supportInitialGuess =
true;
376 if (InitialGuessIsZero || supportInitialGuess) {
379 prec_->apply(tpB, tpX);
381 typedef Teuchos::ScalarTraits<Scalar> TST;
383 RCP<MultiVector> Correction = MultiVectorFactory::Build(
A_->getDomainMap(), X.getNumVectors());
388 prec_->apply(tpB, tpX);
390 X.update(TST::one(), *Correction, TST::one());
394 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
401 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
403 std::ostringstream out;
405 out <<
prec_->description();
408 out <<
"{type = " <<
type_ <<
"}";
413 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
418 out0 <<
"Prec. type: " <<
type_ << std::endl;
421 out0 <<
"Parameter list: " << std::endl;
422 Teuchos::OSTab tab2(out);
424 out0 <<
"Overlap: " <<
overlap_ << std::endl;
428 if (
prec_ != Teuchos::null) {
429 Teuchos::OSTab tab2(out);
430 out << *
prec_ << std::endl;
433 if (verbLevel &
Debug) {
436 <<
"RCP<prec_>: " <<
prec_ << std::endl;
442 #endif // HAVE_MUELU_TPETRA && HAVE_MUELU_IFPACK2 443 #endif // MUELU_IFPACK2SMOOTHER_DEF_HPP void DeclareInput(Level ¤tLevel) const
Input.
Important warning messages (one line)
Exception indicating invalid cast attempted.
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.
Print external lib objects.
friend class Ifpack2Smoother
Constructor.
Timer to be used in factories. Similar to Monitor but with additional timers.
RCP< Matrix > A_
matrix, used in apply if solving residual equation
Print additional debugging information.
One-liner description of what is happening.
void SetParameterList(const Teuchos::ParameterList ¶mList)
Namespace for MueLu classes and methods.
virtual const Teuchos::ParameterList & GetParameterList() const
LO overlap_
overlap when using the smoother in additive Schwarz mode
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
std::string type_
ifpack2-specific key phrase that denote smoother type
Print statistics that do not involve significant additional computation.
RCP< ParameterList > RemoveFactoriesFromList(const ParameterList &list) 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)
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
virtual void SetParameterList(const ParameterList ¶mList)
Set parameters from a parameter list and return with default values.
RCP< SmootherPrototype > Copy() const
Class that holds all level-specific information.
void Setup(Level ¤tLevel)
Set up the smoother.
std::string description() const
Return a simple one-line description of this object.
bool IsSetup() const
Get the state of a smoother prototype.
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)
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
Teuchos::FancyOStream & GetOStream(MsgType type, int thisProcRankOnly=0) const
Get an output stream for outputting the input message type.
void SetPrecParameters(const Teuchos::ParameterList &list=Teuchos::ParameterList()) const
MatrixType::scalar_type getLambdaMaxForApply() const
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the preconditioner.
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.
static RCP< Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraRow(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
void Input(Level &level, const std::string &varName) const
RCP< Ifpack2::Preconditioner< Scalar, LocalOrdinal, GlobalOrdinal, Node > > prec_
pointer to Ifpack2 preconditioner object
virtual std::string description() const
Return a simple one-line description of this object.