10 #ifndef MUELU_IFPACK2SMOOTHER_DEF_HPP
11 #define MUELU_IFPACK2SMOOTHER_DEF_HPP
16 #if defined(HAVE_MUELU_IFPACK2)
20 #include <Tpetra_RowMatrix.hpp>
22 #include <Ifpack2_Chebyshev.hpp>
23 #include <Ifpack2_Hiptmair.hpp>
24 #include <Ifpack2_RILUK.hpp>
25 #include <Ifpack2_Relaxation.hpp>
26 #include <Ifpack2_ILUT.hpp>
27 #include <Ifpack2_BlockRelaxation.hpp>
28 #include <Ifpack2_Factory.hpp>
29 #include <Ifpack2_Parameters.hpp>
33 #include <Xpetra_CrsMatrixWrap.hpp>
34 #include <Xpetra_TpetraBlockCrsMatrix.hpp>
37 #include <Xpetra_MultiVectorFactory.hpp>
38 #include <Xpetra_TpetraMultiVector.hpp>
40 #include <Tpetra_BlockCrsMatrix_Helpers.hpp>
44 #include "MueLu_Utilities.hpp"
46 #include "MueLu_Aggregates.hpp"
48 #ifdef HAVE_MUELU_INTREPID2
51 #include "Intrepid2_Basis.hpp"
52 #include "Kokkos_DynRankView.hpp"
57 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
62 bool isSupported = Ifpack2::Factory::isSupported<tRowMatrix>(
type_) || (
type_ ==
"LINESMOOTHING_TRIDI_RELAXATION" ||
63 type_ ==
"LINESMOOTHING_TRIDI RELAXATION" ||
64 type_ ==
"LINESMOOTHING_TRIDIRELAXATION" ||
65 type_ ==
"LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
66 type_ ==
"LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
67 type_ ==
"LINESMOOTHING_TRIDIAGONALRELAXATION" ||
68 type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
69 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
70 type_ ==
"LINESMOOTHING_BANDEDRELAXATION" ||
71 type_ ==
"LINESMOOTHING_BLOCK_RELAXATION" ||
72 type_ ==
"LINESMOOTHING_BLOCK RELAXATION" ||
73 type_ ==
"LINESMOOTHING_BLOCKRELAXATION" ||
74 type_ ==
"TOPOLOGICAL" ||
75 type_ ==
"AGGREGATE");
81 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
92 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
94 std::string prefix = this->ShortClassName() +
": SetPrecParameters";
100 RCP<ParameterList> precList = this->RemoveFactoriesFromList(this->GetParameterList());
103 precList->
set(
"timer for apply", this->IsPrint(
Timings));
105 if (!precList.
is_null() && precList->
isParameter(
"partitioner: type") && precList->
get<std::string>(
"partitioner: type") ==
"linear" &&
106 !precList->
isParameter(
"partitioner: local parts")) {
107 LO matrixBlockSize = 1;
108 int lclSize = A_->getRangeMap()->getLocalNumElements();
111 lclSize = matA->getLocalNumRows();
112 matrixBlockSize = matA->GetFixedBlockSize();
114 precList->
set(
"partitioner: local parts", lclSize / matrixBlockSize);
122 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
124 this->Input(currentLevel,
"A");
126 if (type_ ==
"LINESMOOTHING_TRIDI_RELAXATION" ||
127 type_ ==
"LINESMOOTHING_TRIDI RELAXATION" ||
128 type_ ==
"LINESMOOTHING_TRIDIRELAXATION" ||
129 type_ ==
"LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
130 type_ ==
"LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
131 type_ ==
"LINESMOOTHING_TRIDIAGONALRELAXATION" ||
132 type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
133 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
134 type_ ==
"LINESMOOTHING_BANDEDRELAXATION" ||
135 type_ ==
"LINESMOOTHING_BLOCK_RELAXATION" ||
136 type_ ==
"LINESMOOTHING_BLOCK RELAXATION" ||
137 type_ ==
"LINESMOOTHING_BLOCKRELAXATION") {
138 this->Input(currentLevel,
"CoarseNumZLayers");
139 this->Input(currentLevel,
"LineDetection_VertLineIds");
140 }
else if (type_ ==
"BLOCK RELAXATION" ||
141 type_ ==
"BLOCK_RELAXATION" ||
142 type_ ==
"BLOCKRELAXATION" ||
144 type_ ==
"BANDED_RELAXATION" ||
145 type_ ==
"BANDED RELAXATION" ||
146 type_ ==
"BANDEDRELAXATION" ||
148 type_ ==
"TRIDI_RELAXATION" ||
149 type_ ==
"TRIDI RELAXATION" ||
150 type_ ==
"TRIDIRELAXATION" ||
151 type_ ==
"TRIDIAGONAL_RELAXATION" ||
152 type_ ==
"TRIDIAGONAL RELAXATION" ||
153 type_ ==
"TRIDIAGONALRELAXATION") {
157 precList.
get<std::string>(
"partitioner: type") ==
"line") {
158 this->Input(currentLevel,
"Coordinates");
160 }
else if (type_ ==
"TOPOLOGICAL") {
162 this->Input(currentLevel,
"pcoarsen: element to node map");
163 }
else if (type_ ==
"AGGREGATE") {
165 this->Input(currentLevel,
"Aggregates");
166 }
else if (type_ ==
"HIPTMAIR") {
168 this->Input(currentLevel,
"NodeMatrix");
169 this->Input(currentLevel,
"D0");
173 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
176 A_ = Factory::Get<RCP<Operator>>(currentLevel,
"A");
181 if (paramList.isParameter(
"smoother: use blockcrsmatrix storage") && paramList.get<
bool>(
"smoother: use blockcrsmatrix storage")) {
185 blocksize = matA->GetFixedBlockSize();
191 throw std::runtime_error(
"Ifpack2Smoother: Cannot convert matrix A to CrsMatrixWrap object.");
194 throw std::runtime_error(
"Ifpack2Smoother: Cannot extract CrsMatrix from matrix A.");
198 throw std::runtime_error(
"Ifpack2Smoother: Cannot extract CrsMatrix or BlockCrsMatrix from matrix A.");
199 this->GetOStream(
Statistics0) <<
"Ifpack2Smoother: Using (native) BlockCrsMatrix storage with blocksize " << blocksize << std::endl;
200 paramList.remove(
"smoother: use blockcrsmatrix storage");
206 this->GetOStream(
Statistics0) <<
"Ifpack2Smoother: Using BlockCrsMatrix storage with blocksize " << blocksize << std::endl;
208 paramList.remove(
"smoother: use blockcrsmatrix storage");
213 if (type_ ==
"SCHWARZ")
214 SetupSchwarz(currentLevel);
216 else if (type_ ==
"LINESMOOTHING_TRIDI_RELAXATION" ||
217 type_ ==
"LINESMOOTHING_TRIDI RELAXATION" ||
218 type_ ==
"LINESMOOTHING_TRIDIRELAXATION" ||
219 type_ ==
"LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
220 type_ ==
"LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
221 type_ ==
"LINESMOOTHING_TRIDIAGONALRELAXATION" ||
222 type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
223 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
224 type_ ==
"LINESMOOTHING_BANDEDRELAXATION" ||
225 type_ ==
"LINESMOOTHING_BLOCK_RELAXATION" ||
226 type_ ==
"LINESMOOTHING_BLOCK RELAXATION" ||
227 type_ ==
"LINESMOOTHING_BLOCKRELAXATION")
228 SetupLineSmoothing(currentLevel);
230 else if (type_ ==
"BLOCK_RELAXATION" ||
231 type_ ==
"BLOCK RELAXATION" ||
232 type_ ==
"BLOCKRELAXATION" ||
234 type_ ==
"BANDED_RELAXATION" ||
235 type_ ==
"BANDED RELAXATION" ||
236 type_ ==
"BANDEDRELAXATION" ||
238 type_ ==
"TRIDI_RELAXATION" ||
239 type_ ==
"TRIDI RELAXATION" ||
240 type_ ==
"TRIDIRELAXATION" ||
241 type_ ==
"TRIDIAGONAL_RELAXATION" ||
242 type_ ==
"TRIDIAGONAL RELAXATION" ||
243 type_ ==
"TRIDIAGONALRELAXATION")
244 SetupBlockRelaxation(currentLevel);
246 else if (type_ ==
"CHEBYSHEV")
247 SetupChebyshev(currentLevel);
249 else if (type_ ==
"TOPOLOGICAL") {
250 #ifdef HAVE_MUELU_INTREPID2
251 SetupTopological(currentLevel);
255 }
else if (type_ ==
"AGGREGATE")
256 SetupAggregate(currentLevel);
258 else if (type_ ==
"HIPTMAIR")
259 SetupHiptmair(currentLevel);
262 SetupGeneric(currentLevel);
267 this->GetOStream(
Statistics1) << description() << std::endl;
270 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
274 bool reusePreconditioner =
false;
275 if (this->IsSetup() ==
true) {
277 this->GetOStream(
Runtime1) <<
"MueLu::Ifpack2Smoother::SetupSchwarz(): Setup() has already been called, assuming reuse" << std::endl;
279 bool isTRowMatrix =
true;
282 tA = Xpetra::toTpetraRowMatrix(A_);
284 isTRowMatrix =
false;
288 if (!prec.
is_null() && isTRowMatrix) {
290 reusePreconditioner =
true;
292 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupSchwarz(): reuse of this type is not available "
293 "(either failed cast to CanChangeMatrix, or to Tpetra Row Matrix), reverting to full construction"
298 if (!reusePreconditioner) {
301 bool isBlockedMatrix =
false;
304 std::string sublistName =
"subdomain solver parameters";
316 std::string partName =
"partitioner: type";
322 if (subList.
isParameter(partName) && subList.
get<std::string>(partName) ==
"vanka user") {
323 isBlockedMatrix =
true;
327 "Matrix A must be of type BlockedCrsMatrix.");
329 size_t numVels = bA->getMatrix(0, 0)->getLocalNumRows();
330 size_t numPres = bA->getMatrix(1, 0)->getLocalNumRows();
331 size_t numRows = rcp_dynamic_cast<Matrix>(A_,
true)->getLocalNumRows();
335 size_t numBlocks = 0;
336 for (
size_t rowOfB = numVels; rowOfB < numVels + numPres; ++rowOfB)
337 blockSeeds[rowOfB] = numBlocks++;
341 "Matrix A must be of type BlockedCrsMatrix.");
343 merged2Mat = bA2->Merge();
347 bool haveBoundary =
false;
348 for (
LO i = 0; i < boundaryNodes.
size(); i++)
349 if (boundaryNodes[i]) {
353 blockSeeds[i] = numBlocks;
359 subList.
set(
"partitioner: type",
"user");
360 subList.
set(
"partitioner: map", blockSeeds);
361 subList.
set(
"partitioner: local parts", as<int>(numBlocks));
366 isBlockedMatrix =
true;
367 merged2Mat = bA->Merge();
373 if (isBlockedMatrix ==
true)
374 tA = Xpetra::toTpetraRowMatrix(merged2Mat);
376 tA = Xpetra::toTpetraRowMatrix(A_);
387 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
391 if (this->IsSetup() ==
true) {
392 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupAggregate(): Setup() has already been called" << std::endl;
393 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupAggregate(): reuse of this type is not available, reverting to full construction" << std::endl;
396 this->GetOStream(
Statistics0) <<
"Ifpack2Smoother: Using Aggregate Smoothing" << std::endl;
398 RCP<Aggregates> aggregates = Factory::Get<RCP<Aggregates>>(currentLevel,
"Aggregates");
401 ArrayRCP<LO> aggregate_ids = rcp_const_cast<LOMultiVector>(vertex2AggId)->getDataNonConst(0);
408 blocksize = matA->GetFixedBlockSize();
410 dof_ids.
resize(aggregate_ids.
size() * blocksize);
411 for (
LO i = 0; i < (
LO)aggregate_ids.
size(); i++) {
412 for (
LO j = 0; j < (
LO)blocksize; j++)
413 dof_ids[i * blocksize + j] = aggregate_ids[i];
416 dof_ids = aggregate_ids;
419 paramList.
set(
"partitioner: map", dof_ids);
420 paramList.
set(
"partitioner: type",
"user");
421 paramList.
set(
"partitioner: overlap", 0);
426 type_ =
"BLOCKRELAXATION";
433 #ifdef HAVE_MUELU_INTREPID2
434 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
444 if (this->IsSetup() ==
true) {
445 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupTopological(): Setup() has already been called" << std::endl;
446 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupTopological(): reuse of this type is not available, reverting to full construction" << std::endl;
449 ParameterList& paramList =
const_cast<ParameterList&
>(this->GetParameterList());
451 typedef typename Node::device_type::execution_space ES;
453 typedef Kokkos::DynRankView<LocalOrdinal, typename Node::device_type> FCO;
459 const Teuchos::RCP<FCO> elemToNode = Factory::Get<Teuchos::RCP<FCO>>(currentLevel,
"pcoarsen: element to node map");
461 string basisString = paramList.
get<
string>(
"pcoarsen: hi basis");
467 auto basis = MueLuIntrepid::BasisFactory<double, ES>(basisString, degree);
469 string topologyTypeString = paramList.get<
string>(
"smoother: neighborhood type");
471 if (topologyTypeString ==
"node")
473 else if (topologyTypeString ==
"edge")
475 else if (topologyTypeString ==
"face")
477 else if (topologyTypeString ==
"cell")
478 dimension = basis->getBaseCellTopology().getDimension();
480 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
"Unrecognized smoother neighborhood type. Supported types are node, edge, face.");
481 RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_,
true);
482 vector<vector<LocalOrdinal>> seeds;
487 int myNodeCount = matA->getRowMap()->getLocalNumElements();
488 ArrayRCP<LocalOrdinal> nodeSeeds(myNodeCount, lo_invalid);
489 int localPartitionNumber = 0;
491 nodeSeeds[seed] = localPartitionNumber++;
494 paramList.remove(
"smoother: neighborhood type");
495 paramList.remove(
"pcoarsen: hi basis");
497 paramList.set(
"partitioner: map", nodeSeeds);
498 paramList.set(
"partitioner: type",
"user");
499 paramList.set(
"partitioner: overlap", 1);
500 paramList.set(
"partitioner: local parts",
int(seeds[dimension].size()));
502 RCP<const Tpetra::RowMatrix<SC, LO, GO, NO>> tA = Xpetra::toTpetraRowMatrix(A_);
504 type_ =
"BLOCKRELAXATION";
512 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
514 if (this->IsSetup() ==
true) {
515 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupLineSmoothing(): Setup() has already been called" << std::endl;
516 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupLineSmoothing(): reuse of this type is not available, reverting to full construction" << std::endl;
521 LO CoarseNumZLayers = Factory::Get<LO>(currentLevel,
"CoarseNumZLayers");
522 if (CoarseNumZLayers > 0) {
523 Teuchos::ArrayRCP<LO> TVertLineIdSmoo = Factory::Get<Teuchos::ArrayRCP<LO>>(currentLevel,
"LineDetection_VertLineIds");
527 for (
size_t k = 0; k < Teuchos::as<size_t>(TVertLineIdSmoo.
size()); k++) {
528 if (maxPart < TVertLineIdSmoo[k]) maxPart = TVertLineIdSmoo[k];
530 RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(A_,
true);
531 size_t numLocalRows = matA->getLocalNumRows();
534 "MueLu::Ifpack2Smoother::Setup(): the number of local nodes is incompatible with the TVertLineIdsSmoo.");
539 int actualDofsPerNode = numLocalRows / TVertLineIdSmoo.
size();
540 LO matrixBlockSize = matA->GetFixedBlockSize();
541 if (matrixBlockSize > 1 && actualDofsPerNode > 1) {
543 "MueLu::Ifpack2Smoother::Setup(): A is a block matrix but its block size and DOFs/node from partitioner disagree");
544 }
else if (matrixBlockSize > 1) {
545 actualDofsPerNode = matrixBlockSize;
547 myparamList.
set(
"partitioner: PDE equations", actualDofsPerNode);
549 if (numLocalRows == Teuchos::as<size_t>(TVertLineIdSmoo.
size())) {
550 myparamList.
set(
"partitioner: type",
"user");
551 myparamList.
set(
"partitioner: map", TVertLineIdSmoo);
552 myparamList.
set(
"partitioner: local parts", maxPart + 1);
554 if (myparamList.
isParameter(
"partitioner: block size") &&
555 myparamList.
get<
int>(
"partitioner: block size") != -1) {
556 int block_size = myparamList.
get<
int>(
"partitioner: block size");
558 "MueLu::Ifpack2Smoother::Setup(): the number of local nodes is incompatible with the specified block size.");
559 numLocalRows /= block_size;
563 size_t numDofsPerNode = numLocalRows / TVertLineIdSmoo.
size();
567 for (
size_t blockRow = 0; blockRow < Teuchos::as<size_t>(TVertLineIdSmoo.
size()); ++blockRow)
568 for (
size_t dof = 0; dof < numDofsPerNode; dof++)
569 partitionerMap[blockRow * numDofsPerNode + dof] = TVertLineIdSmoo[blockRow];
570 myparamList.
set(
"partitioner: type",
"user");
571 myparamList.
set(
"partitioner: map", partitionerMap);
572 myparamList.
set(
"partitioner: local parts", maxPart + 1);
575 if (type_ ==
"LINESMOOTHING_BANDED_RELAXATION" ||
576 type_ ==
"LINESMOOTHING_BANDED RELAXATION" ||
577 type_ ==
"LINESMOOTHING_BANDEDRELAXATION")
578 type_ =
"BANDEDRELAXATION";
579 else if (type_ ==
"LINESMOOTHING_TRIDI_RELAXATION" ||
580 type_ ==
"LINESMOOTHING_TRIDI RELAXATION" ||
581 type_ ==
"LINESMOOTHING_TRIDIRELAXATION" ||
582 type_ ==
"LINESMOOTHING_TRIDIAGONAL_RELAXATION" ||
583 type_ ==
"LINESMOOTHING_TRIDIAGONAL RELAXATION" ||
584 type_ ==
"LINESMOOTHING_TRIDIAGONALRELAXATION")
585 type_ =
"TRIDIAGONALRELAXATION";
587 type_ =
"BLOCKRELAXATION";
590 this->GetOStream(
Runtime0) <<
"Line detection failed: fall back to point-wise relaxation" << std::endl;
591 myparamList.
remove(
"partitioner: type",
false);
592 myparamList.
remove(
"partitioner: map",
false);
593 myparamList.
remove(
"partitioner: local parts",
false);
594 type_ =
"RELAXATION";
605 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
615 bool reusePreconditioner =
false;
616 if (this->IsSetup() ==
true) {
618 this->GetOStream(
Runtime1) <<
"MueLu::Ifpack2Smoother::SetupBlockRelaxation(): Setup() has already been called, assuming reuse" << std::endl;
623 reusePreconditioner =
true;
625 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupBlockRelaxation(): reuse of this type is not available (failed cast to CanChangeMatrix), "
626 "reverting to full construction"
631 if (!reusePreconditioner) {
634 if (myparamList.
isParameter(
"partitioner: type") &&
635 myparamList.
get<std::string>(
"partitioner: type") ==
"line") {
637 Factory::Get<Teuchos::RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LO, GO,
NO>>>(currentLevel,
"Coordinates");
641 size_t lclSize = A_->getRangeMap()->getLocalNumElements();
643 lclSize = matA->getLocalNumRows();
644 size_t numDofsPerNode = lclSize / xCoordinates->getMap()->getLocalNumElements();
645 myparamList.
set(
"partitioner: coordinates", coordinates);
646 myparamList.
set(
"partitioner: PDE equations", (
int)numDofsPerNode);
657 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
667 bool reusePreconditioner =
false;
669 if (this->IsSetup() ==
true) {
671 this->GetOStream(
Runtime1) <<
"MueLu::Ifpack2Smoother::SetupChebyshev(): Setup() has already been called, assuming reuse" << std::endl;
676 reusePreconditioner =
true;
678 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupChebyshev(): reuse of this type is not available (failed cast to CanChangeMatrix), "
679 "reverting to full construction"
686 SC negone = -STS::one();
687 SC lambdaMax = SetupChebyshevEigenvalues(currentLevel,
"A",
"", paramList);
689 if (!reusePreconditioner) {
704 if (lambdaMax == negone) {
708 if (chebyPrec != Teuchos::null) {
709 lambdaMax = chebyPrec->getLambdaMaxForApply();
712 matA->SetMaxEigenvalueEstimate(lambdaMax);
713 this->GetOStream(
Statistics1) <<
"chebyshev: max eigenvalue (calculated by Ifpack2)"
714 <<
" = " << lambdaMax << std::endl;
720 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
729 bool reusePreconditioner =
false;
730 if (this->IsSetup() ==
true) {
732 this->GetOStream(
Runtime1) <<
"MueLu::Ifpack2Smoother::SetupHiptmair(): Setup() has already been called, assuming reuse" << std::endl;
737 reusePreconditioner =
true;
739 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupHiptmair(): reuse of this type is not available (failed cast to CanChangeMatrix), "
740 "reverting to full construction"
748 std::string smoother1 = paramList.
get(
"hiptmair: smoother type 1",
"CHEBYSHEV");
749 std::string smoother2 = paramList.
get(
"hiptmair: smoother type 2",
"CHEBYSHEV");
750 SC lambdaMax11 = negone;
752 if (smoother1 ==
"CHEBYSHEV") {
754 lambdaMax11 = SetupChebyshevEigenvalues(currentLevel,
"A",
"EdgeMatrix ", list1);
756 if (smoother2 ==
"CHEBYSHEV") {
758 SetupChebyshevEigenvalues(currentLevel,
"NodeMatrix",
"NodeMatrix ", list2);
772 newlist.
set(
"P", tD0);
773 newlist.
set(
"PtAP", tNodeMatrix);
774 if (!reusePreconditioner) {
776 SetPrecParameters(newlist);
783 if (smoother1 ==
"CHEBYSHEV" && lambdaMax11 == negone) {
784 using Teuchos::rcp_dynamic_cast;
787 if (hiptmairPrec != Teuchos::null) {
789 if (chebyPrec != Teuchos::null) {
793 matA->SetMaxEigenvalueEstimate(lambdaMax11);
794 this->GetOStream(
Statistics1) <<
"chebyshev: max eigenvalue (calculated by Ifpack2)"
795 <<
" = " << lambdaMax11 << std::endl;
802 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
806 SC negone = -STS::one();
808 RCP<Matrix> matA = rcp_dynamic_cast<Matrix>(currentA);
809 SC lambdaMax = negone;
811 std::string maxEigString =
"chebyshev: max eigenvalue";
812 std::string eigRatioString =
"chebyshev: ratio eigenvalue";
816 if (paramList.
isType<
double>(maxEigString))
817 lambdaMax = paramList.
get<
double>(maxEigString);
819 lambdaMax = paramList.
get<
SC>(maxEigString);
820 this->GetOStream(
Statistics1) << label << maxEigString <<
" (cached with smoother parameter list) = " << lambdaMax << std::endl;
822 matA->SetMaxEigenvalueEstimate(lambdaMax);
825 lambdaMax = matA->GetMaxEigenvalueEstimate();
826 if (lambdaMax != negone) {
827 this->GetOStream(
Statistics1) << label << maxEigString <<
" (cached with matrix) = " << lambdaMax << std::endl;
828 paramList.
set(maxEigString, lambdaMax);
833 const SC defaultEigRatio = 20;
835 SC ratio = defaultEigRatio;
837 if (paramList.
isType<
double>(eigRatioString))
838 ratio = paramList.
get<
double>(eigRatioString);
840 ratio = paramList.
get<
SC>(eigRatioString);
848 size_t nRowsFine = fineA->getDomainMap()->getGlobalNumElements();
849 size_t nRowsCoarse = currentA->getDomainMap()->getGlobalNumElements();
851 SC levelRatio = as<SC>(as<double>(nRowsFine) / nRowsCoarse);
852 if (STS::magnitude(levelRatio) > STS::magnitude(ratio))
856 this->GetOStream(
Statistics1) << label << eigRatioString <<
" (computed) = " << ratio << std::endl;
857 paramList.
set(eigRatioString, ratio);
859 if (paramList.
isParameter(
"chebyshev: use rowsumabs diagonal scaling")) {
860 this->GetOStream(
Runtime1) <<
"chebyshev: using rowsumabs diagonal scaling" << std::endl;
861 bool doScale =
false;
862 doScale = paramList.
get<
bool>(
"chebyshev: use rowsumabs diagonal scaling");
863 paramList.
remove(
"chebyshev: use rowsumabs diagonal scaling");
865 std::string paramName =
"chebyshev: rowsumabs diagonal replacement tolerance";
867 chebyReplaceTol = paramList.
get<
double>(paramName);
868 paramList.
remove(paramName);
871 paramName =
"chebyshev: rowsumabs diagonal replacement value";
873 chebyReplaceVal = paramList.
get<
double>(paramName);
874 paramList.
remove(paramName);
876 bool chebyReplaceSingleEntryRowWithZero =
false;
877 paramName =
"chebyshev: rowsumabs replace single entry row with zero";
879 chebyReplaceSingleEntryRowWithZero = paramList.
get<
bool>(paramName);
880 paramList.
remove(paramName);
882 bool useAverageAbsDiagVal =
false;
883 paramName =
"chebyshev: rowsumabs use automatic diagonal tolerance";
885 useAverageAbsDiagVal = paramList.
get<
bool>(paramName);
886 paramList.
remove(paramName);
890 const bool doReciprocal =
true;
900 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
909 bool reusePreconditioner =
false;
910 if (this->IsSetup() ==
true) {
912 this->GetOStream(
Runtime1) <<
"MueLu::Ifpack2Smoother::SetupGeneric(): Setup() has already been called, assuming reuse" << std::endl;
917 reusePreconditioner =
true;
919 this->GetOStream(
Warnings0) <<
"MueLu::Ifpack2Smoother::SetupGeneric(): reuse of this type is not available (failed cast to CanChangeMatrix), "
920 "reverting to full construction"
925 if (!reusePreconditioner) {
934 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
944 bool supportInitialGuess =
false;
946 if (prec_->supportsZeroStartingSolution()) {
947 prec_->setZeroStartingSolution(InitialGuessIsZero);
948 supportInitialGuess =
true;
949 }
else if (type_ ==
"SCHWARZ") {
951 paramList.
set(
"schwarz: zero starting solution", InitialGuessIsZero);
957 supportInitialGuess =
true;
967 if (InitialGuessIsZero || supportInitialGuess) {
970 prec_->apply(tpB, tpX);
976 std::string prefix = this->ShortClassName() +
": Apply: ";
981 RCP<MultiVector> Correction = MultiVectorFactory::Build(A_->getDomainMap(), X.getNumVectors());
986 prec_->apply(tpB, tpX);
988 X.update(TST::one(), *Correction, TST::one());
992 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
999 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1001 std::ostringstream out;
1003 out << prec_->description();
1006 out <<
"{type = " << type_ <<
"}";
1011 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1016 out0 <<
"Prec. type: " << type_ << std::endl;
1019 out0 <<
"Parameter list: " << std::endl;
1021 out << this->GetParameterList();
1022 out0 <<
"Overlap: " << overlap_ << std::endl;
1026 if (prec_ != Teuchos::null) {
1028 out << *prec_ << std::endl;
1031 if (verbLevel &
Debug) {
1034 <<
"RCP<prec_>: " << prec_ << std::endl;
1038 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1043 if (!pr.
is_null())
return pr->getNodeSmootherComplexity();
1046 if (!pc.
is_null())
return pc->getNodeSmootherComplexity();
1049 if (!pb.
is_null())
return pb->getNodeSmootherComplexity();
1052 if (!pi.
is_null())
return pi->getNodeSmootherComplexity();
1055 if (!pk.
is_null())
return pk->getNodeSmootherComplexity();
1062 #endif // HAVE_MUELU_IFPACK2
1063 #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.
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)
KOKKOS_INLINE_FUNCTION LO GetNumAggregates() const
static magnitudeType eps()
Print external lib objects.
T & get(const std::string &name, T def_value)
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.
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
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)
RCP< Level > & GetPreviousLevel()
Previous level.
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.
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