47 #ifndef MUELU_HIERARCHY_DEF_HPP
48 #define MUELU_HIERARCHY_DEF_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();
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.");
536 TEUCHOS_TEST_FOR_EXCEPTION(!Levels_[startLevel]->IsAvailable(
"A"), Exceptions::RuntimeError,
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);
583 TEUCHOS_TEST_FOR_EXCEPTION(iLevel != Levels_.size() - 1, Exceptions::RuntimeError,
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) +
")";
620 using namespace Teuchos;
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());
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());
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) +
")";
845 else if (!useStackedTimer)
848 std::string iterateLevelTimeLabel = prefix +
"Solve" + levelSuffix;
851 bool zeroGuess = InitialGuessIsZero;
854 using namespace Teuchos;
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);
956 if (!useStackedTimer)
960 if (Fine->IsAvailable(
"PreSmoother")) {
962 preSmoo->Apply(X, B, zeroGuess);
969 if (!useStackedTimer)
973 residual = residual_[startLevel];
977 if (Coarse->IsAvailable(
"Pbar"))
985 if (!useStackedTimer)
988 coarseRhs = coarseRhs_[startLevel];
990 if (implicitTranspose_) {
1000 if (Coarse->IsAvailable(
"Importer"))
1003 coarseX = coarseX_[startLevel];
1004 if (!doPRrebalance_ && !importer.
is_null()) {
1006 if (!useStackedTimer)
1013 coarseRhs.
swap(coarseTmp);
1022 coarseRhs->replaceMap(Ac->getRangeMap());
1023 coarseX ->replaceMap(Ac->getDomainMap());
1026 iterateLevelTime = Teuchos::null;
1028 Iterate(*coarseRhs, *coarseX, 1,
true, startLevel+1);
1031 Iterate(*coarseRhs, *coarseX, 1,
false, startLevel+1);
1034 iterateLevelTime =
rcp(
new TimeMonitor(*
this, iterateLevelTimeLabel));
1036 coarseX->replaceMap(origXMap);
1037 coarseRhs->replaceMap(origRhsMap);
1040 if (!doPRrebalance_ && !importer.
is_null()) {
1042 if (!useStackedTimer)
1049 coarseX.
swap(coarseTmp);
1060 if (!useStackedTimer)
1065 X.update(scalingFactor_, *correction, one);
1070 if (!useStackedTimer)
1074 if (Fine->IsAvailable(
"PostSmoother")) {
1076 postSmoo->Apply(X, B,
false);
1082 if (startLevel == 0 && !isPreconditioner_ &&
1091 rate_ = as<MagnitudeType>(curNorm / prevNorm);
1095 << std::setiosflags(std::ios::left)
1096 << std::setprecision(3) << i
1098 << std::setprecision(10) << rn
1103 for (
LO k = 0; k < rn.
size(); k++)
1117 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1119 LO startLevel = (start != -1 ? start : 0);
1120 LO endLevel = (end != -1 ? end : Levels_.size()-1);
1123 "MueLu::Hierarchy::Write : startLevel must be <= endLevel");
1126 "MueLu::Hierarchy::Write bad start or end level");
1128 for (
LO i = startLevel; i < endLevel + 1; i++) {
1129 RCP<Matrix> A = rcp_dynamic_cast<Matrix>(Levels_[i]->
template Get< RCP< Operator> >(
"A")), P, R;
1131 P = rcp_dynamic_cast<Matrix>(Levels_[i]->
template Get< RCP< Operator> >(
"P"));
1132 if (!implicitTranspose_)
1133 R = rcp_dynamic_cast<Matrix>(Levels_[i]->
template Get< RCP< Operator> >(
"R"));
1146 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1148 for (
Array<
RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1149 (*it)->Keep(ename, factory);
1152 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1154 for (
Array<
RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1155 (*it)->Delete(ename, factory);
1158 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1160 for (
Array<
RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1161 (*it)->AddKeepFlag(ename, factory, keep);
1164 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1166 for (
Array<
RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1167 (*it)->RemoveKeepFlag(ename, factory, keep);
1170 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1172 if ( description_.empty() )
1174 std::ostringstream out;
1176 out <<
"{#levels = " << GetGlobalNumLevels() <<
", complexity = " << GetOperatorComplexity() <<
"}";
1177 description_ = out.str();
1179 return description_;
1182 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1187 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1189 RCP<Operator> A0 = Levels_[0]->template Get<RCP<Operator> >(
"A");
1192 int numLevels = GetNumLevels();
1193 RCP<Operator> Ac = Levels_[numLevels-1]->template Get<RCP<Operator> >(
"A");
1200 int root = comm->getRank();
1203 int smartData = numLevels*comm->getSize() + comm->getRank(), maxSmartData;
1205 root = maxSmartData % comm->getSize();
1209 double smoother_comp =-1.0;
1211 smoother_comp = GetSmootherComplexity();
1215 std::vector<Xpetra::global_size_t> nnzPerLevel;
1216 std::vector<Xpetra::global_size_t> rowsPerLevel;
1217 std::vector<int> numProcsPerLevel;
1218 bool aborted =
false;
1219 for (
int i = 0; i < numLevels; i++) {
1221 "Operator A is not available on level " << i);
1223 RCP<Operator> A = Levels_[i]->template Get<RCP<Operator> >(
"A");
1225 "Operator A on level " << i <<
" is null.");
1229 GetOStream(
Warnings0) <<
"Some level operators are not matrices, statistics calculation aborted" << std::endl;
1235 nnzPerLevel .push_back(nnz);
1236 rowsPerLevel .push_back(Am->getGlobalNumRows());
1237 numProcsPerLevel.push_back(Am->getRowMap()->getComm()->getSize());
1241 std::string label = Levels_[0]->getObjectLabel();
1242 std::ostringstream oss;
1243 oss << std::setfill(
' ');
1244 oss <<
"\n--------------------------------------------------------------------------------\n";
1245 oss <<
"--- Multigrid Summary " << std::setw(28) << std::left << label <<
"---\n";
1246 oss <<
"--------------------------------------------------------------------------------" << std::endl;
1247 oss <<
"Number of levels = " << numLevels << std::endl;
1248 oss <<
"Operator complexity = " << std::setprecision(2) << std::setiosflags(std::ios::fixed)
1249 << GetOperatorComplexity() << std::endl;
1251 if(smoother_comp!=-1.0) {
1252 oss <<
"Smoother complexity = " << std::setprecision(2) << std::setiosflags(std::ios::fixed)
1253 << smoother_comp << std::endl;
1258 oss <<
"Cycle type = V" << std::endl;
1261 oss <<
"Cycle type = W" << std::endl;
1269 int rowspacer = 2;
while (tt != 0) { tt /= 10; rowspacer++; }
1270 tt = nnzPerLevel[0];
1271 int nnzspacer = 2;
while (tt != 0) { tt /= 10; nnzspacer++; }
1272 tt = numProcsPerLevel[0];
1273 int npspacer = 2;
while (tt != 0) { tt /= 10; npspacer++; }
1274 oss <<
"level " << std::setw(rowspacer) <<
" rows " << std::setw(nnzspacer) <<
" nnz " <<
" nnz/row" << std::setw(npspacer) <<
" c ratio" <<
" procs" << std::endl;
1275 for (
size_t i = 0; i < nnzPerLevel.size(); ++i) {
1276 oss <<
" " << i <<
" ";
1277 oss << std::setw(rowspacer) << rowsPerLevel[i];
1278 oss << std::setw(nnzspacer) << nnzPerLevel[i];
1279 oss << std::setprecision(2) << std::setiosflags(std::ios::fixed);
1280 oss << std::setw(9) << as<double>(nnzPerLevel[i]) / rowsPerLevel[i];
1281 if (i) oss << std::setw(9) << as<double>(rowsPerLevel[i-1])/rowsPerLevel[i];
1282 else oss << std::setw(9) <<
" ";
1283 oss <<
" " << std::setw(npspacer) << numProcsPerLevel[i] << std::endl;
1286 for (
int i = 0; i < GetNumLevels(); ++i) {
1288 if (Levels_[i]->IsAvailable(
"PreSmoother"))
1290 if (Levels_[i]->IsAvailable(
"PostSmoother"))
1293 if (preSmoo != null && preSmoo == postSmoo)
1294 oss <<
"Smoother (level " << i <<
") both : " << preSmoo->
description() << std::endl;
1296 oss <<
"Smoother (level " << i <<
") pre : "
1297 << (preSmoo != null ? preSmoo->
description() :
"no smoother") << std::endl;
1298 oss <<
"Smoother (level " << i <<
") post : "
1299 << (postSmoo != null ? postSmoo->
description() :
"no smoother") << std::endl;
1311 MPI_Comm rawComm = (*mpiComm->getRawMpiComm())();
1313 int strLength = outstr.size();
1314 MPI_Bcast(&strLength, 1, MPI_INT, root, rawComm);
1315 if (comm->getRank() != root)
1316 outstr.resize(strLength);
1317 MPI_Bcast(&outstr[0], strLength, MPI_CHAR, root, rawComm);
1324 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1327 for (
int i = 0; i < GetNumLevels(); ++i)
1328 Levels_[i]->print(out, verbLevel);
1331 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1333 isPreconditioner_ = flag;
1336 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1338 if (GetProcRankVerbose() != 0)
1340 #if defined(HAVE_MUELU_BOOST) && defined(HAVE_MUELU_BOOST_FOR_REAL) && defined(BOOST_VERSION) && (BOOST_VERSION >= 104400)
1344 dp.property(
"label", boost::get(boost::vertex_name, graph));
1345 dp.property(
"id", boost::get(boost::vertex_index, graph));
1346 dp.property(
"label", boost::get(boost::edge_name, graph));
1347 dp.property(
"color", boost::get(boost::edge_color, graph));
1350 std::map<const FactoryBase*, BoostVertex> vindices;
1351 typedef std::map<std::pair<BoostVertex,BoostVertex>, std::string> emap; emap edges;
1353 static int call_id=0;
1355 RCP<Operator> A = Levels_[0]->template Get<RCP<Operator> >(
"A");
1356 int rank = A->getDomainMap()->getComm()->getRank();
1359 for (
int i = dumpLevel_; i <= dumpLevel_+1 && i < GetNumLevels(); i++) {
1361 Levels_[i]->UpdateGraph(vindices, edges, dp, graph);
1363 for (emap::const_iterator eit = edges.begin(); eit != edges.end(); eit++) {
1364 std::pair<BoostEdge, bool> boost_edge = boost::add_edge(eit->first.first, eit->first.second, graph);
1367 if(eit->second==std::string(
"Graph")) boost::put(
"label", dp, boost_edge.first, std::string(
"Graph_"));
1368 else boost::put(
"label", dp, boost_edge.first, eit->second);
1369 if (i == dumpLevel_)
1370 boost::put(
"color", dp, boost_edge.first, std::string(
"red"));
1372 boost::put(
"color", dp, boost_edge.first, std::string(
"blue"));
1377 std::ostringstream legend;
1378 legend <<
"< <TABLE BORDER=\"0\" CELLBORDER=\"1\" CELLSPACING=\"0\" CELLPADDING=\"4\"> \
1379 <TR><TD COLSPAN=\"2\">Legend</TD></TR> \
1380 <TR><TD><FONT color=\"red\">Level " << dumpLevel_ <<
"</FONT></TD><TD><FONT color=\"blue\">Level " << dumpLevel_+1 <<
"</FONT></TD></TR> \
1382 BoostVertex boost_vertex = boost::add_vertex(graph);
1383 boost::put(
"label", dp, boost_vertex, legend.str());
1386 std::ofstream out(dumpFile_.c_str() +std::to_string(call_id)+std::string(
"_")+ std::to_string(rank) + std::string(
".dot"));
1387 boost::write_graphviz_dp(out, graph, dp, std::string(
"id"));
1391 GetOStream(
Errors) <<
"Dependency graph output requires boost and MueLu_ENABLE_Boost_for_real" << std::endl;
1396 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1401 GetOStream(
Warnings1) <<
"Hierarchy::ReplaceCoordinateMap: operator is not a matrix, skipping..." << std::endl;
1404 if(Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A) != Teuchos::null) {
1405 GetOStream(
Warnings1) <<
"Hierarchy::ReplaceCoordinateMap: operator is a BlockedCrsMatrix, skipping..." << std::endl;
1413 if (A->getRowMap()->isSameAs(*(coords->getMap()))) {
1414 GetOStream(
Warnings1) <<
"Hierarchy::ReplaceCoordinateMap: matrix and coordinates maps are same, skipping..." << std::endl;
1418 if (A->IsView(
"stridedMaps") && rcp_dynamic_cast<
const StridedMap>(A->getRowMap(
"stridedMaps")) != Teuchos::null) {
1419 RCP<const StridedMap> stridedRowMap = rcp_dynamic_cast<
const StridedMap>(A->getRowMap(
"stridedMaps"));
1423 Exceptions::RuntimeError,
"Hierarchy::ReplaceCoordinateMap: nontrivial maps (block id = " << stridedRowMap->getStridedBlockId()
1424 <<
", offset = " << stridedRowMap->getOffset() <<
")");
1427 GetOStream(
Runtime1) <<
"Replacing coordinate map" << std::endl;
1429 size_t blkSize = A->GetFixedBlockSize();
1435 GO indexBase = dofMap->getIndexBase();
1436 size_t numLocalDOFs = dofMap->getNodeNumElements();
1438 "Hierarchy::ReplaceCoordinateMap: block size (" << blkSize <<
") is incompatible with the number of local dofs in a row map (" << numLocalDOFs);
1441 Array<GO> nodeGIDs(numLocalDOFs/blkSize);
1442 for (
size_t i = 0; i < numLocalDOFs; i += blkSize)
1443 nodeGIDs[i/blkSize] = (GIDs[i] - indexBase)/blkSize + indexBase;
1446 nodeMap = MapFactory::Build(dofMap->lib(), INVALID, nodeGIDs(), indexBase, dofMap->getComm());
1452 if(coords->getLocalLength() != A->getRowMap()->getNodeNumElements()) {
1453 GetOStream(
Warnings) <<
"Coordinate vector does not match row map of matrix A!" << std::endl;
1459 std::vector<ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> > coordData;
1460 for (
size_t i = 0; i < coords->getNumVectors(); i++) {
1461 coordData.
push_back(coords->getData(i));
1462 coordDataView.
push_back(coordData[i]());
1466 level.
Set(
"Coordinates", newCoords);
1469 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1471 int N = Levels_.size();
1472 if( (sizeOfAllocatedLevelMultiVectors_ == numvecs && residual_.size() == N) || numvecs<=0 )
return;
1475 if(residual_.size() != N) {
1476 DeleteLevelMultiVectors();
1478 residual_.resize(N);
1479 coarseRhs_.resize(N);
1481 coarseImport_.resize(N);
1482 coarseExport_.resize(N);
1483 correction_.resize(N);
1486 for(
int i=0; i<N; i++) {
1487 RCP<Operator> A = Levels_[i]->template Get< RCP<Operator> >(
"A");
1494 Adm = A_as_blocked->getFullDomainMap();
1498 residual_[i] = MultiVectorFactory::Build(Arm, numvecs,
true);
1499 correction_[i] = MultiVectorFactory::Build(Adm, numvecs,
false);
1504 if(implicitTranspose_) {
1505 RCP<Operator> P = Levels_[i+1]->template Get< RCP<Operator> >(
"P");
1506 if(!P.
is_null()) coarseRhs_[i] = MultiVectorFactory::Build(P->getDomainMap(),numvecs,
true);
1508 RCP<Operator> R = Levels_[i+1]->template Get< RCP<Operator> >(
"R");
1509 if(!R.
is_null()) coarseRhs_[i] = MultiVectorFactory::Build(R->getRangeMap(),numvecs,
true);
1514 if(Levels_[i+1]->IsAvailable(
"Importer"))
1516 if (doPRrebalance_ || importer.
is_null())
1517 coarseX_[i] = MultiVectorFactory::Build(coarseRhs_[i]->getMap(),numvecs,
false);
1519 coarseImport_[i] = MultiVectorFactory::Build(importer->getTargetMap(), numvecs,
false);
1520 coarseExport_[i] = MultiVectorFactory::Build(importer->getSourceMap(), numvecs,
false);
1521 coarseX_[i] = MultiVectorFactory::Build(importer->getTargetMap(),numvecs,
false);
1525 sizeOfAllocatedLevelMultiVectors_ = numvecs;
1529 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1531 if(sizeOfAllocatedLevelMultiVectors_==0)
return;
1532 residual_.resize(0);
1533 coarseRhs_.resize(0);
1535 coarseImport_.resize(0);
1536 coarseExport_.resize(0);
1537 correction_.resize(0);
1538 sizeOfAllocatedLevelMultiVectors_ = 0;
1544 #endif // MUELU_HIERARCHY_DEF_HPP
Important warning messages (one line)
void Build(Level &level) const
Build an object with this factory.
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)
virtual std::string getObjectLabel() 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.
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 AddLevel(const RCP< Level > &level)
Add a level at the end of the hierarchy.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
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.
virtual void Clean() const
void DumpCurrentGraph() const
void Clear(int startLevel=0)
Clear impermanent data from previous setup.
void SetPreviousLevel(const RCP< Level > &previousLevel)
Exception throws to report incompatible objects (like maps).
Integrates Teuchos::TimeMonitor with MueLu verbosity system.
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)
void print(std::ostream &out=std::cout, const VerbLevel verbLevel=(MueLu::Parameters|MueLu::Statistics0)) const
double GetSmootherComplexity() 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)
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.
std::string description() const
Return a simple one-line description of this object.
Class that provides default factories within Needs class.
void RemoveKeepFlag(const std::string &ename, const FactoryBase *factory, KeepType keep=MueLu::All)
Call Level::RemoveKeepFlag for each level of the Hierarchy.
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)
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())
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
int GetGlobalNumLevels() 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.
static const Teuchos::RCP< Teuchos::StackedTimer > & getStackedTimer()
Array< RCP< Level > > Levels_
Container for Level objects.
int GetLevelID() const
Return level number.
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)
void describe(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the Hierarchy with some verbosity level to a FancyOStream object.
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.
double GetOperatorComplexity() const
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.
RCP< const Teuchos::Comm< int > > GetComm() 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)
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
void Request(const FactoryBase &factory)
Increment the storage counter for all the inputs of a factory.