46 #ifndef MUELU_HIERARCHY_DEF_HPP 47 #define MUELU_HIERARCHY_DEF_HPP 52 #include <Xpetra_Matrix.hpp> 53 #include <Xpetra_MultiVectorFactory.hpp> 54 #include <Xpetra_Operator.hpp> 55 #include <Xpetra_IO.hpp> 60 #include "MueLu_FactoryManager.hpp" 61 #include "MueLu_HierarchyHelpers.hpp" 64 #include "MueLu_PFactory.hpp" 66 #include "MueLu_SmootherFactory.hpp" 71 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
73 : maxCoarseSize_(GetDefaultMaxCoarseSize()), implicitTranspose_(GetDefaultImplicitTranspose()),
74 doPRrebalance_(GetDefaultPRrebalance()), isPreconditioner_(true), Cycle_(GetDefaultCycle()),
75 lib_(Xpetra::UseTpetra), isDumpingEnabled_(false), dumpLevel_(-1), rate_(-1)
80 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
86 lib_ = A->getDomainMap()->lib();
88 RCP<Level> Finest = rcp(
new Level);
94 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
98 if (level->GetLevelID() != -1 && (level->GetLevelID() != levelID))
100 " have been added at the end of the hierarchy\n but its ID have been redefined" <<
101 " because last level ID of the hierarchy was " <<
LastLevelID() <<
"." << std::endl;
104 level->SetLevelID(levelID);
107 level->SetPreviousLevel( (levelID == 0) ? Teuchos::null :
Levels_[
LastLevelID() - 1] );
110 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
113 newLevel->setlib(
lib_);
117 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
120 "MueLu::Hierarchy::GetLevel(): invalid input parameter value: LevelID = " << levelID);
124 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
129 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
131 RCP<Operator> A =
Levels_[0]->template Get<RCP<Operator> >(
"A");
132 RCP<const Teuchos::Comm<int> > comm = A->getDomainMap()->getComm();
136 Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, numLevels, Teuchos::ptr(&numGlobalLevels));
138 return numGlobalLevels;
141 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
143 double totalNnz = 0, lev0Nnz = 1;
146 "Operator complexity cannot be calculated because A is unavailable on level " << i);
147 RCP<Operator> A =
Levels_[i]->template Get<RCP<Operator> >(
"A");
151 RCP<Matrix> Am = rcp_dynamic_cast<Matrix>(A);
153 GetOStream(
Warnings0) <<
"Some level operators are not matrices, operator complexity calculation aborted" << std::endl;
157 totalNnz += as<double>(Am->getGlobalNumEntries());
161 return totalNnz / lev0Nnz;
165 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
168 "MueLu::Hierarchy::CheckLevel(): wrong underlying linear algebra library.");
169 TEUCHOS_TEST_FOR_EXCEPTION(level.
GetLevelID() != levelID, Exceptions::RuntimeError,
170 "MueLu::Hierarchy::CheckLevel(): wrong level ID");
171 TEUCHOS_TEST_FOR_EXCEPTION(levelID != 0 && level.
GetPreviousLevel() !=
Levels_[levelID-1], Exceptions::RuntimeError,
172 "MueLu::Hierarchy::Setup(): wrong level parent");
177 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
179 const RCP<const FactoryManagerBase> fineLevelManager,
180 const RCP<const FactoryManagerBase> coarseLevelManager,
181 const RCP<const FactoryManagerBase> nextLevelManager) {
189 "MueLu::Hierarchy::Setup(): argument coarseLevelManager cannot be null");
195 "MueLu::Hierarchy:Setup(): level " << coarseLevelID <<
" (specified by coarseLevelID argument) " 196 "must be built before calling this function.");
204 bool isFinestLevel = (fineLevelManager.is_null());
205 bool isLastLevel = (nextLevelManager.is_null());
209 RCP<Operator> A = level.
Get< RCP<Operator> >(
"A");
210 RCP<const Map> domainMap = A->getDomainMap();
211 RCP<const Teuchos::Comm<int> > comm = domainMap->getComm();
222 lib_ = domainMap->lib();
236 RCP<SetFactoryManager> SFMFine;
240 if (isFinestLevel &&
Levels_[coarseLevelID]->IsAvailable(
"Coordinates"))
249 RCP<TopSmootherFactory> coarseFact = rcp(
new TopSmootherFactory(coarseLevelManager,
"CoarseSolver"));
250 RCP<TopSmootherFactory> smootherFact = rcp(
new TopSmootherFactory(coarseLevelManager,
"Smoother"));
252 int nextLevelID = coarseLevelID + 1;
254 RCP<SetFactoryManager> SFMNext;
255 if (isLastLevel ==
false) {
263 Levels_[nextLevelID]->Request(TopRAPFactory(coarseLevelManager, nextLevelManager));
298 RCP<Operator> Ac = Teuchos::null;
299 TopRAPFactory coarseRAPFactory(fineLevelManager, coarseLevelManager);
302 Ac = level.
Get<RCP<Operator> >(
"A");
303 }
else if (!isFinestLevel) {
309 Ac = level.
Get<RCP<Operator> >(
"A");
310 RCP<Matrix> Acm = rcp_dynamic_cast<Matrix>(Ac);
314 level.
SetComm(Ac->getDomainMap()->getComm());
317 bool isOrigLastLevel = isLastLevel;
322 }
else if (Ac.is_null()) {
329 if (!Acm.is_null() && Acm->getGlobalNumRows() <=
maxCoarseSize_) {
336 if (!Ac.is_null() && !isFinestLevel) {
337 RCP<Operator> A =
Levels_[coarseLevelID-1]->template Get< RCP<Operator> >(
"A");
338 RCP<Matrix> Am = rcp_dynamic_cast<Matrix>(A);
340 const double maxCoarse2FineRatio = 0.8;
341 if (!Acm.is_null() && !Am.is_null() && Acm->getGlobalNumRows() > maxCoarse2FineRatio * Am->getGlobalNumRows()) {
349 GetOStream(
Warnings0) <<
"Aggregation stagnated. Please check your matrix and/or adjust your configuration file." 350 <<
"Possible fixes:\n" 351 <<
" - reduce the maximum number of levels\n" 352 <<
" - enable repartitioning\n" 353 <<
" - increase the minimum coarse size." << std::endl;
359 if (!isOrigLastLevel) {
369 coarseFact->Build(level);
380 smootherFact->Build(level);
385 if (isLastLevel ==
true) {
386 if (isOrigLastLevel ==
false) {
389 Levels_[nextLevelID]->Release(TopRAPFactory(coarseLevelManager, nextLevelManager));
398 if (!isFinestLevel) {
402 level.
Release(coarseRAPFactory);
411 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
413 int numLevels =
Levels_.size();
415 "Hierarchy::SetupRe: " <<
Levels_.size() <<
" levels, but " <<
levelManagers_.size() <<
" level factory managers");
417 #ifdef HAVE_MUELU_DEBUG 419 for (
int i = 0; i < numLevels; i++)
425 for (levelID = 0; levelID < numLevels;) {
426 bool r =
Setup(levelID,
429 (levelID+1 != numLevels ?
levelManagers_[levelID+1] : Teuchos::null));
441 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
450 "MueLu::Hierarchy::Setup(): fine level (" << startLevel <<
") does not exist");
452 TEUCHOS_TEST_FOR_EXCEPTION(numDesiredLevels <= 0, Exceptions::RuntimeError,
453 "Constructing non-positive (" << numDesiredLevels <<
") number of levels does not make sense.");
456 TEUCHOS_TEST_FOR_EXCEPTION(!
Levels_[startLevel]->IsAvailable(
"A"), Exceptions::RuntimeError,
457 "MueLu::Hierarchy::Setup(): fine level (" << startLevel <<
") has no matrix A! " 458 "Set fine level matrix A using Level.Set()");
460 RCP<Operator> A =
Levels_[startLevel]->template Get<RCP<Operator> >(
"A");
461 lib_ = A->getDomainMap()->lib();
463 RCP<const FactoryManagerBase> rcpmanager = rcpFromRef(manager);
465 const int lastLevel = startLevel + numDesiredLevels - 1;
466 GetOStream(
Runtime0) <<
"Setup loop: startLevel = " << startLevel <<
", lastLevel = " << lastLevel
467 <<
" (stop if numLevels = " << numDesiredLevels <<
" or Ac.size() < " <<
maxCoarseSize_ <<
")" << std::endl;
471 if (numDesiredLevels == 1) {
473 Setup(startLevel, Teuchos::null, rcpmanager, Teuchos::null);
476 bool bIsLastLevel =
Setup(startLevel, Teuchos::null, rcpmanager, rcpmanager);
477 if (bIsLastLevel ==
false) {
478 for (iLevel = startLevel + 1; iLevel < lastLevel; iLevel++) {
479 bIsLastLevel =
Setup(iLevel, rcpmanager, rcpmanager, rcpmanager);
480 if (bIsLastLevel ==
true)
483 if (bIsLastLevel ==
false)
484 Setup(lastLevel, rcpmanager, rcpmanager, Teuchos::null);
489 TEUCHOS_TEST_FOR_EXCEPTION(iLevel !=
Levels_.size() - 1, Exceptions::RuntimeError,
490 "MueLu::Hierarchy::Setup(): number of level");
499 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
504 for (
int iLevel = startLevel; iLevel <
GetNumLevels(); iLevel++)
508 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
517 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
519 bool InitialGuessIsZero, LO startLevel) {
535 std::string levelSuffix =
" (level=" +
toString(startLevel) +
")";
536 RCP<Monitor> iterateTime;
537 RCP<TimeMonitor> iterateTime1;
543 std::string iterateLevelTimeLabel = prefix +
"Solve" + levelSuffix;
544 RCP<TimeMonitor> iterateLevelTime = rcp(
new TimeMonitor(*
this, iterateLevelTimeLabel,
Timings0));
546 bool zeroGuess = InitialGuessIsZero;
548 RCP<Level> Fine =
Levels_[startLevel];
549 RCP<Operator> A = Fine->Get< RCP<Operator> >(
"A");
559 MagnitudeType prevNorm = STS::magnitude(STS::one()), curNorm = STS::magnitude(STS::one());
565 Teuchos::Array<MagnitudeType> rn;
570 for (LO k = 0; k < rn.size(); k++)
580 << std::setiosflags(std::ios::left)
581 << std::setprecision(3) << 0
583 << std::setprecision(10) << rn
587 SC one = STS::one(), zero = STS::zero();
588 for (LO i = 1; i <= nIts; i++) {
589 #ifdef HAVE_MUELU_DEBUG 590 if (A->getDomainMap()->isCompatible(*(X.getMap())) ==
false) {
591 std::ostringstream ss;
592 ss <<
"Level " << startLevel <<
": level A's domain map is not compatible with X";
596 if (A->getRangeMap()->isCompatible(*(B.getMap())) ==
false) {
597 std::ostringstream ss;
598 ss <<
"Level " << startLevel <<
": level A's range map is not compatible with B";
603 if (startLevel == as<LO>(
Levels_.size())-1) {
605 RCP<TimeMonitor> CLevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : coarse" + levelSuffix,
Timings0));
607 bool emptySolve =
true;
610 if (Fine->IsAvailable(
"PreSmoother")) {
611 RCP<SmootherBase> preSmoo = Fine->Get< RCP<SmootherBase> >(
"PreSmoother");
612 preSmoo->Apply(X, B, zeroGuess);
616 if (Fine->IsAvailable(
"PostSmoother")) {
617 RCP<SmootherBase> postSmoo = Fine->Get< RCP<SmootherBase> >(
"PostSmoother");
618 postSmoo->Apply(X, B, zeroGuess);
621 if (emptySolve ==
true)
626 RCP<Level> Coarse =
Levels_[startLevel+1];
630 RCP<TimeMonitor> STime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : smoothing (total)" ,
Timings0));
631 RCP<TimeMonitor> SLevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : smoothing" + levelSuffix,
Timings0));
633 if (Fine->IsAvailable(
"PreSmoother")) {
634 RCP<SmootherBase> preSmoo = Fine->Get< RCP<SmootherBase> >(
"PreSmoother");
635 preSmoo->Apply(X, B, zeroGuess);
641 RCP<MultiVector> residual;
643 RCP<TimeMonitor> ATime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : residual calculation (total)" ,
Timings0));
644 RCP<TimeMonitor> ALevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : residual calculation" + levelSuffix,
Timings0));
648 RCP<Operator> P = Coarse->Get< RCP<Operator> >(
"P");
649 RCP<MultiVector> coarseRhs, coarseX;
650 const bool initializeWithZeros =
true;
653 RCP<TimeMonitor> RTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : restriction (total)" ,
Timings0));
654 RCP<TimeMonitor> RLevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : restriction" + levelSuffix,
Timings0));
657 coarseRhs = MultiVectorFactory::Build(P->getDomainMap(), X.getNumVectors(), !initializeWithZeros);
658 P->apply(*residual, *coarseRhs, Teuchos::TRANS, one, zero);
661 RCP<Operator> R = Coarse->Get< RCP<Operator> >(
"R");
662 coarseRhs = MultiVectorFactory::Build(R->getRangeMap(), X.getNumVectors(), !initializeWithZeros);
663 R->apply(*residual, *coarseRhs, Teuchos::NO_TRANS, one, zero);
667 RCP<const Import> importer;
668 if (Coarse->IsAvailable(
"Importer"))
669 importer = Coarse->Get< RCP<const Import> >(
"Importer");
672 coarseX = MultiVectorFactory::Build(coarseRhs->getMap(), X.getNumVectors(), initializeWithZeros);
675 RCP<TimeMonitor> ITime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : import (total)" ,
Timings0));
676 RCP<TimeMonitor> ILevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : import" + levelSuffix,
Timings0));
679 RCP<MultiVector> coarseTmp = MultiVectorFactory::Build(importer->getTargetMap(), coarseRhs->getNumVectors(),
false);
680 coarseTmp->doImport(*coarseRhs, *importer, Xpetra::INSERT);
681 coarseRhs.swap(coarseTmp);
683 coarseX = MultiVectorFactory::Build(importer->getTargetMap(), X.getNumVectors(), initializeWithZeros);
686 RCP<Operator> Ac = Coarse->Get< RCP<Operator> >(
"A");
688 RCP<const Map> origXMap = coarseX->getMap();
691 coarseRhs->replaceMap(Ac->getRangeMap());
692 coarseX ->replaceMap(Ac->getDomainMap());
695 iterateLevelTime = Teuchos::null;
697 Iterate(*coarseRhs, *coarseX, 1,
true, startLevel+1);
700 Iterate(*coarseRhs, *coarseX, 1,
false, startLevel+1);
703 iterateLevelTime = rcp(
new TimeMonitor(*
this, iterateLevelTimeLabel));
705 coarseX->replaceMap(origXMap);
709 RCP<TimeMonitor> ITime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : export (total)" ,
Timings0));
710 RCP<TimeMonitor> ILevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : export" + levelSuffix,
Timings0));
713 RCP<MultiVector> coarseTmp = MultiVectorFactory::Build(importer->getSourceMap(), coarseX->getNumVectors(),
false);
714 coarseTmp->doExport(*coarseX, *importer, Xpetra::INSERT);
715 coarseX.swap(coarseTmp);
722 RCP<MultiVector> correction = MultiVectorFactory::Build(P->getRangeMap(), X.getNumVectors(),
false);
725 RCP<TimeMonitor> PTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : prolongation (total)" ,
Timings0));
726 RCP<TimeMonitor> PLevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : prolongation" + levelSuffix,
Timings0));
727 P->apply(*coarseX, *correction, Teuchos::NO_TRANS, one, zero);
729 X.update(one, *correction, one);
733 RCP<TimeMonitor> STime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : smoothing (total)" ,
Timings0));
734 RCP<TimeMonitor> SLevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : smoothing" + levelSuffix,
Timings0));
736 if (Fine->IsAvailable(
"PostSmoother")) {
737 RCP<SmootherBase> postSmoo = Fine->Get< RCP<SmootherBase> >(
"PostSmoother");
738 postSmoo->Apply(X, B,
false);
751 Teuchos::Array<MagnitudeType> rn;
756 rate_ = as<MagnitudeType>(curNorm / prevNorm);
760 for (LO k = 0; k < rn.size(); k++)
770 << std::setiosflags(std::ios::left)
771 << std::setprecision(3) << i
773 << std::setprecision(10) << rn
780 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
782 LO startLevel = (start != -1 ? start : 0);
783 LO endLevel = (end != -1 ? end :
Levels_.size()-1);
786 "MueLu::Hierarchy::Write : startLevel must be <= endLevel");
789 "MueLu::Hierarchy::Write bad start or end level");
791 for (LO i = startLevel; i < endLevel + 1; i++) {
792 RCP<Matrix> A = rcp_dynamic_cast<Matrix>(
Levels_[i]->
template Get< RCP< Operator> >(
"A")), P, R;
794 P = rcp_dynamic_cast<Matrix>(
Levels_[i]->
template Get< RCP< Operator> >(
"P"));
796 R = rcp_dynamic_cast<Matrix>(
Levels_[i]->
template Get< RCP< Operator> >(
"R"));
799 if (!A.is_null()) Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(
"A_" +
toString(i) +
".m", *A);
800 if (!P.is_null()) Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(
"P_" +
toString(i) +
".m", *P);
801 if (!R.is_null()) Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(
"R_" +
toString(i) +
".m", *R);
805 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
807 for (Array<RCP<Level> >::iterator it =
Levels_.begin(); it !=
Levels_.end(); ++it)
808 (*it)->Keep(ename, factory);
811 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
813 for (Array<RCP<Level> >::iterator it =
Levels_.begin(); it !=
Levels_.end(); ++it)
814 (*it)->Delete(ename, factory);
817 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
819 for (Array<RCP<Level> >::iterator it =
Levels_.begin(); it !=
Levels_.end(); ++it)
820 (*it)->AddKeepFlag(ename, factory, keep);
823 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
825 for (Array<RCP<Level> >::iterator it =
Levels_.begin(); it !=
Levels_.end(); ++it)
826 (*it)->RemoveKeepFlag(ename, factory, keep);
829 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
831 std::ostringstream out;
837 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
842 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
844 RCP<Operator> A0 =
Levels_[0]->template Get<RCP<Operator> >(
"A");
845 RCP<const Teuchos::Comm<int> > comm = A0->getDomainMap()->getComm();
848 RCP<Operator> Ac =
Levels_[numLevels-1]->template Get<RCP<Operator> >(
"A");
855 int root = comm->getRank();
858 int smartData = numLevels*comm->getSize() + comm->getRank(), maxSmartData;
859 reduceAll(*comm, Teuchos::REDUCE_MAX, smartData, Teuchos::ptr(&maxSmartData));
860 root = maxSmartData % comm->getSize();
865 std::vector<Xpetra::global_size_t> nnzPerLevel;
866 std::vector<Xpetra::global_size_t> rowsPerLevel;
867 std::vector<int> numProcsPerLevel;
868 bool aborted =
false;
869 for (
int i = 0; i < numLevels; i++) {
871 "Operator A is not available on level " << i);
873 RCP<Operator> A =
Levels_[i]->template Get<RCP<Operator> >(
"A");
875 "Operator A on level " << i <<
" is null.");
877 RCP<Matrix> Am = rcp_dynamic_cast<Matrix>(A);
879 GetOStream(
Warnings0) <<
"Some level operators are not matrices, statistics calculation aborted" << std::endl;
884 Xpetra::global_size_t nnz = Am->getGlobalNumEntries();
885 nnzPerLevel .push_back(nnz);
886 rowsPerLevel .push_back(Am->getGlobalNumRows());
887 numProcsPerLevel.push_back(Am->getRowMap()->getComm()->getSize());
891 std::ostringstream oss;
892 oss <<
"\n--------------------------------------------------------------------------------\n" <<
893 "--- Multigrid Summary ---\n" 894 "--------------------------------------------------------------------------------" << std::endl;
895 oss <<
"Number of levels = " << numLevels << std::endl;
896 oss <<
"Operator complexity = " << std::setprecision(2) << std::setiosflags(std::ios::fixed)
900 Xpetra::global_size_t tt = rowsPerLevel[0];
901 int rowspacer = 2;
while (tt != 0) { tt /= 10; rowspacer++; }
903 int nnzspacer = 2;
while (tt != 0) { tt /= 10; nnzspacer++; }
904 tt = numProcsPerLevel[0];
905 int npspacer = 2;
while (tt != 0) { tt /= 10; npspacer++; }
906 oss <<
"level " << std::setw(rowspacer) <<
" rows " << std::setw(nnzspacer) <<
" nnz " <<
" nnz/row" << std::setw(npspacer) <<
" c ratio" <<
" procs" << std::endl;
907 for (
size_t i = 0; i < nnzPerLevel.size(); ++i) {
908 oss <<
" " << i <<
" ";
909 oss << std::setw(rowspacer) << rowsPerLevel[i];
910 oss << std::setw(nnzspacer) << nnzPerLevel[i];
911 oss << std::setprecision(2) << std::setiosflags(std::ios::fixed);
912 oss << std::setw(9) << as<double>(nnzPerLevel[i]) / rowsPerLevel[i];
913 if (i) oss << std::setw(9) << as<double>(rowsPerLevel[i-1])/rowsPerLevel[i];
914 else oss << std::setw(9) <<
" ";
915 oss <<
" " << std::setw(npspacer) << numProcsPerLevel[i] << std::endl;
919 RCP<SmootherBase> preSmoo, postSmoo;
920 if (
Levels_[i]->IsAvailable(
"PreSmoother"))
921 preSmoo =
Levels_[i]->
template Get< RCP<SmootherBase> >(
"PreSmoother");
922 if (
Levels_[i]->IsAvailable(
"PostSmoother"))
923 postSmoo =
Levels_[i]->
template Get< RCP<SmootherBase> >(
"PostSmoother");
925 if (preSmoo != null && preSmoo == postSmoo)
926 oss <<
"Smoother (level " << i <<
") both : " << preSmoo->description() << std::endl;
928 oss <<
"Smoother (level " << i <<
") pre : " 929 << (preSmoo != null ? preSmoo->description() :
"no smoother") << std::endl;
930 oss <<
"Smoother (level " << i <<
") post : " 931 << (postSmoo != null ? postSmoo->description() :
"no smoother") << std::endl;
942 RCP<const Teuchos::MpiComm<int> > mpiComm = rcp_dynamic_cast<
const Teuchos::MpiComm<int> >(comm);
943 MPI_Comm rawComm = (*mpiComm->getRawMpiComm())();
945 int strLength = outstr.size();
946 MPI_Bcast(&strLength, 1, MPI_INT, root, rawComm);
947 if (comm->getRank() != root)
948 outstr.resize(strLength);
949 MPI_Bcast(&outstr[0], strLength, MPI_CHAR, root, rawComm);
956 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
958 Teuchos::OSTab tab2(out);
963 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
968 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
972 #if defined(HAVE_MUELU_BOOST) && defined(HAVE_MUELU_BOOST_FOR_REAL) && defined(BOOST_VERSION) && (BOOST_VERSION >= 104400) 976 dp.property(
"label", boost::get(boost::vertex_name, graph));
977 dp.property(
"id", boost::get(boost::vertex_index, graph));
978 dp.property(
"label", boost::get(boost::edge_name, graph));
979 dp.property(
"color", boost::get(boost::edge_color, graph));
982 std::map<const FactoryBase*, BoostVertex> vindices;
983 typedef std::map<std::pair<BoostVertex,BoostVertex>, std::string> emap; emap edges;
987 Levels_[i]->UpdateGraph(vindices, edges, dp, graph);
989 for (emap::const_iterator eit = edges.begin(); eit != edges.end(); eit++) {
990 std::pair<BoostEdge, bool> boost_edge = boost::add_edge(eit->first.first, eit->first.second, graph);
991 boost::put(
"label", dp, boost_edge.first, eit->second);
993 boost::put(
"color", dp, boost_edge.first, std::string(
"red"));
995 boost::put(
"color", dp, boost_edge.first, std::string(
"blue"));
1000 std::ostringstream legend;
1001 legend <<
"< <TABLE BORDER=\"0\" CELLBORDER=\"1\" CELLSPACING=\"0\" CELLPADDING=\"4\"> \ 1002 <TR><TD COLSPAN=\"2\">Legend</TD></TR> \ 1003 <TR><TD><FONT color=\"red\">Level " <<
dumpLevel_ <<
"</FONT></TD><TD><FONT color=\"blue\">Level " <<
dumpLevel_+1 <<
"</FONT></TD></TR> \ 1005 BoostVertex boost_vertex = boost::add_vertex(graph);
1006 boost::put(
"label", dp, boost_vertex, legend.str());
1009 boost::write_graphviz_dp(out, graph, dp, std::string(
"id"));
1011 GetOStream(
Errors) <<
"Dependency graph output requires boost and MueLu_ENABLE_Boost_for_real" << std::endl;
1016 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1018 RCP<Operator> Ao = level.
Get<RCP<Operator> >(
"A");
1019 RCP<Matrix> A = rcp_dynamic_cast<Matrix>(Ao);
1021 GetOStream(
Warnings0) <<
"Hierarchy::ReplaceCoordinateMap: operator is not a matrix, skipping..." << std::endl;
1025 typedef Xpetra::MultiVector<double,LO,GO,NO> xdMV;
1027 RCP<xdMV> coords = level.
Get<RCP<xdMV> >(
"Coordinates");
1029 if (A->getRowMap()->isSameAs(*(coords->getMap()))) {
1030 GetOStream(
Warnings0) <<
"Hierarchy::ReplaceCoordinateMap: matrix and coordinates maps are same, skipping..." << std::endl;
1034 if (A->IsView(
"stridedMaps") && rcp_dynamic_cast<
const StridedMap>(A->getRowMap(
"stridedMaps")) != Teuchos::null) {
1035 RCP<const StridedMap> stridedRowMap = rcp_dynamic_cast<
const StridedMap>(A->getRowMap(
"stridedMaps"));
1038 TEUCHOS_TEST_FOR_EXCEPTION(stridedRowMap->getStridedBlockId() != -1 || stridedRowMap->getOffset() != 0,
1039 Exceptions::RuntimeError,
"Hierarchy::ReplaceCoordinateMap: nontrivial maps (block id = " << stridedRowMap->getStridedBlockId()
1040 <<
", offset = " << stridedRowMap->getOffset() <<
")");
1045 size_t blkSize = A->GetFixedBlockSize();
1047 RCP<const Map> nodeMap = A->getRowMap();
1050 RCP<const Map> dofMap = A->getRowMap();
1051 GO indexBase = dofMap->getIndexBase();
1052 size_t numLocalDOFs = dofMap->getNodeNumElements();
1054 "Hierarchy::ReplaceCoordinateMap: block size (" << blkSize <<
") is incompatible with the number of local dofs in a row map (" << numLocalDOFs);
1055 ArrayView<const GO> GIDs = dofMap->getNodeElementList();
1057 Array<GO> nodeGIDs(numLocalDOFs/blkSize);
1058 for (
size_t i = 0; i < numLocalDOFs; i += blkSize)
1059 nodeGIDs[i/blkSize] = (GIDs[i] - indexBase)/blkSize + indexBase;
1061 Xpetra::global_size_t INVALID = Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid();
1062 nodeMap = MapFactory::Build(dofMap->lib(), INVALID, nodeGIDs(), indexBase, dofMap->getComm());
1065 Array<ArrayView<const double> > coordDataView;
1066 std::vector<ArrayRCP<const double> > coordData;
1067 for (
size_t i = 0; i < coords->getNumVectors(); i++) {
1068 coordData.push_back(coords->getData(i));
1069 coordDataView.push_back(coordData[i]());
1072 RCP<xdMV> newCoords = Xpetra::MultiVectorFactory<double,LO,GO,NO>::Build(nodeMap, coordDataView(), coords->getNumVectors());
1073 level.
Set(
"Coordinates", newCoords);
1078 #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 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.
static Xpetra::global_size_t GetDefaultMaxCoarseSize()
virtual void Clean() const
void DumpCurrentGraph() const
void Clear(int startLevel=0)
Clear impermanent data from previous setup.
Namespace for MueLu classes and methods.
void Write(const LO &start=-1, const LO &end=-1)
Print matrices in the multigrid hierarchy to file.
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.
MagnitudeType rate_
Convergece rate.
virtual std::string ShortClassName() const
Return the class name of the object, without template parameters and without namespace.
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 SetProcRankVerbose(int procRank) const
Set proc rank used for printing.
void print(std::ostream &out=std::cout, const VerbLevel verbLevel=(MueLu::Parameters|MueLu::Statistics0)) const
VerbLevel GetVerbLevel() const
Get the verbosity level.
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)
static CycleType GetDefaultCycle()
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()
static bool GetDefaultImplicitTranspose()
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
Teuchos::FancyOStream & GetOStream(MsgType type, int thisProcRankOnly=0) const
Get an output stream for outputting the input message type.
bool IsPrint(MsgType type, int thisProcRankOnly=-1) const
Find out whether we need to print out information for a specific message type.
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.
Array< RCP< Level > > Levels_
Container for Level objects.
static bool GetDefaultPRrebalance()
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.
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.
bool isDumpingEnabled_
Graph dumping.
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.
Array< RCP< const FactoryManagerBase > > levelManagers_
Xpetra::global_size_t maxCoarseSize_
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
int GetProcRankVerbose() const
Get proc rank used for printing. Do not use this information for any other purpose.
void Request(const FactoryBase &factory)
Increment the storage counter for all the inputs of a factory.