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_(-2), 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_(-2), 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 || dumpLevel_ == -1) && 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_) || dumpLevel_ == -1 ) )
510 DumpCurrentGraph(coarseLevelID);
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_;
557 AllocateLevelMultiVectors(sizeOfVecs,
true);
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;
1305 oss <<
"Scalar = " << Teuchos::ScalarTraits<Scalar>::name() << std::endl;
1306 oss <<
"Number of levels = " << numLevels << std::endl;
1307 oss <<
"Operator complexity = " << std::setprecision(2) << std::setiosflags(std::ios::fixed)
1308 << GetOperatorComplexity() << std::endl;
1310 if(smoother_comp!=-1.0) {
1311 oss <<
"Smoother complexity = " << std::setprecision(2) << std::setiosflags(std::ios::fixed)
1312 << smoother_comp << std::endl;
1317 oss <<
"Cycle type = V" << std::endl;
1320 oss <<
"Cycle type = W" << std::endl;
1321 if (WCycleStartLevel_ > 0)
1322 oss <<
"Cycle start level = " << WCycleStartLevel_ << std::endl;
1329 Xpetra::global_size_t tt = rowsPerLevel[0];
1330 int rowspacer = 2;
while (tt != 0) { tt /= 10; rowspacer++; }
1331 tt = nnzPerLevel[0];
1332 int nnzspacer = 2;
while (tt != 0) { tt /= 10; nnzspacer++; }
1333 tt = numProcsPerLevel[0];
1334 int npspacer = 2;
while (tt != 0) { tt /= 10; npspacer++; }
1335 oss <<
"level " << std::setw(rowspacer) <<
" rows " << std::setw(nnzspacer) <<
" nnz " <<
" nnz/row" << std::setw(npspacer) <<
" c ratio" <<
" procs" << std::endl;
1336 for (
size_t i = 0; i < nnzPerLevel.size(); ++i) {
1337 oss <<
" " << i <<
" ";
1338 oss << std::setw(rowspacer) << rowsPerLevel[i];
1339 oss << std::setw(nnzspacer) << nnzPerLevel[i];
1340 oss << std::setprecision(2) << std::setiosflags(std::ios::fixed);
1341 oss << std::setw(9) << as<double>(nnzPerLevel[i]) / rowsPerLevel[i];
1342 if (i) oss << std::setw(9) << as<double>(rowsPerLevel[i-1])/rowsPerLevel[i];
1343 else oss << std::setw(9) <<
" ";
1344 oss <<
" " << std::setw(npspacer) << numProcsPerLevel[i] << std::endl;
1347 for (
int i = 0; i < GetNumLevels(); ++i) {
1348 RCP<SmootherBase> preSmoo, postSmoo;
1349 if (Levels_[i]->IsAvailable(
"PreSmoother"))
1350 preSmoo = Levels_[i]->
template Get< RCP<SmootherBase> >(
"PreSmoother");
1351 if (Levels_[i]->IsAvailable(
"PostSmoother"))
1352 postSmoo = Levels_[i]->
template Get< RCP<SmootherBase> >(
"PostSmoother");
1354 if (preSmoo != null && preSmoo == postSmoo)
1355 oss <<
"Smoother (level " << i <<
") both : " << preSmoo->description() << std::endl;
1357 oss <<
"Smoother (level " << i <<
") pre : " 1358 << (preSmoo != null ? preSmoo->description() :
"no smoother") << std::endl;
1359 oss <<
"Smoother (level " << i <<
") post : " 1360 << (postSmoo != null ? postSmoo->description() :
"no smoother") << std::endl;
1371 RCP<const Teuchos::MpiComm<int> > mpiComm = rcp_dynamic_cast<
const Teuchos::MpiComm<int> >(comm);
1372 MPI_Comm rawComm = (*mpiComm->getRawMpiComm())();
1374 int strLength = outstr.size();
1375 MPI_Bcast(&strLength, 1, MPI_INT, root, rawComm);
1376 if (comm->getRank() != root)
1377 outstr.resize(strLength);
1378 MPI_Bcast(&outstr[0], strLength, MPI_CHAR, root, rawComm);
1385 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1387 Teuchos::OSTab tab2(out);
1388 for (
int i = 0; i < GetNumLevels(); ++i)
1389 Levels_[i]->print(out, verbLevel);
1392 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1394 isPreconditioner_ = flag;
1397 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1399 if (GetProcRankVerbose() != 0)
1401 #if defined(HAVE_MUELU_BOOST) && defined(HAVE_MUELU_BOOST_FOR_REAL) && defined(BOOST_VERSION) && (BOOST_VERSION >= 104400) 1406 dp.property(
"label", boost::get(boost::vertex_name, graph));
1407 dp.property(
"id", boost::get(boost::vertex_index, graph));
1408 dp.property(
"label", boost::get(boost::edge_name, graph));
1409 dp.property(
"color", boost::get(boost::edge_color, graph));
1412 std::map<const FactoryBase*, BoostVertex> vindices;
1413 typedef std::map<std::pair<BoostVertex,BoostVertex>, std::string> emap; emap edges;
1415 static int call_id=0;
1417 RCP<Operator> A = Levels_[0]->template Get<RCP<Operator> >(
"A");
1418 int rank = A->getDomainMap()->getComm()->getRank();
1421 for (
int i = currLevel; i <= currLevel+1 && i < GetNumLevels(); i++) {
1423 Levels_[i]->UpdateGraph(vindices, edges, dp, graph);
1425 for (emap::const_iterator eit = edges.begin(); eit != edges.end(); eit++) {
1426 std::pair<BoostEdge, bool> boost_edge = boost::add_edge(eit->first.first, eit->first.second, graph);
1429 if(eit->second==std::string(
"Graph")) boost::put(
"label", dp, boost_edge.first, std::string(
"Graph_"));
1430 else boost::put(
"label", dp, boost_edge.first, eit->second);
1432 boost::put(
"color", dp, boost_edge.first, std::string(
"red"));
1434 boost::put(
"color", dp, boost_edge.first, std::string(
"blue"));
1438 std::ofstream out(dumpFile_.c_str()+std::string(
"_")+std::to_string(currLevel)+std::string(
"_")+std::to_string(call_id)+std::string(
"_")+ std::to_string(rank) + std::string(
".dot"));
1439 boost::write_graphviz_dp(out, graph, dp, std::string(
"id"));
1443 GetOStream(
Errors) <<
"Dependency graph output requires boost and MueLu_ENABLE_Boost_for_real" << std::endl;
1448 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1450 RCP<Operator> Ao = level.
Get<RCP<Operator> >(
"A");
1451 RCP<Matrix> A = rcp_dynamic_cast<Matrix>(Ao);
1453 GetOStream(
Warnings1) <<
"Hierarchy::ReplaceCoordinateMap: operator is not a matrix, skipping..." << std::endl;
1456 if(Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A) != Teuchos::null) {
1457 GetOStream(
Warnings1) <<
"Hierarchy::ReplaceCoordinateMap: operator is a BlockedCrsMatrix, skipping..." << std::endl;
1461 typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO> xdMV;
1463 RCP<xdMV> coords = level.
Get<RCP<xdMV> >(
"Coordinates");
1465 if (A->getRowMap()->isSameAs(*(coords->getMap()))) {
1466 GetOStream(
Warnings1) <<
"Hierarchy::ReplaceCoordinateMap: matrix and coordinates maps are same, skipping..." << std::endl;
1470 if (A->IsView(
"stridedMaps") && rcp_dynamic_cast<
const StridedMap>(A->getRowMap(
"stridedMaps")) != Teuchos::null) {
1471 RCP<const StridedMap> stridedRowMap = rcp_dynamic_cast<
const StridedMap>(A->getRowMap(
"stridedMaps"));
1474 TEUCHOS_TEST_FOR_EXCEPTION(stridedRowMap->getStridedBlockId() != -1 || stridedRowMap->getOffset() != 0,
1475 Exceptions::RuntimeError,
"Hierarchy::ReplaceCoordinateMap: nontrivial maps (block id = " << stridedRowMap->getStridedBlockId()
1476 <<
", offset = " << stridedRowMap->getOffset() <<
")");
1479 GetOStream(
Runtime1) <<
"Replacing coordinate map" << std::endl;
1481 size_t blkSize = A->GetFixedBlockSize();
1483 RCP<const Map> nodeMap = A->getRowMap();
1486 RCP<const Map> dofMap = A->getRowMap();
1487 GO indexBase = dofMap->getIndexBase();
1488 size_t numLocalDOFs = dofMap->getNodeNumElements();
1490 "Hierarchy::ReplaceCoordinateMap: block size (" << blkSize <<
") is incompatible with the number of local dofs in a row map (" << numLocalDOFs);
1491 ArrayView<const GO> GIDs = dofMap->getNodeElementList();
1493 Array<GO> nodeGIDs(numLocalDOFs/blkSize);
1494 for (
size_t i = 0; i < numLocalDOFs; i += blkSize)
1495 nodeGIDs[i/blkSize] = (GIDs[i] - indexBase)/blkSize + indexBase;
1497 Xpetra::global_size_t INVALID = Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid();
1498 nodeMap = MapFactory::Build(dofMap->lib(), INVALID, nodeGIDs(), indexBase, dofMap->getComm());
1504 if(coords->getLocalLength() != A->getRowMap()->getNodeNumElements()) {
1505 GetOStream(
Warnings) <<
"Coordinate vector does not match row map of matrix A!" << std::endl;
1510 Array<ArrayView<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> > coordDataView;
1511 std::vector<ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> > coordData;
1512 for (
size_t i = 0; i < coords->getNumVectors(); i++) {
1513 coordData.push_back(coords->getData(i));
1514 coordDataView.push_back(coordData[i]());
1517 RCP<xdMV> newCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO>::Build(nodeMap, coordDataView(), coords->getNumVectors());
1518 level.
Set(
"Coordinates", newCoords);
1521 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1523 int N = Levels_.size();
1524 if( ( (sizeOfAllocatedLevelMultiVectors_ == numvecs && residual_.size() == N) || numvecs<=0 ) && !forceMapCheck)
return;
1527 if(residual_.size() != N) {
1528 DeleteLevelMultiVectors();
1530 residual_.resize(N);
1531 coarseRhs_.resize(N);
1533 coarseImport_.resize(N);
1534 coarseExport_.resize(N);
1535 correction_.resize(N);
1538 for(
int i=0; i<N; i++) {
1539 RCP<Operator> A = Levels_[i]->template Get< RCP<Operator> >(
"A");
1542 RCP<const BlockedCrsMatrix> A_as_blocked = Teuchos::rcp_dynamic_cast<
const BlockedCrsMatrix>(A);
1543 RCP<const Map> Arm = A->getRangeMap();
1544 RCP<const Map> Adm = A->getDomainMap();
1545 if(!A_as_blocked.is_null()) {
1546 Adm = A_as_blocked->getFullDomainMap();
1549 if (residual_[i].is_null() || !residual_[i]->getMap()->isSameAs(*Arm))
1551 residual_[i] = MultiVectorFactory::Build(Arm, numvecs,
true);
1552 if (correction_[i].is_null() || !correction_[i]->getMap()->isSameAs(*Adm))
1553 correction_[i] = MultiVectorFactory::Build(Adm, numvecs,
false);
1558 if(implicitTranspose_) {
1559 RCP<Operator> P = Levels_[i+1]->template Get< RCP<Operator> >(
"P");
1561 RCP<const Map> map = P->getDomainMap();
1562 if (coarseRhs_[i].is_null() || !coarseRhs_[i]->getMap()->isSameAs(*map))
1563 coarseRhs_[i] = MultiVectorFactory::Build(map, numvecs,
true);
1566 RCP<Operator> R = Levels_[i+1]->template Get< RCP<Operator> >(
"R");
1568 RCP<const Map> map = R->getRangeMap();
1569 if (coarseRhs_[i].is_null() || !coarseRhs_[i]->getMap()->isSameAs(*map))
1570 coarseRhs_[i] = MultiVectorFactory::Build(map, numvecs,
true);
1575 RCP<const Import> importer;
1576 if(Levels_[i+1]->IsAvailable(
"Importer"))
1577 importer = Levels_[i+1]->template Get< RCP<const Import> >(
"Importer");
1578 if (doPRrebalance_ || importer.is_null()) {
1579 RCP<const Map> map = coarseRhs_[i]->getMap();
1580 if (coarseX_[i].is_null() || !coarseX_[i]->getMap()->isSameAs(*map))
1581 coarseX_[i] = MultiVectorFactory::Build(map, numvecs,
true);
1584 map = importer->getTargetMap();
1585 if (coarseImport_[i].is_null() || !coarseImport_[i]->getMap()->isSameAs(*map)) {
1586 coarseImport_[i] = MultiVectorFactory::Build(map, numvecs,
false);
1587 coarseX_[i] = MultiVectorFactory::Build(map,numvecs,
false);
1589 map = importer->getSourceMap();
1590 if (coarseExport_[i].is_null() || !coarseExport_[i]->getMap()->isSameAs(*map))
1591 coarseExport_[i] = MultiVectorFactory::Build(map, numvecs,
false);
1595 sizeOfAllocatedLevelMultiVectors_ = numvecs;
1599 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1601 if(sizeOfAllocatedLevelMultiVectors_==0)
return;
1602 residual_.resize(0);
1603 coarseRhs_.resize(0);
1605 coarseImport_.resize(0);
1606 coarseExport_.resize(0);
1607 correction_.resize(0);
1608 sizeOfAllocatedLevelMultiVectors_ = 0;
1614 #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)
Translate Teuchos verbosity level to MueLu verbosity level.
One-liner description of what is happening.
void Clear(int startLevel=0)
Clear impermanent data from previous setup.
void DumpCurrentGraph(int level) const
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 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.
void AllocateLevelMultiVectors(int numvecs, bool forceMapCheck=false)
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.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
Print class parameters (more parameters, more verbose)
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.