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_PFactory.hpp" 71 #include "MueLu_SmootherFactory.hpp" 74 #include "Teuchos_TimeMonitor.hpp" 78 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
80 : maxCoarseSize_(GetDefaultMaxCoarseSize()), implicitTranspose_(GetDefaultImplicitTranspose()),
81 doPRrebalance_(GetDefaultPRrebalance()), isPreconditioner_(true), Cycle_(GetDefaultCycle()),
82 lib_(Xpetra::UseTpetra), isDumpingEnabled_(false), dumpLevel_(-1), rate_(-1)
87 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
89 : maxCoarseSize_(GetDefaultMaxCoarseSize()), implicitTranspose_(GetDefaultImplicitTranspose()),
90 doPRrebalance_(GetDefaultPRrebalance()), isPreconditioner_(true), Cycle_(GetDefaultCycle()),
91 isDumpingEnabled_(false), dumpLevel_(-1), rate_(-1)
93 lib_ = A->getDomainMap()->lib();
95 RCP<Level> Finest = rcp(
new Level);
101 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
103 int levelID = LastLevelID() + 1;
105 if (level->GetLevelID() != -1 && (level->GetLevelID() != levelID))
106 GetOStream(
Warnings1) <<
"Hierarchy::AddLevel(): Level with ID=" << level->GetLevelID() <<
107 " have been added at the end of the hierarchy\n but its ID have been redefined" <<
108 " because last level ID of the hierarchy was " << LastLevelID() <<
"." << std::endl;
110 Levels_.push_back(level);
111 level->SetLevelID(levelID);
114 level->SetPreviousLevel( (levelID == 0) ? Teuchos::null : Levels_[LastLevelID() - 1] );
117 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
119 RCP<Level> newLevel = Levels_[LastLevelID()]->Build();
120 newLevel->setlib(lib_);
121 this->AddLevel(newLevel);
124 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
127 "MueLu::Hierarchy::GetLevel(): invalid input parameter value: LevelID = " << levelID);
128 return Levels_[levelID];
131 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
133 return Levels_.size();
136 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
138 RCP<Operator> A = Levels_[0]->template Get<RCP<Operator> >(
"A");
139 RCP<const Teuchos::Comm<int> > comm = A->getDomainMap()->getComm();
141 int numLevels = GetNumLevels();
143 Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, numLevels, Teuchos::ptr(&numGlobalLevels));
145 return numGlobalLevels;
148 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
150 double totalNnz = 0, lev0Nnz = 1;
151 for (
int i = 0; i < GetNumLevels(); ++i) {
153 "Operator complexity cannot be calculated because A is unavailable on level " << i);
154 RCP<Operator> A = Levels_[i]->template Get<RCP<Operator> >(
"A");
158 RCP<Matrix> Am = rcp_dynamic_cast<Matrix>(A);
160 GetOStream(
Warnings0) <<
"Some level operators are not matrices, operator complexity calculation aborted" << std::endl;
164 totalNnz += as<double>(Am->getGlobalNumEntries());
168 return totalNnz / lev0Nnz;
172 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
175 "MueLu::Hierarchy::CheckLevel(): wrong underlying linear algebra library.");
177 "MueLu::Hierarchy::CheckLevel(): wrong level ID");
179 "MueLu::Hierarchy::Setup(): wrong level parent");
184 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
186 const RCP<const FactoryManagerBase> fineLevelManager,
187 const RCP<const FactoryManagerBase> coarseLevelManager,
188 const RCP<const FactoryManagerBase> nextLevelManager) {
191 TimeMonitor m1(*
this, this->ShortClassName() +
": " +
"Setup (total)");
196 "MueLu::Hierarchy::Setup(): argument coarseLevelManager cannot be null");
202 "MueLu::Hierarchy:Setup(): level " << coarseLevelID <<
" (specified by coarseLevelID argument) " 203 "must be built before calling this function.");
205 Level& level = *Levels_[coarseLevelID];
207 if (levelManagers_.size() < coarseLevelID+1)
208 levelManagers_.resize(coarseLevelID+1);
209 levelManagers_[coarseLevelID] = coarseLevelManager;
211 bool isFinestLevel = (fineLevelManager.is_null());
212 bool isLastLevel = (nextLevelManager.is_null());
216 RCP<Operator> A = level.
Get< RCP<Operator> >(
"A");
217 RCP<const Map> domainMap = A->getDomainMap();
218 RCP<const Teuchos::Comm<int> > comm = domainMap->getComm();
225 oldRank = SetProcRankVerbose(comm->getRank());
229 lib_ = domainMap->lib();
236 Level& prevLevel = *Levels_[coarseLevelID-1];
237 oldRank = SetProcRankVerbose(prevLevel.
GetComm()->getRank());
240 CheckLevel(level, coarseLevelID);
243 RCP<SetFactoryManager> SFMFine;
245 SFMFine = rcp(
new SetFactoryManager(Levels_[coarseLevelID-1], fineLevelManager));
247 if (isFinestLevel && Levels_[coarseLevelID]->IsAvailable(
"Coordinates"))
248 ReplaceCoordinateMap(*Levels_[coarseLevelID]);
253 if (isDumpingEnabled_ && dumpLevel_ == 0 && coarseLevelID == 1)
256 RCP<TopSmootherFactory> coarseFact = rcp(
new TopSmootherFactory(coarseLevelManager,
"CoarseSolver"));
257 RCP<TopSmootherFactory> smootherFact = rcp(
new TopSmootherFactory(coarseLevelManager,
"Smoother"));
259 int nextLevelID = coarseLevelID + 1;
261 RCP<SetFactoryManager> SFMNext;
262 if (isLastLevel ==
false) {
264 if (nextLevelID > LastLevelID())
266 CheckLevel(*Levels_[nextLevelID], nextLevelID);
270 Levels_[nextLevelID]->Request(
TopRAPFactory(coarseLevelManager, nextLevelManager));
305 RCP<Operator> Ac = Teuchos::null;
306 TopRAPFactory coarseRAPFactory(fineLevelManager, coarseLevelManager);
309 Ac = level.
Get<RCP<Operator> >(
"A");
310 }
else if (!isFinestLevel) {
316 Ac = level.
Get<RCP<Operator> >(
"A");
317 RCP<Matrix> Acm = rcp_dynamic_cast<Matrix>(Ac);
321 level.
SetComm(Ac->getDomainMap()->getComm());
324 bool isOrigLastLevel = isLastLevel;
329 }
else if (Ac.is_null()) {
336 if (!Acm.is_null() && Acm->getGlobalNumRows() <= maxCoarseSize_) {
338 GetOStream(
Runtime0) <<
"Max coarse size (<= " << maxCoarseSize_ <<
") achieved" << std::endl;
343 if (!Ac.is_null() && !isFinestLevel) {
344 RCP<Operator> A = Levels_[coarseLevelID-1]->template Get< RCP<Operator> >(
"A");
345 RCP<Matrix> Am = rcp_dynamic_cast<Matrix>(A);
347 const double maxCoarse2FineRatio = 0.8;
348 if (!Acm.is_null() && !Am.is_null() && Acm->getGlobalNumRows() > maxCoarse2FineRatio * Am->getGlobalNumRows()) {
356 GetOStream(
Warnings0) <<
"Aggregation stagnated. Please check your matrix and/or adjust your configuration file." 357 <<
"Possible fixes:\n" 358 <<
" - reduce the maximum number of levels\n" 359 <<
" - enable repartitioning\n" 360 <<
" - increase the minimum coarse size." << std::endl;
366 if (!isOrigLastLevel) {
376 coarseFact->Build(level);
387 smootherFact->Build(level);
392 if (isLastLevel ==
true) {
393 if (isOrigLastLevel ==
false) {
396 Levels_[nextLevelID]->Release(
TopRAPFactory(coarseLevelManager, nextLevelManager));
398 Levels_.resize(nextLevelID);
402 if (isDumpingEnabled_ && dumpLevel_ > 0 && coarseLevelID == dumpLevel_)
405 if (!isFinestLevel) {
409 level.
Release(coarseRAPFactory);
413 SetProcRankVerbose(oldRank);
418 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
420 int numLevels = Levels_.size();
422 "Hierarchy::SetupRe: " << Levels_.size() <<
" levels, but " << levelManagers_.size() <<
" level factory managers");
424 const int startLevel = 0;
427 #ifdef HAVE_MUELU_DEBUG 429 for (
int i = 0; i < numLevels; i++)
430 levelManagers_[i]->ResetDebugData();
435 for (levelID = startLevel; levelID < numLevels;) {
436 bool r = Setup(levelID,
437 (levelID != 0 ? levelManagers_[levelID-1] : Teuchos::null),
438 levelManagers_[levelID],
439 (levelID+1 != numLevels ? levelManagers_[levelID+1] : Teuchos::null));
445 Levels_ .resize(levelID);
446 levelManagers_.resize(levelID);
451 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
460 "MueLu::Hierarchy::Setup(): fine level (" << startLevel <<
") does not exist");
463 "Constructing non-positive (" << numDesiredLevels <<
") number of levels does not make sense.");
467 "MueLu::Hierarchy::Setup(): fine level (" << startLevel <<
") has no matrix A! " 468 "Set fine level matrix A using Level.Set()");
470 RCP<Operator> A = Levels_[startLevel]->template Get<RCP<Operator> >(
"A");
471 lib_ = A->getDomainMap()->lib();
473 RCP<const FactoryManagerBase> rcpmanager = rcpFromRef(manager);
475 const int lastLevel = startLevel + numDesiredLevels - 1;
476 GetOStream(
Runtime0) <<
"Setup loop: startLevel = " << startLevel <<
", lastLevel = " << lastLevel
477 <<
" (stop if numLevels = " << numDesiredLevels <<
" or Ac.size() < " << maxCoarseSize_ <<
")" << std::endl;
481 if (numDesiredLevels == 1) {
483 Setup(startLevel, Teuchos::null, rcpmanager, Teuchos::null);
486 bool bIsLastLevel = Setup(startLevel, Teuchos::null, rcpmanager, rcpmanager);
487 if (bIsLastLevel ==
false) {
488 for (iLevel = startLevel + 1; iLevel < lastLevel; iLevel++) {
489 bIsLastLevel = Setup(iLevel, rcpmanager, rcpmanager, rcpmanager);
490 if (bIsLastLevel ==
true)
493 if (bIsLastLevel ==
false)
494 Setup(lastLevel, rcpmanager, rcpmanager, Teuchos::null);
500 "MueLu::Hierarchy::Setup(): number of level");
509 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
511 if (startLevel < GetNumLevels())
512 GetOStream(
Runtime0) <<
"Clearing old data (if any)" << std::endl;
514 for (
int iLevel = startLevel; iLevel < GetNumLevels(); iLevel++)
515 Levels_[iLevel]->Clear();
518 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
520 GetOStream(
Runtime0) <<
"Clearing old data (expert)" << std::endl;
521 for (
int iLevel = 0; iLevel < GetNumLevels(); iLevel++)
522 Levels_[iLevel]->ExpertClear();
525 #if defined(HAVE_MUELU_EXPERIMENTAL) && defined(HAVE_MUELU_ADDITIVE_VARIANT) 526 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
528 bool InitialGuessIsZero, LO startLevel) {
529 LO nIts = conv.maxIts_;
530 MagnitudeType tol = conv.tol_;
532 std::string prefix = this->ShortClassName() +
": ";
533 std::string levelSuffix =
" (level=" +
toString(startLevel) +
")";
536 RCP<Time> CompTime = Teuchos::TimeMonitor::getNewCounter(prefix +
"Computational Time (total)");
537 RCP<Time> Concurrent = Teuchos::TimeMonitor::getNewCounter(prefix +
"Concurrent portion");
538 RCP<Time> ApplyR = Teuchos::TimeMonitor::getNewCounter(prefix +
"R: Computational Time");
539 RCP<Time> ApplyPbar = Teuchos::TimeMonitor::getNewCounter(prefix +
"Pbar: Computational Time");
540 RCP<Time> CompFine = Teuchos::TimeMonitor::getNewCounter(prefix +
"Fine: Computational Time");
541 RCP<Time> CompCoarse = Teuchos::TimeMonitor::getNewCounter(prefix +
"Coarse: Computational Time");
542 RCP<Time> ApplySum = Teuchos::TimeMonitor::getNewCounter(prefix +
"Sum: Computational Time");
543 RCP<Time> Synchronize_beginning = Teuchos::TimeMonitor::getNewCounter(prefix +
"Synchronize_beginning");
544 RCP<Time> Synchronize_center = Teuchos::TimeMonitor::getNewCounter(prefix +
"Synchronize_center");
545 RCP<Time> Synchronize_end = Teuchos::TimeMonitor::getNewCounter(prefix +
"Synchronize_end");
547 RCP<Level> Fine = Levels_[0];
550 RCP<Operator> A = Fine->Get< RCP<Operator> >(
"A");
551 Teuchos::RCP< const Teuchos::Comm< int > > communicator = A->getDomainMap()->getComm();
559 SC one = STS::one(), zero = STS::zero();
561 bool zeroGuess = InitialGuessIsZero;
567 RCP< Operator > Pbar;
569 RCP< MultiVector > coarseRhs, coarseX;
571 RCP<SmootherBase> preSmoo_coarse, postSmoo_coarse;
572 bool emptyCoarseSolve =
true;
573 RCP<MultiVector> coarseX_prolonged = MultiVectorFactory::Build(X.getMap(), X.getNumVectors(),
true);
575 RCP<const Import> importer;
577 if( Levels_.size()>1 )
580 if (Coarse->IsAvailable(
"Importer"))
581 importer = Coarse->Get< RCP<const Import> >(
"Importer");
583 R = Coarse->Get< RCP<Operator> >(
"R");
584 P = Coarse->Get< RCP<Operator> >(
"P");
587 Pbar = Coarse->Get< RCP<Operator> >(
"Pbar");
589 coarseRhs = MultiVectorFactory::Build(R->getRangeMap(), B.getNumVectors(),
true);
591 Ac = Coarse->Get< RCP< Operator > >(
"A");
594 R->apply(B, *coarseRhs, Teuchos::NO_TRANS, one, zero);
598 if (doPRrebalance_ || importer.is_null()) {
600 coarseX = MultiVectorFactory::Build(coarseRhs->getMap(), X.getNumVectors(),
true);
604 RCP<TimeMonitor> ITime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : import (total)" ,
Timings0));
605 RCP<TimeMonitor> ILevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : import" + levelSuffix,
Timings0));
608 RCP<MultiVector> coarseTmp = MultiVectorFactory::Build(importer->getTargetMap(), coarseRhs->getNumVectors());
609 coarseTmp->doImport(*coarseRhs, *importer, Xpetra::INSERT);
610 coarseRhs.swap(coarseTmp);
612 coarseX = MultiVectorFactory::Build(importer->getTargetMap(), X.getNumVectors(),
true);
615 if (Coarse->IsAvailable(
"PreSmoother"))
616 preSmoo_coarse = Coarse->Get< RCP<SmootherBase> >(
"PreSmoother");
617 if (Coarse->IsAvailable(
"PostSmoother"))
618 postSmoo_coarse = Coarse->Get< RCP<SmootherBase> >(
"PostSmoother");
625 MagnitudeType prevNorm = STS::magnitude(STS::one()), curNorm = STS::magnitude(STS::one());
628 for (LO i = 1; i <= nIts; i++) {
629 #ifdef HAVE_MUELU_DEBUG 630 if (A->getDomainMap()->isCompatible(*(X.getMap())) ==
false) {
631 std::ostringstream ss;
632 ss <<
"Level " << startLevel <<
": level A's domain map is not compatible with X";
633 throw Exceptions::Incompatible(ss.str());
636 if (A->getRangeMap()->isCompatible(*(B.getMap())) ==
false) {
637 std::ostringstream ss;
638 ss <<
"Level " << startLevel <<
": level A's range map is not compatible with B";
639 throw Exceptions::Incompatible(ss.str());
644 bool emptyFineSolve =
true;
646 RCP< MultiVector > fineX;
647 fineX = MultiVectorFactory::Build(X.getMap(), X.getNumVectors(),
true);
656 if (Fine->IsAvailable(
"PreSmoother")) {
657 RCP<SmootherBase> preSmoo = Fine->Get< RCP<SmootherBase> >(
"PreSmoother");
659 preSmoo->Apply(*fineX, B, zeroGuess);
661 emptyFineSolve =
false;
663 if (Fine->IsAvailable(
"PostSmoother")) {
664 RCP<SmootherBase> postSmoo = Fine->Get< RCP<SmootherBase> >(
"PostSmoother");
667 postSmoo->Apply(*fineX, B, zeroGuess);
670 emptyFineSolve =
false;
672 if (emptyFineSolve ==
true) {
673 GetOStream(
Warnings1) <<
"No fine grid smoother" << std::endl;
675 fineX->update(one, B, zero);
678 if( Levels_.size()>1 )
681 if (Coarse->IsAvailable(
"PreSmoother")) {
683 preSmoo_coarse->Apply(*coarseX, *coarseRhs, zeroGuess);
685 emptyCoarseSolve =
false;
688 if (Coarse->IsAvailable(
"PostSmoother")) {
690 postSmoo_coarse->Apply(*coarseX, *coarseRhs, zeroGuess);
692 emptyCoarseSolve =
false;
695 if (emptyCoarseSolve ==
true) {
696 GetOStream(
Warnings1) <<
"No coarse grid solver" << std::endl;
698 coarseX->update(one, *coarseRhs, zero);
706 if (!doPRrebalance_ && !importer.is_null()) {
707 RCP<TimeMonitor> ITime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : export (total)" ,
Timings0));
708 RCP<TimeMonitor> ILevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : export" + levelSuffix,
Timings0));
711 RCP<MultiVector> coarseTmp = MultiVectorFactory::Build(importer->getSourceMap(), coarseX->getNumVectors());
712 coarseTmp->doExport(*coarseX, *importer, Xpetra::INSERT);
713 coarseX.swap(coarseTmp);
717 Pbar->apply(*coarseX, *coarseX_prolonged, Teuchos::NO_TRANS, one, zero);
722 X.update(1.0, *fineX, 1.0, *coarseX_prolonged, 0.0);
733 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
735 bool InitialGuessIsZero, LO startLevel) {
750 std::string prefix = this->ShortClassName() +
": ";
751 std::string levelSuffix =
" (level=" +
toString(startLevel) +
")";
752 RCP<Monitor> iterateTime;
753 RCP<TimeMonitor> iterateTime1;
759 std::string iterateLevelTimeLabel = prefix +
"Solve" + levelSuffix;
760 RCP<TimeMonitor> iterateLevelTime = rcp(
new TimeMonitor(*
this, iterateLevelTimeLabel,
Timings0));
762 bool zeroGuess = InitialGuessIsZero;
764 RCP<Level> Fine = Levels_[startLevel];
765 RCP<Operator> A = Fine->Get< RCP<Operator> >(
"A");
767 RCP<Time> CompCoarse = Teuchos::TimeMonitor::getNewCounter(prefix +
"Coarse: Computational Time");
777 MagnitudeType prevNorm = STS::magnitude(STS::one()), curNorm = STS::magnitude(STS::one());
779 if (startLevel == 0 && !isPreconditioner_ &&
783 Teuchos::Array<MagnitudeType> rn;
788 for (LO k = 0; k < rn.size(); k++)
798 << std::setiosflags(std::ios::left)
799 << std::setprecision(3) << 0
801 << std::setprecision(10) << rn
805 SC one = STS::one(), zero = STS::zero();
806 for (LO i = 1; i <= nIts; i++) {
807 #ifdef HAVE_MUELU_DEBUG 808 if (A->getDomainMap()->isCompatible(*(X.getMap())) ==
false) {
809 std::ostringstream ss;
810 ss <<
"Level " << startLevel <<
": level A's domain map is not compatible with X";
814 if (A->getRangeMap()->isCompatible(*(B.getMap())) ==
false) {
815 std::ostringstream ss;
816 ss <<
"Level " << startLevel <<
": level A's range map is not compatible with B";
821 if (startLevel == as<LO>(Levels_.size())-1) {
823 RCP<TimeMonitor> CLevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : coarse" + levelSuffix,
Timings0));
825 bool emptySolve =
true;
828 if (Fine->IsAvailable(
"PreSmoother")) {
829 RCP<SmootherBase> preSmoo = Fine->Get< RCP<SmootherBase> >(
"PreSmoother");
831 preSmoo->Apply(X, B, zeroGuess);
836 if (Fine->IsAvailable(
"PostSmoother")) {
837 RCP<SmootherBase> postSmoo = Fine->Get< RCP<SmootherBase> >(
"PostSmoother");
839 postSmoo->Apply(X, B, zeroGuess);
843 if (emptySolve ==
true) {
844 GetOStream(
Warnings1) <<
"No coarse grid solver" << std::endl;
846 X.update(one, B, zero);
851 RCP<Level> Coarse = Levels_[startLevel+1];
854 RCP<TimeMonitor> STime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : smoothing (total)" ,
Timings0));
855 RCP<TimeMonitor> SLevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : smoothing" + levelSuffix,
Timings0));
857 if (Fine->IsAvailable(
"PreSmoother")) {
858 RCP<SmootherBase> preSmoo = Fine->Get< RCP<SmootherBase> >(
"PreSmoother");
859 preSmoo->Apply(X, B, zeroGuess);
861 GetOStream(
Warnings1) <<
"Level " << startLevel <<
": No PreSmoother!" << std::endl;
865 RCP<MultiVector> residual;
867 RCP<TimeMonitor> ATime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : residual calculation (total)" ,
Timings0));
868 RCP<TimeMonitor> ALevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : residual calculation" + levelSuffix,
Timings0));
872 RCP<Operator> P = Coarse->Get< RCP<Operator> >(
"P");
873 if (Coarse->IsAvailable(
"Pbar"))
874 P = Coarse->Get< RCP<Operator> >(
"Pbar");
876 RCP<MultiVector> coarseRhs, coarseX;
877 const bool initializeWithZeros =
true;
880 RCP<TimeMonitor> RTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : restriction (total)" ,
Timings0));
881 RCP<TimeMonitor> RLevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : restriction" + levelSuffix,
Timings0));
883 if (implicitTranspose_) {
884 coarseRhs = MultiVectorFactory::Build(P->getDomainMap(), X.getNumVectors(), !initializeWithZeros);
885 P->apply(*residual, *coarseRhs, Teuchos::TRANS, one, zero);
888 RCP<Operator> R = Coarse->Get< RCP<Operator> >(
"R");
889 coarseRhs = MultiVectorFactory::Build(R->getRangeMap(), X.getNumVectors(), !initializeWithZeros);
890 R->apply(*residual, *coarseRhs, Teuchos::NO_TRANS, one, zero);
894 RCP<const Import> importer;
895 if (Coarse->IsAvailable(
"Importer"))
896 importer = Coarse->Get< RCP<const Import> >(
"Importer");
898 if (doPRrebalance_ || importer.is_null()) {
900 coarseX = MultiVectorFactory::Build(coarseRhs->getMap(), X.getNumVectors(), initializeWithZeros);
904 RCP<TimeMonitor> ITime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : import (total)" ,
Timings0));
905 RCP<TimeMonitor> ILevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : import" + levelSuffix,
Timings0));
908 RCP<MultiVector> coarseTmp = MultiVectorFactory::Build(importer->getTargetMap(), coarseRhs->getNumVectors());
909 coarseTmp->doImport(*coarseRhs, *importer, Xpetra::INSERT);
910 coarseRhs.swap(coarseTmp);
912 coarseX = MultiVectorFactory::Build(importer->getTargetMap(), X.getNumVectors(), initializeWithZeros);
915 RCP<Operator> Ac = Coarse->Get< RCP<Operator> >(
"A");
917 RCP<const Map> origXMap = coarseX->getMap();
920 coarseRhs->replaceMap(Ac->getRangeMap());
921 coarseX ->replaceMap(Ac->getDomainMap());
924 iterateLevelTime = Teuchos::null;
926 Iterate(*coarseRhs, *coarseX, 1,
true, startLevel+1);
929 Iterate(*coarseRhs, *coarseX, 1,
false, startLevel+1);
932 iterateLevelTime = rcp(
new TimeMonitor(*
this, iterateLevelTimeLabel));
934 coarseX->replaceMap(origXMap);
937 if (!doPRrebalance_ && !importer.is_null()) {
938 RCP<TimeMonitor> ITime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : export (total)" ,
Timings0));
939 RCP<TimeMonitor> ILevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : export" + levelSuffix,
Timings0));
942 RCP<MultiVector> coarseTmp = MultiVectorFactory::Build(importer->getSourceMap(), coarseX->getNumVectors());
943 coarseTmp->doExport(*coarseX, *importer, Xpetra::INSERT);
944 coarseX.swap(coarseTmp);
951 RCP<MultiVector> correction = MultiVectorFactory::Build(P->getRangeMap(), X.getNumVectors(),
false);
954 RCP<TimeMonitor> PTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : prolongation (total)" ,
Timings0));
955 RCP<TimeMonitor> PLevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : prolongation" + levelSuffix,
Timings0));
956 P->apply(*coarseX, *correction, Teuchos::NO_TRANS, one, zero);
958 X.update(one, *correction, one);
962 RCP<TimeMonitor> STime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : smoothing (total)" ,
Timings0));
963 RCP<TimeMonitor> SLevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : smoothing" + levelSuffix,
Timings0));
965 if (Fine->IsAvailable(
"PostSmoother")) {
966 RCP<SmootherBase> postSmoo = Fine->Get< RCP<SmootherBase> >(
"PostSmoother");
967 postSmoo->Apply(X, B,
false);
970 GetOStream(
Warnings1) <<
"Level " << startLevel <<
": No PostSmoother!" << std::endl;
976 if (startLevel == 0 && !isPreconditioner_ &&
980 Teuchos::Array<MagnitudeType> rn;
985 rate_ = as<MagnitudeType>(curNorm / prevNorm);
989 for (LO k = 0; k < rn.size(); k++)
999 << std::setiosflags(std::ios::left)
1000 << std::setprecision(3) << i
1002 << std::setprecision(10) << rn
1011 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1013 LO startLevel = (start != -1 ? start : 0);
1014 LO endLevel = (end != -1 ? end : Levels_.size()-1);
1017 "MueLu::Hierarchy::Write : startLevel must be <= endLevel");
1020 "MueLu::Hierarchy::Write bad start or end level");
1022 for (LO i = startLevel; i < endLevel + 1; i++) {
1023 RCP<Matrix> A = rcp_dynamic_cast<Matrix>(Levels_[i]->
template Get< RCP< Operator> >(
"A")), P, R;
1025 P = rcp_dynamic_cast<Matrix>(Levels_[i]->
template Get< RCP< Operator> >(
"P"));
1026 if (!implicitTranspose_)
1027 R = rcp_dynamic_cast<Matrix>(Levels_[i]->
template Get< RCP< Operator> >(
"R"));
1030 if (!A.is_null()) Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(
"A_" +
toString(i) + suffix +
".m", *A);
1031 if (!P.is_null()) Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(
"P_" +
toString(i) + suffix +
".m", *P);
1032 if (!R.is_null()) Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(
"R_" +
toString(i) + suffix +
".m", *R);
1036 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1038 for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1039 (*it)->Keep(ename, factory);
1042 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1044 for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1045 (*it)->Delete(ename, factory);
1048 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1050 for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1051 (*it)->AddKeepFlag(ename, factory, keep);
1054 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1056 for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1057 (*it)->RemoveKeepFlag(ename, factory, keep);
1060 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1062 std::ostringstream out;
1064 out <<
"{#levels = " << GetGlobalNumLevels() <<
", complexity = " << GetOperatorComplexity() <<
"}";
1068 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1073 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1075 RCP<Operator> A0 = Levels_[0]->template Get<RCP<Operator> >(
"A");
1076 RCP<const Teuchos::Comm<int> > comm = A0->getDomainMap()->getComm();
1078 int numLevels = GetNumLevels();
1079 RCP<Operator> Ac = Levels_[numLevels-1]->template Get<RCP<Operator> >(
"A");
1086 int root = comm->getRank();
1089 int smartData = numLevels*comm->getSize() + comm->getRank(), maxSmartData;
1090 reduceAll(*comm, Teuchos::REDUCE_MAX, smartData, Teuchos::ptr(&maxSmartData));
1091 root = maxSmartData % comm->getSize();
1096 std::vector<Xpetra::global_size_t> nnzPerLevel;
1097 std::vector<Xpetra::global_size_t> rowsPerLevel;
1098 std::vector<int> numProcsPerLevel;
1099 bool aborted =
false;
1100 for (
int i = 0; i < numLevels; i++) {
1102 "Operator A is not available on level " << i);
1104 RCP<Operator> A = Levels_[i]->template Get<RCP<Operator> >(
"A");
1106 "Operator A on level " << i <<
" is null.");
1108 RCP<Matrix> Am = rcp_dynamic_cast<Matrix>(A);
1110 GetOStream(
Warnings0) <<
"Some level operators are not matrices, statistics calculation aborted" << std::endl;
1115 Xpetra::global_size_t nnz = Am->getGlobalNumEntries();
1116 nnzPerLevel .push_back(nnz);
1117 rowsPerLevel .push_back(Am->getGlobalNumRows());
1118 numProcsPerLevel.push_back(Am->getRowMap()->getComm()->getSize());
1122 std::ostringstream oss;
1123 oss <<
"\n--------------------------------------------------------------------------------\n" <<
1124 "--- Multigrid Summary ---\n" 1125 "--------------------------------------------------------------------------------" << std::endl;
1126 oss <<
"Number of levels = " << numLevels << std::endl;
1127 oss <<
"Operator complexity = " << std::setprecision(2) << std::setiosflags(std::ios::fixed)
1128 << GetOperatorComplexity() << std::endl;
1131 oss <<
"Cycle type = V" << std::endl;
1134 oss <<
"Cycle type = W" << std::endl;
1141 Xpetra::global_size_t tt = rowsPerLevel[0];
1142 int rowspacer = 2;
while (tt != 0) { tt /= 10; rowspacer++; }
1143 tt = nnzPerLevel[0];
1144 int nnzspacer = 2;
while (tt != 0) { tt /= 10; nnzspacer++; }
1145 tt = numProcsPerLevel[0];
1146 int npspacer = 2;
while (tt != 0) { tt /= 10; npspacer++; }
1147 oss <<
"level " << std::setw(rowspacer) <<
" rows " << std::setw(nnzspacer) <<
" nnz " <<
" nnz/row" << std::setw(npspacer) <<
" c ratio" <<
" procs" << std::endl;
1148 for (
size_t i = 0; i < nnzPerLevel.size(); ++i) {
1149 oss <<
" " << i <<
" ";
1150 oss << std::setw(rowspacer) << rowsPerLevel[i];
1151 oss << std::setw(nnzspacer) << nnzPerLevel[i];
1152 oss << std::setprecision(2) << std::setiosflags(std::ios::fixed);
1153 oss << std::setw(9) << as<double>(nnzPerLevel[i]) / rowsPerLevel[i];
1154 if (i) oss << std::setw(9) << as<double>(rowsPerLevel[i-1])/rowsPerLevel[i];
1155 else oss << std::setw(9) <<
" ";
1156 oss <<
" " << std::setw(npspacer) << numProcsPerLevel[i] << std::endl;
1159 for (
int i = 0; i < GetNumLevels(); ++i) {
1160 RCP<SmootherBase> preSmoo, postSmoo;
1161 if (Levels_[i]->IsAvailable(
"PreSmoother"))
1162 preSmoo = Levels_[i]->
template Get< RCP<SmootherBase> >(
"PreSmoother");
1163 if (Levels_[i]->IsAvailable(
"PostSmoother"))
1164 postSmoo = Levels_[i]->
template Get< RCP<SmootherBase> >(
"PostSmoother");
1166 if (preSmoo != null && preSmoo == postSmoo)
1167 oss <<
"Smoother (level " << i <<
") both : " << preSmoo->description() << std::endl;
1169 oss <<
"Smoother (level " << i <<
") pre : " 1170 << (preSmoo != null ? preSmoo->description() :
"no smoother") << std::endl;
1171 oss <<
"Smoother (level " << i <<
") post : " 1172 << (postSmoo != null ? postSmoo->description() :
"no smoother") << std::endl;
1183 RCP<const Teuchos::MpiComm<int> > mpiComm = rcp_dynamic_cast<
const Teuchos::MpiComm<int> >(comm);
1184 MPI_Comm rawComm = (*mpiComm->getRawMpiComm())();
1186 int strLength = outstr.size();
1187 MPI_Bcast(&strLength, 1, MPI_INT, root, rawComm);
1188 if (comm->getRank() != root)
1189 outstr.resize(strLength);
1190 MPI_Bcast(&outstr[0], strLength, MPI_CHAR, root, rawComm);
1197 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1199 Teuchos::OSTab tab2(out);
1200 for (
int i = 0; i < GetNumLevels(); ++i)
1201 Levels_[i]->print(out, verbLevel);
1204 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1206 isPreconditioner_ = flag;
1209 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1211 if (GetProcRankVerbose() != 0)
1213 #if defined(HAVE_MUELU_BOOST) && defined(HAVE_MUELU_BOOST_FOR_REAL) && defined(BOOST_VERSION) && (BOOST_VERSION >= 104400) 1217 dp.property(
"label", boost::get(boost::vertex_name, graph));
1218 dp.property(
"id", boost::get(boost::vertex_index, graph));
1219 dp.property(
"label", boost::get(boost::edge_name, graph));
1220 dp.property(
"color", boost::get(boost::edge_color, graph));
1223 std::map<const FactoryBase*, BoostVertex> vindices;
1224 typedef std::map<std::pair<BoostVertex,BoostVertex>, std::string> emap; emap edges;
1226 for (
int i = dumpLevel_; i <= dumpLevel_+1 && i < GetNumLevels(); i++) {
1228 Levels_[i]->UpdateGraph(vindices, edges, dp, graph);
1230 for (emap::const_iterator eit = edges.begin(); eit != edges.end(); eit++) {
1231 std::pair<BoostEdge, bool> boost_edge = boost::add_edge(eit->first.first, eit->first.second, graph);
1232 boost::put(
"label", dp, boost_edge.first, eit->second);
1233 if (i == dumpLevel_)
1234 boost::put(
"color", dp, boost_edge.first, std::string(
"red"));
1236 boost::put(
"color", dp, boost_edge.first, std::string(
"blue"));
1241 std::ostringstream legend;
1242 legend <<
"< <TABLE BORDER=\"0\" CELLBORDER=\"1\" CELLSPACING=\"0\" CELLPADDING=\"4\"> \ 1243 <TR><TD COLSPAN=\"2\">Legend</TD></TR> \ 1244 <TR><TD><FONT color=\"red\">Level " << dumpLevel_ <<
"</FONT></TD><TD><FONT color=\"blue\">Level " << dumpLevel_+1 <<
"</FONT></TD></TR> \ 1246 BoostVertex boost_vertex = boost::add_vertex(graph);
1247 boost::put(
"label", dp, boost_vertex, legend.str());
1249 std::ofstream out(dumpFile_.c_str());
1250 boost::write_graphviz_dp(out, graph, dp, std::string(
"id"));
1252 GetOStream(
Errors) <<
"Dependency graph output requires boost and MueLu_ENABLE_Boost_for_real" << std::endl;
1257 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1259 RCP<Operator> Ao = level.
Get<RCP<Operator> >(
"A");
1260 RCP<Matrix> A = rcp_dynamic_cast<Matrix>(Ao);
1262 GetOStream(
Warnings1) <<
"Hierarchy::ReplaceCoordinateMap: operator is not a matrix, skipping..." << std::endl;
1266 typedef Xpetra::MultiVector<double,LO,GO,NO> xdMV;
1268 RCP<xdMV> coords = level.
Get<RCP<xdMV> >(
"Coordinates");
1270 if (A->getRowMap()->isSameAs(*(coords->getMap()))) {
1271 GetOStream(
Warnings1) <<
"Hierarchy::ReplaceCoordinateMap: matrix and coordinates maps are same, skipping..." << std::endl;
1275 if (A->IsView(
"stridedMaps") && rcp_dynamic_cast<
const StridedMap>(A->getRowMap(
"stridedMaps")) != Teuchos::null) {
1276 RCP<const StridedMap> stridedRowMap = rcp_dynamic_cast<
const StridedMap>(A->getRowMap(
"stridedMaps"));
1279 TEUCHOS_TEST_FOR_EXCEPTION(stridedRowMap->getStridedBlockId() != -1 || stridedRowMap->getOffset() != 0,
1280 Exceptions::RuntimeError,
"Hierarchy::ReplaceCoordinateMap: nontrivial maps (block id = " << stridedRowMap->getStridedBlockId()
1281 <<
", offset = " << stridedRowMap->getOffset() <<
")");
1284 GetOStream(
Runtime1) <<
"Replacing coordinate map" << std::endl;
1286 size_t blkSize = A->GetFixedBlockSize();
1288 RCP<const Map> nodeMap = A->getRowMap();
1291 RCP<const Map> dofMap = A->getRowMap();
1292 GO indexBase = dofMap->getIndexBase();
1293 size_t numLocalDOFs = dofMap->getNodeNumElements();
1295 "Hierarchy::ReplaceCoordinateMap: block size (" << blkSize <<
") is incompatible with the number of local dofs in a row map (" << numLocalDOFs);
1296 ArrayView<const GO> GIDs = dofMap->getNodeElementList();
1298 Array<GO> nodeGIDs(numLocalDOFs/blkSize);
1299 for (
size_t i = 0; i < numLocalDOFs; i += blkSize)
1300 nodeGIDs[i/blkSize] = (GIDs[i] - indexBase)/blkSize + indexBase;
1302 Xpetra::global_size_t INVALID = Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid();
1303 nodeMap = MapFactory::Build(dofMap->lib(), INVALID, nodeGIDs(), indexBase, dofMap->getComm());
1306 Array<ArrayView<const double> > coordDataView;
1307 std::vector<ArrayRCP<const double> > coordData;
1308 for (
size_t i = 0; i < coords->getNumVectors(); i++) {
1309 coordData.push_back(coords->getData(i));
1310 coordDataView.push_back(coordData[i]());
1313 RCP<xdMV> newCoords = Xpetra::MultiVectorFactory<double,LO,GO,NO>::Build(nodeMap, coordDataView(), coords->getNumVectors());
1314 level.
Set(
"Coordinates", newCoords);
1319 #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)
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.
void Release(const FactoryBase &factory)
Decrement the storage counter for all the inputs of a factory.
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.
void AddNewLevel()
Add a new level at the end of the hierarchy.
VerbLevel toMueLuVerbLevel(const Teuchos::EVerbosityLevel verbLevel)
One-liner description of what is happening.
void Clear(int startLevel=0)
Clear impermanent data from previous setup.
Namespace for MueLu classes and methods.
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 AddKeepFlag(const std::string &ename, const FactoryBase *factory=NoFactory::get(), KeepType keep=MueLu::Keep)
Call Level::AddKeepFlag for each level of the Hierarchy.
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 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.
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()
void describe(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the Hierarchy with some verbosity level to a FancyOStream object.
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
double GetOperatorComplexity() const
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)
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.
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)
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.
void Request(const FactoryBase &factory)
Increment the storage counter for all the inputs of a factory.