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_MultiVectorFactory.hpp> 71 #include <Xpetra_TpetraMultiVector.hpp> 73 #include <Tpetra_BlockCrsMatrix_Helpers.hpp> 78 #include "MueLu_Utilities.hpp" 80 #include "MueLu_Aggregates.hpp" 83 #ifdef HAVE_MUELU_INTREPID2 86 #include "Intrepid2_Basis.hpp" 93 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
95 : type_(type), overlap_(overlap)
97 typedef Tpetra::RowMatrix<SC,LO,GO,NO> tRowMatrix;
98 bool isSupported = Ifpack2::Factory::isSupported<tRowMatrix>(
type_) || (
type_ ==
"LINESMOOTHING_TRIDI_RELAXATION" ||
99 type_ ==
"LINESMOOTHING_TRIDI RELAXATION" ||
100 type_ ==
"LINESMOOTHING_TRIDIRELAXATION" ||
101 type_ ==
"LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
102 type_ ==
"LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
103 type_ ==
"LINESMOOTHING_TRIDIAGONALRELAXATION" ||
104 type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
105 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
106 type_ ==
"LINESMOOTHING_BANDEDRELAXATION" ||
107 type_ ==
"LINESMOOTHING_BLOCK_RELAXATION" ||
108 type_ ==
"LINESMOOTHING_BLOCK RELAXATION" ||
109 type_ ==
"LINESMOOTHING_BLOCKRELAXATION" ||
110 type_ ==
"TOPOLOGICAL" ||
111 type_ ==
"AGGREGATE");
117 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
128 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
130 std::string prefix = this->ShortClassName() +
": SetPrecParameters";
132 ParameterList& paramList =
const_cast<ParameterList&
>(this->GetParameterList());
134 paramList.setParameters(list);
136 RCP<ParameterList> precList = this->RemoveFactoriesFromList(this->GetParameterList());
138 if(!precList.is_null() && precList->isParameter(
"partitioner: type") && precList->get<std::string>(
"partitioner: type") ==
"linear" &&
139 !precList->isParameter(
"partitioner: local parts")) {
140 precList->set(
"partitioner: local parts", (
int)A_->getNodeNumRows() / A_->GetFixedBlockSize());
143 prec_->setParameters(*precList);
145 paramList.setParameters(*precList);
148 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
150 this->Input(currentLevel,
"A");
152 if (type_ ==
"LINESMOOTHING_TRIDI_RELAXATION" ||
153 type_ ==
"LINESMOOTHING_TRIDI RELAXATION" ||
154 type_ ==
"LINESMOOTHING_TRIDIRELAXATION" ||
155 type_ ==
"LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
156 type_ ==
"LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
157 type_ ==
"LINESMOOTHING_TRIDIAGONALRELAXATION" ||
158 type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
159 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
160 type_ ==
"LINESMOOTHING_BANDEDRELAXATION" ||
161 type_ ==
"LINESMOOTHING_BLOCK_RELAXATION" ||
162 type_ ==
"LINESMOOTHING_BLOCK RELAXATION" ||
163 type_ ==
"LINESMOOTHING_BLOCKRELAXATION") {
164 this->Input(currentLevel,
"CoarseNumZLayers");
165 this->Input(currentLevel,
"LineDetection_VertLineIds");
167 else if (type_ ==
"BLOCK RELAXATION" ||
168 type_ ==
"BLOCK_RELAXATION" ||
169 type_ ==
"BLOCKRELAXATION" ||
171 type_ ==
"BANDED_RELAXATION" ||
172 type_ ==
"BANDED RELAXATION" ||
173 type_ ==
"BANDEDRELAXATION" ||
175 type_ ==
"TRIDI_RELAXATION" ||
176 type_ ==
"TRIDI RELAXATION" ||
177 type_ ==
"TRIDIRELAXATION" ||
178 type_ ==
"TRIDIAGONAL_RELAXATION" ||
179 type_ ==
"TRIDIAGONAL RELAXATION" ||
180 type_ ==
"TRIDIAGONALRELAXATION")
183 ParameterList precList = this->GetParameterList();
184 if(precList.isParameter(
"partitioner: type") &&
185 precList.get<std::string>(
"partitioner: type") ==
"line") {
186 this->Input(currentLevel,
"Coordinates");
189 else if (type_ ==
"TOPOLOGICAL")
192 this->Input(currentLevel,
"pcoarsen: element to node map");
194 else if (type_ ==
"AGGREGATE")
197 this->Input(currentLevel,
"Aggregates");
199 else if (type_ ==
"HIPTMAIR") {
201 this->Input(currentLevel,
"NodeMatrix");
202 this->Input(currentLevel,
"D0");
207 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
210 A_ = Factory::Get< RCP<Matrix> >(currentLevel,
"A");
211 ParameterList& paramList =
const_cast<ParameterList&
>(this->GetParameterList());
214 if(paramList.isParameter(
"smoother: use blockcrsmatrix storage") && paramList.get<
bool>(
"smoother: use blockcrsmatrix storage") && A_->GetFixedBlockSize()) {
216 int blocksize = A_->GetFixedBlockSize();
217 using TpetraBlockCrsMatrix = Xpetra::TpetraBlockCrsMatrix<SC,LO,GO,NO>;
218 RCP<CrsMatrixWrap> AcrsWrap = rcp_dynamic_cast<CrsMatrixWrap>(A_);
219 if(AcrsWrap.is_null())
220 throw std::runtime_error(
"Ifpack2Smoother: Cannot convert matrix A to CrsMatrixWrap object.");
221 RCP<CrsMatrix> Acrs = AcrsWrap->getCrsMatrix();
223 throw std::runtime_error(
"Ifpack2Smoother: Cannot extract CrsMatrix from matrix A.");
224 RCP<TpetraCrsMatrix> At = rcp_dynamic_cast<TpetraCrsMatrix>(Acrs);
226 throw std::runtime_error(
"Ifpack2Smoother: Cannot extract TpetraCrsMatrix from matrix A.");
228 RCP<Tpetra::BlockCrsMatrix<Scalar, LO, GO, Node> > blockCrs = Tpetra::convertToBlockCrsMatrix(*At->getTpetra_CrsMatrix(),blocksize);
229 RCP<CrsMatrix> blockCrs_as_crs = rcp(
new TpetraBlockCrsMatrix(blockCrs));
230 RCP<CrsMatrixWrap> blockWrap = rcp(
new CrsMatrixWrap(blockCrs_as_crs));
232 this->GetOStream(
Statistics0) <<
"Ifpack2Smoother: Using BlockCrsMatrix storage with blocksize "<<blocksize<<std::endl;
234 paramList.remove(
"smoother: use blockcrsmatrix storage");
237 if (type_ ==
"SCHWARZ")
238 SetupSchwarz(currentLevel);
240 else if (type_ ==
"LINESMOOTHING_TRIDI_RELAXATION" ||
241 type_ ==
"LINESMOOTHING_TRIDI RELAXATION" ||
242 type_ ==
"LINESMOOTHING_TRIDIRELAXATION" ||
243 type_ ==
"LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
244 type_ ==
"LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
245 type_ ==
"LINESMOOTHING_TRIDIAGONALRELAXATION" ||
246 type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
247 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
248 type_ ==
"LINESMOOTHING_BANDEDRELAXATION" ||
249 type_ ==
"LINESMOOTHING_BLOCK_RELAXATION" ||
250 type_ ==
"LINESMOOTHING_BLOCK RELAXATION" ||
251 type_ ==
"LINESMOOTHING_BLOCKRELAXATION")
252 SetupLineSmoothing(currentLevel);
254 else if (type_ ==
"BLOCK_RELAXATION" ||
255 type_ ==
"BLOCK RELAXATION" ||
256 type_ ==
"BLOCKRELAXATION" ||
258 type_ ==
"BANDED_RELAXATION" ||
259 type_ ==
"BANDED RELAXATION" ||
260 type_ ==
"BANDEDRELAXATION" ||
262 type_ ==
"TRIDI_RELAXATION" ||
263 type_ ==
"TRIDI RELAXATION" ||
264 type_ ==
"TRIDIRELAXATION" ||
265 type_ ==
"TRIDIAGONAL_RELAXATION" ||
266 type_ ==
"TRIDIAGONAL RELAXATION" ||
267 type_ ==
"TRIDIAGONALRELAXATION")
268 SetupBlockRelaxation(currentLevel);
270 else if (type_ ==
"CHEBYSHEV")
271 SetupChebyshev(currentLevel);
273 else if (type_ ==
"TOPOLOGICAL")
275 #ifdef HAVE_MUELU_INTREPID2 276 SetupTopological(currentLevel);
278 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"'TOPOLOGICAL' smoother choice requires Intrepid2");
281 else if (type_ ==
"AGGREGATE")
282 SetupAggregate(currentLevel);
284 else if (type_ ==
"HIPTMAIR")
285 SetupHiptmair(currentLevel);
289 SetupGeneric(currentLevel);
294 this->GetOStream(
Statistics1) << description() << std::endl;
297 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
299 typedef Tpetra::RowMatrix<SC,LO,GO,NO> tRowMatrix;
301 bool reusePreconditioner =
false;
302 if (this->IsSetup() ==
true) {
304 this->GetOStream(
Runtime1) <<
"MueLu::Ifpack2Smoother::SetupSchwarz(): Setup() has already been called, assuming reuse" << std::endl;
306 bool isTRowMatrix =
true;
307 RCP<const tRowMatrix> tA;
311 isTRowMatrix =
false;
315 if (!prec.is_null() && isTRowMatrix) {
317 reusePreconditioner =
true;
319 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupSchwarz(): reuse of this type is not available " 320 "(either failed cast to CanChangeMatrix, or to Tpetra Row Matrix), reverting to full construction" << std::endl;
324 if (!reusePreconditioner) {
325 ParameterList& paramList =
const_cast<ParameterList&
>(this->GetParameterList());
327 bool isBlockedMatrix =
false;
328 RCP<Matrix> merged2Mat;
330 std::string sublistName =
"subdomain solver parameters";
331 if (paramList.isSublist(sublistName)) {
340 ParameterList& subList = paramList.sublist(sublistName);
342 std::string partName =
"partitioner: type";
343 if (subList.isParameter(partName) && subList.get<std::string>(partName) ==
"user") {
344 isBlockedMatrix =
true;
346 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
348 "Matrix A must be of type BlockedCrsMatrix.");
350 size_t numVels = bA->getMatrix(0,0)->getNodeNumRows();
351 size_t numPres = bA->getMatrix(1,0)->getNodeNumRows();
352 size_t numRows = A_->getNodeNumRows();
354 ArrayRCP<LocalOrdinal> blockSeeds(numRows, Teuchos::OrdinalTraits<LocalOrdinal>::invalid());
356 size_t numBlocks = 0;
357 for (
size_t rowOfB = numVels; rowOfB < numVels+numPres; ++rowOfB)
358 blockSeeds[rowOfB] = numBlocks++;
360 RCP<BlockedCrsMatrix> bA2 = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
362 "Matrix A must be of type BlockedCrsMatrix.");
364 merged2Mat = bA2->Merge();
368 bool haveBoundary =
false;
369 for (LO i = 0; i < boundaryNodes.size(); i++)
370 if (boundaryNodes[i]) {
374 blockSeeds[i] = numBlocks;
380 subList.set(
"partitioner: map", blockSeeds);
381 subList.set(
"partitioner: local parts", as<int>(numBlocks));
384 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
386 isBlockedMatrix =
true;
387 merged2Mat = bA->Merge();
392 RCP<const tRowMatrix> tA;
407 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
410 ParameterList& paramList =
const_cast<ParameterList&
>(this->GetParameterList());
412 if (this->IsSetup() ==
true) {
413 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupAggregate(): Setup() has already been called" << std::endl;
414 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupAggregate(): reuse of this type is not available, reverting to full construction" << std::endl;
417 this->GetOStream(
Statistics0) <<
"Ifpack2Smoother: Using Aggregate Smoothing"<<std::endl;
419 RCP<Aggregates> aggregates = Factory::Get<RCP<Aggregates> >(currentLevel,
"Aggregates");
421 RCP<const LOMultiVector> vertex2AggId = aggregates->GetVertex2AggId();
422 ArrayRCP<LO> aggregate_ids = rcp_const_cast<LOMultiVector>(vertex2AggId)->getDataNonConst(0);
423 ArrayRCP<LO> dof_ids;
426 if(A_->GetFixedBlockSize() > 1) {
427 LO blocksize = (LO) A_->GetFixedBlockSize();
428 dof_ids.resize(aggregate_ids.size() * blocksize);
429 for(LO i=0; i<(LO)aggregate_ids.size(); i++) {
430 for(LO j=0; j<(LO)blocksize; j++)
431 dof_ids[i*blocksize+j] = aggregate_ids[i];
435 dof_ids = aggregate_ids;
439 paramList.set(
"partitioner: map", dof_ids);
440 paramList.set(
"partitioner: type",
"user");
441 paramList.set(
"partitioner: overlap", 0);
442 paramList.set(
"partitioner: local parts", (
int)aggregates->GetNumAggregates());
446 type_ =
"BLOCKRELAXATION";
453 #ifdef HAVE_MUELU_INTREPID2 454 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
464 if (this->IsSetup() ==
true) {
465 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupTopological(): Setup() has already been called" << std::endl;
466 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupTopological(): reuse of this type is not available, reverting to full construction" << std::endl;
469 ParameterList& paramList =
const_cast<ParameterList&
>(this->GetParameterList());
471 typedef typename Node::device_type::execution_space ES;
473 typedef Kokkos::DynRankView<LocalOrdinal,typename Node::device_type> FCO;
475 LocalOrdinal lo_invalid = Teuchos::OrdinalTraits<LO>::invalid();
479 const Teuchos::RCP<FCO> elemToNode = Factory::Get<Teuchos::RCP<FCO> >(currentLevel,
"pcoarsen: element to node map");
481 string basisString = paramList.get<
string>(
"pcoarsen: hi basis");
487 auto basis = MueLuIntrepid::BasisFactory<double,ES>(basisString, degree);
489 string topologyTypeString = paramList.get<
string>(
"smoother: neighborhood type");
491 if (topologyTypeString ==
"node")
493 else if (topologyTypeString ==
"edge")
495 else if (topologyTypeString ==
"face")
497 else if (topologyTypeString ==
"cell")
498 dimension = basis->getBaseCellTopology().getDimension();
500 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::invalid_argument,
"Unrecognized smoother neighborhood type. Supported types are node, edge, face.");
501 vector<vector<LocalOrdinal>> seeds;
506 int myNodeCount = A_->getRowMap()->getNodeNumElements();
507 ArrayRCP<LocalOrdinal> nodeSeeds(myNodeCount,lo_invalid);
508 int localPartitionNumber = 0;
511 nodeSeeds[seed] = localPartitionNumber++;
514 paramList.remove(
"smoother: neighborhood type");
515 paramList.remove(
"pcoarsen: hi basis");
517 paramList.set(
"partitioner: map", nodeSeeds);
518 paramList.set(
"partitioner: type",
"user");
519 paramList.set(
"partitioner: overlap", 1);
520 paramList.set(
"partitioner: local parts",
int(seeds[dimension].size()));
524 type_ =
"BLOCKRELAXATION";
532 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
534 if (this->IsSetup() ==
true) {
535 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupLineSmoothing(): Setup() has already been called" << std::endl;
536 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupLineSmoothing(): reuse of this type is not available, reverting to full construction" << std::endl;
539 ParameterList& myparamList =
const_cast<ParameterList&
>(this->GetParameterList());
541 LO CoarseNumZLayers = Factory::Get<LO>(currentLevel,
"CoarseNumZLayers");
542 if (CoarseNumZLayers > 0) {
543 Teuchos::ArrayRCP<LO> TVertLineIdSmoo = Factory::Get< Teuchos::ArrayRCP<LO> >(currentLevel,
"LineDetection_VertLineIds");
547 for(
size_t k = 0; k < Teuchos::as<size_t>(TVertLineIdSmoo.size()); k++) {
548 if(maxPart < TVertLineIdSmoo[k]) maxPart = TVertLineIdSmoo[k];
550 size_t numLocalRows = A_->getNodeNumRows();
553 "MueLu::Ifpack2Smoother::Setup(): the number of local nodes is incompatible with the TVertLineIdsSmoo.");
558 int actualDofsPerNode = numLocalRows / TVertLineIdSmoo.size();
559 LO matrixBlockSize = A_->GetFixedBlockSize();
560 if(matrixBlockSize > 1 && actualDofsPerNode > 1)
563 "MueLu::Ifpack2Smoother::Setup(): A is a block matrix but its block size and DOFs/node from partitioner disagree");
565 else if(matrixBlockSize > 1)
567 actualDofsPerNode = A_->GetFixedBlockSize();
569 myparamList.set(
"partitioner: PDE equations", actualDofsPerNode);
571 if (numLocalRows == Teuchos::as<size_t>(TVertLineIdSmoo.size())) {
572 myparamList.set(
"partitioner: type",
"user");
573 myparamList.set(
"partitioner: map",TVertLineIdSmoo);
574 myparamList.set(
"partitioner: local parts",maxPart+1);
577 size_t numDofsPerNode = numLocalRows / TVertLineIdSmoo.size();
580 Teuchos::ArrayRCP<LO> partitionerMap(numLocalRows, Teuchos::OrdinalTraits<LocalOrdinal>::invalid());
581 for (
size_t blockRow = 0; blockRow < Teuchos::as<size_t>(TVertLineIdSmoo.size()); ++blockRow)
582 for (
size_t dof = 0; dof < numDofsPerNode; dof++)
583 partitionerMap[blockRow * numDofsPerNode + dof] = TVertLineIdSmoo[blockRow];
584 myparamList.set(
"partitioner: type",
"user");
585 myparamList.set(
"partitioner: map",partitionerMap);
586 myparamList.set(
"partitioner: local parts",maxPart + 1);
589 if (type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
590 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
591 type_ ==
"LINESMOOTHING_BANDEDRELAXATION")
592 type_ =
"BANDEDRELAXATION";
593 else if (type_ ==
"LINESMOOTHING_TRIDI_RELAXATION" ||
594 type_ ==
"LINESMOOTHING_TRIDI RELAXATION" ||
595 type_ ==
"LINESMOOTHING_TRIDIRELAXATION" ||
596 type_ ==
"LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
597 type_ ==
"LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
598 type_ ==
"LINESMOOTHING_TRIDIAGONALRELAXATION")
599 type_ =
"TRIDIAGONALRELAXATION";
601 type_ =
"BLOCKRELAXATION";
604 this->GetOStream(
Runtime0) <<
"Line detection failed: fall back to point-wise relaxation" << std::endl;
605 myparamList.remove(
"partitioner: type",
false);
606 myparamList.remove(
"partitioner: map",
false);
607 myparamList.remove(
"partitioner: local parts",
false);
608 type_ =
"RELAXATION";
619 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
621 typedef Tpetra::RowMatrix<SC,LO,GO,NO> tRowMatrix;
623 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
629 bool reusePreconditioner =
false;
630 if (this->IsSetup() ==
true) {
632 this->GetOStream(
Runtime1) <<
"MueLu::Ifpack2Smoother::SetupBlockRelaxation(): Setup() has already been called, assuming reuse" << std::endl;
635 if (!prec.is_null()) {
637 reusePreconditioner =
true;
639 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupBlockRelaxation(): reuse of this type is not available (failed cast to CanChangeMatrix), " 640 "reverting to full construction" << std::endl;
644 if (!reusePreconditioner) {
645 ParameterList& myparamList =
const_cast<ParameterList&
>(this->GetParameterList());
647 if(myparamList.isParameter(
"partitioner: type") &&
648 myparamList.get<std::string>(
"partitioner: type") ==
"line") {
649 Teuchos::RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> > xCoordinates =
650 Factory::Get<Teuchos::RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> > >(currentLevel,
"Coordinates");
651 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));
653 size_t numDofsPerNode = A_->getNodeNumRows() / xCoordinates->getMap()->getNodeNumElements();
654 myparamList.set(
"partitioner: coordinates", coordinates);
655 myparamList.set(
"partitioner: PDE equations", (
int) numDofsPerNode);
666 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
668 typedef Tpetra::RowMatrix<SC,LO,GO,NO> tRowMatrix;
669 using STS = Teuchos::ScalarTraits<SC>;
670 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
676 bool reusePreconditioner =
false;
678 if (this->IsSetup() ==
true) {
680 this->GetOStream(
Runtime1) <<
"MueLu::Ifpack2Smoother::SetupChebyshev(): Setup() has already been called, assuming reuse" << std::endl;
683 if (!prec.is_null()) {
685 reusePreconditioner =
true;
687 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupChebyshev(): reuse of this type is not available (failed cast to CanChangeMatrix), " 688 "reverting to full construction" << std::endl;
693 ParameterList& paramList =
const_cast<ParameterList&
>(this->GetParameterList());
694 SC negone = -STS::one();
695 SC lambdaMax = SetupChebyshevEigenvalues(currentLevel,
"A",
"",paramList);
698 if (!reusePreconditioner) {
713 if (lambdaMax == negone) {
714 typedef Tpetra::RowMatrix<SC, LO, GO, NO> MatrixType;
717 if (chebyPrec != Teuchos::null) {
719 A_->SetMaxEigenvalueEstimate(lambdaMax);
720 this->GetOStream(
Statistics1) <<
"chebyshev: max eigenvalue (calculated by Ifpack2)" <<
" = " << lambdaMax << std::endl;
722 TEUCHOS_TEST_FOR_EXCEPTION(lambdaMax == negone,
Exceptions::RuntimeError,
"MueLu::Ifpack2Smoother::Setup(): no maximum eigenvalue estimate");
727 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
729 typedef Tpetra::RowMatrix<SC,LO,GO,NO> tRowMatrix;
730 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
736 bool reusePreconditioner =
false;
737 if (this->IsSetup() ==
true) {
739 this->GetOStream(
Runtime1) <<
"MueLu::Ifpack2Smoother::SetupHiptmair(): Setup() has already been called, assuming reuse" << std::endl;
742 if (!prec.is_null()) {
744 reusePreconditioner =
true;
746 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupHiptmair(): reuse of this type is not available (failed cast to CanChangeMatrix), " 747 "reverting to full construction" << std::endl;
752 ParameterList& paramList =
const_cast<ParameterList&
>(this->GetParameterList());
753 std::string smoother1 = paramList.get(
"hiptmair: smoother type 1",
"CHEBYSHEV");
754 std::string smoother2 = paramList.get(
"hiptmair: smoother type 2",
"CHEBYSHEV");
757 if(smoother1 ==
"CHEBYSHEV") {
758 ParameterList & list1 = paramList.sublist(
"hiptmair: smoother list 1");
760 SetupChebyshevEigenvalues(currentLevel,
"A",
"EdgeMatrix ",list1);
762 if(smoother2 ==
"CHEBYSHEV") {
763 ParameterList & list2 = paramList.sublist(
"hiptmair: smoother list 2");
765 SetupChebyshevEigenvalues(currentLevel,
"A",
"EdgeMatrix ",list2);
772 RCP<Matrix> NodeMatrix = currentLevel.
Get<RCP<Matrix> >(
"NodeMatrix");
773 RCP<Matrix> D0 = currentLevel.
Get<RCP<Matrix> >(
"D0");
779 Teuchos::ParameterList newlist;
780 newlist.set(
"P",tD0);
781 newlist.set(
"PtAP",tNodeMatrix);
782 if (!reusePreconditioner) {
784 SetPrecParameters(newlist);
793 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
796 typedef Teuchos::ScalarTraits<SC> STS;
797 SC negone = -STS::one();
798 RCP<const Matrix> currentA = currentLevel.
Get<RCP<Matrix> >(matrixName);
799 SC lambdaMax = negone;
801 std::string maxEigString =
"chebyshev: max eigenvalue";
802 std::string eigRatioString =
"chebyshev: ratio eigenvalue";
805 if (paramList.isParameter(maxEigString)) {
806 if (paramList.isType<
double>(maxEigString))
807 lambdaMax = paramList.get<
double>(maxEigString);
809 lambdaMax = paramList.get<SC>(maxEigString);
810 this->GetOStream(
Statistics1) << label << maxEigString <<
" (cached with smoother parameter list) = " << lambdaMax << std::endl;
813 lambdaMax = currentA->GetMaxEigenvalueEstimate();
814 if (lambdaMax != negone) {
815 this->GetOStream(
Statistics1) << label << maxEigString <<
" (cached with matrix) = " << lambdaMax << std::endl;
816 paramList.set(maxEigString, lambdaMax);
821 const SC defaultEigRatio = 20;
823 SC ratio = defaultEigRatio;
824 if (paramList.isParameter(eigRatioString)) {
825 if (paramList.isType<
double>(eigRatioString))
826 ratio = paramList.get<
double>(eigRatioString);
828 ratio = paramList.get<SC>(eigRatioString);
835 RCP<const Matrix> fineA = currentLevel.
GetPreviousLevel()->Get<RCP<Matrix> >(matrixName);
836 size_t nRowsFine = fineA->getGlobalNumRows();
837 size_t nRowsCoarse = currentA->getGlobalNumRows();
839 SC levelRatio = as<SC>(as<float>(nRowsFine)/nRowsCoarse);
840 if (STS::magnitude(levelRatio) > STS::magnitude(ratio))
844 this->GetOStream(
Statistics1) << label << eigRatioString <<
" (computed) = " << ratio << std::endl;
845 paramList.set(eigRatioString, ratio);
847 if (paramList.isParameter(
"chebyshev: use rowsumabs diagonal scaling")) {
848 this->GetOStream(
Runtime1) <<
"chebyshev: using rowsumabs diagonal scaling" << std::endl;
849 bool doScale =
false;
850 doScale = paramList.get<
bool>(
"chebyshev: use rowsumabs diagonal scaling");
851 paramList.remove(
"chebyshev: use rowsumabs diagonal scaling");
852 double chebyReplaceTol = Teuchos::ScalarTraits<Scalar>::eps()*100;
853 if (paramList.isParameter(
"chebyshev: rowsumabs diagonal replacement tolerance")) {
854 paramList.get<
double>(
"chebyshev: rowsumabs diagonal replacement tolerance",chebyReplaceTol);
855 paramList.remove(
"chebyshev: rowsumabs diagonal replacement tolerance");
857 double chebyReplaceVal = Teuchos::ScalarTraits<double>::zero();
858 if (paramList.isParameter(
"chebyshev: rowsumabs diagonal replacement value")) {
859 paramList.get<
double>(
"chebyshev: rowsumabs diagonal replacement value",chebyReplaceVal);
860 paramList.remove(
"chebyshev: rowsumabs diagonal replacement value");
864 const Xpetra::TpetraVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmpVec =
dynamic_cast<const Xpetra::TpetraVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>&
>(*lumpedDiagonal);
865 paramList.set(
"chebyshev: operator inv diagonal",tmpVec.getTpetra_Vector());
875 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
877 typedef Tpetra::RowMatrix<SC,LO,GO,NO> tRowMatrix;
878 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
884 bool reusePreconditioner =
false;
885 if (this->IsSetup() ==
true) {
887 this->GetOStream(
Runtime1) <<
"MueLu::Ifpack2Smoother::SetupGeneric(): Setup() has already been called, assuming reuse" << std::endl;
890 if (!prec.is_null()) {
892 reusePreconditioner =
true;
894 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupGeneric(): reuse of this type is not available (failed cast to CanChangeMatrix), " 895 "reverting to full construction" << std::endl;
899 if (!reusePreconditioner) {
908 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
918 Teuchos::ParameterList paramList;
919 bool supportInitialGuess =
false;
920 const Teuchos::ParameterList params = this->GetParameterList();
922 if (prec_->supportsZeroStartingSolution()) {
923 prec_->setZeroStartingSolution(InitialGuessIsZero);
924 supportInitialGuess =
true;
925 }
else if (type_ ==
"SCHWARZ") {
926 paramList.set(
"schwarz: zero starting solution", InitialGuessIsZero);
931 prec_->setParameters(paramList);
932 supportInitialGuess =
true;
942 if (InitialGuessIsZero || supportInitialGuess) {
945 prec_->apply(tpB, tpX);
947 typedef Teuchos::ScalarTraits<Scalar> TST;
949 RCP<MultiVector> Residual;
951 std::string prefix = this->ShortClassName() +
": Apply: ";
952 RCP<TimeMonitor> tM = rcp(
new TimeMonitor(*
this, prefix +
"residual calculation",
Timings0));
956 RCP<MultiVector> Correction = MultiVectorFactory::Build(A_->getDomainMap(), X.getNumVectors());
961 prec_->apply(tpB, tpX);
963 X.update(TST::one(), *Correction, TST::one());
967 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
970 smoother->SetParameterList(this->GetParameterList());
974 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
976 std::ostringstream out;
978 out << prec_->description();
981 out <<
"{type = " << type_ <<
"}";
986 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
991 out0 <<
"Prec. type: " << type_ << std::endl;
994 out0 <<
"Parameter list: " << std::endl;
995 Teuchos::OSTab tab2(out);
996 out << this->GetParameterList();
997 out0 <<
"Overlap: " << overlap_ << std::endl;
1001 if (prec_ != Teuchos::null) {
1002 Teuchos::OSTab tab2(out);
1003 out << *prec_ << std::endl;
1006 if (verbLevel &
Debug) {
1009 <<
"RCP<prec_>: " << prec_ << std::endl;
1013 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1015 typedef Tpetra::RowMatrix<SC,LO,GO,NO> MatrixType;
1033 return Teuchos::OrdinalTraits<size_t>::invalid();
1039 #endif // HAVE_MUELU_TPETRA && HAVE_MUELU_IFPACK2 1040 #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.
Scalar SetupChebyshevEigenvalues(Level ¤tLevel, const std::string &matrixName, const std::string &label, Teuchos::ParameterList ¶mList) const
High level timing information (use Teuchos::TimeMonitor::summarize() to print)
void SetPrecParameters(const Teuchos::ParameterList &list=Teuchos::ParameterList()) const
MueLu::DefaultLocalOrdinal LocalOrdinal
size_t getNodeSmootherComplexity() const
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level->Get< RCP<Matrix> >("A", factory) if factory == NULL => use default factory.
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)
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.
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)
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)
Set parameters from a parameter list and return with default values.
Namespace for MueLu classes and methods.
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)
Integrates Teuchos::TimeMonitor with MueLu verbosity system.
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
int GetLevelID() const
Return level number.
std::string type_
ifpack2-specific key phrase that denote smoother type
void SetupBlockRelaxation(Level ¤tLevel)
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
void SetupSchwarz(Level ¤tLevel)
void SetupLineSmoothing(Level ¤tLevel)
Print statistics that do not involve significant additional computation.
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 Teuchos::ParameterList ¶mList)
Set parameters from a parameter list and return with default values.
MueLu::DefaultScalar Scalar
void declareConstructionOutcome(bool fail, std::string msg)
Class that holds all level-specific information.
void Setup(Level ¤tLevel)
Set up the smoother.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
size_t getNodeSmootherComplexity() const
void SetupHiptmair(Level ¤tLevel)
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the preconditioner.
void SetupTopological(Level ¤tLevel)
size_t getNodeSmootherComplexity() const
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.
void SetupAggregate(Level ¤tLevel)
size_t getNodeSmootherComplexity() const
Print class parameters (more parameters, more verbose)
size_t getNodeSmootherComplexity() const
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.
virtual void setMatrix(const Teuchos::RCP< const RowMatrixType > &A)=0