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)");
305 TimeMonitor m2(*
this, label + this->ShortClassName() +
": " +
"Setup" +
" (total, level=" + Teuchos::toString(coarseLevelID) +
")");
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;
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));
406 if (coarseFact.is_null())
415 RCP<Operator> Ac = Teuchos::null;
416 TopRAPFactory coarseRAPFactory(fineLevelManager, coarseLevelManager);
419 Ac = level.
Get<RCP<Operator> >(
"A");
420 }
else if (!isFinestLevel) {
425 bool setLastLevelviaMaxCoarseSize =
false;
427 Ac = level.
Get<RCP<Operator> >(
"A");
428 RCP<Matrix> Acm = rcp_dynamic_cast<Matrix>(Ac);
432 level.
SetComm(Ac->getDomainMap()->getComm());
435 bool isOrigLastLevel = isLastLevel;
440 }
else if (Ac.is_null()) {
447 if (!Acm.is_null() && Acm->getGlobalNumRows() <= maxCoarseSize_) {
449 GetOStream(
Runtime0) <<
"Max coarse size (<= " << maxCoarseSize_ <<
") achieved" << std::endl;
451 if (Acm->getGlobalNumRows() != 0) setLastLevelviaMaxCoarseSize =
true;
455 if (!Ac.is_null() && !isFinestLevel) {
456 RCP<Operator> A = Levels_[coarseLevelID-1]->template Get< RCP<Operator> >(
"A");
457 RCP<Matrix> Am = rcp_dynamic_cast<Matrix>(A);
459 const double maxCoarse2FineRatio = 0.8;
460 if (!Acm.is_null() && !Am.is_null() && Acm->getGlobalNumRows() > maxCoarse2FineRatio * Am->getGlobalNumRows()) {
468 GetOStream(
Warnings0) <<
"Aggregation stagnated. Please check your matrix and/or adjust your configuration file."
469 <<
"Possible fixes:\n"
470 <<
" - reduce the maximum number of levels\n"
471 <<
" - enable repartitioning\n"
472 <<
" - increase the minimum coarse size." << std::endl;
478 if (!isOrigLastLevel) {
482 if (coarseFact.is_null())
490 coarseFact->Build(level);
501 smootherFact->Build(level);
506 if (isLastLevel ==
true) {
507 int actualNumLevels = nextLevelID;
508 if (isOrigLastLevel ==
false) {
511 Levels_[nextLevelID]->Release(
TopRAPFactory(coarseLevelManager, nextLevelManager));
518 if (!setLastLevelviaMaxCoarseSize) {
519 if (Levels_[nextLevelID-1]->IsAvailable(
"P")) {
520 if (Levels_[nextLevelID-1]->
template Get<RCP<Matrix> >(
"P") == Teuchos::null) actualNumLevels = nextLevelID-1;
522 else actualNumLevels = nextLevelID-1;
525 if (actualNumLevels == nextLevelID-1) {
527 Levels_[nextLevelID-2]->Release(*smootherFact);
529 if (Levels_[nextLevelID-2]->IsAvailable(
"PreSmoother") ) Levels_[nextLevelID-2]->RemoveKeepFlag(
"PreSmoother" ,
NoFactory::get());
530 if (Levels_[nextLevelID-2]->IsAvailable(
"PostSmoother")) Levels_[nextLevelID-2]->RemoveKeepFlag(
"PostSmoother",
NoFactory::get());
531 if (coarseFact.is_null())
533 Levels_[nextLevelID-2]->Request(*coarseFact);
534 if ( !(Levels_[nextLevelID-2]->
template Get<RCP<Matrix> >(
"A").is_null() ))
535 coarseFact->Build( *(Levels_[nextLevelID-2]));
536 Levels_[nextLevelID-2]->Release(*coarseFact);
538 Levels_.resize(actualNumLevels);
542 if (isDumpingEnabled_ && ( (dumpLevel_ > 0 && coarseLevelID == dumpLevel_) || dumpLevel_ == -1 ) )
543 DumpCurrentGraph(coarseLevelID);
545 if (!isFinestLevel) {
549 level.
Release(coarseRAPFactory);
553 SetProcRankVerbose(oldRank);
558 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
560 int numLevels = Levels_.size();
562 "Hierarchy::SetupRe: " << Levels_.size() <<
" levels, but " << levelManagers_.size() <<
" level factory managers");
564 const int startLevel = 0;
567#ifdef HAVE_MUELU_DEBUG
569 for (
int i = 0; i < numLevels; i++)
570 levelManagers_[i]->ResetDebugData();
575 for (levelID = startLevel; levelID < numLevels;) {
576 bool r = Setup(levelID,
577 (levelID != 0 ? levelManagers_[levelID-1] : Teuchos::null),
578 levelManagers_[levelID],
579 (levelID+1 != numLevels ? levelManagers_[levelID+1] : Teuchos::null));
585 Levels_ .resize(levelID);
586 levelManagers_.resize(levelID);
588 int sizeOfVecs = sizeOfAllocatedLevelMultiVectors_;
590 AllocateLevelMultiVectors(sizeOfVecs,
true);
598 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
607 "MueLu::Hierarchy::Setup(): fine level (" << startLevel <<
") does not exist");
610 "Constructing non-positive (" << numDesiredLevels <<
") number of levels does not make sense.");
614 "MueLu::Hierarchy::Setup(): fine level (" << startLevel <<
") has no matrix A! "
615 "Set fine level matrix A using Level.Set()");
617 RCP<Operator> A = Levels_[startLevel]->template Get<RCP<Operator> >(
"A");
618 lib_ = A->getDomainMap()->lib();
621 RCP<Matrix> Amat = rcp_dynamic_cast<Matrix>(A);
623 if (!Amat.is_null()) {
624 RCP<ParameterList> params = rcp(
new ParameterList());
625 params->set(
"printLoadBalancingInfo",
true);
626 params->set(
"printCommInfo",
true);
630 GetOStream(
Warnings1) <<
"Fine level operator is not a matrix, statistics are not available" << std::endl;
634 RCP<const FactoryManagerBase> rcpmanager = rcpFromRef(manager);
636 const int lastLevel = startLevel + numDesiredLevels - 1;
637 GetOStream(
Runtime0) <<
"Setup loop: startLevel = " << startLevel <<
", lastLevel = " << lastLevel
638 <<
" (stop if numLevels = " << numDesiredLevels <<
" or Ac.size() < " << maxCoarseSize_ <<
")" << std::endl;
642 if (numDesiredLevels == 1) {
644 Setup(startLevel, Teuchos::null, rcpmanager, Teuchos::null);
647 bool bIsLastLevel = Setup(startLevel, Teuchos::null, rcpmanager, rcpmanager);
648 if (bIsLastLevel ==
false) {
649 for (iLevel = startLevel + 1; iLevel < lastLevel; iLevel++) {
650 bIsLastLevel = Setup(iLevel, rcpmanager, rcpmanager, rcpmanager);
651 if (bIsLastLevel ==
true)
654 if (bIsLastLevel ==
false)
655 Setup(lastLevel, rcpmanager, rcpmanager, Teuchos::null);
661 "MueLu::Hierarchy::Setup(): number of level");
670 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
672 if (startLevel < GetNumLevels())
673 GetOStream(
Runtime0) <<
"Clearing old data (if any)" << std::endl;
675 for (
int iLevel = startLevel; iLevel < GetNumLevels(); iLevel++)
676 Levels_[iLevel]->Clear();
679 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
681 GetOStream(
Runtime0) <<
"Clearing old data (expert)" << std::endl;
682 for (
int iLevel = 0; iLevel < GetNumLevels(); iLevel++)
683 Levels_[iLevel]->ExpertClear();
686#if defined(HAVE_MUELU_EXPERIMENTAL) && defined(HAVE_MUELU_ADDITIVE_VARIANT)
687 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
689 bool InitialGuessIsZero, LO startLevel) {
690 LO nIts = conv.maxIts_;
691 MagnitudeType tol = conv.tol_;
693 std::string prefix = this->ShortClassName() +
": ";
694 std::string levelSuffix =
" (level=" +
toString(startLevel) +
")";
695 std::string levelSuffix1 =
" (level=" +
toString(startLevel+1) +
")";
698 RCP<Time> CompTime = Teuchos::TimeMonitor::getNewCounter(prefix +
"Computational Time (total)");
699 RCP<Time> Concurrent = Teuchos::TimeMonitor::getNewCounter(prefix +
"Concurrent portion");
700 RCP<Time> ApplyR = Teuchos::TimeMonitor::getNewCounter(prefix +
"R: Computational Time");
701 RCP<Time> ApplyPbar = Teuchos::TimeMonitor::getNewCounter(prefix +
"Pbar: Computational Time");
702 RCP<Time> CompFine = Teuchos::TimeMonitor::getNewCounter(prefix +
"Fine: Computational Time");
703 RCP<Time> CompCoarse = Teuchos::TimeMonitor::getNewCounter(prefix +
"Coarse: Computational Time");
704 RCP<Time> ApplySum = Teuchos::TimeMonitor::getNewCounter(prefix +
"Sum: Computational Time");
705 RCP<Time> Synchronize_beginning = Teuchos::TimeMonitor::getNewCounter(prefix +
"Synchronize_beginning");
706 RCP<Time> Synchronize_center = Teuchos::TimeMonitor::getNewCounter(prefix +
"Synchronize_center");
707 RCP<Time> Synchronize_end = Teuchos::TimeMonitor::getNewCounter(prefix +
"Synchronize_end");
709 RCP<Level> Fine = Levels_[0];
712 RCP<Operator> A = Fine->Get< RCP<Operator> >(
"A");
713 Teuchos::RCP< const Teuchos::Comm< int > > communicator = A->getDomainMap()->getComm();
722 SC one = STS::one(), zero = STS::zero();
724 bool zeroGuess = InitialGuessIsZero;
730 RCP< Operator > Pbar;
732 RCP< MultiVector > coarseRhs, coarseX;
734 RCP<SmootherBase> preSmoo_coarse, postSmoo_coarse;
735 bool emptyCoarseSolve =
true;
736 RCP<MultiVector> coarseX_prolonged = MultiVectorFactory::Build(X.getMap(), X.getNumVectors(),
true);
738 RCP<const Import> importer;
740 if (Levels_.size() > 1) {
742 if (Coarse->IsAvailable(
"Importer"))
743 importer = Coarse->Get< RCP<const Import> >(
"Importer");
745 R = Coarse->Get< RCP<Operator> >(
"R");
746 P = Coarse->Get< RCP<Operator> >(
"P");
750 Pbar = Coarse->Get< RCP<Operator> >(
"Pbar");
752 coarseRhs = MultiVectorFactory::Build(R->getRangeMap(), B.getNumVectors(),
true);
754 Ac = Coarse->Get< RCP< Operator > >(
"A");
757 R->apply(B, *coarseRhs, Teuchos::NO_TRANS, one, zero);
761 if (doPRrebalance_ || importer.is_null()) {
762 coarseX = MultiVectorFactory::Build(coarseRhs->getMap(), X.getNumVectors(),
true);
766 RCP<TimeMonitor> ITime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : import (total)" ,
Timings0));
767 RCP<TimeMonitor> ILevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : import" + levelSuffix1,
Timings0));
770 RCP<MultiVector> coarseTmp = MultiVectorFactory::Build(importer->getTargetMap(), coarseRhs->getNumVectors());
771 coarseTmp->doImport(*coarseRhs, *importer, Xpetra::INSERT);
772 coarseRhs.swap(coarseTmp);
774 coarseX = MultiVectorFactory::Build(importer->getTargetMap(), X.getNumVectors(),
true);
777 if (Coarse->IsAvailable(
"PreSmoother"))
778 preSmoo_coarse = Coarse->Get< RCP<SmootherBase> >(
"PreSmoother");
779 if (Coarse->IsAvailable(
"PostSmoother"))
780 postSmoo_coarse = Coarse->Get< RCP<SmootherBase> >(
"PostSmoother");
786 MagnitudeType prevNorm = STS::magnitude(STS::one()), curNorm = STS::magnitude(STS::one());
789 for (LO i = 1; i <= nIts; i++) {
790#ifdef HAVE_MUELU_DEBUG
791 if (A->getDomainMap()->isCompatible(*(X.getMap())) ==
false) {
792 std::ostringstream ss;
793 ss <<
"Level " << startLevel <<
": level A's domain map is not compatible with X";
794 throw Exceptions::Incompatible(ss.str());
797 if (A->getRangeMap()->isCompatible(*(B.getMap())) ==
false) {
798 std::ostringstream ss;
799 ss <<
"Level " << startLevel <<
": level A's range map is not compatible with B";
800 throw Exceptions::Incompatible(ss.str());
805 bool emptyFineSolve =
true;
807 RCP< MultiVector > fineX;
808 fineX = MultiVectorFactory::Build(X.getMap(), X.getNumVectors(),
true);
817 if (Fine->IsAvailable(
"PreSmoother")) {
818 RCP<SmootherBase> preSmoo = Fine->Get< RCP<SmootherBase> >(
"PreSmoother");
820 preSmoo->Apply(*fineX, B, zeroGuess);
822 emptyFineSolve =
false;
824 if (Fine->IsAvailable(
"PostSmoother")) {
825 RCP<SmootherBase> postSmoo = Fine->Get< RCP<SmootherBase> >(
"PostSmoother");
827 postSmoo->Apply(*fineX, B, zeroGuess);
830 emptyFineSolve =
false;
832 if (emptyFineSolve ==
true) {
833 GetOStream(
Warnings1) <<
"No fine grid smoother" << std::endl;
835 fineX->update(one, B, zero);
838 if (Levels_.size() > 1) {
840 if (Coarse->IsAvailable(
"PreSmoother")) {
842 preSmoo_coarse->Apply(*coarseX, *coarseRhs, zeroGuess);
844 emptyCoarseSolve =
false;
847 if (Coarse->IsAvailable(
"PostSmoother")) {
849 postSmoo_coarse->Apply(*coarseX, *coarseRhs, zeroGuess);
851 emptyCoarseSolve =
false;
854 if (emptyCoarseSolve ==
true) {
855 GetOStream(
Warnings1) <<
"No coarse grid solver" << std::endl;
857 coarseX->update(one, *coarseRhs, zero);
864 if (!doPRrebalance_ && !importer.is_null()) {
865 RCP<TimeMonitor> ITime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : export (total)" ,
Timings0));
866 RCP<TimeMonitor> ILevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : export" + levelSuffix1,
Timings0));
869 RCP<MultiVector> coarseTmp = MultiVectorFactory::Build(importer->getSourceMap(), coarseX->getNumVectors());
870 coarseTmp->doExport(*coarseX, *importer, Xpetra::INSERT);
871 coarseX.swap(coarseTmp);
875 Pbar->apply(*coarseX, *coarseX_prolonged, Teuchos::NO_TRANS, one, zero);
880 X.update(1.0, *fineX, 1.0, *coarseX_prolonged, 0.0);
891 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
893 bool InitialGuessIsZero, LO startLevel) {
909 RCP<Level> Fine = Levels_[startLevel];
912 std::string prefix = label + this->ShortClassName() +
": ";
913 std::string levelSuffix =
" (level=" +
toString(startLevel) +
")";
914 std::string levelSuffix1 =
" (level=" +
toString(startLevel+1) +
")";
916 bool useStackedTimer = !Teuchos::TimeMonitor::getStackedTimer().is_null();
918 RCP<Monitor> iterateTime;
919 RCP<TimeMonitor> iterateTime1;
922 else if (!useStackedTimer)
925 std::string iterateLevelTimeLabel = prefix +
"Solve" + levelSuffix;
926 RCP<TimeMonitor> iterateLevelTime = rcp(
new TimeMonitor(*
this, iterateLevelTimeLabel,
Timings0));
928 bool zeroGuess = InitialGuessIsZero;
930 RCP<Operator> A = Fine->Get< RCP<Operator> >(
"A");
932 RCP<Time> CompCoarse = Teuchos::TimeMonitor::getNewCounter(prefix +
"Coarse: Computational Time");
945 const BlockedMultiVector * Bblocked =
dynamic_cast<const BlockedMultiVector*
>(&B);
946 if(residual_.size() > startLevel &&
947 ( ( Bblocked && !Bblocked->isSameSize(*residual_[startLevel])) ||
948 (!Bblocked && !residual_[startLevel]->isSameSize(B))))
949 DeleteLevelMultiVectors();
950 AllocateLevelMultiVectors(X.getNumVectors());
953 typedef Teuchos::ScalarTraits<typename STS::magnitudeType> STM;
956 if (IsCalculationOfResidualRequired(startLevel, conv)) {
957 ConvergenceStatus convergenceStatus = ComputeResidualAndPrintHistory(*A, X, B, Teuchos::ScalarTraits<LO>::zero(), startLevel, conv, prevNorm);
959 return convergenceStatus;
962 SC one = STS::one(), zero = STS::zero();
963 for (LO iteration = 1; iteration <= nIts; iteration++) {
964#ifdef HAVE_MUELU_DEBUG
966 if (A->getDomainMap()->isCompatible(*(X.getMap())) ==
false) {
967 std::ostringstream ss;
968 ss <<
"Level " << startLevel <<
": level A's domain map is not compatible with X";
972 if (A->getRangeMap()->isCompatible(*(B.getMap())) ==
false) {
973 std::ostringstream ss;
974 ss <<
"Level " << startLevel <<
": level A's range map is not compatible with B";
980 if (startLevel == as<LO>(Levels_.size())-1) {
982 RCP<TimeMonitor> CLevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : coarse" + levelSuffix,
Timings0));
984 bool emptySolve =
true;
987 if (Fine->IsAvailable(
"PreSmoother")) {
988 RCP<SmootherBase> preSmoo = Fine->Get< RCP<SmootherBase> >(
"PreSmoother");
990 preSmoo->Apply(X, B, zeroGuess);
995 if (Fine->IsAvailable(
"PostSmoother")) {
996 RCP<SmootherBase> postSmoo = Fine->Get< RCP<SmootherBase> >(
"PostSmoother");
998 postSmoo->Apply(X, B, zeroGuess);
1003 if (emptySolve ==
true) {
1004 GetOStream(
Warnings1) <<
"No coarse grid solver" << std::endl;
1006 X.update(one, B, zero);
1011 RCP<Level> Coarse = Levels_[startLevel+1];
1014 RCP<TimeMonitor> STime;
1015 if (!useStackedTimer)
1017 RCP<TimeMonitor> SLevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : smoothing" + levelSuffix,
Timings0));
1019 if (Fine->IsAvailable(
"PreSmoother")) {
1020 RCP<SmootherBase> preSmoo = Fine->Get< RCP<SmootherBase> >(
"PreSmoother");
1021 preSmoo->Apply(X, B, zeroGuess);
1026 RCP<MultiVector> residual;
1028 RCP<TimeMonitor> ATime;
1029 if (!useStackedTimer)
1030 ATime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : residual calculation (total)" ,
Timings0));
1031 RCP<TimeMonitor> ALevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : residual calculation" + levelSuffix,
Timings0));
1039 residual = residual_[startLevel];
1042 RCP<Operator> P = Coarse->Get< RCP<Operator> >(
"P");
1043 if (Coarse->IsAvailable(
"Pbar"))
1044 P = Coarse->Get< RCP<Operator> >(
"Pbar");
1046 RCP<MultiVector> coarseRhs, coarseX;
1050 RCP<TimeMonitor> RTime;
1051 if (!useStackedTimer)
1053 RCP<TimeMonitor> RLevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : restriction" + levelSuffix,
Timings0));
1054 coarseRhs = coarseRhs_[startLevel];
1056 if (implicitTranspose_) {
1057 P->apply(*residual, *coarseRhs, Teuchos::TRANS, one, zero);
1060 RCP<Operator> R = Coarse->Get< RCP<Operator> >(
"R");
1061 R->apply(*residual, *coarseRhs, Teuchos::NO_TRANS, one, zero);
1065 RCP<const Import> importer;
1066 if (Coarse->IsAvailable(
"Importer"))
1067 importer = Coarse->Get< RCP<const Import> >(
"Importer");
1069 coarseX = coarseX_[startLevel];
1070 if (!doPRrebalance_ && !importer.is_null()) {
1071 RCP<TimeMonitor> ITime;
1072 if (!useStackedTimer)
1074 RCP<TimeMonitor> ILevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : import" + levelSuffix1,
Timings0));
1077 RCP<MultiVector> coarseTmp = coarseImport_[startLevel];
1078 coarseTmp->doImport(*coarseRhs, *importer, Xpetra::INSERT);
1079 coarseRhs.swap(coarseTmp);
1082 RCP<Operator> Ac = Coarse->Get< RCP<Operator> >(
"A");
1083 if (!Ac.is_null()) {
1084 RCP<const Map> origXMap = coarseX->getMap();
1085 RCP<const Map> origRhsMap = coarseRhs->getMap();
1088 coarseRhs->replaceMap(Ac->getRangeMap());
1089 coarseX ->replaceMap(Ac->getDomainMap());
1092 iterateLevelTime = Teuchos::null;
1094 Iterate(*coarseRhs, *coarseX, 1,
true, startLevel+1);
1096 if (Cycle_ ==
WCYCLE && WCycleStartLevel_ >= startLevel)
1097 Iterate(*coarseRhs, *coarseX, 1,
false, startLevel+1);
1100 iterateLevelTime = rcp(
new TimeMonitor(*
this, iterateLevelTimeLabel));
1102 coarseX->replaceMap(origXMap);
1103 coarseRhs->replaceMap(origRhsMap);
1106 if (!doPRrebalance_ && !importer.is_null()) {
1107 RCP<TimeMonitor> ITime;
1108 if (!useStackedTimer)
1110 RCP<TimeMonitor> ILevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : export" + levelSuffix1,
Timings0));
1113 RCP<MultiVector> coarseTmp = coarseExport_[startLevel];
1114 coarseTmp->doExport(*coarseX, *importer, Xpetra::INSERT);
1115 coarseX.swap(coarseTmp);
1120 RCP<TimeMonitor> PTime;
1121 if (!useStackedTimer)
1123 RCP<TimeMonitor> PLevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : prolongation" + levelSuffix,
Timings0));
1128 if (fuseProlongationAndUpdate_) {
1129 P->apply(*coarseX, X, Teuchos::NO_TRANS, scalingFactor_, one);
1131 RCP<MultiVector> correction = correction_[startLevel];
1132 P->apply(*coarseX, *correction, Teuchos::NO_TRANS, one, zero);
1133 X.update(scalingFactor_, *correction, one);
1139 RCP<TimeMonitor> STime;
1140 if (!useStackedTimer)
1142 RCP<TimeMonitor> SLevelTime = rcp(
new TimeMonitor(*
this, prefix +
"Solve : smoothing" + levelSuffix,
Timings0));
1144 if (Fine->IsAvailable(
"PostSmoother")) {
1145 RCP<SmootherBase> postSmoo = Fine->Get< RCP<SmootherBase> >(
"PostSmoother");
1146 postSmoo->Apply(X, B,
false);
1153 if (IsCalculationOfResidualRequired(startLevel, conv)) {
1154 ConvergenceStatus convergenceStatus = ComputeResidualAndPrintHistory(*A, X, B, iteration, startLevel, conv, prevNorm);
1156 return convergenceStatus;
1164 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1166 LO startLevel = (start != -1 ? start : 0);
1167 LO endLevel = (end != -1 ? end : Levels_.size()-1);
1170 "MueLu::Hierarchy::Write : startLevel must be <= endLevel");
1173 "MueLu::Hierarchy::Write bad start or end level");
1175 for (LO i = startLevel; i < endLevel + 1; i++) {
1176 RCP<Matrix> A = rcp_dynamic_cast<Matrix>(Levels_[i]->
template Get< RCP< Operator> >(
"A")), P, R;
1178 P = rcp_dynamic_cast<Matrix>(Levels_[i]->
template Get< RCP< Operator> >(
"P"));
1179 if (!implicitTranspose_)
1180 R = rcp_dynamic_cast<Matrix>(Levels_[i]->
template Get< RCP< Operator> >(
"R"));
1183 if (!A.is_null()) Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(
"A_" +
toString(i) + suffix +
".m", *A);
1185 Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(
"P_" +
toString(i) + suffix +
".m", *P);
1188 Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(
"R_" +
toString(i) + suffix +
".m", *R);
1193 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1195 for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1196 (*it)->Keep(ename, factory);
1199 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1201 for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1202 (*it)->Delete(ename, factory);
1205 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1207 for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1208 (*it)->AddKeepFlag(ename, factory, keep);
1211 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1213 for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1214 (*it)->RemoveKeepFlag(ename, factory, keep);
1217 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1219 if ( description_.empty() )
1221 std::ostringstream out;
1223 out <<
"{#levels = " << GetGlobalNumLevels() <<
", complexity = " << GetOperatorComplexity() <<
"}";
1224 description_ = out.str();
1226 return description_;
1229 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1234 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1236 RCP<Operator> A0 = Levels_[0]->template Get<RCP<Operator> >(
"A");
1237 RCP<const Teuchos::Comm<int> > comm = A0->getDomainMap()->getComm();
1239 int numLevels = GetNumLevels();
1240 RCP<Operator> Ac = Levels_[numLevels-1]->template Get<RCP<Operator> >(
"A");
1247 int root = comm->getRank();
1250 int smartData = numLevels*comm->getSize() + comm->getRank(), maxSmartData;
1251 reduceAll(*comm, Teuchos::REDUCE_MAX, smartData, Teuchos::ptr(&maxSmartData));
1252 root = maxSmartData % comm->getSize();
1256 double smoother_comp =-1.0;
1258 smoother_comp = GetSmootherComplexity();
1262 std::vector<Xpetra::global_size_t> nnzPerLevel;
1263 std::vector<Xpetra::global_size_t> rowsPerLevel;
1264 std::vector<int> numProcsPerLevel;
1265 bool someOpsNotMatrices =
false;
1266 const Xpetra::global_size_t INVALID = Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid();
1267 for (
int i = 0; i < numLevels; i++) {
1269 "Operator A is not available on level " << i);
1271 RCP<Operator> A = Levels_[i]->template Get<RCP<Operator> >(
"A");
1273 "Operator A on level " << i <<
" is null.");
1275 RCP<Matrix> Am = rcp_dynamic_cast<Matrix>(A);
1277 someOpsNotMatrices =
true;
1278 nnzPerLevel .push_back(INVALID);
1279 rowsPerLevel .push_back(A->getDomainMap()->getGlobalNumElements());
1280 numProcsPerLevel.push_back(A->getDomainMap()->getComm()->getSize());
1282 LO storageblocksize=Am->GetStorageBlockSize();
1283 Xpetra::global_size_t nnz = Am->getGlobalNumEntries()*storageblocksize*storageblocksize;
1284 nnzPerLevel .push_back(nnz);
1285 rowsPerLevel .push_back(Am->getGlobalNumRows()*storageblocksize);
1286 numProcsPerLevel.push_back(Am->getRowMap()->getComm()->getSize());
1289 if (someOpsNotMatrices)
1290 GetOStream(
Warnings0) <<
"Some level operators are not matrices, statistics calculation are incomplete" << std::endl;
1293 std::string label = Levels_[0]->getObjectLabel();
1294 std::ostringstream oss;
1295 oss << std::setfill(
' ');
1296 oss <<
"\n--------------------------------------------------------------------------------\n";
1297 oss <<
"--- Multigrid Summary " << std::setw(28) << std::left << label <<
"---\n";
1298 oss <<
"--------------------------------------------------------------------------------" << std::endl;
1300 oss <<
"Scalar = " << Teuchos::ScalarTraits<Scalar>::name() << std::endl;
1301 oss <<
"Number of levels = " << numLevels << std::endl;
1302 oss <<
"Operator complexity = " << std::setprecision(2) << std::setiosflags(std::ios::fixed);
1303 if (!someOpsNotMatrices)
1304 oss << GetOperatorComplexity() << std::endl;
1306 oss <<
"not available (Some operators in hierarchy are not matrices.)" << 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 for (
size_t i = 0; i < nnzPerLevel.size(); ++i) {
1330 tt = nnzPerLevel[i];
1335 int nnzspacer = 2;
while (tt != 0) { tt /= 10; nnzspacer++; }
1336 tt = numProcsPerLevel[0];
1337 int npspacer = 2;
while (tt != 0) { tt /= 10; npspacer++; }
1338 oss <<
"level " << std::setw(rowspacer) <<
" rows " << std::setw(nnzspacer) <<
" nnz " <<
" nnz/row" << std::setw(npspacer) <<
" c ratio" <<
" procs" << std::endl;
1339 for (
size_t i = 0; i < nnzPerLevel.size(); ++i) {
1340 oss <<
" " << i <<
" ";
1341 oss << std::setw(rowspacer) << rowsPerLevel[i];
1342 if (nnzPerLevel[i] != INVALID) {
1343 oss << std::setw(nnzspacer) << nnzPerLevel[i];
1344 oss << std::setprecision(2) << std::setiosflags(std::ios::fixed);
1345 oss << std::setw(9) << as<double>(nnzPerLevel[i]) / rowsPerLevel[i];
1347 oss << std::setw(nnzspacer) <<
"Operator";
1348 oss << std::setprecision(2) << std::setiosflags(std::ios::fixed);
1349 oss << std::setw(9) <<
" ";
1351 if (i) oss << std::setw(9) << as<double>(rowsPerLevel[i-1])/rowsPerLevel[i];
1352 else oss << std::setw(9) <<
" ";
1353 oss <<
" " << std::setw(npspacer) << numProcsPerLevel[i] << std::endl;
1356 for (
int i = 0; i < GetNumLevels(); ++i) {
1357 RCP<SmootherBase> preSmoo, postSmoo;
1358 if (Levels_[i]->IsAvailable(
"PreSmoother"))
1359 preSmoo = Levels_[i]->
template Get< RCP<SmootherBase> >(
"PreSmoother");
1360 if (Levels_[i]->IsAvailable(
"PostSmoother"))
1361 postSmoo = Levels_[i]->
template Get< RCP<SmootherBase> >(
"PostSmoother");
1363 if (preSmoo != null && preSmoo == postSmoo)
1364 oss <<
"Smoother (level " << i <<
") both : " << preSmoo->description() << std::endl;
1366 oss <<
"Smoother (level " << i <<
") pre : "
1367 << (preSmoo != null ? preSmoo->description() :
"no smoother") << std::endl;
1368 oss <<
"Smoother (level " << i <<
") post : "
1369 << (postSmoo != null ? postSmoo->description() :
"no smoother") << std::endl;
1380 RCP<const Teuchos::MpiComm<int> > mpiComm = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(comm);
1381 MPI_Comm rawComm = (*mpiComm->getRawMpiComm())();
1383 int strLength = outstr.size();
1384 MPI_Bcast(&strLength, 1, MPI_INT, root, rawComm);
1385 if (comm->getRank() != root)
1386 outstr.resize(strLength);
1387 MPI_Bcast(&outstr[0], strLength, MPI_CHAR, root, rawComm);
1394 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1396 Teuchos::OSTab tab2(out);
1397 for (
int i = 0; i < GetNumLevels(); ++i)
1398 Levels_[i]->print(out, verbLevel);
1401 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1403 isPreconditioner_ = flag;
1406 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1408 if (GetProcRankVerbose() != 0)
1410#if defined(HAVE_MUELU_BOOST) && defined(HAVE_MUELU_BOOST_FOR_REAL) && defined(BOOST_VERSION) && (BOOST_VERSION >= 104400)
1415 dp.property(
"label", boost::get(boost::vertex_name, graph));
1416 dp.property(
"id", boost::get(boost::vertex_index, graph));
1417 dp.property(
"label", boost::get(boost::edge_name, graph));
1418 dp.property(
"color", boost::get(boost::edge_color, graph));
1421 std::map<const FactoryBase*, BoostVertex> vindices;
1422 typedef std::map<std::pair<BoostVertex,BoostVertex>, std::string> emap; emap edges;
1424 static int call_id=0;
1426 RCP<Operator> A = Levels_[0]->template Get<RCP<Operator> >(
"A");
1427 int rank = A->getDomainMap()->getComm()->getRank();
1430 for (
int i = currLevel; i <= currLevel+1 && i < GetNumLevels(); i++) {
1432 Levels_[i]->UpdateGraph(vindices, edges, dp, graph);
1434 for (emap::const_iterator eit = edges.begin(); eit != edges.end(); eit++) {
1435 std::pair<BoostEdge, bool> boost_edge = boost::add_edge(eit->first.first, eit->first.second, graph);
1438 if(eit->second==std::string(
"Graph")) boost::put(
"label", dp, boost_edge.first, std::string(
"Graph_"));
1439 else boost::put(
"label", dp, boost_edge.first, eit->second);
1441 boost::put(
"color", dp, boost_edge.first, std::string(
"red"));
1443 boost::put(
"color", dp, boost_edge.first, std::string(
"blue"));
1447 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"));
1448 boost::write_graphviz_dp(out, graph, dp, std::string(
"id"));
1452 GetOStream(
Errors) <<
"Dependency graph output requires boost and MueLu_ENABLE_Boost_for_real" << std::endl;
1457 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1459 RCP<Operator> Ao = level.
Get<RCP<Operator> >(
"A");
1460 RCP<Matrix> A = rcp_dynamic_cast<Matrix>(Ao);
1462 GetOStream(
Runtime1) <<
"Hierarchy::ReplaceCoordinateMap: operator is not a matrix, skipping..." << std::endl;
1465 if(Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A) != Teuchos::null) {
1466 GetOStream(
Runtime1) <<
"Hierarchy::ReplaceCoordinateMap: operator is a BlockedCrsMatrix, skipping..." << std::endl;
1470 typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO> xdMV;
1472 RCP<xdMV> coords = level.
Get<RCP<xdMV> >(
"Coordinates");
1474 if (A->getRowMap()->isSameAs(*(coords->getMap()))) {
1475 GetOStream(
Runtime1) <<
"Hierarchy::ReplaceCoordinateMap: matrix and coordinates maps are same, skipping..." << std::endl;
1479 if (A->IsView(
"stridedMaps") && rcp_dynamic_cast<const StridedMap>(A->getRowMap(
"stridedMaps")) != Teuchos::null) {
1480 RCP<const StridedMap> stridedRowMap = rcp_dynamic_cast<const StridedMap>(A->getRowMap(
"stridedMaps"));
1483 TEUCHOS_TEST_FOR_EXCEPTION(stridedRowMap->getStridedBlockId() != -1 || stridedRowMap->getOffset() != 0,
1484 Exceptions::RuntimeError,
"Hierarchy::ReplaceCoordinateMap: nontrivial maps (block id = " << stridedRowMap->getStridedBlockId()
1485 <<
", offset = " << stridedRowMap->getOffset() <<
")");
1488 GetOStream(
Runtime1) <<
"Replacing coordinate map" << std::endl;
1489 TEUCHOS_TEST_FOR_EXCEPTION(A->GetFixedBlockSize() % A->GetStorageBlockSize() != 0,
Exceptions::RuntimeError,
"Hierarchy::ReplaceCoordinateMap: Storage block size does not evenly divide fixed block size");
1491 size_t blkSize = A->GetFixedBlockSize() / A->GetStorageBlockSize();
1493 RCP<const Map> nodeMap = A->getRowMap();
1496 RCP<const Map> dofMap = A->getRowMap();
1497 GO indexBase = dofMap->getIndexBase();
1498 size_t numLocalDOFs = dofMap->getLocalNumElements();
1500 "Hierarchy::ReplaceCoordinateMap: block size (" << blkSize <<
") is incompatible with the number of local dofs in a row map (" << numLocalDOFs);
1501 ArrayView<const GO> GIDs = dofMap->getLocalElementList();
1503 Array<GO> nodeGIDs(numLocalDOFs/blkSize);
1504 for (
size_t i = 0; i < numLocalDOFs; i += blkSize)
1505 nodeGIDs[i/blkSize] = (GIDs[i] - indexBase)/blkSize + indexBase;
1507 Xpetra::global_size_t INVALID = Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid();
1508 nodeMap = MapFactory::Build(dofMap->lib(), INVALID, nodeGIDs(), indexBase, dofMap->getComm());
1514 if(coords->getLocalLength() != A->getRowMap()->getLocalNumElements()) {
1515 GetOStream(
Warnings) <<
"Coordinate vector does not match row map of matrix A!" << std::endl;
1520 Array<ArrayView<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> > coordDataView;
1521 std::vector<ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> > coordData;
1522 for (
size_t i = 0; i < coords->getNumVectors(); i++) {
1523 coordData.push_back(coords->getData(i));
1524 coordDataView.push_back(coordData[i]());
1527 RCP<xdMV> newCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO>::Build(nodeMap, coordDataView(), coords->getNumVectors());
1528 level.
Set(
"Coordinates", newCoords);
1531 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1533 int N = Levels_.size();
1534 if( ( (sizeOfAllocatedLevelMultiVectors_ == numvecs && residual_.size() == N) || numvecs<=0 ) && !forceMapCheck)
return;
1537 if(residual_.size() != N) {
1538 DeleteLevelMultiVectors();
1540 residual_.resize(N);
1541 coarseRhs_.resize(N);
1543 coarseImport_.resize(N);
1544 coarseExport_.resize(N);
1545 correction_.resize(N);
1548 for(
int i=0; i<N; i++) {
1549 RCP<Operator> A = Levels_[i]->template Get< RCP<Operator> >(
"A");
1552 RCP<const BlockedCrsMatrix> A_as_blocked = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(A);
1553 RCP<const Map> Arm = A->getRangeMap();
1554 RCP<const Map> Adm = A->getDomainMap();
1555 if(!A_as_blocked.is_null()) {
1556 Adm = A_as_blocked->getFullDomainMap();
1559 if (residual_[i].is_null() || !residual_[i]->getMap()->isSameAs(*Arm))
1561 residual_[i] = MultiVectorFactory::Build(Arm, numvecs,
true);
1562 if (correction_[i].is_null() || !correction_[i]->getMap()->isSameAs(*Adm))
1563 correction_[i] = MultiVectorFactory::Build(Adm, numvecs,
false);
1568 if(implicitTranspose_) {
1569 RCP<Operator> P = Levels_[i+1]->template Get< RCP<Operator> >(
"P");
1571 RCP<const Map> map = P->getDomainMap();
1572 if (coarseRhs_[i].is_null() || !coarseRhs_[i]->getMap()->isSameAs(*map))
1573 coarseRhs_[i] = MultiVectorFactory::Build(map, numvecs,
true);
1576 RCP<Operator> R = Levels_[i+1]->template Get< RCP<Operator> >(
"R");
1578 RCP<const Map> map = R->getRangeMap();
1579 if (coarseRhs_[i].is_null() || !coarseRhs_[i]->getMap()->isSameAs(*map))
1580 coarseRhs_[i] = MultiVectorFactory::Build(map, numvecs,
true);
1585 RCP<const Import> importer;
1586 if(Levels_[i+1]->IsAvailable(
"Importer"))
1587 importer = Levels_[i+1]->
template Get< RCP<const Import> >(
"Importer");
1588 if (doPRrebalance_ || importer.is_null()) {
1589 RCP<const Map> map = coarseRhs_[i]->getMap();
1590 if (coarseX_[i].is_null() || !coarseX_[i]->getMap()->isSameAs(*map))
1591 coarseX_[i] = MultiVectorFactory::Build(map, numvecs,
true);
1594 map = importer->getTargetMap();
1595 if (coarseImport_[i].is_null() || !coarseImport_[i]->getMap()->isSameAs(*map)) {
1596 coarseImport_[i] = MultiVectorFactory::Build(map, numvecs,
false);
1597 coarseX_[i] = MultiVectorFactory::Build(map,numvecs,
false);
1599 map = importer->getSourceMap();
1600 if (coarseExport_[i].is_null() || !coarseExport_[i]->getMap()->isSameAs(*map))
1601 coarseExport_[i] = MultiVectorFactory::Build(map, numvecs,
false);
1605 sizeOfAllocatedLevelMultiVectors_ = numvecs;
1609template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1611 if(sizeOfAllocatedLevelMultiVectors_==0)
return;
1612 residual_.resize(0);
1613 coarseRhs_.resize(0);
1615 coarseImport_.resize(0);
1616 coarseExport_.resize(0);
1617 correction_.resize(0);
1618 sizeOfAllocatedLevelMultiVectors_ = 0;
1622template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1624 const LO startLevel,
const ConvData& conv)
const
1626 return (startLevel == 0 && !isPreconditioner_ && (IsPrint(
Statistics1) || conv.
tol_ > 0));
1630template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1632 const Teuchos::Array<MagnitudeType>& residualNorm,
const MagnitudeType convergenceTolerance)
const
1636 if (convergenceTolerance > Teuchos::ScalarTraits<MagnitudeType>::zero())
1639 for (LO k = 0; k < residualNorm.size(); k++)
1640 if (residualNorm[k] >= convergenceTolerance)
1649 return convergenceStatus;
1653template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1655 const LO iteration,
const Teuchos::Array<MagnitudeType>& residualNorm)
const
1658 << std::setiosflags(std::ios::left)
1659 << std::setprecision(3) << std::setw(4) << iteration
1661 << std::setprecision(10) << residualNorm
1665template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1667 const Operator& A,
const MultiVector& X,
const MultiVector& B,
const LO iteration,
1670 Teuchos::Array<MagnitudeType> residualNorm;
1674 rate_ = currentResidualNorm / previousResidualNorm;
1675 previousResidualNorm = currentResidualNorm;
1678 PrintResidualHistory(iteration, residualNorm);
1680 return IsConverged(residualNorm, conv.
tol_);
virtual std::string description() const
Return a simple one-line description of this object.
Exception throws to report incompatible objects (like maps).
Exception throws to report errors in the internal logical of the program.
Base class for factories (e.g., R, P, and A_coarse).
Class that provides default factories within Needs class.
virtual void Clean() const
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
void AddLevel(const RCP< Level > &level)
Add a level at the end of the hierarchy.
double GetSmootherComplexity() const
void Write(const LO &start=-1, const LO &end=-1, const std::string &suffix="")
Print matrices in the multigrid hierarchy to file.
RCP< Level > & GetLevel(const int levelID=0)
Retrieve a certain level from hierarchy.
void CheckLevel(Level &level, int levelID)
Helper function.
std::string description() const
Return a simple one-line description of this object.
void IsPreconditioner(const bool flag)
Array< RCP< Level > > Levels_
Container for Level objects.
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.
STS::magnitudeType MagnitudeType
void describe(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the Hierarchy with some verbosity level to a FancyOStream object.
ConvergenceStatus IsConverged(const Teuchos::Array< MagnitudeType > &residualNorm, const MagnitudeType convergenceTolerance) const
Decide if the multigrid iteration is converged.
void DeleteLevelMultiVectors()
ConvergenceStatus Iterate(const MultiVector &B, MultiVector &X, ConvData conv=ConvData(), bool InitialGuessIsZero=false, LO startLevel=0)
Apply the multigrid preconditioner.
void DumpCurrentGraph(int level) const
void SetMatvecParams(RCP< ParameterList > matvecParams)
Xpetra::UnderlyingLib lib_
Epetra/Tpetra mode.
void Clear(int startLevel=0)
Clear impermanent data from previous setup.
bool IsCalculationOfResidualRequired(const LO startLevel, const ConvData &conv) const
Decide if the residual needs to be computed.
ConvergenceStatus ComputeResidualAndPrintHistory(const Operator &A, const MultiVector &X, const MultiVector &B, const LO iteration, const LO startLevel, const ConvData &conv, MagnitudeType &previousResidualNorm)
Compute the residual norm and print it depending on the verbosity level.
double GetOperatorComplexity() const
void PrintResidualHistory(const LO iteration, const Teuchos::Array< MagnitudeType > &residualNorm) const
Print residualNorm for this iteration to the screen.
void AllocateLevelMultiVectors(int numvecs, bool forceMapCheck=false)
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.
void Delete(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Call Level::Delete(ename, factory) for each level of the Hierarchy.
int GetGlobalNumLevels() const
void Keep(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Call Level::Keep(ename, factory) for each level of the Hierarchy.
void AddKeepFlag(const std::string &ename, const FactoryBase *factory=NoFactory::get(), KeepType keep=MueLu::Keep)
Call Level::AddKeepFlag for each level of the Hierarchy.
void AddNewLevel()
Add a new level at the end of the hierarchy.
void RemoveKeepFlag(const std::string &ename, const FactoryBase *factory, KeepType keep=MueLu::All)
Call Level::RemoveKeepFlag for each level of the Hierarchy.
void ReplaceCoordinateMap(Level &level)
Class that holds all level-specific information.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
void SetComm(RCP< const Teuchos::Comm< int > > const &comm)
void setlib(Xpetra::UnderlyingLib lib2)
RCP< Level > & GetPreviousLevel()
Previous level.
void Release(const FactoryBase &factory)
Decrement the storage counter for all the inputs of a factory.
RCP< const Teuchos::Comm< int > > GetComm() const
int GetLevelID() const
Return level number.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
void Request(const FactoryBase &factory)
Increment the storage counter for all the inputs of a factory.
Xpetra::UnderlyingLib lib()
Timer to be used in non-factories.
static const NoFactory * get()
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
An exception safe way to call the method 'Level::SetFactoryManager()'.
Integrates Teuchos::TimeMonitor with MueLu verbosity system.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
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)
static RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &RHS)
static void SetRandomSeed(const Teuchos::Comm< int > &comm)
Namespace for MueLu classes and methods.
@ Warnings0
Important warning messages (one line)
@ Statistics2
Print even more statistics.
@ Warnings
Print all warning messages.
@ Statistics1
Print more statistics.
@ Timings0
High level timing information (use Teuchos::TimeMonitor::summarize() to print)
@ Runtime0
One-liner description of what is happening.
@ Runtime1
Description of what is happening (more verbose)
@ Warnings1
Additional warnings.
@ Statistics0
Print statistics that do not involve significant additional computation.
@ Parameters1
Print class parameters (more parameters, more verbose)
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
VerbLevel toMueLuVerbLevel(const Teuchos::EVerbosityLevel verbLevel)
Translate Teuchos verbosity level to MueLu verbosity level.
Data struct for defining stopping criteria of multigrid iteration.