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_RILUK.hpp>
59#include <Ifpack2_Relaxation.hpp>
60#include <Ifpack2_ILUT.hpp>
61#include <Ifpack2_BlockRelaxation.hpp>
62#include <Ifpack2_Factory.hpp>
63#include <Ifpack2_Parameters.hpp>
65#include <Xpetra_BlockedCrsMatrix.hpp>
66#include <Xpetra_CrsMatrix.hpp>
67#include <Xpetra_CrsMatrixWrap.hpp>
68#include <Xpetra_TpetraBlockCrsMatrix.hpp>
69#include <Xpetra_Matrix.hpp>
70#include <Xpetra_MatrixMatrix.hpp>
71#include <Xpetra_MultiVectorFactory.hpp>
72#include <Xpetra_TpetraMultiVector.hpp>
74#include <Tpetra_BlockCrsMatrix_Helpers.hpp>
79#include "MueLu_Utilities.hpp"
81#include "MueLu_Aggregates.hpp"
84#ifdef HAVE_MUELU_INTREPID2
87#include "Intrepid2_Basis.hpp"
88#include "Kokkos_DynRankView.hpp"
94 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
96 : type_(type), overlap_(overlap)
98 typedef Tpetra::RowMatrix<SC,LO,GO,NO> tRowMatrix;
99 bool isSupported = Ifpack2::Factory::isSupported<tRowMatrix>(
type_) || (
type_ ==
"LINESMOOTHING_TRIDI_RELAXATION" ||
100 type_ ==
"LINESMOOTHING_TRIDI RELAXATION" ||
101 type_ ==
"LINESMOOTHING_TRIDIRELAXATION" ||
102 type_ ==
"LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
103 type_ ==
"LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
104 type_ ==
"LINESMOOTHING_TRIDIAGONALRELAXATION" ||
105 type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
106 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
107 type_ ==
"LINESMOOTHING_BANDEDRELAXATION" ||
108 type_ ==
"LINESMOOTHING_BLOCK_RELAXATION" ||
109 type_ ==
"LINESMOOTHING_BLOCK RELAXATION" ||
110 type_ ==
"LINESMOOTHING_BLOCKRELAXATION" ||
111 type_ ==
"TOPOLOGICAL" ||
112 type_ ==
"AGGREGATE");
118 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
129 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
131 std::string prefix = this->ShortClassName() +
": SetPrecParameters";
133 ParameterList& paramList =
const_cast<ParameterList&
>(this->GetParameterList());
135 paramList.setParameters(list);
137 RCP<ParameterList> precList = this->RemoveFactoriesFromList(this->GetParameterList());
139 if(!precList.is_null() && precList->isParameter(
"partitioner: type") && precList->get<std::string>(
"partitioner: type") ==
"linear" &&
140 !precList->isParameter(
"partitioner: local parts")) {
141 LO matrixBlockSize = 1;
142 int lclSize = A_->getRangeMap()->getLocalNumElements();
143 RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_);
144 if (!matA.is_null()) {
145 lclSize = matA->getLocalNumRows();
146 matrixBlockSize = matA->GetFixedBlockSize();
148 precList->set(
"partitioner: local parts", lclSize / matrixBlockSize);
151 prec_->setParameters(*precList);
153 paramList.setParameters(*precList);
156 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
158 this->Input(currentLevel,
"A");
160 if (type_ ==
"LINESMOOTHING_TRIDI_RELAXATION" ||
161 type_ ==
"LINESMOOTHING_TRIDI RELAXATION" ||
162 type_ ==
"LINESMOOTHING_TRIDIRELAXATION" ||
163 type_ ==
"LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
164 type_ ==
"LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
165 type_ ==
"LINESMOOTHING_TRIDIAGONALRELAXATION" ||
166 type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
167 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
168 type_ ==
"LINESMOOTHING_BANDEDRELAXATION" ||
169 type_ ==
"LINESMOOTHING_BLOCK_RELAXATION" ||
170 type_ ==
"LINESMOOTHING_BLOCK RELAXATION" ||
171 type_ ==
"LINESMOOTHING_BLOCKRELAXATION") {
172 this->Input(currentLevel,
"CoarseNumZLayers");
173 this->Input(currentLevel,
"LineDetection_VertLineIds");
175 else if (type_ ==
"BLOCK RELAXATION" ||
176 type_ ==
"BLOCK_RELAXATION" ||
177 type_ ==
"BLOCKRELAXATION" ||
179 type_ ==
"BANDED_RELAXATION" ||
180 type_ ==
"BANDED RELAXATION" ||
181 type_ ==
"BANDEDRELAXATION" ||
183 type_ ==
"TRIDI_RELAXATION" ||
184 type_ ==
"TRIDI RELAXATION" ||
185 type_ ==
"TRIDIRELAXATION" ||
186 type_ ==
"TRIDIAGONAL_RELAXATION" ||
187 type_ ==
"TRIDIAGONAL RELAXATION" ||
188 type_ ==
"TRIDIAGONALRELAXATION")
191 ParameterList precList = this->GetParameterList();
192 if(precList.isParameter(
"partitioner: type") &&
193 precList.get<std::string>(
"partitioner: type") ==
"line") {
194 this->Input(currentLevel,
"Coordinates");
197 else if (type_ ==
"TOPOLOGICAL")
200 this->Input(currentLevel,
"pcoarsen: element to node map");
202 else if (type_ ==
"AGGREGATE")
205 this->Input(currentLevel,
"Aggregates");
207 else if (type_ ==
"HIPTMAIR") {
209 this->Input(currentLevel,
"NodeMatrix");
210 this->Input(currentLevel,
"D0");
215 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
218 A_ = Factory::Get< RCP<Operator> >(currentLevel,
"A");
219 ParameterList& paramList =
const_cast<ParameterList&
>(this->GetParameterList());
223 if(paramList.isParameter(
"smoother: use blockcrsmatrix storage") && paramList.get<
bool>(
"smoother: use blockcrsmatrix storage")) {
225 RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_);
227 blocksize = matA->GetFixedBlockSize();
231 RCP<CrsMatrixWrap> AcrsWrap = rcp_dynamic_cast<CrsMatrixWrap>(A_);
232 if(AcrsWrap.is_null())
233 throw std::runtime_error(
"Ifpack2Smoother: Cannot convert matrix A to CrsMatrixWrap object.");
234 RCP<CrsMatrix> Acrs = AcrsWrap->getCrsMatrix();
236 throw std::runtime_error(
"Ifpack2Smoother: Cannot extract CrsMatrix from matrix A.");
237 RCP<TpetraCrsMatrix> At = rcp_dynamic_cast<TpetraCrsMatrix>(Acrs);
239 if(!Xpetra::Helpers<Scalar,LO,GO,Node>::isTpetraBlockCrs(matA))
240 throw std::runtime_error(
"Ifpack2Smoother: Cannot extract CrsMatrix or BlockCrsMatrix from matrix A.");
241 this->GetOStream(
Statistics0) <<
"Ifpack2Smoother: Using (native) BlockCrsMatrix storage with blocksize "<<blocksize<<std::endl;
242 paramList.remove(
"smoother: use blockcrsmatrix storage");
245 RCP<Tpetra::BlockCrsMatrix<Scalar, LO, GO, Node> > blockCrs = Tpetra::convertToBlockCrsMatrix(*At->getTpetra_CrsMatrix(),blocksize);
246 RCP<CrsMatrix> blockCrs_as_crs = rcp(
new TpetraBlockCrsMatrix(blockCrs));
247 RCP<CrsMatrixWrap> blockWrap = rcp(
new CrsMatrixWrap(blockCrs_as_crs));
249 this->GetOStream(
Statistics0) <<
"Ifpack2Smoother: Using BlockCrsMatrix storage with blocksize "<<blocksize<<std::endl;
251 paramList.remove(
"smoother: use blockcrsmatrix storage");
256 if (type_ ==
"SCHWARZ")
257 SetupSchwarz(currentLevel);
259 else if (type_ ==
"LINESMOOTHING_TRIDI_RELAXATION" ||
260 type_ ==
"LINESMOOTHING_TRIDI RELAXATION" ||
261 type_ ==
"LINESMOOTHING_TRIDIRELAXATION" ||
262 type_ ==
"LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
263 type_ ==
"LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
264 type_ ==
"LINESMOOTHING_TRIDIAGONALRELAXATION" ||
265 type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
266 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
267 type_ ==
"LINESMOOTHING_BANDEDRELAXATION" ||
268 type_ ==
"LINESMOOTHING_BLOCK_RELAXATION" ||
269 type_ ==
"LINESMOOTHING_BLOCK RELAXATION" ||
270 type_ ==
"LINESMOOTHING_BLOCKRELAXATION")
271 SetupLineSmoothing(currentLevel);
273 else if (type_ ==
"BLOCK_RELAXATION" ||
274 type_ ==
"BLOCK RELAXATION" ||
275 type_ ==
"BLOCKRELAXATION" ||
277 type_ ==
"BANDED_RELAXATION" ||
278 type_ ==
"BANDED RELAXATION" ||
279 type_ ==
"BANDEDRELAXATION" ||
281 type_ ==
"TRIDI_RELAXATION" ||
282 type_ ==
"TRIDI RELAXATION" ||
283 type_ ==
"TRIDIRELAXATION" ||
284 type_ ==
"TRIDIAGONAL_RELAXATION" ||
285 type_ ==
"TRIDIAGONAL RELAXATION" ||
286 type_ ==
"TRIDIAGONALRELAXATION")
287 SetupBlockRelaxation(currentLevel);
289 else if (type_ ==
"CHEBYSHEV")
290 SetupChebyshev(currentLevel);
292 else if (type_ ==
"TOPOLOGICAL")
294#ifdef HAVE_MUELU_INTREPID2
295 SetupTopological(currentLevel);
297 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"'TOPOLOGICAL' smoother choice requires Intrepid2");
300 else if (type_ ==
"AGGREGATE")
301 SetupAggregate(currentLevel);
303 else if (type_ ==
"HIPTMAIR")
304 SetupHiptmair(currentLevel);
308 SetupGeneric(currentLevel);
313 this->GetOStream(
Statistics1) << description() << std::endl;
316 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
318 typedef Tpetra::RowMatrix<SC,LO,GO,NO> tRowMatrix;
320 bool reusePreconditioner =
false;
321 if (this->IsSetup() ==
true) {
323 this->GetOStream(
Runtime1) <<
"MueLu::Ifpack2Smoother::SetupSchwarz(): Setup() has already been called, assuming reuse" << std::endl;
325 bool isTRowMatrix =
true;
326 RCP<const tRowMatrix> tA;
330 isTRowMatrix =
false;
333 RCP<Ifpack2::Details::CanChangeMatrix<tRowMatrix> > prec = rcp_dynamic_cast<Ifpack2::Details::CanChangeMatrix<tRowMatrix> >(prec_);
334 if (!prec.is_null() && isTRowMatrix) {
336 reusePreconditioner =
true;
338 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupSchwarz(): reuse of this type is not available "
339 "(either failed cast to CanChangeMatrix, or to Tpetra Row Matrix), reverting to full construction" << std::endl;
343 if (!reusePreconditioner) {
344 ParameterList& paramList =
const_cast<ParameterList&
>(this->GetParameterList());
346 bool isBlockedMatrix =
false;
347 RCP<Matrix> merged2Mat;
349 std::string sublistName =
"subdomain solver parameters";
350 if (paramList.isSublist(sublistName)) {
359 ParameterList& subList = paramList.sublist(sublistName);
361 std::string partName =
"partitioner: type";
367 if (subList.isParameter(partName) && subList.get<std::string>(partName) ==
"vanka user") {
368 isBlockedMatrix =
true;
370 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
372 "Matrix A must be of type BlockedCrsMatrix.");
374 size_t numVels = bA->getMatrix(0,0)->getLocalNumRows();
375 size_t numPres = bA->getMatrix(1,0)->getLocalNumRows();
376 size_t numRows = rcp_dynamic_cast<Matrix>(A_,
true)->getLocalNumRows();
378 ArrayRCP<LocalOrdinal> blockSeeds(numRows, Teuchos::OrdinalTraits<LocalOrdinal>::invalid());
380 size_t numBlocks = 0;
381 for (
size_t rowOfB = numVels; rowOfB < numVels+numPres; ++rowOfB)
382 blockSeeds[rowOfB] = numBlocks++;
384 RCP<BlockedCrsMatrix> bA2 = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
386 "Matrix A must be of type BlockedCrsMatrix.");
388 merged2Mat = bA2->Merge();
392 bool haveBoundary =
false;
393 for (LO i = 0; i < boundaryNodes.size(); i++)
394 if (boundaryNodes[i]) {
398 blockSeeds[i] = numBlocks;
404 subList.set(
"partitioner: type",
"user");
405 subList.set(
"partitioner: map", blockSeeds);
406 subList.set(
"partitioner: local parts", as<int>(numBlocks));
409 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
411 isBlockedMatrix =
true;
412 merged2Mat = bA->Merge();
417 RCP<const tRowMatrix> tA;
432 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
435 ParameterList& paramList =
const_cast<ParameterList&
>(this->GetParameterList());
437 if (this->IsSetup() ==
true) {
438 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupAggregate(): Setup() has already been called" << std::endl;
439 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupAggregate(): reuse of this type is not available, reverting to full construction" << std::endl;
442 this->GetOStream(
Statistics0) <<
"Ifpack2Smoother: Using Aggregate Smoothing"<<std::endl;
444 RCP<Aggregates> aggregates = Factory::Get<RCP<Aggregates> >(currentLevel,
"Aggregates");
446 RCP<const LOMultiVector> vertex2AggId = aggregates->GetVertex2AggId();
447 ArrayRCP<LO> aggregate_ids = rcp_const_cast<LOMultiVector>(vertex2AggId)->getDataNonConst(0);
448 ArrayRCP<LO> dof_ids;
452 RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_);
454 blocksize = matA->GetFixedBlockSize();
456 dof_ids.resize(aggregate_ids.size() * blocksize);
457 for(LO i=0; i<(LO)aggregate_ids.size(); i++) {
458 for(LO j=0; j<(LO)blocksize; j++)
459 dof_ids[i*blocksize+j] = aggregate_ids[i];
463 dof_ids = aggregate_ids;
467 paramList.set(
"partitioner: map", dof_ids);
468 paramList.set(
"partitioner: type",
"user");
469 paramList.set(
"partitioner: overlap", 0);
470 paramList.set(
"partitioner: local parts", (
int)aggregates->GetNumAggregates());
474 type_ =
"BLOCKRELAXATION";
481#ifdef HAVE_MUELU_INTREPID2
482 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
492 if (this->IsSetup() ==
true) {
493 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupTopological(): Setup() has already been called" << std::endl;
494 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupTopological(): reuse of this type is not available, reverting to full construction" << std::endl;
497 ParameterList& paramList =
const_cast<ParameterList&
>(this->GetParameterList());
499 typedef typename Node::device_type::execution_space ES;
501 typedef Kokkos::DynRankView<LocalOrdinal,typename Node::device_type> FCO;
503 LocalOrdinal lo_invalid = Teuchos::OrdinalTraits<LO>::invalid();
507 const Teuchos::RCP<FCO> elemToNode = Factory::Get<Teuchos::RCP<FCO> >(currentLevel,
"pcoarsen: element to node map");
509 string basisString = paramList.get<
string>(
"pcoarsen: hi basis");
515 auto basis = MueLuIntrepid::BasisFactory<double,ES>(basisString, degree);
517 string topologyTypeString = paramList.get<
string>(
"smoother: neighborhood type");
519 if (topologyTypeString ==
"node")
521 else if (topologyTypeString ==
"edge")
523 else if (topologyTypeString ==
"face")
525 else if (topologyTypeString ==
"cell")
526 dimension = basis->getBaseCellTopology().getDimension();
528 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::invalid_argument,
"Unrecognized smoother neighborhood type. Supported types are node, edge, face.");
529 RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_,
true);
530 vector<vector<LocalOrdinal>> seeds;
535 int myNodeCount = matA->getRowMap()->getLocalNumElements();
536 ArrayRCP<LocalOrdinal> nodeSeeds(myNodeCount,lo_invalid);
537 int localPartitionNumber = 0;
540 nodeSeeds[seed] = localPartitionNumber++;
543 paramList.remove(
"smoother: neighborhood type");
544 paramList.remove(
"pcoarsen: hi basis");
546 paramList.set(
"partitioner: map", nodeSeeds);
547 paramList.set(
"partitioner: type",
"user");
548 paramList.set(
"partitioner: overlap", 1);
549 paramList.set(
"partitioner: local parts",
int(seeds[dimension].size()));
553 type_ =
"BLOCKRELAXATION";
561 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
563 if (this->IsSetup() ==
true) {
564 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupLineSmoothing(): Setup() has already been called" << std::endl;
565 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupLineSmoothing(): reuse of this type is not available, reverting to full construction" << std::endl;
568 ParameterList& myparamList =
const_cast<ParameterList&
>(this->GetParameterList());
570 LO CoarseNumZLayers = Factory::Get<LO>(currentLevel,
"CoarseNumZLayers");
571 if (CoarseNumZLayers > 0) {
572 Teuchos::ArrayRCP<LO> TVertLineIdSmoo = Factory::Get< Teuchos::ArrayRCP<LO> >(currentLevel,
"LineDetection_VertLineIds");
576 for(
size_t k = 0; k < Teuchos::as<size_t>(TVertLineIdSmoo.size()); k++) {
577 if(maxPart < TVertLineIdSmoo[k]) maxPart = TVertLineIdSmoo[k];
579 RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_,
true);
580 size_t numLocalRows = matA->getLocalNumRows();
583 "MueLu::Ifpack2Smoother::Setup(): the number of local nodes is incompatible with the TVertLineIdsSmoo.");
588 int actualDofsPerNode = numLocalRows / TVertLineIdSmoo.size();
589 LO matrixBlockSize = matA->GetFixedBlockSize();
590 if(matrixBlockSize > 1 && actualDofsPerNode > 1)
593 "MueLu::Ifpack2Smoother::Setup(): A is a block matrix but its block size and DOFs/node from partitioner disagree");
595 else if(matrixBlockSize > 1)
597 actualDofsPerNode = matrixBlockSize;
599 myparamList.set(
"partitioner: PDE equations", actualDofsPerNode);
601 if (numLocalRows == Teuchos::as<size_t>(TVertLineIdSmoo.size())) {
602 myparamList.set(
"partitioner: type",
"user");
603 myparamList.set(
"partitioner: map",TVertLineIdSmoo);
604 myparamList.set(
"partitioner: local parts",maxPart+1);
607 size_t numDofsPerNode = numLocalRows / TVertLineIdSmoo.size();
610 Teuchos::ArrayRCP<LO> partitionerMap(numLocalRows, Teuchos::OrdinalTraits<LocalOrdinal>::invalid());
611 for (
size_t blockRow = 0; blockRow < Teuchos::as<size_t>(TVertLineIdSmoo.size()); ++blockRow)
612 for (
size_t dof = 0; dof < numDofsPerNode; dof++)
613 partitionerMap[blockRow * numDofsPerNode + dof] = TVertLineIdSmoo[blockRow];
614 myparamList.set(
"partitioner: type",
"user");
615 myparamList.set(
"partitioner: map",partitionerMap);
616 myparamList.set(
"partitioner: local parts",maxPart + 1);
619 if (type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
620 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
621 type_ ==
"LINESMOOTHING_BANDEDRELAXATION")
622 type_ =
"BANDEDRELAXATION";
623 else if (type_ ==
"LINESMOOTHING_TRIDI_RELAXATION" ||
624 type_ ==
"LINESMOOTHING_TRIDI RELAXATION" ||
625 type_ ==
"LINESMOOTHING_TRIDIRELAXATION" ||
626 type_ ==
"LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
627 type_ ==
"LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
628 type_ ==
"LINESMOOTHING_TRIDIAGONALRELAXATION")
629 type_ =
"TRIDIAGONALRELAXATION";
631 type_ =
"BLOCKRELAXATION";
634 this->GetOStream(
Runtime0) <<
"Line detection failed: fall back to point-wise relaxation" << std::endl;
635 myparamList.remove(
"partitioner: type",
false);
636 myparamList.remove(
"partitioner: map",
false);
637 myparamList.remove(
"partitioner: local parts",
false);
638 type_ =
"RELAXATION";
649 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
651 typedef Tpetra::RowMatrix<SC,LO,GO,NO> tRowMatrix;
653 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
659 bool reusePreconditioner =
false;
660 if (this->IsSetup() ==
true) {
662 this->GetOStream(
Runtime1) <<
"MueLu::Ifpack2Smoother::SetupBlockRelaxation(): Setup() has already been called, assuming reuse" << std::endl;
664 RCP<Ifpack2::Details::CanChangeMatrix<tRowMatrix> > prec = rcp_dynamic_cast<Ifpack2::Details::CanChangeMatrix<tRowMatrix> >(prec_);
665 if (!prec.is_null()) {
667 reusePreconditioner =
true;
669 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupBlockRelaxation(): reuse of this type is not available (failed cast to CanChangeMatrix), "
670 "reverting to full construction" << std::endl;
674 if (!reusePreconditioner) {
675 ParameterList& myparamList =
const_cast<ParameterList&
>(this->GetParameterList());
677 if(myparamList.isParameter(
"partitioner: type") &&
678 myparamList.get<std::string>(
"partitioner: type") ==
"line") {
679 Teuchos::RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> > xCoordinates =
680 Factory::Get<Teuchos::RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> > >(currentLevel,
"Coordinates");
681 Teuchos::RCP<Tpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> > coordinates = Teuchos::rcpFromRef(Xpetra::toTpetra<
typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>(*xCoordinates));
683 RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_);
684 size_t lclSize = A_->getRangeMap()->getLocalNumElements();
686 lclSize = matA->getLocalNumRows();
687 size_t numDofsPerNode = lclSize / xCoordinates->getMap()->getLocalNumElements();
688 myparamList.set(
"partitioner: coordinates", coordinates);
689 myparamList.set(
"partitioner: PDE equations", (
int) numDofsPerNode);
700 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
702 typedef Tpetra::RowMatrix<SC,LO,GO,NO> tRowMatrix;
703 using STS = Teuchos::ScalarTraits<SC>;
704 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
710 bool reusePreconditioner =
false;
712 if (this->IsSetup() ==
true) {
714 this->GetOStream(
Runtime1) <<
"MueLu::Ifpack2Smoother::SetupChebyshev(): Setup() has already been called, assuming reuse" << std::endl;
716 RCP<Ifpack2::Details::CanChangeMatrix<tRowMatrix> > prec = rcp_dynamic_cast<Ifpack2::Details::CanChangeMatrix<tRowMatrix> >(prec_);
717 if (!prec.is_null()) {
719 reusePreconditioner =
true;
721 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupChebyshev(): reuse of this type is not available (failed cast to CanChangeMatrix), "
722 "reverting to full construction" << std::endl;
727 ParameterList& paramList =
const_cast<ParameterList&
>(this->GetParameterList());
728 SC negone = -STS::one();
729 SC lambdaMax = SetupChebyshevEigenvalues(currentLevel,
"A",
"",paramList);
732 if (!reusePreconditioner) {
747 if (lambdaMax == negone) {
748 typedef Tpetra::RowMatrix<SC, LO, GO, NO> MatrixType;
750 Teuchos::RCP<Ifpack2::Chebyshev<MatrixType> > chebyPrec = rcp_dynamic_cast<Ifpack2::Chebyshev<MatrixType> >(prec_);
751 if (chebyPrec != Teuchos::null) {
752 lambdaMax = chebyPrec->getLambdaMaxForApply();
753 RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_);
755 matA->SetMaxEigenvalueEstimate(lambdaMax);
756 this->GetOStream(
Statistics1) <<
"chebyshev: max eigenvalue (calculated by Ifpack2)" <<
" = " << lambdaMax << std::endl;
758 TEUCHOS_TEST_FOR_EXCEPTION(lambdaMax == negone,
Exceptions::RuntimeError,
"MueLu::Ifpack2Smoother::Setup(): no maximum eigenvalue estimate");
763 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
765 typedef Tpetra::RowMatrix<SC,LO,GO,NO> tRowMatrix;
766 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
772 bool reusePreconditioner =
false;
773 if (this->IsSetup() ==
true) {
775 this->GetOStream(
Runtime1) <<
"MueLu::Ifpack2Smoother::SetupHiptmair(): Setup() has already been called, assuming reuse" << std::endl;
777 RCP<Ifpack2::Details::CanChangeMatrix<tRowMatrix> > prec = rcp_dynamic_cast<Ifpack2::Details::CanChangeMatrix<tRowMatrix> >(prec_);
778 if (!prec.is_null()) {
780 reusePreconditioner =
true;
782 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupHiptmair(): reuse of this type is not available (failed cast to CanChangeMatrix), "
783 "reverting to full construction" << std::endl;
788 ParameterList& paramList =
const_cast<ParameterList&
>(this->GetParameterList());
789 std::string smoother1 = paramList.get(
"hiptmair: smoother type 1",
"CHEBYSHEV");
790 std::string smoother2 = paramList.get(
"hiptmair: smoother type 2",
"CHEBYSHEV");
793 if(smoother1 ==
"CHEBYSHEV") {
794 ParameterList & list1 = paramList.sublist(
"hiptmair: smoother list 1");
796 SetupChebyshevEigenvalues(currentLevel,
"A",
"EdgeMatrix ",list1);
798 if(smoother2 ==
"CHEBYSHEV") {
799 ParameterList & list2 = paramList.sublist(
"hiptmair: smoother list 2");
801 SetupChebyshevEigenvalues(currentLevel,
"A",
"EdgeMatrix ",list2);
808 RCP<Operator> NodeMatrix = currentLevel.
Get<RCP<Operator> >(
"NodeMatrix");
809 RCP<Operator> D0 = currentLevel.
Get<RCP<Operator> >(
"D0");
815 Teuchos::ParameterList newlist;
816 newlist.set(
"P",tD0);
817 newlist.set(
"PtAP",tNodeMatrix);
818 if (!reusePreconditioner) {
820 SetPrecParameters(newlist);
829 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
832 typedef Teuchos::ScalarTraits<SC> STS;
833 SC negone = -STS::one();
834 RCP<const Matrix> currentA = currentLevel.
Get<RCP<Matrix> >(matrixName);
835 SC lambdaMax = negone;
837 std::string maxEigString =
"chebyshev: max eigenvalue";
838 std::string eigRatioString =
"chebyshev: ratio eigenvalue";
841 if (paramList.isParameter(maxEigString)) {
842 if (paramList.isType<
double>(maxEigString))
843 lambdaMax = paramList.get<
double>(maxEigString);
845 lambdaMax = paramList.get<SC>(maxEigString);
846 this->GetOStream(
Statistics1) << label << maxEigString <<
" (cached with smoother parameter list) = " << lambdaMax << std::endl;
847 RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_);
849 matA->SetMaxEigenvalueEstimate(lambdaMax);
852 lambdaMax = currentA->GetMaxEigenvalueEstimate();
853 if (lambdaMax != negone) {
854 this->GetOStream(
Statistics1) << label << maxEigString <<
" (cached with matrix) = " << lambdaMax << std::endl;
855 paramList.set(maxEigString, lambdaMax);
860 const SC defaultEigRatio = 20;
862 SC ratio = defaultEigRatio;
863 if (paramList.isParameter(eigRatioString)) {
864 if (paramList.isType<
double>(eigRatioString))
865 ratio = paramList.get<
double>(eigRatioString);
867 ratio = paramList.get<SC>(eigRatioString);
874 RCP<const Matrix> fineA = currentLevel.
GetPreviousLevel()->Get<RCP<Matrix> >(matrixName);
875 size_t nRowsFine = fineA->getGlobalNumRows();
876 size_t nRowsCoarse = currentA->getGlobalNumRows();
878 SC levelRatio = as<SC>(as<float>(nRowsFine)/nRowsCoarse);
879 if (STS::magnitude(levelRatio) > STS::magnitude(ratio))
883 this->GetOStream(
Statistics1) << label << eigRatioString <<
" (computed) = " << ratio << std::endl;
884 paramList.set(eigRatioString, ratio);
886 if (paramList.isParameter(
"chebyshev: use rowsumabs diagonal scaling")) {
887 this->GetOStream(
Runtime1) <<
"chebyshev: using rowsumabs diagonal scaling" << std::endl;
888 bool doScale =
false;
889 doScale = paramList.get<
bool>(
"chebyshev: use rowsumabs diagonal scaling");
890 paramList.remove(
"chebyshev: use rowsumabs diagonal scaling");
891 double chebyReplaceTol = Teuchos::ScalarTraits<Scalar>::eps()*100;
892 std::string paramName =
"chebyshev: rowsumabs diagonal replacement tolerance";
893 if (paramList.isParameter(paramName)) {
894 chebyReplaceTol = paramList.get<
double>(paramName);
895 paramList.remove(paramName);
897 double chebyReplaceVal = Teuchos::ScalarTraits<double>::zero();
898 paramName =
"chebyshev: rowsumabs diagonal replacement value";
899 if (paramList.isParameter(paramName)) {
900 chebyReplaceVal = paramList.get<
double>(paramName);
901 paramList.remove(paramName);
903 bool chebyReplaceSingleEntryRowWithZero =
false;
904 paramName =
"chebyshev: rowsumabs replace single entry row with zero";
905 if (paramList.isParameter(paramName)) {
906 chebyReplaceSingleEntryRowWithZero = paramList.get<
bool>(paramName);
907 paramList.remove(paramName);
909 bool useAverageAbsDiagVal =
false;
910 paramName =
"chebyshev: rowsumabs use automatic diagonal tolerance";
911 if (paramList.isParameter(paramName)) {
912 useAverageAbsDiagVal = paramList.get<
bool>(paramName);
913 paramList.remove(paramName);
916 const bool doReciprocal =
true;
917 RCP<Vector> lumpedDiagonal =
Utilities::GetLumpedMatrixDiagonal(*currentA, doReciprocal, chebyReplaceTol, chebyReplaceVal, chebyReplaceSingleEntryRowWithZero, useAverageAbsDiagVal);
918 const Xpetra::TpetraVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmpVec =
dynamic_cast<const Xpetra::TpetraVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>&
>(*lumpedDiagonal);
919 paramList.set(
"chebyshev: operator inv diagonal",tmpVec.getTpetra_Vector());
929 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
931 typedef Tpetra::RowMatrix<SC,LO,GO,NO> tRowMatrix;
932 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
938 bool reusePreconditioner =
false;
939 if (this->IsSetup() ==
true) {
941 this->GetOStream(
Runtime1) <<
"MueLu::Ifpack2Smoother::SetupGeneric(): Setup() has already been called, assuming reuse" << std::endl;
943 RCP<Ifpack2::Details::CanChangeMatrix<tRowMatrix> > prec = rcp_dynamic_cast<Ifpack2::Details::CanChangeMatrix<tRowMatrix> >(prec_);
944 if (!prec.is_null()) {
946 reusePreconditioner =
true;
948 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupGeneric(): reuse of this type is not available (failed cast to CanChangeMatrix), "
949 "reverting to full construction" << std::endl;
953 if (!reusePreconditioner) {
962 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
972 Teuchos::ParameterList paramList;
973 bool supportInitialGuess =
false;
974 const Teuchos::ParameterList params = this->GetParameterList();
976 if (prec_->supportsZeroStartingSolution()) {
977 prec_->setZeroStartingSolution(InitialGuessIsZero);
978 supportInitialGuess =
true;
979 }
else if (type_ ==
"SCHWARZ") {
980 paramList.set(
"schwarz: zero starting solution", InitialGuessIsZero);
985 prec_->setParameters(paramList);
986 supportInitialGuess =
true;
996 if (InitialGuessIsZero || supportInitialGuess) {
999 prec_->apply(tpB, tpX);
1001 typedef Teuchos::ScalarTraits<Scalar> TST;
1003 RCP<MultiVector> Residual;
1005 std::string prefix = this->ShortClassName() +
": Apply: ";
1006 RCP<TimeMonitor> tM = rcp(
new TimeMonitor(*
this, prefix +
"residual calculation",
Timings0));
1010 RCP<MultiVector> Correction = MultiVectorFactory::Build(A_->getDomainMap(), X.getNumVectors());
1015 prec_->apply(tpB, tpX);
1017 X.update(TST::one(), *Correction, TST::one());
1021 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1024 smoother->SetParameterList(this->GetParameterList());
1028 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1030 std::ostringstream out;
1032 out << prec_->description();
1035 out <<
"{type = " << type_ <<
"}";
1040 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1045 out0 <<
"Prec. type: " << type_ << std::endl;
1048 out0 <<
"Parameter list: " << std::endl;
1049 Teuchos::OSTab tab2(out);
1050 out << this->GetParameterList();
1051 out0 <<
"Overlap: " << overlap_ << std::endl;
1055 if (prec_ != Teuchos::null) {
1056 Teuchos::OSTab tab2(out);
1057 out << *prec_ << std::endl;
1060 if (verbLevel &
Debug) {
1063 <<
"RCP<prec_>: " << prec_ << std::endl;
1067 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1069 typedef Tpetra::RowMatrix<SC,LO,GO,NO> MatrixType;
1071 RCP<Ifpack2::Relaxation<MatrixType> > pr = rcp_dynamic_cast<Ifpack2::Relaxation<MatrixType> >(prec_);
1072 if(!pr.is_null())
return pr->getNodeSmootherComplexity();
1074 RCP<Ifpack2::Chebyshev<MatrixType> > pc = rcp_dynamic_cast<Ifpack2::Chebyshev<MatrixType> >(prec_);
1075 if(!pc.is_null())
return pc->getNodeSmootherComplexity();
1077 RCP<Ifpack2::BlockRelaxation<MatrixType> > pb = rcp_dynamic_cast<Ifpack2::BlockRelaxation<MatrixType> >(prec_);
1078 if(!pb.is_null())
return pb->getNodeSmootherComplexity();
1080 RCP<Ifpack2::ILUT<MatrixType> > pi = rcp_dynamic_cast<Ifpack2::ILUT<MatrixType> >(prec_);
1081 if(!pi.is_null())
return pi->getNodeSmootherComplexity();
1083 RCP<Ifpack2::RILUK<MatrixType> > pk = rcp_dynamic_cast<Ifpack2::RILUK<MatrixType> >(prec_);
1084 if(!pk.is_null())
return pk->getNodeSmootherComplexity();
1087 return Teuchos::OrdinalTraits<size_t>::invalid();
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultScalar Scalar
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)
virtual std::string description() const
Return a simple one-line description of this object.
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that encapsulates Ifpack2 smoothers.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
void SetupSchwarz(Level ¤tLevel)
void Setup(Level ¤tLevel)
Set up the smoother.
Scalar SetupChebyshevEigenvalues(Level ¤tLevel, const std::string &matrixName, const std::string &label, Teuchos::ParameterList ¶mList) const
void SetupAggregate(Level ¤tLevel)
void SetupBlockRelaxation(Level ¤tLevel)
void SetupChebyshev(Level ¤tLevel)
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
void SetParameterList(const Teuchos::ParameterList ¶mList)
Set parameters from a parameter list and return with default values.
RCP< SmootherPrototype > Copy() const
std::string description() const
Return a simple one-line description of this object.
void SetupHiptmair(Level ¤tLevel)
void SetPrecParameters(const Teuchos::ParameterList &list=Teuchos::ParameterList()) const
void SetupTopological(Level ¤tLevel)
std::string type_
ifpack2-specific key phrase that denote smoother type
friend class Ifpack2Smoother
Constructor.
void DeclareInput(Level ¤tLevel) const
Input.
void SetupLineSmoothing(Level ¤tLevel)
void SetupGeneric(Level ¤tLevel)
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the preconditioner.
Class that holds all level-specific information.
RCP< Level > & GetPreviousLevel()
Previous level.
int GetLevelID() const
Return level number.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
virtual void SetParameterList(const Teuchos::ParameterList ¶mList)=0
Set parameters from a parameter list and return with default values.
void declareConstructionOutcome(bool fail, std::string msg)
bool IsSetup() const
Get the state of a smoother prototype.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Integrates Teuchos::TimeMonitor with MueLu verbosity system.
static RCP< Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraRow(RCP< Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static Teuchos::RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetLumpedMatrixDiagonal(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > const &A, const bool doReciprocal=false, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100, Scalar tolReplacement=Teuchos::ScalarTraits< Scalar >::zero(), const bool replaceSingleEntryRowWithZero=false, const bool useAverageAbsDiagVal=false)
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.
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
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)
static Teuchos::ArrayRCP< const bool > DetectDirichletRows(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::magnitude(0.), const bool count_twos_as_dirichlet=false)
void FindGeometricSeedOrdinals(Teuchos::RCP< Basis > basis, const LOFieldContainer &elementToNodeMap, std::vector< std::vector< LocalOrdinal > > &seeds, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &rowMap, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &columnMap)
Namespace for MueLu classes and methods.
@ Warnings0
Important warning messages (one line)
@ Debug
Print additional debugging information.
@ Statistics1
Print more statistics.
@ External
Print external lib objects.
@ Timings0
High level timing information (use Teuchos::TimeMonitor::summarize() to print)
@ Runtime0
One-liner description of what is happening.
@ Runtime1
Description of what is happening (more verbose)
@ Parameters0
Print class parameters.
@ Statistics0
Print statistics that do not involve significant additional computation.
@ Parameters1
Print class parameters (more parameters, more verbose)