47 #ifndef MUELU_HIERARCHY_DEF_HPP 48 #define MUELU_HIERARCHY_DEF_HPP 55 #include <Xpetra_Matrix.hpp> 56 #include <Xpetra_MultiVectorFactory.hpp> 57 #include <Xpetra_Operator.hpp> 58 #include <Xpetra_IO.hpp> 63 #include "MueLu_FactoryManager.hpp" 64 #include "MueLu_HierarchyUtils.hpp" 65 #include "MueLu_TopRAPFactory.hpp" 66 #include "MueLu_TopSmootherFactory.hpp" 69 #include "MueLu_PerfUtils.hpp" 70 #include "MueLu_PFactory.hpp" 72 #include "MueLu_SmootherFactory.hpp" 75 #include "Teuchos_TimeMonitor.hpp" 81 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
83 : maxCoarseSize_(GetDefaultMaxCoarseSize()), implicitTranspose_(GetDefaultImplicitTranspose()),
84 fuseProlongationAndUpdate_(GetDefaultFuseProlongationAndUpdate()),
85 doPRrebalance_(GetDefaultPRrebalance()), isPreconditioner_(true), Cycle_(GetDefaultCycle()), WCycleStartLevel_(0),
86 scalingFactor_(
Teuchos::ScalarTraits<double>::one()), lib_(Xpetra::UseTpetra), isDumpingEnabled_(false), dumpLevel_(-1), rate_(-1),
87 sizeOfAllocatedLevelMultiVectors_(0)
92 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
96 setObjectLabel(label);
97 Levels_[0]->setObjectLabel(label);
100 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
102 : maxCoarseSize_(GetDefaultMaxCoarseSize()), implicitTranspose_(GetDefaultImplicitTranspose()),
103 fuseProlongationAndUpdate_(GetDefaultFuseProlongationAndUpdate()),
104 doPRrebalance_(GetDefaultPRrebalance()), isPreconditioner_(true), Cycle_(GetDefaultCycle()), WCycleStartLevel_(0),
105 scalingFactor_(
Teuchos::ScalarTraits<double>::one()), isDumpingEnabled_(false), dumpLevel_(-1), rate_(-1),
106 sizeOfAllocatedLevelMultiVectors_(0)
108 lib_ = A->getDomainMap()->lib();
110 RCP<Level> Finest = rcp(
new Level);
116 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
120 setObjectLabel(label);
121 Levels_[0]->setObjectLabel(label);
124 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
126 int levelID = LastLevelID() + 1;
128 if (level->GetLevelID() != -1 && (level->GetLevelID() != levelID))
129 GetOStream(
Warnings1) <<
"Hierarchy::AddLevel(): Level with ID=" << level->GetLevelID() <<
130 " have been added at the end of the hierarchy\n but its ID have been redefined" <<
131 " because last level ID of the hierarchy was " << LastLevelID() <<
"." << std::endl;
133 Levels_.push_back(level);
134 level->SetLevelID(levelID);
137 level->SetPreviousLevel( (levelID == 0) ? Teuchos::null : Levels_[LastLevelID() - 1] );
138 level->setObjectLabel(this->getObjectLabel());
141 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
143 RCP<Level> newLevel = Levels_[LastLevelID()]->Build();
144 newLevel->setlib(lib_);
145 this->AddLevel(newLevel);
148 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
151 "MueLu::Hierarchy::GetLevel(): invalid input parameter value: LevelID = " << levelID);
152 return Levels_[levelID];
155 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
157 return Levels_.size();
160 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
162 RCP<Operator> A = Levels_[0]->template Get<RCP<Operator> >(
"A");
163 RCP<const Teuchos::Comm<int> > comm = A->getDomainMap()->getComm();
165 int numLevels = GetNumLevels();
167 Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, numLevels, Teuchos::ptr(&numGlobalLevels));
169 return numGlobalLevels;
172 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
174 double totalNnz = 0, lev0Nnz = 1;
175 for (
int i = 0; i < GetNumLevels(); ++i) {
177 "Operator complexity cannot be calculated because A is unavailable on level " << i);
178 RCP<Operator> A = Levels_[i]->template Get<RCP<Operator> >(
"A");
182 RCP<Matrix> Am = rcp_dynamic_cast<Matrix>(A);
184 GetOStream(
Warnings0) <<
"Some level operators are not matrices, operator complexity calculation aborted" << std::endl;
188 totalNnz += as<double>(Am->getGlobalNumEntries());
192 return totalNnz / lev0Nnz;
195 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
197 double node_sc = 0, global_sc=0;
199 const size_t INVALID = Teuchos::OrdinalTraits<size_t>::invalid();
201 if (GetNumLevels() <= 0)
return -1.0;
202 if (!Levels_[0]->IsAvailable(
"A"))
return -1.0;
204 RCP<Operator> A = Levels_[0]->template Get<RCP<Operator> >(
"A");
205 if (A.is_null())
return -1.0;
206 RCP<Matrix> Am = rcp_dynamic_cast<Matrix>(A);
207 if(Am.is_null())
return -1.0;
208 a0_nnz = as<double>(Am->getGlobalNumEntries());
211 for (
int i = 0; i < GetNumLevels(); ++i) {
213 if(!Levels_[i]->IsAvailable(
"PreSmoother"))
continue;
214 RCP<SmootherBase> S = Levels_[i]->template Get<RCP<SmootherBase> >(
"PreSmoother");
215 if (S.is_null())
continue;
216 level_sc = S->getNodeSmootherComplexity();
217 if(level_sc == INVALID) {global_sc=-1.0;
break;}
219 node_sc += as<double>(level_sc);
223 RCP<const Teuchos::Comm<int> > comm =A->getDomainMap()->getComm();
224 Teuchos::reduceAll(*comm,Teuchos::REDUCE_SUM,node_sc,Teuchos::ptr(&global_sc));
225 Teuchos::reduceAll(*comm,Teuchos::REDUCE_MIN,node_sc,Teuchos::ptr(&min_sc));
227 if(min_sc < 0.0)
return -1.0;
228 else return global_sc / a0_nnz;
235 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
238 "MueLu::Hierarchy::CheckLevel(): wrong underlying linear algebra library.");
240 "MueLu::Hierarchy::CheckLevel(): wrong level ID");
242 "MueLu::Hierarchy::Setup(): wrong level parent");
245 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
247 for (
int i = 0; i < GetNumLevels(); ++i) {
248 RCP<Level> level = Levels_[i];
249 if (level->IsAvailable(
"A")) {
250 RCP<Operator> Aop = level->Get<RCP<Operator> >(
"A");
251 RCP<Matrix> A = rcp_dynamic_cast<Matrix>(Aop);
253 RCP<const Import> xpImporter = A->getCrsGraph()->getImporter();
254 if (!xpImporter.is_null())
255 xpImporter->setDistributorParameters(matvecParams);
256 RCP<const Export> xpExporter = A->getCrsGraph()->getExporter();
257 if (!xpExporter.is_null())
258 xpExporter->setDistributorParameters(matvecParams);
261 if (level->IsAvailable(
"P")) {
262 RCP<Matrix> P = level->Get<RCP<Matrix> >(
"P");
263 RCP<const Import> xpImporter = P->getCrsGraph()->getImporter();
264 if (!xpImporter.is_null())
265 xpImporter->setDistributorParameters(matvecParams);
266 RCP<const Export> xpExporter = P->getCrsGraph()->getExporter();
267 if (!xpExporter.is_null())
268 xpExporter->setDistributorParameters(matvecParams);
270 if (level->IsAvailable(
"R")) {
271 RCP<Matrix> R = level->Get<RCP<Matrix> >(
"R");
272 RCP<const Import> xpImporter = R->getCrsGraph()->getImporter();
273 if (!xpImporter.is_null())
274 xpImporter->setDistributorParameters(matvecParams);
275 RCP<const Export> xpExporter = R->getCrsGraph()->getExporter();
276 if (!xpExporter.is_null())
277 xpExporter->setDistributorParameters(matvecParams);
279 if (level->IsAvailable(
"Importer")) {
280 RCP<const Import> xpImporter = level->Get< RCP<const Import> >(
"Importer");
281 if (!xpImporter.is_null())
282 xpImporter->setDistributorParameters(matvecParams);
289 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
291 const RCP<const FactoryManagerBase> fineLevelManager,
292 const RCP<const FactoryManagerBase> coarseLevelManager,
293 const RCP<const FactoryManagerBase> nextLevelManager) {
298 "MueLu::Hierarchy:Setup(): level " << coarseLevelID <<
" (specified by coarseLevelID argument) " 299 "must be built before calling this function.");
301 Level& level = *Levels_[coarseLevelID];
304 TimeMonitor m1(*
this, label + this->ShortClassName() +
": " +
"Setup (total)");
309 "MueLu::Hierarchy::Setup(): argument coarseLevelManager cannot be null");
314 if (levelManagers_.size() < coarseLevelID+1)
315 levelManagers_.resize(coarseLevelID+1);
316 levelManagers_[coarseLevelID] = coarseLevelManager;
318 bool isFinestLevel = (fineLevelManager.is_null());
319 bool isLastLevel = (nextLevelManager.is_null());
323 RCP<Operator> A = level.
Get< RCP<Operator> >(
"A");
324 RCP<const Map> domainMap = A->getDomainMap();
325 RCP<const Teuchos::Comm<int> > comm = domainMap->getComm();
332 oldRank = SetProcRankVerbose(comm->getRank());
336 lib_ = domainMap->lib();
343 Level& prevLevel = *Levels_[coarseLevelID-1];
344 oldRank = SetProcRankVerbose(prevLevel.
GetComm()->getRank());
347 CheckLevel(level, coarseLevelID);
350 RCP<SetFactoryManager> SFMFine;
352 SFMFine = rcp(
new SetFactoryManager(Levels_[coarseLevelID-1], fineLevelManager));
354 if (isFinestLevel && Levels_[coarseLevelID]->IsAvailable(
"Coordinates"))
355 ReplaceCoordinateMap(*Levels_[coarseLevelID]);
360 if (isDumpingEnabled_ && dumpLevel_ == 0 && coarseLevelID == 1)
363 RCP<TopSmootherFactory> coarseFact = rcp(
new TopSmootherFactory(coarseLevelManager,
"CoarseSolver"));
364 RCP<TopSmootherFactory> smootherFact = rcp(
new TopSmootherFactory(coarseLevelManager,
"Smoother"));
366 int nextLevelID = coarseLevelID + 1;
368 RCP<SetFactoryManager> SFMNext;
369 if (isLastLevel ==
false) {
371 if (nextLevelID > LastLevelID())
373 CheckLevel(*Levels_[nextLevelID], nextLevelID);
377 Levels_[nextLevelID]->Request(
TopRAPFactory(coarseLevelManager, nextLevelManager));
412 RCP<Operator> Ac = Teuchos::null;
413 TopRAPFactory coarseRAPFactory(fineLevelManager, coarseLevelManager);
416 Ac = level.
Get<RCP<Operator> >(
"A");
417 }
else if (!isFinestLevel) {
423 Ac = level.
Get<RCP<Operator> >(
"A");
424 RCP<Matrix> Acm = rcp_dynamic_cast<Matrix>(Ac);
428 level.
SetComm(Ac->getDomainMap()->getComm());
431 bool isOrigLastLevel = isLastLevel;
436 }
else if (Ac.is_null()) {
443 if (!Acm.is_null() && Acm->getGlobalNumRows() <= maxCoarseSize_) {
445 GetOStream(
Runtime0) <<
"Max coarse size (<= " << maxCoarseSize_ <<
") achieved" << std::endl;
450 if (!Ac.is_null() && !isFinestLevel) {
451 RCP<Operator> A = Levels_[coarseLevelID-1]->template Get< RCP<Operator> >(
"A");
452 RCP<Matrix> Am = rcp_dynamic_cast<Matrix>(A);
454 const double maxCoarse2FineRatio = 0.8;
455 if (!Acm.is_null() && !Am.is_null() && Acm->getGlobalNumRows() > maxCoarse2FineRatio * Am->getGlobalNumRows()) {
463 GetOStream(
Warnings0) <<
"Aggregation stagnated. Please check your matrix and/or adjust your configuration file." 464 <<
"Possible fixes:\n" 465 <<
" - reduce the maximum number of levels\n" 466 <<
" - enable repartitioning\n" 467 <<
" - increase the minimum coarse size." << std::endl;
473 if (!isOrigLastLevel) {
483 coarseFact->Build(level);
494 smootherFact->Build(level);
499 if (isLastLevel ==
true) {
500 if (isOrigLastLevel ==
false) {
503 Levels_[nextLevelID]->Release(
TopRAPFactory(coarseLevelManager, nextLevelManager));
505 Levels_.resize(nextLevelID);
509 if (isDumpingEnabled_ && dumpLevel_ > 0 && coarseLevelID == dumpLevel_)
512 if (!isFinestLevel) {
516 level.
Release(coarseRAPFactory);
520 SetProcRankVerbose(oldRank);
525 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
527 int numLevels = Levels_.size();
529 "Hierarchy::SetupRe: " << Levels_.size() <<
" levels, but " << levelManagers_.size() <<
" level factory managers");
531 const int startLevel = 0;
534 #ifdef HAVE_MUELU_DEBUG 536 for (
int i = 0; i < numLevels; i++)
537 levelManagers_[i]->ResetDebugData();
542 for (levelID = startLevel; levelID < numLevels;) {
543 bool r = Setup(levelID,
544 (levelID != 0 ? levelManagers_[levelID-1] : Teuchos::null),
545 levelManagers_[levelID],
546 (levelID+1 != numLevels ? levelManagers_[levelID+1] : Teuchos::null));
552 Levels_ .resize(levelID);
553 levelManagers_.resize(levelID);
555 int sizeOfVecs = sizeOfAllocatedLevelMultiVectors_;
556 DeleteLevelMultiVectors();
557 AllocateLevelMultiVectors(sizeOfVecs);
565 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
574 "MueLu::Hierarchy::Setup(): fine level (" << startLevel <<
") does not exist");
577 "Constructing non-positive (" << numDesiredLevels <<
") number of levels does not make sense.");
581 "MueLu::Hierarchy::Setup(): fine level (" << startLevel <<
") has no matrix A! " 582 "Set fine level matrix A using Level.Set()");
584 RCP<Operator> A = Levels_[startLevel]->template Get<RCP<Operator> >(
"A");
585 lib_ = A->getDomainMap()->lib();
588 RCP<Matrix> Amat = rcp_dynamic_cast<Matrix>(A);
590 if (!Amat.is_null()) {
591 RCP<ParameterList> params = rcp(
new ParameterList());
592 params->set(
"printLoadBalancingInfo",
true);
593 params->set(
"printCommInfo",
true);
597 GetOStream(
Warnings1) <<
"Fine level operator is not a matrix, statistics are not available" << std::endl;
601 RCP<const FactoryManagerBase> rcpmanager = rcpFromRef(manager);
603 const int lastLevel = startLevel + numDesiredLevels - 1;
604 GetOStream(
Runtime0) <<
"Setup loop: startLevel = " << startLevel <<
", lastLevel = " << lastLevel
605 <<
" (stop if numLevels = " << numDesiredLevels <<
" or Ac.size() < " << maxCoarseSize_ <<
")" << std::endl;
609 if (numDesiredLevels == 1) {
611 Setup(startLevel, Teuchos::null, rcpmanager, Teuchos::null);
614 bool bIsLastLevel = Setup(startLevel, Teuchos::null, rcpmanager, rcpmanager);
615 if (bIsLastLevel ==
false) {
616 for (iLevel = startLevel + 1; iLevel < lastLevel; iLevel++) {
617 bIsLastLevel = Setup(iLevel, rcpmanager, rcpmanager, rcpmanager);
618 if (bIsLastLevel ==
true)
621 if (bIsLastLevel ==
false)
622 Setup(lastLevel, rcpmanager, rcpmanager, Teuchos::null);
628 "MueLu::Hierarchy::Setup(): number of level");
637 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
639 if (startLevel < GetNumLevels())
640 GetOStream(
Runtime0) <<
"Clearing old data (if any)" << std::endl;
642 for (
int iLevel = startLevel; iLevel < GetNumLevels(); iLevel++)
643 Levels_[iLevel]->Clear();
646 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
648 GetOStream(
Runtime0) <<
"Clearing old data (expert)" << std::endl;
649 for (
int iLevel = 0; iLevel < GetNumLevels(); iLevel++)
650 Levels_[iLevel]->ExpertClear();
653 #if defined(HAVE_MUELU_EXPERIMENTAL) && defined(HAVE_MUELU_ADDITIVE_VARIANT) 654 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
656 bool InitialGuessIsZero, LO startLevel) {
657 LO nIts = conv.maxIts_;
658 MagnitudeType tol = conv.tol_;
660 std::string prefix = this->ShortClassName() +
": ";
661 std::string levelSuffix =
" (level=" +
toString(startLevel) +
")";
662 std::string levelSuffix1 =
" (level=" +
toString(startLevel+1) +
")";
665 RCP<Time> CompTime = Teuchos::TimeMonitor::getNewCounter(prefix +
"Computational Time (total)");
666 RCP<Time> Concurrent = Teuchos::TimeMonitor::getNewCounter(prefix +
"Concurrent portion");
667 RCP<Time> ApplyR = Teuchos::TimeMonitor::getNewCounter(prefix +
"R: Computational Time");
668 RCP<Time> ApplyPbar = Teuchos::TimeMonitor::getNewCounter(prefix +
"Pbar: Computational Time");
669 RCP<Time> CompFine = Teuchos::TimeMonitor::getNewCounter(prefix +
"Fine: Computational Time");
670 RCP<Time> CompCoarse = Teuchos::TimeMonitor::getNewCounter(prefix +
"Coarse: Computational Time");
671 RCP<Time> ApplySum = Teuchos::TimeMonitor::getNewCounter(prefix +
"Sum: Computational Time");
672 RCP<Time> Synchronize_beginning = Teuchos::TimeMonitor::getNewCounter(prefix +
"Synchronize_beginning");
673 RCP<Time> Synchronize_center = Teuchos::TimeMonitor::getNewCounter(prefix +
"Synchronize_center");
674 RCP<Time> Synchronize_end = Teuchos::TimeMonitor::getNewCounter(prefix +
"Synchronize_end");
676 RCP<Level> Fine = Levels_[0];
679 RCP<Operator> A = Fine->Get< RCP<Operator> >(
"A");
680 Teuchos::RCP< const Teuchos::Comm< int > > communicator = A->getDomainMap()->getComm();
689 SC one = STS::one(), zero = STS::zero();
691 bool zeroGuess = InitialGuessIsZero;
697 RCP< Operator > Pbar;
699 RCP< MultiVector > coarseRhs, coarseX;
701 RCP<SmootherBase> preSmoo_coarse, postSmoo_coarse;
702 bool emptyCoarseSolve =
true;
703 RCP<MultiVector> coarseX_prolonged = MultiVectorFactory::Build(X.getMap(), X.getNumVectors(),
true);
705 RCP<const Import> importer;
707 if (Levels_.size() > 1) {
709 if (Coarse->IsAvailable(
"Importer"))
710 importer = Coarse->Get< RCP<const Import> >(
"Importer");
712 R = Coarse->Get< RCP<Operator> >(
"R");
713 P = Coarse->Get< RCP<Operator> >(
"P");
717 Pbar = Coarse->Get< RCP<Operator> >(
"Pbar");
719 coarseRhs = MultiVectorFactory::Build(R->getRangeMap(), B.getNumVectors(),
true);
721 Ac = Coarse->Get< RCP< Operator > >(
"A");
724 R->apply(B, *coarseRhs, Teuchos::NO_TRANS, one, zero);
728 if (doPRrebalance_ || importer.is_null()) {
729 coarseX = MultiVectorFactory::Build(coarseRhs->getMap(), X.getNumVectors(),
true);
733 RCP<TimeMonitor> ITime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : import (total)" ,
Timings0));
734 RCP<TimeMonitor> ILevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : import" + levelSuffix1,
Timings0));
737 RCP<MultiVector> coarseTmp = MultiVectorFactory::Build(importer->getTargetMap(), coarseRhs->getNumVectors());
738 coarseTmp->doImport(*coarseRhs, *importer, Xpetra::INSERT);
739 coarseRhs.swap(coarseTmp);
741 coarseX = MultiVectorFactory::Build(importer->getTargetMap(), X.getNumVectors(),
true);
744 if (Coarse->IsAvailable(
"PreSmoother"))
745 preSmoo_coarse = Coarse->Get< RCP<SmootherBase> >(
"PreSmoother");
746 if (Coarse->IsAvailable(
"PostSmoother"))
747 postSmoo_coarse = Coarse->Get< RCP<SmootherBase> >(
"PostSmoother");
753 MagnitudeType prevNorm = STS::magnitude(STS::one()), curNorm = STS::magnitude(STS::one());
756 for (LO i = 1; i <= nIts; i++) {
757 #ifdef HAVE_MUELU_DEBUG 758 if (A->getDomainMap()->isCompatible(*(X.getMap())) ==
false) {
759 std::ostringstream ss;
760 ss <<
"Level " << startLevel <<
": level A's domain map is not compatible with X";
761 throw Exceptions::Incompatible(ss.str());
764 if (A->getRangeMap()->isCompatible(*(B.getMap())) ==
false) {
765 std::ostringstream ss;
766 ss <<
"Level " << startLevel <<
": level A's range map is not compatible with B";
767 throw Exceptions::Incompatible(ss.str());
772 bool emptyFineSolve =
true;
774 RCP< MultiVector > fineX;
775 fineX = MultiVectorFactory::Build(X.getMap(), X.getNumVectors(),
true);
784 if (Fine->IsAvailable(
"PreSmoother")) {
785 RCP<SmootherBase> preSmoo = Fine->Get< RCP<SmootherBase> >(
"PreSmoother");
787 preSmoo->Apply(*fineX, B, zeroGuess);
789 emptyFineSolve =
false;
791 if (Fine->IsAvailable(
"PostSmoother")) {
792 RCP<SmootherBase> postSmoo = Fine->Get< RCP<SmootherBase> >(
"PostSmoother");
794 postSmoo->Apply(*fineX, B, zeroGuess);
797 emptyFineSolve =
false;
799 if (emptyFineSolve ==
true) {
800 GetOStream(
Warnings1) <<
"No fine grid smoother" << std::endl;
802 fineX->update(one, B, zero);
805 if (Levels_.size() > 1) {
807 if (Coarse->IsAvailable(
"PreSmoother")) {
809 preSmoo_coarse->Apply(*coarseX, *coarseRhs, zeroGuess);
811 emptyCoarseSolve =
false;
814 if (Coarse->IsAvailable(
"PostSmoother")) {
816 postSmoo_coarse->Apply(*coarseX, *coarseRhs, zeroGuess);
818 emptyCoarseSolve =
false;
821 if (emptyCoarseSolve ==
true) {
822 GetOStream(
Warnings1) <<
"No coarse grid solver" << std::endl;
824 coarseX->update(one, *coarseRhs, zero);
831 if (!doPRrebalance_ && !importer.is_null()) {
832 RCP<TimeMonitor> ITime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : export (total)" ,
Timings0));
833 RCP<TimeMonitor> ILevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : export" + levelSuffix1,
Timings0));
836 RCP<MultiVector> coarseTmp = MultiVectorFactory::Build(importer->getSourceMap(), coarseX->getNumVectors());
837 coarseTmp->doExport(*coarseX, *importer, Xpetra::INSERT);
838 coarseX.swap(coarseTmp);
842 Pbar->apply(*coarseX, *coarseX_prolonged, Teuchos::NO_TRANS, one, zero);
847 X.update(1.0, *fineX, 1.0, *coarseX_prolonged, 0.0);
858 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
860 bool InitialGuessIsZero, LO startLevel) {
876 RCP<Level> Fine = Levels_[startLevel];
879 std::string prefix = label + this->ShortClassName() +
": ";
880 std::string levelSuffix =
" (level=" +
toString(startLevel) +
")";
881 std::string levelSuffix1 =
" (level=" +
toString(startLevel+1) +
")";
883 bool useStackedTimer = !Teuchos::TimeMonitor::getStackedTimer().is_null();
885 RCP<Monitor> iterateTime;
886 RCP<TimeMonitor> iterateTime1;
889 else if (!useStackedTimer)
892 std::string iterateLevelTimeLabel = prefix +
"Solve" + levelSuffix;
893 RCP<TimeMonitor> iterateLevelTime = rcp(
new TimeMonitor(*
this, iterateLevelTimeLabel,
Timings0));
895 bool zeroGuess = InitialGuessIsZero;
897 RCP<Operator> A = Fine->Get< RCP<Operator> >(
"A");
899 RCP<Time> CompCoarse = Teuchos::TimeMonitor::getNewCounter(prefix +
"Coarse: Computational Time");
912 const BlockedMultiVector * Bblocked =
dynamic_cast<const BlockedMultiVector*
>(&B);
913 if(residual_.size() > startLevel &&
914 ( ( Bblocked && !Bblocked->isSameSize(*residual_[startLevel])) ||
915 (!Bblocked && !residual_[startLevel]->isSameSize(B))))
916 DeleteLevelMultiVectors();
917 AllocateLevelMultiVectors(X.getNumVectors());
920 typedef Teuchos::ScalarTraits<typename STS::magnitudeType> STM;
923 if (startLevel == 0 && !isPreconditioner_ &&
927 Teuchos::Array<MagnitudeType> rn;
932 for (LO k = 0; k < rn.size(); k++)
942 << std::setiosflags(std::ios::left)
943 << std::setprecision(3) << 0
945 << std::setprecision(10) << rn
949 SC one = STS::one(), zero = STS::zero();
950 for (LO i = 1; i <= nIts; i++) {
951 #ifdef HAVE_MUELU_DEBUG 953 if (A->getDomainMap()->isCompatible(*(X.getMap())) ==
false) {
954 std::ostringstream ss;
955 ss <<
"Level " << startLevel <<
": level A's domain map is not compatible with X";
959 if (A->getRangeMap()->isCompatible(*(B.getMap())) ==
false) {
960 std::ostringstream ss;
961 ss <<
"Level " << startLevel <<
": level A's range map is not compatible with B";
967 if (startLevel == as<LO>(Levels_.size())-1) {
969 RCP<TimeMonitor> CLevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : coarse" + levelSuffix,
Timings0));
971 bool emptySolve =
true;
974 if (Fine->IsAvailable(
"PreSmoother")) {
975 RCP<SmootherBase> preSmoo = Fine->Get< RCP<SmootherBase> >(
"PreSmoother");
977 preSmoo->Apply(X, B, zeroGuess);
982 if (Fine->IsAvailable(
"PostSmoother")) {
983 RCP<SmootherBase> postSmoo = Fine->Get< RCP<SmootherBase> >(
"PostSmoother");
985 postSmoo->Apply(X, B, zeroGuess);
990 if (emptySolve ==
true) {
991 GetOStream(
Warnings1) <<
"No coarse grid solver" << std::endl;
993 X.update(one, B, zero);
998 RCP<Level> Coarse = Levels_[startLevel+1];
1001 RCP<TimeMonitor> STime;
1002 if (!useStackedTimer)
1004 RCP<TimeMonitor> SLevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : smoothing" + levelSuffix,
Timings0));
1006 if (Fine->IsAvailable(
"PreSmoother")) {
1007 RCP<SmootherBase> preSmoo = Fine->Get< RCP<SmootherBase> >(
"PreSmoother");
1008 preSmoo->Apply(X, B, zeroGuess);
1013 RCP<MultiVector> residual;
1015 RCP<TimeMonitor> ATime;
1016 if (!useStackedTimer)
1017 ATime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : residual calculation (total)" ,
Timings0));
1018 RCP<TimeMonitor> ALevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : residual calculation" + levelSuffix,
Timings0));
1026 residual = residual_[startLevel];
1029 RCP<Operator> P = Coarse->Get< RCP<Operator> >(
"P");
1030 if (Coarse->IsAvailable(
"Pbar"))
1031 P = Coarse->Get< RCP<Operator> >(
"Pbar");
1033 RCP<MultiVector> coarseRhs, coarseX;
1037 RCP<TimeMonitor> RTime;
1038 if (!useStackedTimer)
1040 RCP<TimeMonitor> RLevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : restriction" + levelSuffix,
Timings0));
1041 coarseRhs = coarseRhs_[startLevel];
1043 if (implicitTranspose_) {
1044 P->apply(*residual, *coarseRhs, Teuchos::TRANS, one, zero);
1047 RCP<Operator> R = Coarse->Get< RCP<Operator> >(
"R");
1048 R->apply(*residual, *coarseRhs, Teuchos::NO_TRANS, one, zero);
1052 RCP<const Import> importer;
1053 if (Coarse->IsAvailable(
"Importer"))
1054 importer = Coarse->Get< RCP<const Import> >(
"Importer");
1056 coarseX = coarseX_[startLevel];
1057 if (!doPRrebalance_ && !importer.is_null()) {
1058 RCP<TimeMonitor> ITime;
1059 if (!useStackedTimer)
1061 RCP<TimeMonitor> ILevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : import" + levelSuffix1,
Timings0));
1064 RCP<MultiVector> coarseTmp = coarseImport_[startLevel];
1065 coarseTmp->doImport(*coarseRhs, *importer, Xpetra::INSERT);
1066 coarseRhs.swap(coarseTmp);
1069 RCP<Operator> Ac = Coarse->Get< RCP<Operator> >(
"A");
1070 if (!Ac.is_null()) {
1071 RCP<const Map> origXMap = coarseX->getMap();
1072 RCP<const Map> origRhsMap = coarseRhs->getMap();
1075 coarseRhs->replaceMap(Ac->getRangeMap());
1076 coarseX ->replaceMap(Ac->getDomainMap());
1079 iterateLevelTime = Teuchos::null;
1081 Iterate(*coarseRhs, *coarseX, 1,
true, startLevel+1);
1083 if (Cycle_ ==
WCYCLE && WCycleStartLevel_ >= startLevel)
1084 Iterate(*coarseRhs, *coarseX, 1,
false, startLevel+1);
1087 iterateLevelTime = rcp(
new TimeMonitor(*
this, iterateLevelTimeLabel));
1089 coarseX->replaceMap(origXMap);
1090 coarseRhs->replaceMap(origRhsMap);
1093 if (!doPRrebalance_ && !importer.is_null()) {
1094 RCP<TimeMonitor> ITime;
1095 if (!useStackedTimer)
1097 RCP<TimeMonitor> ILevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : export" + levelSuffix1,
Timings0));
1100 RCP<MultiVector> coarseTmp = coarseExport_[startLevel];
1101 coarseTmp->doExport(*coarseX, *importer, Xpetra::INSERT);
1102 coarseX.swap(coarseTmp);
1107 RCP<TimeMonitor> PTime;
1108 if (!useStackedTimer)
1110 RCP<TimeMonitor> PLevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : prolongation" + levelSuffix,
Timings0));
1115 if (fuseProlongationAndUpdate_) {
1116 P->apply(*coarseX, X, Teuchos::NO_TRANS, scalingFactor_, one);
1118 RCP<MultiVector> correction = correction_[startLevel];
1119 P->apply(*coarseX, *correction, Teuchos::NO_TRANS, one, zero);
1120 X.update(scalingFactor_, *correction, one);
1126 RCP<TimeMonitor> STime;
1127 if (!useStackedTimer)
1129 RCP<TimeMonitor> SLevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : smoothing" + levelSuffix,
Timings0));
1131 if (Fine->IsAvailable(
"PostSmoother")) {
1132 RCP<SmootherBase> postSmoo = Fine->Get< RCP<SmootherBase> >(
"PostSmoother");
1133 postSmoo->Apply(X, B,
false);
1139 if (startLevel == 0 && !isPreconditioner_ &&
1143 Teuchos::Array<MagnitudeType> rn;
1148 rate_ = as<MagnitudeType>(curNorm / prevNorm);
1152 << std::setiosflags(std::ios::left)
1153 << std::setprecision(3) << i
1155 << std::setprecision(10) << rn
1160 for (LO k = 0; k < rn.size(); k++)
1174 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1176 LO startLevel = (start != -1 ? start : 0);
1177 LO endLevel = (end != -1 ? end : Levels_.size()-1);
1180 "MueLu::Hierarchy::Write : startLevel must be <= endLevel");
1183 "MueLu::Hierarchy::Write bad start or end level");
1185 for (LO i = startLevel; i < endLevel + 1; i++) {
1186 RCP<Matrix> A = rcp_dynamic_cast<Matrix>(Levels_[i]->
template Get< RCP< Operator> >(
"A")), P, R;
1188 P = rcp_dynamic_cast<Matrix>(Levels_[i]->
template Get< RCP< Operator> >(
"P"));
1189 if (!implicitTranspose_)
1190 R = rcp_dynamic_cast<Matrix>(Levels_[i]->
template Get< RCP< Operator> >(
"R"));
1193 if (!A.is_null()) Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(
"A_" +
toString(i) + suffix +
".m", *A);
1195 Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(
"P_" +
toString(i) + suffix +
".m", *P);
1198 Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(
"R_" +
toString(i) + suffix +
".m", *R);
1203 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1205 for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1206 (*it)->Keep(ename, factory);
1209 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1211 for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1212 (*it)->Delete(ename, factory);
1215 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1217 for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1218 (*it)->AddKeepFlag(ename, factory, keep);
1221 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1223 for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1224 (*it)->RemoveKeepFlag(ename, factory, keep);
1227 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1229 if ( description_.empty() )
1231 std::ostringstream out;
1233 out <<
"{#levels = " << GetGlobalNumLevels() <<
", complexity = " << GetOperatorComplexity() <<
"}";
1234 description_ = out.str();
1236 return description_;
1239 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1244 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1246 RCP<Operator> A0 = Levels_[0]->template Get<RCP<Operator> >(
"A");
1247 RCP<const Teuchos::Comm<int> > comm = A0->getDomainMap()->getComm();
1249 int numLevels = GetNumLevels();
1250 RCP<Operator> Ac = Levels_[numLevels-1]->template Get<RCP<Operator> >(
"A");
1257 int root = comm->getRank();
1260 int smartData = numLevels*comm->getSize() + comm->getRank(), maxSmartData;
1261 reduceAll(*comm, Teuchos::REDUCE_MAX, smartData, Teuchos::ptr(&maxSmartData));
1262 root = maxSmartData % comm->getSize();
1266 double smoother_comp =-1.0;
1268 smoother_comp = GetSmootherComplexity();
1272 std::vector<Xpetra::global_size_t> nnzPerLevel;
1273 std::vector<Xpetra::global_size_t> rowsPerLevel;
1274 std::vector<int> numProcsPerLevel;
1275 bool aborted =
false;
1276 for (
int i = 0; i < numLevels; i++) {
1278 "Operator A is not available on level " << i);
1280 RCP<Operator> A = Levels_[i]->template Get<RCP<Operator> >(
"A");
1282 "Operator A on level " << i <<
" is null.");
1284 RCP<Matrix> Am = rcp_dynamic_cast<Matrix>(A);
1286 GetOStream(
Warnings0) <<
"Some level operators are not matrices, statistics calculation aborted" << std::endl;
1291 Xpetra::global_size_t nnz = Am->getGlobalNumEntries();
1292 nnzPerLevel .push_back(nnz);
1293 rowsPerLevel .push_back(Am->getGlobalNumRows());
1294 numProcsPerLevel.push_back(Am->getRowMap()->getComm()->getSize());
1298 std::string label = Levels_[0]->getObjectLabel();
1299 std::ostringstream oss;
1300 oss << std::setfill(
' ');
1301 oss <<
"\n--------------------------------------------------------------------------------\n";
1302 oss <<
"--- Multigrid Summary " << std::setw(28) << std::left << label <<
"---\n";
1303 oss <<
"--------------------------------------------------------------------------------" << std::endl;
1304 oss <<
"Number of levels = " << numLevels << std::endl;
1305 oss <<
"Operator complexity = " << std::setprecision(2) << std::setiosflags(std::ios::fixed)
1306 << GetOperatorComplexity() << std::endl;
1308 if(smoother_comp!=-1.0) {
1309 oss <<
"Smoother complexity = " << std::setprecision(2) << std::setiosflags(std::ios::fixed)
1310 << smoother_comp << std::endl;
1315 oss <<
"Cycle type = V" << std::endl;
1318 oss <<
"Cycle type = W" << std::endl;
1319 if (WCycleStartLevel_ > 0)
1320 oss <<
"Cycle start level = " << WCycleStartLevel_ << std::endl;
1327 Xpetra::global_size_t tt = rowsPerLevel[0];
1328 int rowspacer = 2;
while (tt != 0) { tt /= 10; rowspacer++; }
1329 tt = nnzPerLevel[0];
1330 int nnzspacer = 2;
while (tt != 0) { tt /= 10; nnzspacer++; }
1331 tt = numProcsPerLevel[0];
1332 int npspacer = 2;
while (tt != 0) { tt /= 10; npspacer++; }
1333 oss <<
"level " << std::setw(rowspacer) <<
" rows " << std::setw(nnzspacer) <<
" nnz " <<
" nnz/row" << std::setw(npspacer) <<
" c ratio" <<
" procs" << std::endl;
1334 for (
size_t i = 0; i < nnzPerLevel.size(); ++i) {
1335 oss <<
" " << i <<
" ";
1336 oss << std::setw(rowspacer) << rowsPerLevel[i];
1337 oss << std::setw(nnzspacer) << nnzPerLevel[i];
1338 oss << std::setprecision(2) << std::setiosflags(std::ios::fixed);
1339 oss << std::setw(9) << as<double>(nnzPerLevel[i]) / rowsPerLevel[i];
1340 if (i) oss << std::setw(9) << as<double>(rowsPerLevel[i-1])/rowsPerLevel[i];
1341 else oss << std::setw(9) <<
" ";
1342 oss <<
" " << std::setw(npspacer) << numProcsPerLevel[i] << std::endl;
1345 for (
int i = 0; i < GetNumLevels(); ++i) {
1346 RCP<SmootherBase> preSmoo, postSmoo;
1347 if (Levels_[i]->IsAvailable(
"PreSmoother"))
1348 preSmoo = Levels_[i]->
template Get< RCP<SmootherBase> >(
"PreSmoother");
1349 if (Levels_[i]->IsAvailable(
"PostSmoother"))
1350 postSmoo = Levels_[i]->
template Get< RCP<SmootherBase> >(
"PostSmoother");
1352 if (preSmoo != null && preSmoo == postSmoo)
1353 oss <<
"Smoother (level " << i <<
") both : " << preSmoo->description() << std::endl;
1355 oss <<
"Smoother (level " << i <<
") pre : " 1356 << (preSmoo != null ? preSmoo->description() :
"no smoother") << std::endl;
1357 oss <<
"Smoother (level " << i <<
") post : " 1358 << (postSmoo != null ? postSmoo->description() :
"no smoother") << std::endl;
1369 RCP<const Teuchos::MpiComm<int> > mpiComm = rcp_dynamic_cast<
const Teuchos::MpiComm<int> >(comm);
1370 MPI_Comm rawComm = (*mpiComm->getRawMpiComm())();
1372 int strLength = outstr.size();
1373 MPI_Bcast(&strLength, 1, MPI_INT, root, rawComm);
1374 if (comm->getRank() != root)
1375 outstr.resize(strLength);
1376 MPI_Bcast(&outstr[0], strLength, MPI_CHAR, root, rawComm);
1383 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1385 Teuchos::OSTab tab2(out);
1386 for (
int i = 0; i < GetNumLevels(); ++i)
1387 Levels_[i]->print(out, verbLevel);
1390 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1392 isPreconditioner_ = flag;
1395 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1397 if (GetProcRankVerbose() != 0)
1399 #if defined(HAVE_MUELU_BOOST) && defined(HAVE_MUELU_BOOST_FOR_REAL) && defined(BOOST_VERSION) && (BOOST_VERSION >= 104400) 1403 dp.property(
"label", boost::get(boost::vertex_name, graph));
1404 dp.property(
"id", boost::get(boost::vertex_index, graph));
1405 dp.property(
"label", boost::get(boost::edge_name, graph));
1406 dp.property(
"color", boost::get(boost::edge_color, graph));
1409 std::map<const FactoryBase*, BoostVertex> vindices;
1410 typedef std::map<std::pair<BoostVertex,BoostVertex>, std::string> emap; emap edges;
1412 static int call_id=0;
1414 RCP<Operator> A = Levels_[0]->template Get<RCP<Operator> >(
"A");
1415 int rank = A->getDomainMap()->getComm()->getRank();
1418 for (
int i = dumpLevel_; i <= dumpLevel_+1 && i < GetNumLevels(); i++) {
1420 Levels_[i]->UpdateGraph(vindices, edges, dp, graph);
1422 for (emap::const_iterator eit = edges.begin(); eit != edges.end(); eit++) {
1423 std::pair<BoostEdge, bool> boost_edge = boost::add_edge(eit->first.first, eit->first.second, graph);
1426 if(eit->second==std::string(
"Graph")) boost::put(
"label", dp, boost_edge.first, std::string(
"Graph_"));
1427 else boost::put(
"label", dp, boost_edge.first, eit->second);
1428 if (i == dumpLevel_)
1429 boost::put(
"color", dp, boost_edge.first, std::string(
"red"));
1431 boost::put(
"color", dp, boost_edge.first, std::string(
"blue"));
1436 std::ostringstream legend;
1437 legend <<
"< <TABLE BORDER=\"0\" CELLBORDER=\"1\" CELLSPACING=\"0\" CELLPADDING=\"4\"> \ 1438 <TR><TD COLSPAN=\"2\">Legend</TD></TR> \ 1439 <TR><TD><FONT color=\"red\">Level " << dumpLevel_ <<
"</FONT></TD><TD><FONT color=\"blue\">Level " << dumpLevel_+1 <<
"</FONT></TD></TR> \ 1441 BoostVertex boost_vertex = boost::add_vertex(graph);
1442 boost::put(
"label", dp, boost_vertex, legend.str());
1445 std::ofstream out(dumpFile_.c_str() +std::to_string(call_id)+std::string(
"_")+ std::to_string(rank) + std::string(
".dot"));
1446 boost::write_graphviz_dp(out, graph, dp, std::string(
"id"));
1450 GetOStream(
Errors) <<
"Dependency graph output requires boost and MueLu_ENABLE_Boost_for_real" << std::endl;
1455 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1457 RCP<Operator> Ao = level.
Get<RCP<Operator> >(
"A");
1458 RCP<Matrix> A = rcp_dynamic_cast<Matrix>(Ao);
1460 GetOStream(
Warnings1) <<
"Hierarchy::ReplaceCoordinateMap: operator is not a matrix, skipping..." << std::endl;
1463 if(Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A) != Teuchos::null) {
1464 GetOStream(
Warnings1) <<
"Hierarchy::ReplaceCoordinateMap: operator is a BlockedCrsMatrix, skipping..." << std::endl;
1468 typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO> xdMV;
1470 RCP<xdMV> coords = level.
Get<RCP<xdMV> >(
"Coordinates");
1472 if (A->getRowMap()->isSameAs(*(coords->getMap()))) {
1473 GetOStream(
Warnings1) <<
"Hierarchy::ReplaceCoordinateMap: matrix and coordinates maps are same, skipping..." << std::endl;
1477 if (A->IsView(
"stridedMaps") && rcp_dynamic_cast<
const StridedMap>(A->getRowMap(
"stridedMaps")) != Teuchos::null) {
1478 RCP<const StridedMap> stridedRowMap = rcp_dynamic_cast<
const StridedMap>(A->getRowMap(
"stridedMaps"));
1481 TEUCHOS_TEST_FOR_EXCEPTION(stridedRowMap->getStridedBlockId() != -1 || stridedRowMap->getOffset() != 0,
1482 Exceptions::RuntimeError,
"Hierarchy::ReplaceCoordinateMap: nontrivial maps (block id = " << stridedRowMap->getStridedBlockId()
1483 <<
", offset = " << stridedRowMap->getOffset() <<
")");
1486 GetOStream(
Runtime1) <<
"Replacing coordinate map" << std::endl;
1488 size_t blkSize = A->GetFixedBlockSize();
1490 RCP<const Map> nodeMap = A->getRowMap();
1493 RCP<const Map> dofMap = A->getRowMap();
1494 GO indexBase = dofMap->getIndexBase();
1495 size_t numLocalDOFs = dofMap->getNodeNumElements();
1497 "Hierarchy::ReplaceCoordinateMap: block size (" << blkSize <<
") is incompatible with the number of local dofs in a row map (" << numLocalDOFs);
1498 ArrayView<const GO> GIDs = dofMap->getNodeElementList();
1500 Array<GO> nodeGIDs(numLocalDOFs/blkSize);
1501 for (
size_t i = 0; i < numLocalDOFs; i += blkSize)
1502 nodeGIDs[i/blkSize] = (GIDs[i] - indexBase)/blkSize + indexBase;
1504 Xpetra::global_size_t INVALID = Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid();
1505 nodeMap = MapFactory::Build(dofMap->lib(), INVALID, nodeGIDs(), indexBase, dofMap->getComm());
1511 if(coords->getLocalLength() != A->getRowMap()->getNodeNumElements()) {
1512 GetOStream(
Warnings) <<
"Coordinate vector does not match row map of matrix A!" << std::endl;
1517 Array<ArrayView<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> > coordDataView;
1518 std::vector<ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> > coordData;
1519 for (
size_t i = 0; i < coords->getNumVectors(); i++) {
1520 coordData.push_back(coords->getData(i));
1521 coordDataView.push_back(coordData[i]());
1524 RCP<xdMV> newCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO>::Build(nodeMap, coordDataView(), coords->getNumVectors());
1525 level.
Set(
"Coordinates", newCoords);
1528 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1530 int N = Levels_.size();
1531 if( (sizeOfAllocatedLevelMultiVectors_ == numvecs && residual_.size() == N) || numvecs<=0 )
return;
1534 if(residual_.size() != N) {
1535 DeleteLevelMultiVectors();
1537 residual_.resize(N);
1538 coarseRhs_.resize(N);
1540 coarseImport_.resize(N);
1541 coarseExport_.resize(N);
1542 correction_.resize(N);
1545 for(
int i=0; i<N; i++) {
1546 RCP<Operator> A = Levels_[i]->template Get< RCP<Operator> >(
"A");
1549 RCP<const BlockedCrsMatrix> A_as_blocked = Teuchos::rcp_dynamic_cast<
const BlockedCrsMatrix>(A);
1550 RCP<const Map> Arm = A->getRangeMap();
1551 RCP<const Map> Adm = A->getDomainMap();
1552 if(!A_as_blocked.is_null()) {
1553 Adm = A_as_blocked->getFullDomainMap();
1557 residual_[i] = MultiVectorFactory::Build(Arm, numvecs,
true);
1558 correction_[i] = MultiVectorFactory::Build(Adm, numvecs,
false);
1563 if(implicitTranspose_) {
1564 RCP<Operator> P = Levels_[i+1]->template Get< RCP<Operator> >(
"P");
1565 if(!P.is_null()) coarseRhs_[i] = MultiVectorFactory::Build(P->getDomainMap(),numvecs,
true);
1567 RCP<Operator> R = Levels_[i+1]->template Get< RCP<Operator> >(
"R");
1568 if(!R.is_null()) coarseRhs_[i] = MultiVectorFactory::Build(R->getRangeMap(),numvecs,
true);
1572 RCP<const Import> importer;
1573 if(Levels_[i+1]->IsAvailable(
"Importer"))
1574 importer = Levels_[i+1]->
template Get< RCP<const Import> >(
"Importer");
1575 if (doPRrebalance_ || importer.is_null())
1576 coarseX_[i] = MultiVectorFactory::Build(coarseRhs_[i]->getMap(),numvecs,
true);
1578 coarseImport_[i] = MultiVectorFactory::Build(importer->getTargetMap(), numvecs,
false);
1579 coarseExport_[i] = MultiVectorFactory::Build(importer->getSourceMap(), numvecs,
false);
1580 coarseX_[i] = MultiVectorFactory::Build(importer->getTargetMap(),numvecs,
false);
1584 sizeOfAllocatedLevelMultiVectors_ = numvecs;
1588 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1590 if(sizeOfAllocatedLevelMultiVectors_==0)
return;
1591 residual_.resize(0);
1592 coarseRhs_.resize(0);
1594 coarseImport_.resize(0);
1595 coarseExport_.resize(0);
1596 correction_.resize(0);
1597 sizeOfAllocatedLevelMultiVectors_ = 0;
1603 #endif // MUELU_HIERARCHY_DEF_HPP Important warning messages (one line)
void IsPreconditioner(const bool flag)
RCP< Level > & GetPreviousLevel()
Previous level.
static void SetRandomSeed(const Teuchos::Comm< int > &comm)
High level timing information (use Teuchos::TimeMonitor::summarize() to print)
double GetSmootherComplexity() const
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level->Get< RCP<Matrix> >("A", factory) if factory == NULL => use default factory.
RCP< Level > & GetLevel(const int levelID=0)
Retrieve a certain level from hierarchy.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Hierarchy()
Default constructor.
void DeleteLevelMultiVectors()
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.
void ReplaceCoordinateMap(Level &level)
void AllocateLevelMultiVectors(int numvecs)
void AddKeepFlag(const std::string &ename, const FactoryBase *factory=NoFactory::get(), KeepType keep=MueLu::Keep)
Call Level::AddKeepFlag for each level of the Hierarchy.
Print even more statistics.
void setlib(Xpetra::UnderlyingLib lib2)
int GetLevelID() const
Return level number.
void print(std::ostream &out=std::cout, const VerbLevel verbLevel=(MueLu::Parameters|MueLu::Statistics0)) const
Hierarchy::print is local hierarchy function, thus the statistics can be different from global ones...
Base class for factories (e.g., R, P, and A_coarse).
Print statistics that do not involve significant additional computation.
Xpetra::UnderlyingLib lib_
Epetra/Tpetra mode.
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.
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 SetMatvecParams(RCP< ParameterList > matvecParams)
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.
Array< RCP< Level > > Levels_
Container for Level objects.
void DumpCurrentGraph() const
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
Exception throws to report errors in the internal logical of the program.
An exception safe way to call the method 'Level::SetFactoryManager()'.
void SetComm(RCP< const Teuchos::Comm< int > > const &comm)
Print all warning messages.
ReturnType Iterate(const MultiVector &B, MultiVector &X, ConvData conv=ConvData(), bool InitialGuessIsZero=false, LO startLevel=0)
Apply the multigrid preconditioner.
Description of what is happening (more verbose)
void Keep(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Call Level::Keep(ename, factory) for each level of the Hierarchy.
std::string description() const
Return a simple one-line description of this object.
bool Setup(int coarseLevelID, const RCP< const FactoryManagerBase > fineLevelManager, const RCP< const FactoryManagerBase > coarseLevelManager, const RCP< const FactoryManagerBase > nextLevelManager=Teuchos::null)
Multi-level setup phase: build a new level of the hierarchy.
virtual void Clean() const
virtual std::string description() const
Return a simple one-line description of this object.
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
void Request(const FactoryBase &factory)
Increment the storage counter for all the inputs of a factory.