46 #ifndef MUELU_IFPACK2SMOOTHER_DEF_HPP
47 #define MUELU_IFPACK2SMOOTHER_DEF_HPP
51 #if defined(HAVE_MUELU_IFPACK2)
55 #include <Tpetra_RowMatrix.hpp>
57 #include <Ifpack2_Chebyshev.hpp>
58 #include <Ifpack2_Hiptmair.hpp>
59 #include <Ifpack2_RILUK.hpp>
60 #include <Ifpack2_Relaxation.hpp>
61 #include <Ifpack2_ILUT.hpp>
62 #include <Ifpack2_BlockRelaxation.hpp>
63 #include <Ifpack2_Factory.hpp>
64 #include <Ifpack2_Parameters.hpp>
68 #include <Xpetra_CrsMatrixWrap.hpp>
69 #include <Xpetra_TpetraBlockCrsMatrix.hpp>
72 #include <Xpetra_MultiVectorFactory.hpp>
73 #include <Xpetra_TpetraMultiVector.hpp>
75 #include <Tpetra_BlockCrsMatrix_Helpers.hpp>
79 #include "MueLu_Utilities.hpp"
81 #include "MueLu_Aggregates.hpp"
83 #ifdef HAVE_MUELU_INTREPID2
86 #include "Intrepid2_Basis.hpp"
87 #include "Kokkos_DynRankView.hpp"
92 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
97 bool isSupported = Ifpack2::Factory::isSupported<tRowMatrix>(
type_) || (
type_ ==
"LINESMOOTHING_TRIDI_RELAXATION" ||
98 type_ ==
"LINESMOOTHING_TRIDI RELAXATION" ||
99 type_ ==
"LINESMOOTHING_TRIDIRELAXATION" ||
100 type_ ==
"LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
101 type_ ==
"LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
102 type_ ==
"LINESMOOTHING_TRIDIAGONALRELAXATION" ||
103 type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
104 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
105 type_ ==
"LINESMOOTHING_BANDEDRELAXATION" ||
106 type_ ==
"LINESMOOTHING_BLOCK_RELAXATION" ||
107 type_ ==
"LINESMOOTHING_BLOCK RELAXATION" ||
108 type_ ==
"LINESMOOTHING_BLOCKRELAXATION" ||
109 type_ ==
"TOPOLOGICAL" ||
110 type_ ==
"AGGREGATE");
116 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
127 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
129 std::string prefix = this->ShortClassName() +
": SetPrecParameters";
135 RCP<ParameterList> precList = this->RemoveFactoriesFromList(this->GetParameterList());
138 precList->
set(
"timer for apply", this->IsPrint(
Timings));
140 if (!precList.
is_null() && precList->
isParameter(
"partitioner: type") && precList->
get<std::string>(
"partitioner: type") ==
"linear" &&
141 !precList->
isParameter(
"partitioner: local parts")) {
142 LO matrixBlockSize = 1;
143 int lclSize = A_->getRangeMap()->getLocalNumElements();
146 lclSize = matA->getLocalNumRows();
147 matrixBlockSize = matA->GetFixedBlockSize();
149 precList->
set(
"partitioner: local parts", lclSize / matrixBlockSize);
157 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
159 this->Input(currentLevel,
"A");
161 if (type_ ==
"LINESMOOTHING_TRIDI_RELAXATION" ||
162 type_ ==
"LINESMOOTHING_TRIDI RELAXATION" ||
163 type_ ==
"LINESMOOTHING_TRIDIRELAXATION" ||
164 type_ ==
"LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
165 type_ ==
"LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
166 type_ ==
"LINESMOOTHING_TRIDIAGONALRELAXATION" ||
167 type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
168 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
169 type_ ==
"LINESMOOTHING_BANDEDRELAXATION" ||
170 type_ ==
"LINESMOOTHING_BLOCK_RELAXATION" ||
171 type_ ==
"LINESMOOTHING_BLOCK RELAXATION" ||
172 type_ ==
"LINESMOOTHING_BLOCKRELAXATION") {
173 this->Input(currentLevel,
"CoarseNumZLayers");
174 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") {
192 precList.
get<std::string>(
"partitioner: type") ==
"line") {
193 this->Input(currentLevel,
"Coordinates");
195 }
else if (type_ ==
"TOPOLOGICAL") {
197 this->Input(currentLevel,
"pcoarsen: element to node map");
198 }
else if (type_ ==
"AGGREGATE") {
200 this->Input(currentLevel,
"Aggregates");
201 }
else if (type_ ==
"HIPTMAIR") {
203 this->Input(currentLevel,
"NodeMatrix");
204 this->Input(currentLevel,
"D0");
208 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
211 A_ = Factory::Get<RCP<Operator>>(currentLevel,
"A");
216 if (paramList.isParameter(
"smoother: use blockcrsmatrix storage") && paramList.get<
bool>(
"smoother: use blockcrsmatrix storage")) {
220 blocksize = matA->GetFixedBlockSize();
226 throw std::runtime_error(
"Ifpack2Smoother: Cannot convert matrix A to CrsMatrixWrap object.");
229 throw std::runtime_error(
"Ifpack2Smoother: Cannot extract CrsMatrix from matrix A.");
233 throw std::runtime_error(
"Ifpack2Smoother: Cannot extract CrsMatrix or BlockCrsMatrix from matrix A.");
234 this->GetOStream(
Statistics0) <<
"Ifpack2Smoother: Using (native) BlockCrsMatrix storage with blocksize " << blocksize << std::endl;
235 paramList.remove(
"smoother: use blockcrsmatrix storage");
241 this->GetOStream(
Statistics0) <<
"Ifpack2Smoother: Using BlockCrsMatrix storage with blocksize " << blocksize << std::endl;
243 paramList.remove(
"smoother: use blockcrsmatrix storage");
248 if (type_ ==
"SCHWARZ")
249 SetupSchwarz(currentLevel);
251 else if (type_ ==
"LINESMOOTHING_TRIDI_RELAXATION" ||
252 type_ ==
"LINESMOOTHING_TRIDI RELAXATION" ||
253 type_ ==
"LINESMOOTHING_TRIDIRELAXATION" ||
254 type_ ==
"LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
255 type_ ==
"LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
256 type_ ==
"LINESMOOTHING_TRIDIAGONALRELAXATION" ||
257 type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
258 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
259 type_ ==
"LINESMOOTHING_BANDEDRELAXATION" ||
260 type_ ==
"LINESMOOTHING_BLOCK_RELAXATION" ||
261 type_ ==
"LINESMOOTHING_BLOCK RELAXATION" ||
262 type_ ==
"LINESMOOTHING_BLOCKRELAXATION")
263 SetupLineSmoothing(currentLevel);
265 else if (type_ ==
"BLOCK_RELAXATION" ||
266 type_ ==
"BLOCK RELAXATION" ||
267 type_ ==
"BLOCKRELAXATION" ||
269 type_ ==
"BANDED_RELAXATION" ||
270 type_ ==
"BANDED RELAXATION" ||
271 type_ ==
"BANDEDRELAXATION" ||
273 type_ ==
"TRIDI_RELAXATION" ||
274 type_ ==
"TRIDI RELAXATION" ||
275 type_ ==
"TRIDIRELAXATION" ||
276 type_ ==
"TRIDIAGONAL_RELAXATION" ||
277 type_ ==
"TRIDIAGONAL RELAXATION" ||
278 type_ ==
"TRIDIAGONALRELAXATION")
279 SetupBlockRelaxation(currentLevel);
281 else if (type_ ==
"CHEBYSHEV")
282 SetupChebyshev(currentLevel);
284 else if (type_ ==
"TOPOLOGICAL") {
285 #ifdef HAVE_MUELU_INTREPID2
286 SetupTopological(currentLevel);
290 }
else if (type_ ==
"AGGREGATE")
291 SetupAggregate(currentLevel);
293 else if (type_ ==
"HIPTMAIR")
294 SetupHiptmair(currentLevel);
297 SetupGeneric(currentLevel);
302 this->GetOStream(
Statistics1) << description() << std::endl;
305 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
309 bool reusePreconditioner =
false;
310 if (this->IsSetup() ==
true) {
312 this->GetOStream(
Runtime1) <<
"MueLu::Ifpack2Smoother::SetupSchwarz(): Setup() has already been called, assuming reuse" << std::endl;
314 bool isTRowMatrix =
true;
319 isTRowMatrix =
false;
323 if (!prec.
is_null() && isTRowMatrix) {
325 reusePreconditioner =
true;
327 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupSchwarz(): reuse of this type is not available "
328 "(either failed cast to CanChangeMatrix, or to Tpetra Row Matrix), reverting to full construction"
333 if (!reusePreconditioner) {
336 bool isBlockedMatrix =
false;
339 std::string sublistName =
"subdomain solver parameters";
351 std::string partName =
"partitioner: type";
357 if (subList.
isParameter(partName) && subList.
get<std::string>(partName) ==
"vanka user") {
358 isBlockedMatrix =
true;
362 "Matrix A must be of type BlockedCrsMatrix.");
364 size_t numVels = bA->getMatrix(0, 0)->getLocalNumRows();
365 size_t numPres = bA->getMatrix(1, 0)->getLocalNumRows();
366 size_t numRows = rcp_dynamic_cast<Matrix>(A_,
true)->getLocalNumRows();
370 size_t numBlocks = 0;
371 for (
size_t rowOfB = numVels; rowOfB < numVels + numPres; ++rowOfB)
372 blockSeeds[rowOfB] = numBlocks++;
376 "Matrix A must be of type BlockedCrsMatrix.");
378 merged2Mat = bA2->Merge();
382 bool haveBoundary =
false;
383 for (
LO i = 0; i < boundaryNodes.
size(); i++)
384 if (boundaryNodes[i]) {
388 blockSeeds[i] = numBlocks;
394 subList.
set(
"partitioner: type",
"user");
395 subList.
set(
"partitioner: map", blockSeeds);
396 subList.
set(
"partitioner: local parts", as<int>(numBlocks));
401 isBlockedMatrix =
true;
402 merged2Mat = bA->Merge();
408 if (isBlockedMatrix ==
true)
422 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
426 if (this->IsSetup() ==
true) {
427 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupAggregate(): Setup() has already been called" << std::endl;
428 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupAggregate(): reuse of this type is not available, reverting to full construction" << std::endl;
431 this->GetOStream(
Statistics0) <<
"Ifpack2Smoother: Using Aggregate Smoothing" << std::endl;
433 RCP<Aggregates> aggregates = Factory::Get<RCP<Aggregates>>(currentLevel,
"Aggregates");
436 ArrayRCP<LO> aggregate_ids = rcp_const_cast<LOMultiVector>(vertex2AggId)->getDataNonConst(0);
443 blocksize = matA->GetFixedBlockSize();
445 dof_ids.
resize(aggregate_ids.
size() * blocksize);
446 for (
LO i = 0; i < (
LO)aggregate_ids.
size(); i++) {
447 for (
LO j = 0; j < (
LO)blocksize; j++)
448 dof_ids[i * blocksize + j] = aggregate_ids[i];
451 dof_ids = aggregate_ids;
454 paramList.
set(
"partitioner: map", dof_ids);
455 paramList.
set(
"partitioner: type",
"user");
456 paramList.
set(
"partitioner: overlap", 0);
461 type_ =
"BLOCKRELAXATION";
468 #ifdef HAVE_MUELU_INTREPID2
469 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
479 if (this->IsSetup() ==
true) {
480 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupTopological(): Setup() has already been called" << std::endl;
481 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupTopological(): reuse of this type is not available, reverting to full construction" << std::endl;
484 ParameterList& paramList =
const_cast<ParameterList&
>(this->GetParameterList());
486 typedef typename Node::device_type::execution_space ES;
488 typedef Kokkos::DynRankView<LocalOrdinal, typename Node::device_type> FCO;
494 const Teuchos::RCP<FCO> elemToNode = Factory::Get<Teuchos::RCP<FCO>>(currentLevel,
"pcoarsen: element to node map");
496 string basisString = paramList.
get<
string>(
"pcoarsen: hi basis");
502 auto basis = MueLuIntrepid::BasisFactory<double, ES>(basisString, degree);
504 string topologyTypeString = paramList.get<
string>(
"smoother: neighborhood type");
506 if (topologyTypeString ==
"node")
508 else if (topologyTypeString ==
"edge")
510 else if (topologyTypeString ==
"face")
512 else if (topologyTypeString ==
"cell")
513 dimension = basis->getBaseCellTopology().getDimension();
515 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unrecognized smoother neighborhood type. Supported types are node, edge, face.");
516 RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_,
true);
517 vector<vector<LocalOrdinal>> seeds;
522 int myNodeCount = matA->getRowMap()->getLocalNumElements();
523 ArrayRCP<LocalOrdinal> nodeSeeds(myNodeCount, lo_invalid);
524 int localPartitionNumber = 0;
526 nodeSeeds[seed] = localPartitionNumber++;
529 paramList.remove(
"smoother: neighborhood type");
530 paramList.remove(
"pcoarsen: hi basis");
532 paramList.set(
"partitioner: map", nodeSeeds);
533 paramList.set(
"partitioner: type",
"user");
534 paramList.set(
"partitioner: overlap", 1);
535 paramList.set(
"partitioner: local parts",
int(seeds[dimension].size()));
539 type_ =
"BLOCKRELAXATION";
547 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
549 if (this->IsSetup() ==
true) {
550 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupLineSmoothing(): Setup() has already been called" << std::endl;
551 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupLineSmoothing(): reuse of this type is not available, reverting to full construction" << std::endl;
556 LO CoarseNumZLayers = Factory::Get<LO>(currentLevel,
"CoarseNumZLayers");
557 if (CoarseNumZLayers > 0) {
558 Teuchos::ArrayRCP<LO> TVertLineIdSmoo = Factory::Get<Teuchos::ArrayRCP<LO>>(currentLevel,
"LineDetection_VertLineIds");
562 for (
size_t k = 0; k < Teuchos::as<size_t>(TVertLineIdSmoo.
size()); k++) {
563 if (maxPart < TVertLineIdSmoo[k]) maxPart = TVertLineIdSmoo[k];
565 RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_,
true);
566 size_t numLocalRows = matA->getLocalNumRows();
569 "MueLu::Ifpack2Smoother::Setup(): the number of local nodes is incompatible with the TVertLineIdsSmoo.");
574 int actualDofsPerNode = numLocalRows / TVertLineIdSmoo.
size();
575 LO matrixBlockSize = matA->GetFixedBlockSize();
576 if (matrixBlockSize > 1 && actualDofsPerNode > 1) {
578 "MueLu::Ifpack2Smoother::Setup(): A is a block matrix but its block size and DOFs/node from partitioner disagree");
579 }
else if (matrixBlockSize > 1) {
580 actualDofsPerNode = matrixBlockSize;
582 myparamList.
set(
"partitioner: PDE equations", actualDofsPerNode);
584 if (numLocalRows == Teuchos::as<size_t>(TVertLineIdSmoo.
size())) {
585 myparamList.
set(
"partitioner: type",
"user");
586 myparamList.
set(
"partitioner: map", TVertLineIdSmoo);
587 myparamList.
set(
"partitioner: local parts", maxPart + 1);
589 if (myparamList.
isParameter(
"partitioner: block size") &&
590 myparamList.
get<
int>(
"partitioner: block size") != -1) {
591 int block_size = myparamList.
get<
int>(
"partitioner: block size");
593 "MueLu::Ifpack2Smoother::Setup(): the number of local nodes is incompatible with the specified block size.");
594 numLocalRows /= block_size;
598 size_t numDofsPerNode = numLocalRows / TVertLineIdSmoo.
size();
602 for (
size_t blockRow = 0; blockRow < Teuchos::as<size_t>(TVertLineIdSmoo.
size()); ++blockRow)
603 for (
size_t dof = 0; dof < numDofsPerNode; dof++)
604 partitionerMap[blockRow * numDofsPerNode + dof] = TVertLineIdSmoo[blockRow];
605 myparamList.
set(
"partitioner: type",
"user");
606 myparamList.
set(
"partitioner: map", partitionerMap);
607 myparamList.
set(
"partitioner: local parts", maxPart + 1);
610 if (type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
611 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
612 type_ ==
"LINESMOOTHING_BANDEDRELAXATION")
613 type_ =
"BANDEDRELAXATION";
614 else if (type_ ==
"LINESMOOTHING_TRIDI_RELAXATION" ||
615 type_ ==
"LINESMOOTHING_TRIDI RELAXATION" ||
616 type_ ==
"LINESMOOTHING_TRIDIRELAXATION" ||
617 type_ ==
"LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
618 type_ ==
"LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
619 type_ ==
"LINESMOOTHING_TRIDIAGONALRELAXATION")
620 type_ =
"TRIDIAGONALRELAXATION";
622 type_ =
"BLOCKRELAXATION";
625 this->GetOStream(
Runtime0) <<
"Line detection failed: fall back to point-wise relaxation" << std::endl;
626 myparamList.
remove(
"partitioner: type",
false);
627 myparamList.
remove(
"partitioner: map",
false);
628 myparamList.
remove(
"partitioner: local parts",
false);
629 type_ =
"RELAXATION";
640 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
650 bool reusePreconditioner =
false;
651 if (this->IsSetup() ==
true) {
653 this->GetOStream(
Runtime1) <<
"MueLu::Ifpack2Smoother::SetupBlockRelaxation(): Setup() has already been called, assuming reuse" << std::endl;
658 reusePreconditioner =
true;
660 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupBlockRelaxation(): reuse of this type is not available (failed cast to CanChangeMatrix), "
661 "reverting to full construction"
666 if (!reusePreconditioner) {
669 if (myparamList.
isParameter(
"partitioner: type") &&
670 myparamList.
get<std::string>(
"partitioner: type") ==
"line") {
672 Factory::Get<Teuchos::RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LO, GO,
NO>>>(currentLevel,
"Coordinates");
676 size_t lclSize = A_->getRangeMap()->getLocalNumElements();
678 lclSize = matA->getLocalNumRows();
679 size_t numDofsPerNode = lclSize / xCoordinates->getMap()->getLocalNumElements();
680 myparamList.
set(
"partitioner: coordinates", coordinates);
681 myparamList.
set(
"partitioner: PDE equations", (
int)numDofsPerNode);
692 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
702 bool reusePreconditioner =
false;
704 if (this->IsSetup() ==
true) {
706 this->GetOStream(
Runtime1) <<
"MueLu::Ifpack2Smoother::SetupChebyshev(): Setup() has already been called, assuming reuse" << std::endl;
711 reusePreconditioner =
true;
713 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupChebyshev(): reuse of this type is not available (failed cast to CanChangeMatrix), "
714 "reverting to full construction"
721 SC negone = -STS::one();
722 SC lambdaMax = SetupChebyshevEigenvalues(currentLevel,
"A",
"", paramList);
724 if (!reusePreconditioner) {
739 if (lambdaMax == negone) {
743 if (chebyPrec != Teuchos::null) {
744 lambdaMax = chebyPrec->getLambdaMaxForApply();
747 matA->SetMaxEigenvalueEstimate(lambdaMax);
748 this->GetOStream(
Statistics1) <<
"chebyshev: max eigenvalue (calculated by Ifpack2)"
749 <<
" = " << lambdaMax << std::endl;
755 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
764 bool reusePreconditioner =
false;
765 if (this->IsSetup() ==
true) {
767 this->GetOStream(
Runtime1) <<
"MueLu::Ifpack2Smoother::SetupHiptmair(): Setup() has already been called, assuming reuse" << std::endl;
772 reusePreconditioner =
true;
774 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupHiptmair(): reuse of this type is not available (failed cast to CanChangeMatrix), "
775 "reverting to full construction"
783 std::string smoother1 = paramList.
get(
"hiptmair: smoother type 1",
"CHEBYSHEV");
784 std::string smoother2 = paramList.
get(
"hiptmair: smoother type 2",
"CHEBYSHEV");
785 SC lambdaMax11 = negone;
787 if (smoother1 ==
"CHEBYSHEV") {
789 lambdaMax11 = SetupChebyshevEigenvalues(currentLevel,
"A",
"EdgeMatrix ", list1);
791 if (smoother2 ==
"CHEBYSHEV") {
793 SetupChebyshevEigenvalues(currentLevel,
"NodeMatrix",
"NodeMatrix ", list2);
807 newlist.
set(
"P", tD0);
808 newlist.
set(
"PtAP", tNodeMatrix);
809 if (!reusePreconditioner) {
811 SetPrecParameters(newlist);
818 if (smoother1 ==
"CHEBYSHEV" && lambdaMax11 == negone) {
819 using Teuchos::rcp_dynamic_cast;
822 if (hiptmairPrec != Teuchos::null) {
824 if (chebyPrec != Teuchos::null) {
828 matA->SetMaxEigenvalueEstimate(lambdaMax11);
829 this->GetOStream(
Statistics1) <<
"chebyshev: max eigenvalue (calculated by Ifpack2)"
830 <<
" = " << lambdaMax11 << std::endl;
837 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
841 SC negone = -STS::one();
843 RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(currentA);
844 SC lambdaMax = negone;
846 std::string maxEigString =
"chebyshev: max eigenvalue";
847 std::string eigRatioString =
"chebyshev: ratio eigenvalue";
851 if (paramList.
isType<
double>(maxEigString))
852 lambdaMax = paramList.
get<
double>(maxEigString);
854 lambdaMax = paramList.
get<
SC>(maxEigString);
855 this->GetOStream(
Statistics1) << label << maxEigString <<
" (cached with smoother parameter list) = " << lambdaMax << std::endl;
857 matA->SetMaxEigenvalueEstimate(lambdaMax);
860 lambdaMax = matA->GetMaxEigenvalueEstimate();
861 if (lambdaMax != negone) {
862 this->GetOStream(
Statistics1) << label << maxEigString <<
" (cached with matrix) = " << lambdaMax << std::endl;
863 paramList.
set(maxEigString, lambdaMax);
868 const SC defaultEigRatio = 20;
870 SC ratio = defaultEigRatio;
872 if (paramList.
isType<
double>(eigRatioString))
873 ratio = paramList.
get<
double>(eigRatioString);
875 ratio = paramList.
get<
SC>(eigRatioString);
883 size_t nRowsFine = fineA->getDomainMap()->getGlobalNumElements();
884 size_t nRowsCoarse = currentA->getDomainMap()->getGlobalNumElements();
886 SC levelRatio = as<SC>(as<double>(nRowsFine) / nRowsCoarse);
887 if (STS::magnitude(levelRatio) > STS::magnitude(ratio))
891 this->GetOStream(
Statistics1) << label << eigRatioString <<
" (computed) = " << ratio << std::endl;
892 paramList.
set(eigRatioString, ratio);
894 if (paramList.
isParameter(
"chebyshev: use rowsumabs diagonal scaling")) {
895 this->GetOStream(
Runtime1) <<
"chebyshev: using rowsumabs diagonal scaling" << std::endl;
896 bool doScale =
false;
897 doScale = paramList.
get<
bool>(
"chebyshev: use rowsumabs diagonal scaling");
898 paramList.
remove(
"chebyshev: use rowsumabs diagonal scaling");
900 std::string paramName =
"chebyshev: rowsumabs diagonal replacement tolerance";
902 chebyReplaceTol = paramList.
get<
double>(paramName);
903 paramList.
remove(paramName);
906 paramName =
"chebyshev: rowsumabs diagonal replacement value";
908 chebyReplaceVal = paramList.
get<
double>(paramName);
909 paramList.
remove(paramName);
911 bool chebyReplaceSingleEntryRowWithZero =
false;
912 paramName =
"chebyshev: rowsumabs replace single entry row with zero";
914 chebyReplaceSingleEntryRowWithZero = paramList.
get<
bool>(paramName);
915 paramList.
remove(paramName);
917 bool useAverageAbsDiagVal =
false;
918 paramName =
"chebyshev: rowsumabs use automatic diagonal tolerance";
920 useAverageAbsDiagVal = paramList.
get<
bool>(paramName);
921 paramList.
remove(paramName);
925 const bool doReciprocal =
true;
935 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
944 bool reusePreconditioner =
false;
945 if (this->IsSetup() ==
true) {
947 this->GetOStream(
Runtime1) <<
"MueLu::Ifpack2Smoother::SetupGeneric(): Setup() has already been called, assuming reuse" << std::endl;
952 reusePreconditioner =
true;
954 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupGeneric(): reuse of this type is not available (failed cast to CanChangeMatrix), "
955 "reverting to full construction"
960 if (!reusePreconditioner) {
969 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
980 bool supportInitialGuess =
false;
983 if (prec_->supportsZeroStartingSolution()) {
984 prec_->setZeroStartingSolution(InitialGuessIsZero);
985 supportInitialGuess =
true;
986 }
else if (type_ ==
"SCHWARZ") {
987 paramList.
set(
"schwarz: zero starting solution", InitialGuessIsZero);
993 supportInitialGuess =
true;
1003 if (InitialGuessIsZero || supportInitialGuess) {
1006 prec_->apply(tpB, tpX);
1012 std::string prefix = this->ShortClassName() +
": Apply: ";
1017 RCP<MultiVector> Correction = MultiVectorFactory::Build(A_->getDomainMap(), X.getNumVectors());
1022 prec_->apply(tpB, tpX);
1024 X.update(TST::one(), *Correction, TST::one());
1028 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1035 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1037 std::ostringstream out;
1039 out << prec_->description();
1042 out <<
"{type = " << type_ <<
"}";
1047 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1052 out0 <<
"Prec. type: " << type_ << std::endl;
1055 out0 <<
"Parameter list: " << std::endl;
1057 out << this->GetParameterList();
1058 out0 <<
"Overlap: " << overlap_ << std::endl;
1062 if (prec_ != Teuchos::null) {
1064 out << *prec_ << std::endl;
1067 if (verbLevel &
Debug) {
1070 <<
"RCP<prec_>: " << prec_ << std::endl;
1074 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1079 if (!pr.
is_null())
return pr->getNodeSmootherComplexity();
1082 if (!pc.
is_null())
return pc->getNodeSmootherComplexity();
1085 if (!pb.
is_null())
return pb->getNodeSmootherComplexity();
1088 if (!pi.
is_null())
return pi->getNodeSmootherComplexity();
1091 if (!pk.
is_null())
return pk->getNodeSmootherComplexity();
1098 #endif // HAVE_MUELU_IFPACK2
1099 #endif // MUELU_IFPACK2SMOOTHER_DEF_HPP
void DeclareInput(Level ¤tLevel) const
Input.
Important warning messages (one line)
void SetupGeneric(Level ¤tLevel)
Exception indicating invalid cast attempted.
RCP< Level > & GetPreviousLevel()
Previous level.
High level timing information (use Teuchos::TimeMonitor::summarize() to print)
MueLu::DefaultLocalOrdinal LocalOrdinal
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.
KOKKOS_INLINE_FUNCTION LO GetNumAggregates() const
static magnitudeType eps()
Print external lib objects.
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
static RCP< Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraRow(RCP< Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node >> Op)
friend class Ifpack2Smoother
Constructor.
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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.
Integrates Teuchos::TimeMonitor with MueLu verbosity system.
std::string type_
ifpack2-specific key phrase that denote smoother type
void SetupBlockRelaxation(Level ¤tLevel)
void resize(const size_type n, const T &val=T())
bool isParameter(const std::string &name) const
void SetupSchwarz(Level ¤tLevel)
void SetupLineSmoothing(Level ¤tLevel)
Print statistics that do not involve significant additional computation.
bool remove(std::string const &name, bool throwIfNotExists=true)
static RCP< MultiVector > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &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 Teuchos::ParameterList ¶mList)
Set parameters from a parameter list and return with default values.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultScalar Scalar
RCP< SmootherPrototype > Copy() const
void declareConstructionOutcome(bool fail, std::string msg)
Class that holds all level-specific information.
bool isSublist(const std::string &name) const
void Setup(Level ¤tLevel)
Set up the smoother.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
static Teuchos::ArrayRCP< const bool > DetectDirichletRows(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude &tol=Teuchos::ScalarTraits< Magnitude >::zero(), bool count_twos_as_dirichlet=false)
Detect Dirichlet rows.
void SetupHiptmair(Level ¤tLevel)
void SetupTopological(Level ¤tLevel)
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)
const RCP< LOMultiVector > & GetVertex2AggId() const
Returns constant vector that maps local node IDs to local aggregates IDs.
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
ParameterList & setParameters(const ParameterList &source)
Print all timing information.
void SetPrecParameters(const Teuchos::ParameterList &list=Teuchos::ParameterList()) const
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > toTpetra(const RCP< const CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &graph)
MatrixType::scalar_type getLambdaMaxForApply() const
void SetupAggregate(Level ¤tLevel)
bool isType(const std::string &name) const
RCP< Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetra_Vector() const
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the preconditioner.
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Teuchos::RCP< BlockCrsMatrix< Scalar, LO, GO, Node > > convertToBlockCrsMatrix(const Tpetra::CrsMatrix< Scalar, LO, GO, Node > &pointMatrix, const LO &blockSize, bool use_local_ID=true)
int GetLevelID() const
Return level number.
Print class parameters (more parameters, more verbose)
Exception throws to report errors in the internal logical of the program.
void SetupChebyshev(Level ¤tLevel)
#define TEUCHOS_ASSERT(assertion_test)
Description of what is happening (more verbose)
Class that encapsulates Ifpack2 smoothers.
static Teuchos::RCP< Vector > GetLumpedMatrixDiagonal(Matrix const &A, const bool doReciprocal=false, Magnitude tol=Teuchos::ScalarTraits< Scalar >::magnitude(Teuchos::ScalarTraits< Scalar >::zero()), Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero(), const bool replaceSingleEntryRowWithZero=false, const bool useAverageAbsDiagVal=false)
Extract Matrix Diagonal of lumped matrix.
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node >> vec)
virtual std::string description() const
Return a simple one-line description of this object.
std::string toString(const T &t)
Scalar SetupChebyshevEigenvalues(Level ¤tLevel, const std::string &matrixName, const std::string &label, Teuchos::ParameterList ¶mList) const