47 #ifndef MUELU_HIERARCHY_DEF_HPP 48 #define MUELU_HIERARCHY_DEF_HPP 55 #include <Xpetra_Matrix.hpp> 56 #include <Xpetra_MultiVectorFactory.hpp> 57 #include <Xpetra_Operator.hpp> 58 #include <Xpetra_IO.hpp> 63 #include "MueLu_FactoryManager.hpp" 64 #include "MueLu_HierarchyUtils.hpp" 65 #include "MueLu_TopRAPFactory.hpp" 66 #include "MueLu_TopSmootherFactory.hpp" 69 #include "MueLu_PerfUtils.hpp" 70 #include "MueLu_PFactory.hpp" 72 #include "MueLu_SmootherFactory.hpp" 81 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
83 : maxCoarseSize_(GetDefaultMaxCoarseSize()), implicitTranspose_(GetDefaultImplicitTranspose()),
84 doPRrebalance_(GetDefaultPRrebalance()), isPreconditioner_(true), Cycle_(GetDefaultCycle()),
85 scalingFactor_(
Teuchos::ScalarTraits<double>::one()), lib_(Xpetra::UseTpetra), isDumpingEnabled_(false), dumpLevel_(-1), rate_(-1),
86 sizeOfAllocatedLevelMultiVectors_(0)
91 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
96 Levels_[0]->setObjectLabel(label);
99 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
101 : maxCoarseSize_(GetDefaultMaxCoarseSize()), implicitTranspose_(GetDefaultImplicitTranspose()),
102 doPRrebalance_(GetDefaultPRrebalance()), isPreconditioner_(true), Cycle_(GetDefaultCycle()),
103 scalingFactor_(
Teuchos::ScalarTraits<double>::one()), isDumpingEnabled_(false), dumpLevel_(-1), rate_(-1),
104 sizeOfAllocatedLevelMultiVectors_(0)
106 lib_ = A->getDomainMap()->lib();
114 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
119 Levels_[0]->setObjectLabel(label);
122 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
124 int levelID = LastLevelID() + 1;
128 " have been added at the end of the hierarchy\n but its ID have been redefined" <<
129 " because last level ID of the hierarchy was " << LastLevelID() <<
"." << std::endl;
131 Levels_.push_back(level);
135 level->
SetPreviousLevel( (levelID == 0) ? Teuchos::null : Levels_[LastLevelID() - 1] );
139 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
141 RCP<Level> newLevel = Levels_[LastLevelID()]->Build();
143 this->AddLevel(newLevel);
146 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
149 "MueLu::Hierarchy::GetLevel(): invalid input parameter value: LevelID = " << levelID);
150 return Levels_[levelID];
153 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
155 return Levels_.size();
158 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
160 RCP<Operator> A = Levels_[0]->template Get<RCP<Operator> >(
"A");
163 int numLevels = GetNumLevels();
165 Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, numLevels, Teuchos::ptr(&numGlobalLevels));
167 return numGlobalLevels;
170 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
172 double totalNnz = 0, lev0Nnz = 1;
173 for (
int i = 0; i < GetNumLevels(); ++i) {
175 "Operator complexity cannot be calculated because A is unavailable on level " << i);
176 RCP<Operator> A = Levels_[i]->template Get<RCP<Operator> >(
"A");
182 GetOStream(
Warnings0) <<
"Some level operators are not matrices, operator complexity calculation aborted" << std::endl;
186 totalNnz += as<double>(Am->getGlobalNumEntries());
190 return totalNnz / lev0Nnz;
193 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
195 double node_sc = 0, global_sc=0;
199 if (GetNumLevels() <= 0)
return -1.0;
200 if (!Levels_[0]->IsAvailable(
"A"))
return -1.0;
202 RCP<Operator> A = Levels_[0]->template Get<RCP<Operator> >(
"A");
206 a0_nnz = as<double>(Am->getGlobalNumEntries());
209 for (
int i = 0; i < GetNumLevels(); ++i) {
211 if(!Levels_[i]->IsAvailable(
"PreSmoother"))
continue;
212 RCP<SmootherBase> S = Levels_[i]->template Get<RCP<SmootherBase> >(
"PreSmoother");
215 if(level_sc == INVALID) {global_sc=-1.0;
break;}
217 node_sc += as<double>(level_sc);
225 if(min_sc < 0.0)
return -1.0;
226 else return global_sc / a0_nnz;
233 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
236 "MueLu::Hierarchy::CheckLevel(): wrong underlying linear algebra library.");
238 "MueLu::Hierarchy::CheckLevel(): wrong level ID");
240 "MueLu::Hierarchy::Setup(): wrong level parent");
245 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
254 "MueLu::Hierarchy:Setup(): level " << coarseLevelID <<
" (specified by coarseLevelID argument) " 255 "must be built before calling this function.");
257 Level& level = *Levels_[coarseLevelID];
260 TimeMonitor m1(*
this, label + this->ShortClassName() +
": " +
"Setup (total)");
265 "MueLu::Hierarchy::Setup(): argument coarseLevelManager cannot be null");
270 if (levelManagers_.size() < coarseLevelID+1)
271 levelManagers_.resize(coarseLevelID+1);
272 levelManagers_[coarseLevelID] = coarseLevelManager;
274 bool isFinestLevel = (fineLevelManager.
is_null());
275 bool isLastLevel = (nextLevelManager.
is_null());
288 oldRank = SetProcRankVerbose(comm->getRank());
292 lib_ = domainMap->lib();
299 Level& prevLevel = *Levels_[coarseLevelID-1];
300 oldRank = SetProcRankVerbose(prevLevel.
GetComm()->getRank());
303 CheckLevel(level, coarseLevelID);
310 if (isFinestLevel && Levels_[coarseLevelID]->IsAvailable(
"Coordinates"))
311 ReplaceCoordinateMap(*Levels_[coarseLevelID]);
316 if (isDumpingEnabled_ && dumpLevel_ == 0 && coarseLevelID == 1)
322 int nextLevelID = coarseLevelID + 1;
325 if (isLastLevel ==
false) {
327 if (nextLevelID > LastLevelID())
329 CheckLevel(*Levels_[nextLevelID], nextLevelID);
333 Levels_[nextLevelID]->Request(
TopRAPFactory(coarseLevelManager, nextLevelManager));
369 TopRAPFactory coarseRAPFactory(fineLevelManager, coarseLevelManager);
373 }
else if (!isFinestLevel) {
384 level.
SetComm(Ac->getDomainMap()->getComm());
387 bool isOrigLastLevel = isLastLevel;
399 if (!Acm.
is_null() && Acm->getGlobalNumRows() <= maxCoarseSize_) {
401 GetOStream(
Runtime0) <<
"Max coarse size (<= " << maxCoarseSize_ <<
") achieved" << std::endl;
406 if (!Ac.
is_null() && !isFinestLevel) {
407 RCP<Operator> A = Levels_[coarseLevelID-1]->template Get< RCP<Operator> >(
"A");
410 const double maxCoarse2FineRatio = 0.8;
411 if (!Acm.
is_null() && !Am.is_null() && Acm->getGlobalNumRows() > maxCoarse2FineRatio * Am->getGlobalNumRows()) {
419 GetOStream(
Warnings0) <<
"Aggregation stagnated. Please check your matrix and/or adjust your configuration file." 420 <<
"Possible fixes:\n" 421 <<
" - reduce the maximum number of levels\n" 422 <<
" - enable repartitioning\n" 423 <<
" - increase the minimum coarse size." << std::endl;
429 if (!isOrigLastLevel) {
439 coarseFact->
Build(level);
450 smootherFact->
Build(level);
455 if (isLastLevel ==
true) {
456 if (isOrigLastLevel ==
false) {
459 Levels_[nextLevelID]->Release(
TopRAPFactory(coarseLevelManager, nextLevelManager));
461 Levels_.resize(nextLevelID);
465 if (isDumpingEnabled_ && dumpLevel_ > 0 && coarseLevelID == dumpLevel_)
468 if (!isFinestLevel) {
472 level.
Release(coarseRAPFactory);
476 SetProcRankVerbose(oldRank);
481 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
483 int numLevels = Levels_.size();
485 "Hierarchy::SetupRe: " << Levels_.size() <<
" levels, but " << levelManagers_.size() <<
" level factory managers");
487 const int startLevel = 0;
490 #ifdef HAVE_MUELU_DEBUG 492 for (
int i = 0; i < numLevels; i++)
493 levelManagers_[i]->ResetDebugData();
498 for (levelID = startLevel; levelID < numLevels;) {
499 bool r = Setup(levelID,
500 (levelID != 0 ? levelManagers_[levelID-1] : Teuchos::null),
501 levelManagers_[levelID],
502 (levelID+1 != numLevels ? levelManagers_[levelID+1] : Teuchos::null));
508 Levels_ .resize(levelID);
509 levelManagers_.resize(levelID);
521 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
530 "MueLu::Hierarchy::Setup(): fine level (" << startLevel <<
") does not exist");
533 "Constructing non-positive (" << numDesiredLevels <<
") number of levels does not make sense.");
537 "MueLu::Hierarchy::Setup(): fine level (" << startLevel <<
") has no matrix A! " 538 "Set fine level matrix A using Level.Set()");
540 RCP<Operator> A = Levels_[startLevel]->template Get<RCP<Operator> >(
"A");
541 lib_ = A->getDomainMap()->lib();
548 params->
set(
"printLoadBalancingInfo",
true);
549 params->
set(
"printCommInfo",
true);
553 GetOStream(
Warnings1) <<
"Fine level operator is not a matrix, statistics are not available" << std::endl;
559 const int lastLevel = startLevel + numDesiredLevels - 1;
560 GetOStream(
Runtime0) <<
"Setup loop: startLevel = " << startLevel <<
", lastLevel = " << lastLevel
561 <<
" (stop if numLevels = " << numDesiredLevels <<
" or Ac.size() < " << maxCoarseSize_ <<
")" << std::endl;
565 if (numDesiredLevels == 1) {
567 Setup(startLevel, Teuchos::null, rcpmanager, Teuchos::null);
570 bool bIsLastLevel = Setup(startLevel, Teuchos::null, rcpmanager, rcpmanager);
571 if (bIsLastLevel ==
false) {
572 for (iLevel = startLevel + 1; iLevel < lastLevel; iLevel++) {
573 bIsLastLevel = Setup(iLevel, rcpmanager, rcpmanager, rcpmanager);
574 if (bIsLastLevel ==
true)
577 if (bIsLastLevel ==
false)
578 Setup(lastLevel, rcpmanager, rcpmanager, Teuchos::null);
584 "MueLu::Hierarchy::Setup(): number of level");
593 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
595 if (startLevel < GetNumLevels())
596 GetOStream(
Runtime0) <<
"Clearing old data (if any)" << std::endl;
598 for (
int iLevel = startLevel; iLevel < GetNumLevels(); iLevel++)
599 Levels_[iLevel]->Clear();
602 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
604 GetOStream(
Runtime0) <<
"Clearing old data (expert)" << std::endl;
605 for (
int iLevel = 0; iLevel < GetNumLevels(); iLevel++)
606 Levels_[iLevel]->ExpertClear();
609 #if defined(HAVE_MUELU_EXPERIMENTAL) && defined(HAVE_MUELU_ADDITIVE_VARIANT) 610 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
612 bool InitialGuessIsZero, LO startLevel) {
613 LO nIts = conv.maxIts_;
614 MagnitudeType tol = conv.tol_;
616 std::string prefix = this->ShortClassName() +
": ";
617 std::string levelSuffix =
" (level=" +
toString(startLevel) +
")";
618 std::string levelSuffix1 =
" (level=" +
toString(startLevel+1) +
")";
645 SC one = STS::one(), zero = STS::zero();
647 bool zeroGuess = InitialGuessIsZero;
658 bool emptyCoarseSolve =
true;
659 RCP<MultiVector> coarseX_prolonged = MultiVectorFactory::Build(X.getMap(), X.getNumVectors(),
true);
663 if (Levels_.size() > 1) {
665 if (Coarse->IsAvailable(
"Importer"))
675 coarseRhs = MultiVectorFactory::Build(R->getRangeMap(), B.getNumVectors(),
true);
684 if (doPRrebalance_ || importer.is_null()) {
685 coarseX = MultiVectorFactory::Build(coarseRhs->getMap(), X.getNumVectors(),
true);
693 RCP<MultiVector> coarseTmp = MultiVectorFactory::Build(importer->getTargetMap(), coarseRhs->getNumVectors());
694 coarseTmp->doImport(*coarseRhs, *importer, Xpetra::INSERT);
695 coarseRhs.
swap(coarseTmp);
697 coarseX = MultiVectorFactory::Build(importer->getTargetMap(), X.getNumVectors(),
true);
700 if (Coarse->IsAvailable(
"PreSmoother"))
702 if (Coarse->IsAvailable(
"PostSmoother"))
709 MagnitudeType prevNorm = STS::magnitude(STS::one()), curNorm = STS::magnitude(STS::one());
712 for (LO i = 1; i <= nIts; i++) {
713 #ifdef HAVE_MUELU_DEBUG 714 if (A->getDomainMap()->isCompatible(*(X.getMap())) ==
false) {
715 std::ostringstream ss;
716 ss <<
"Level " << startLevel <<
": level A's domain map is not compatible with X";
717 throw Exceptions::Incompatible(ss.str());
720 if (A->getRangeMap()->isCompatible(*(B.getMap())) ==
false) {
721 std::ostringstream ss;
722 ss <<
"Level " << startLevel <<
": level A's range map is not compatible with B";
723 throw Exceptions::Incompatible(ss.str());
728 bool emptyFineSolve =
true;
731 fineX = MultiVectorFactory::Build(X.getMap(), X.getNumVectors(),
true);
740 if (Fine->IsAvailable(
"PreSmoother")) {
743 preSmoo->Apply(*fineX, B, zeroGuess);
745 emptyFineSolve =
false;
747 if (Fine->IsAvailable(
"PostSmoother")) {
750 postSmoo->Apply(*fineX, B, zeroGuess);
753 emptyFineSolve =
false;
755 if (emptyFineSolve ==
true) {
756 GetOStream(
Warnings1) <<
"No fine grid smoother" << std::endl;
758 fineX->update(one, B, zero);
761 if (Levels_.size() > 1) {
763 if (Coarse->IsAvailable(
"PreSmoother")) {
765 preSmoo_coarse->Apply(*coarseX, *coarseRhs, zeroGuess);
767 emptyCoarseSolve =
false;
770 if (Coarse->IsAvailable(
"PostSmoother")) {
772 postSmoo_coarse->Apply(*coarseX, *coarseRhs, zeroGuess);
774 emptyCoarseSolve =
false;
777 if (emptyCoarseSolve ==
true) {
778 GetOStream(
Warnings1) <<
"No coarse grid solver" << std::endl;
780 coarseX->update(one, *coarseRhs, zero);
787 if (!doPRrebalance_ && !importer.is_null()) {
792 RCP<MultiVector> coarseTmp = MultiVectorFactory::Build(importer->getSourceMap(), coarseX->getNumVectors());
793 coarseTmp->doExport(*coarseX, *importer, Xpetra::INSERT);
794 coarseX.
swap(coarseTmp);
803 X.update(1.0, *fineX, 1.0, *coarseX_prolonged, 0.0);
814 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
816 bool InitialGuessIsZero, LO startLevel) {
835 std::string prefix = label + this->ShortClassName() +
": ";
836 std::string levelSuffix =
" (level=" +
toString(startLevel) +
")";
837 std::string levelSuffix1 =
" (level=" +
toString(startLevel+1) +
")";
846 std::string iterateLevelTimeLabel = prefix +
"Solve" + levelSuffix;
849 bool zeroGuess = InitialGuessIsZero;
868 const BlockedMultiVector * Bblocked =
dynamic_cast<const BlockedMultiVector*
>(&B);
869 if(residual_.size() > startLevel &&
870 ( ( Bblocked && !Bblocked->isSameSize(*residual_[startLevel])) ||
871 (!Bblocked && !residual_[startLevel]->isSameSize(B))))
872 DeleteLevelMultiVectors();
873 AllocateLevelMultiVectors(X.getNumVectors());
876 MagnitudeType prevNorm = STS::magnitude(STS::one()), curNorm = STS::magnitude(STS::one());
878 if (startLevel == 0 && !isPreconditioner_ &&
887 for (LO k = 0; k < rn.
size(); k++)
897 << std::setiosflags(std::ios::left)
898 << std::setprecision(3) << 0
900 << std::setprecision(10) << rn
904 SC one = STS::one(), zero = STS::zero();
905 for (LO i = 1; i <= nIts; i++) {
906 #ifdef HAVE_MUELU_DEBUG 908 if (A->getDomainMap()->isCompatible(*(X.getMap())) ==
false) {
909 std::ostringstream ss;
910 ss <<
"Level " << startLevel <<
": level A's domain map is not compatible with X";
914 if (A->getRangeMap()->isCompatible(*(B.getMap())) ==
false) {
915 std::ostringstream ss;
916 ss <<
"Level " << startLevel <<
": level A's range map is not compatible with B";
922 if (startLevel == as<LO>(Levels_.size())-1) {
926 bool emptySolve =
true;
929 if (Fine->IsAvailable(
"PreSmoother")) {
932 preSmoo->Apply(X, B, zeroGuess);
937 if (Fine->IsAvailable(
"PostSmoother")) {
940 postSmoo->Apply(X, B, zeroGuess);
944 if (emptySolve ==
true) {
945 GetOStream(
Warnings1) <<
"No coarse grid solver" << std::endl;
947 X.update(one, B, zero);
958 if (Fine->IsAvailable(
"PreSmoother")) {
960 preSmoo->Apply(X, B, zeroGuess);
969 residual = residual_[startLevel];
973 if (Coarse->IsAvailable(
"Pbar"))
982 coarseRhs = coarseRhs_[startLevel];
984 if (implicitTranspose_) {
994 if (Coarse->IsAvailable(
"Importer"))
997 coarseX = coarseX_[startLevel]; coarseX->putScalar(SC_ZERO);
998 if (!doPRrebalance_ && !importer.
is_null()) {
1004 coarseTmp->doImport(*coarseRhs, *importer, Xpetra::INSERT);
1005 coarseRhs.
swap(coarseTmp);
1014 coarseRhs->replaceMap(Ac->getRangeMap());
1015 coarseX ->replaceMap(Ac->getDomainMap());
1018 iterateLevelTime = Teuchos::null;
1020 Iterate(*coarseRhs, *coarseX, 1,
true, startLevel+1);
1023 Iterate(*coarseRhs, *coarseX, 1,
false, startLevel+1);
1026 iterateLevelTime =
rcp(
new TimeMonitor(*
this, iterateLevelTimeLabel));
1028 coarseX->replaceMap(origXMap);
1029 coarseRhs->replaceMap(origRhsMap);
1032 if (!doPRrebalance_ && !importer.
is_null()) {
1038 coarseTmp->doExport(*coarseX, *importer, Xpetra::INSERT);
1039 coarseX.
swap(coarseTmp);
1053 X.update(scalingFactor_, *correction, one);
1060 if (Fine->IsAvailable(
"PostSmoother")) {
1062 postSmoo->Apply(X, B,
false);
1068 if (startLevel == 0 && !isPreconditioner_ &&
1077 rate_ = as<MagnitudeType>(curNorm / prevNorm);
1081 for (LO k = 0; k < rn.
size(); k++)
1091 << std::setiosflags(std::ios::left)
1092 << std::setprecision(3) << i
1094 << std::setprecision(10) << rn
1103 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1105 LO startLevel = (start != -1 ? start : 0);
1106 LO endLevel = (end != -1 ? end : Levels_.size()-1);
1109 "MueLu::Hierarchy::Write : startLevel must be <= endLevel");
1112 "MueLu::Hierarchy::Write bad start or end level");
1114 for (LO i = startLevel; i < endLevel + 1; i++) {
1115 RCP<Matrix> A = rcp_dynamic_cast<Matrix>(Levels_[i]->
template Get< RCP< Operator> >(
"A")), P, R;
1117 P = rcp_dynamic_cast<Matrix>(Levels_[i]->
template Get< RCP< Operator> >(
"P"));
1118 if (!implicitTranspose_)
1119 R = rcp_dynamic_cast<Matrix>(Levels_[i]->
template Get< RCP< Operator> >(
"R"));
1122 if (!A.
is_null()) Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(
"A_" +
toString(i) + suffix +
".m", *A);
1124 Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(
"P_" +
toString(i) + suffix +
".m", *P);
1127 Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(
"R_" +
toString(i) + suffix +
".m", *R);
1132 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1134 for (
Array<
RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1135 (*it)->Keep(ename, factory);
1138 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1140 for (
Array<
RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1141 (*it)->Delete(ename, factory);
1144 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1146 for (
Array<
RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1147 (*it)->AddKeepFlag(ename, factory, keep);
1150 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1152 for (
Array<
RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1153 (*it)->RemoveKeepFlag(ename, factory, keep);
1156 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1158 if ( description_.empty() )
1160 std::ostringstream out;
1162 out <<
"{#levels = " << GetGlobalNumLevels() <<
", complexity = " << GetOperatorComplexity() <<
"}";
1163 description_ = out.str();
1165 return description_;
1168 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1173 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1175 RCP<Operator> A0 = Levels_[0]->template Get<RCP<Operator> >(
"A");
1178 int numLevels = GetNumLevels();
1179 RCP<Operator> Ac = Levels_[numLevels-1]->template Get<RCP<Operator> >(
"A");
1186 int root = comm->getRank();
1189 int smartData = numLevels*comm->getSize() + comm->getRank(), maxSmartData;
1190 reduceAll(*comm, Teuchos::REDUCE_MAX, smartData, Teuchos::ptr(&maxSmartData));
1191 root = maxSmartData % comm->getSize();
1195 double smoother_comp =-1.0;
1197 smoother_comp = GetSmootherComplexity();
1201 std::vector<Xpetra::global_size_t> nnzPerLevel;
1202 std::vector<Xpetra::global_size_t> rowsPerLevel;
1203 std::vector<int> numProcsPerLevel;
1204 bool aborted =
false;
1205 for (
int i = 0; i < numLevels; i++) {
1207 "Operator A is not available on level " << i);
1209 RCP<Operator> A = Levels_[i]->template Get<RCP<Operator> >(
"A");
1211 "Operator A on level " << i <<
" is null.");
1215 GetOStream(
Warnings0) <<
"Some level operators are not matrices, statistics calculation aborted" << std::endl;
1220 Xpetra::global_size_t nnz = Am->getGlobalNumEntries();
1221 nnzPerLevel .push_back(nnz);
1222 rowsPerLevel .push_back(Am->getGlobalNumRows());
1223 numProcsPerLevel.push_back(Am->getRowMap()->getComm()->getSize());
1227 std::string label = Levels_[0]->getObjectLabel();
1228 std::ostringstream oss;
1229 oss << std::setfill(
' ');
1230 oss <<
"\n--------------------------------------------------------------------------------\n";
1231 oss <<
"--- Multigrid Summary " << std::setw(28) << std::left << label <<
"---\n";
1232 oss <<
"--------------------------------------------------------------------------------" << std::endl;
1233 oss <<
"Number of levels = " << numLevels << std::endl;
1234 oss <<
"Operator complexity = " << std::setprecision(2) << std::setiosflags(std::ios::fixed)
1235 << GetOperatorComplexity() << std::endl;
1237 if(smoother_comp!=-1.0) {
1238 oss <<
"Smoother complexity = " << std::setprecision(2) << std::setiosflags(std::ios::fixed)
1239 << smoother_comp << std::endl;
1244 oss <<
"Cycle type = V" << std::endl;
1247 oss <<
"Cycle type = W" << std::endl;
1254 Xpetra::global_size_t tt = rowsPerLevel[0];
1255 int rowspacer = 2;
while (tt != 0) { tt /= 10; rowspacer++; }
1256 tt = nnzPerLevel[0];
1257 int nnzspacer = 2;
while (tt != 0) { tt /= 10; nnzspacer++; }
1258 tt = numProcsPerLevel[0];
1259 int npspacer = 2;
while (tt != 0) { tt /= 10; npspacer++; }
1260 oss <<
"level " << std::setw(rowspacer) <<
" rows " << std::setw(nnzspacer) <<
" nnz " <<
" nnz/row" << std::setw(npspacer) <<
" c ratio" <<
" procs" << std::endl;
1261 for (
size_t i = 0; i < nnzPerLevel.size(); ++i) {
1262 oss <<
" " << i <<
" ";
1263 oss << std::setw(rowspacer) << rowsPerLevel[i];
1264 oss << std::setw(nnzspacer) << nnzPerLevel[i];
1265 oss << std::setprecision(2) << std::setiosflags(std::ios::fixed);
1266 oss << std::setw(9) << as<double>(nnzPerLevel[i]) / rowsPerLevel[i];
1267 if (i) oss << std::setw(9) << as<double>(rowsPerLevel[i-1])/rowsPerLevel[i];
1268 else oss << std::setw(9) <<
" ";
1269 oss <<
" " << std::setw(npspacer) << numProcsPerLevel[i] << std::endl;
1272 for (
int i = 0; i < GetNumLevels(); ++i) {
1274 if (Levels_[i]->IsAvailable(
"PreSmoother"))
1276 if (Levels_[i]->IsAvailable(
"PostSmoother"))
1279 if (preSmoo != null && preSmoo == postSmoo)
1280 oss <<
"Smoother (level " << i <<
") both : " << preSmoo->
description() << std::endl;
1282 oss <<
"Smoother (level " << i <<
") pre : " 1283 << (preSmoo != null ? preSmoo->
description() :
"no smoother") << std::endl;
1284 oss <<
"Smoother (level " << i <<
") post : " 1285 << (postSmoo != null ? postSmoo->
description() :
"no smoother") << std::endl;
1297 MPI_Comm rawComm = (*mpiComm->getRawMpiComm())();
1299 int strLength = outstr.size();
1300 MPI_Bcast(&strLength, 1, MPI_INT, root, rawComm);
1301 if (comm->getRank() != root)
1302 outstr.resize(strLength);
1303 MPI_Bcast(&outstr[0], strLength, MPI_CHAR, root, rawComm);
1310 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1313 for (
int i = 0; i < GetNumLevels(); ++i)
1314 Levels_[i]->print(out, verbLevel);
1317 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1319 isPreconditioner_ = flag;
1322 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1324 if (GetProcRankVerbose() != 0)
1326 #if defined(HAVE_MUELU_BOOST) && defined(HAVE_MUELU_BOOST_FOR_REAL) && defined(BOOST_VERSION) && (BOOST_VERSION >= 104400) 1330 dp.property(
"label", boost::get(boost::vertex_name, graph));
1331 dp.property(
"id", boost::get(boost::vertex_index, graph));
1332 dp.property(
"label", boost::get(boost::edge_name, graph));
1333 dp.property(
"color", boost::get(boost::edge_color, graph));
1336 std::map<const FactoryBase*, BoostVertex> vindices;
1337 typedef std::map<std::pair<BoostVertex,BoostVertex>, std::string> emap; emap edges;
1339 static int call_id=0;
1341 RCP<Operator> A = Levels_[0]->template Get<RCP<Operator> >(
"A");
1342 int rank = A->getDomainMap()->getComm()->getRank();
1345 for (
int i = dumpLevel_; i <= dumpLevel_+1 && i < GetNumLevels(); i++) {
1347 Levels_[i]->UpdateGraph(vindices, edges, dp, graph);
1349 for (emap::const_iterator eit = edges.begin(); eit != edges.end(); eit++) {
1350 std::pair<BoostEdge, bool> boost_edge = boost::add_edge(eit->first.first, eit->first.second, graph);
1353 if(eit->second==std::string(
"Graph")) boost::put(
"label", dp, boost_edge.first, std::string(
"Graph_"));
1354 else boost::put(
"label", dp, boost_edge.first, eit->second);
1355 if (i == dumpLevel_)
1356 boost::put(
"color", dp, boost_edge.first, std::string(
"red"));
1358 boost::put(
"color", dp, boost_edge.first, std::string(
"blue"));
1363 std::ostringstream legend;
1364 legend <<
"< <TABLE BORDER=\"0\" CELLBORDER=\"1\" CELLSPACING=\"0\" CELLPADDING=\"4\"> \ 1365 <TR><TD COLSPAN=\"2\">Legend</TD></TR> \ 1366 <TR><TD><FONT color=\"red\">Level " << dumpLevel_ <<
"</FONT></TD><TD><FONT color=\"blue\">Level " << dumpLevel_+1 <<
"</FONT></TD></TR> \ 1368 BoostVertex boost_vertex = boost::add_vertex(graph);
1369 boost::put(
"label", dp, boost_vertex, legend.str());
1372 std::ofstream out(dumpFile_.c_str() +std::to_string(call_id)+std::string(
"_")+ std::to_string(rank) + std::string(
".dot"));
1373 boost::write_graphviz_dp(out, graph, dp, std::string(
"id"));
1377 GetOStream(
Errors) <<
"Dependency graph output requires boost and MueLu_ENABLE_Boost_for_real" << std::endl;
1382 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1387 GetOStream(
Warnings1) <<
"Hierarchy::ReplaceCoordinateMap: operator is not a matrix, skipping..." << std::endl;
1390 if(Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A) != Teuchos::null) {
1391 GetOStream(
Warnings1) <<
"Hierarchy::ReplaceCoordinateMap: operator is a BlockedCrsMatrix, skipping..." << std::endl;
1395 typedef Xpetra::MultiVector<double,LO,GO,NO> xdMV;
1399 if (A->getRowMap()->isSameAs(*(coords->getMap()))) {
1400 GetOStream(
Warnings1) <<
"Hierarchy::ReplaceCoordinateMap: matrix and coordinates maps are same, skipping..." << std::endl;
1404 if (A->IsView(
"stridedMaps") && rcp_dynamic_cast<
const StridedMap>(A->getRowMap(
"stridedMaps")) != Teuchos::null) {
1405 RCP<const StridedMap> stridedRowMap = rcp_dynamic_cast<
const StridedMap>(A->getRowMap(
"stridedMaps"));
1409 Exceptions::RuntimeError,
"Hierarchy::ReplaceCoordinateMap: nontrivial maps (block id = " << stridedRowMap->getStridedBlockId()
1410 <<
", offset = " << stridedRowMap->getOffset() <<
")");
1413 GetOStream(
Runtime1) <<
"Replacing coordinate map" << std::endl;
1415 size_t blkSize = A->GetFixedBlockSize();
1421 GO indexBase = dofMap->getIndexBase();
1422 size_t numLocalDOFs = dofMap->getNodeNumElements();
1424 "Hierarchy::ReplaceCoordinateMap: block size (" << blkSize <<
") is incompatible with the number of local dofs in a row map (" << numLocalDOFs);
1427 Array<GO> nodeGIDs(numLocalDOFs/blkSize);
1428 for (
size_t i = 0; i < numLocalDOFs; i += blkSize)
1429 nodeGIDs[i/blkSize] = (GIDs[i] - indexBase)/blkSize + indexBase;
1432 nodeMap = MapFactory::Build(dofMap->lib(), INVALID, nodeGIDs(), indexBase, dofMap->getComm());
1438 if(coords->getLocalLength() != A->getRowMap()->getNodeNumElements()) {
1439 GetOStream(
Warnings) <<
"Coordinate vector does not match row map of matrix A!" << std::endl;
1445 std::vector<ArrayRCP<const double> > coordData;
1446 for (
size_t i = 0; i < coords->getNumVectors(); i++) {
1447 coordData.
push_back(coords->getData(i));
1448 coordDataView.
push_back(coordData[i]());
1451 RCP<xdMV> newCoords = Xpetra::MultiVectorFactory<double,LO,GO,NO>::Build(nodeMap, coordDataView(), coords->getNumVectors());
1452 level.
Set(
"Coordinates", newCoords);
1455 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1457 int N = Levels_.size();
1458 if( (sizeOfAllocatedLevelMultiVectors_ == numvecs && residual_.size() == N) || numvecs<=0 )
return;
1461 if(residual_.size() != N) {
1462 DeleteLevelMultiVectors();
1464 residual_.resize(N);
1465 coarseRhs_.resize(N);
1467 coarseImport_.resize(N);
1468 coarseExport_.resize(N);
1469 correction_.resize(N);
1472 for(
int i=0; i<N; i++) {
1473 RCP<Operator> A = Levels_[i]->template Get< RCP<Operator> >(
"A");
1480 Adm = A_as_blocked->getFullDomainMap();
1484 residual_[i] = MultiVectorFactory::Build(Arm, numvecs,
true);
1485 correction_[i] = MultiVectorFactory::Build(Adm, numvecs,
false);
1490 if(implicitTranspose_) {
1491 RCP<Operator> P = Levels_[i+1]->template Get< RCP<Operator> >(
"P");
1492 if(!P.
is_null()) coarseRhs_[i] = MultiVectorFactory::Build(P->getDomainMap(),numvecs,
true);
1494 RCP<Operator> R = Levels_[i+1]->template Get< RCP<Operator> >(
"R");
1495 if(!R.
is_null()) coarseRhs_[i] = MultiVectorFactory::Build(R->getRangeMap(),numvecs,
true);
1500 if(Levels_[i+1]->IsAvailable(
"Importer"))
1502 if (doPRrebalance_ || importer.
is_null())
1503 coarseX_[i] = MultiVectorFactory::Build(coarseRhs_[i]->getMap(),numvecs,
false);
1505 coarseImport_[i] = MultiVectorFactory::Build(importer->getTargetMap(), numvecs,
false);
1506 coarseExport_[i] = MultiVectorFactory::Build(importer->getSourceMap(), numvecs,
false);
1507 coarseX_[i] = MultiVectorFactory::Build(importer->getTargetMap(),numvecs,
false);
1511 sizeOfAllocatedLevelMultiVectors_ = numvecs;
1515 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1517 if(sizeOfAllocatedLevelMultiVectors_==0)
return;
1518 residual_.resize(0);
1519 coarseRhs_.resize(0);
1521 coarseImport_.resize(0);
1522 coarseExport_.resize(0);
1523 correction_.resize(0);
1524 sizeOfAllocatedLevelMultiVectors_ = 0;
1530 #endif // MUELU_HIERARCHY_DEF_HPP Important warning messages (one line)
void IsPreconditioner(const bool flag)
RCP< Level > & GetPreviousLevel()
Previous level.
static void SetRandomSeed(const Teuchos::Comm< int > &comm)
High level timing information (use Teuchos::TimeMonitor::summarize() to print)
double GetSmootherComplexity() 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.
RCP< Level > & GetLevel(const int levelID=0)
Retrieve a certain level from hierarchy.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
virtual std::string getObjectLabel() const
Hierarchy()
Default constructor.
virtual size_t getNodeSmootherComplexity() const =0
Compute a rough estimate of the cost to apply this smoother on this MPI rank. Return Teuchos::Ordinal...
void DeleteLevelMultiVectors()
void Release(const FactoryBase &factory)
Decrement the storage counter for all the inputs of a factory.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
void AddLevel(const RCP< Level > &level)
Add a level at the end of the hierarchy.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void AddNewLevel()
Add a new level at the end of the hierarchy.
VerbLevel toMueLuVerbLevel(const Teuchos::EVerbosityLevel verbLevel)
void swap(RCP< T > &r_ptr)
One-liner description of what is happening.
void Clear(int startLevel=0)
Clear impermanent data from previous setup.
Namespace for MueLu classes and methods.
void SetPreviousLevel(const RCP< Level > &previousLevel)
Exception throws to report incompatible objects (like maps).
Integrates Teuchos::TimeMonitor with MueLu verbosity system.
Print skeleton for the run, i.e. factory calls and used parameters.
void ReplaceCoordinateMap(Level &level)
void AllocateLevelMultiVectors(int numvecs)
void AddKeepFlag(const std::string &ename, const FactoryBase *factory=NoFactory::get(), KeepType keep=MueLu::Keep)
Call Level::AddKeepFlag for each level of the Hierarchy.
Print even more statistics.
void setlib(Xpetra::UnderlyingLib lib2)
int GetLevelID() const
Return level number.
void print(std::ostream &out=std::cout, const VerbLevel verbLevel=(MueLu::Parameters|MueLu::Statistics0)) const
Base class for factories (e.g., R, P, and A_coarse).
Print statistics that do not involve significant additional computation.
Xpetra::UnderlyingLib lib_
static RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &RHS)
void start(bool reset=false)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void Write(const LO &start=-1, const LO &end=-1, const std::string &suffix="")
Print matrices in the multigrid hierarchy to file.
Class that holds all level-specific information.
Class that provides default factories within Needs class.
void Build(Level &level) const
Build an object with this factory.
RCP< const Teuchos::Comm< int > > GetComm() const
void RemoveKeepFlag(const std::string &ename, const FactoryBase *factory, KeepType keep=MueLu::All)
Call Level::RemoveKeepFlag for each level of the Hierarchy.
int GetGlobalNumLevels() const
void CheckLevel(Level &level, int levelID)
Helper function.
Xpetra::UnderlyingLib lib()
virtual void setObjectLabel(const std::string &objectLabel)
TEUCHOS_DEPRECATED void reduceAll(const Comm< Ordinal > &comm, const EReductionType reductType, const Packet &send, Packet *globalReduct)
void describe(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the Hierarchy with some verbosity level to a FancyOStream object.
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
double GetOperatorComplexity() const
void push_back(const value_type &x)
Data struct for defining stopping criteria of multigrid iteration.
static Teuchos::Array< Magnitude > ResidualNorm(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &RHS)
void SetLevelID(int levelID)
Set level number.
STS::magnitudeType MagnitudeType
void Delete(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Call Level::Delete(ename, factory) for each level of the Hierarchy.
Timer to be used in non-factories.
Classes::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > MultiVector
Array< RCP< Level > > Levels_
Container for Level objects.
void DumpCurrentGraph() const
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
Exception throws to report errors in the internal logical of the program.
An exception safe way to call the method 'Level::SetFactoryManager()'.
void SetComm(RCP< const Teuchos::Comm< int > > const &comm)
Print all warning messages.
ReturnType Iterate(const MultiVector &B, MultiVector &X, ConvData conv=ConvData(), bool InitialGuessIsZero=false, LO startLevel=0)
Apply the multigrid preconditioner.
Description of what is happening (more verbose)
void Keep(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Call Level::Keep(ename, factory) for each level of the Hierarchy.
std::string description() const
Return a simple one-line description of this object.
bool Setup(int coarseLevelID, const RCP< const FactoryManagerBase > fineLevelManager, const RCP< const FactoryManagerBase > coarseLevelManager, const RCP< const FactoryManagerBase > nextLevelManager=Teuchos::null)
Multi-level setup phase: build a new level of the hierarchy.
virtual void Clean() const
virtual std::string description() const
Return a simple one-line description of this object.
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
std::string toString(const T &t)
void Request(const FactoryBase &factory)
Increment the storage counter for all the inputs of a factory.