MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_Hierarchy_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // MueLu: A package for multigrid based preconditioning
4 //
5 // Copyright 2012 NTESS and the MueLu contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef MUELU_HIERARCHY_DEF_HPP
11 #define MUELU_HIERARCHY_DEF_HPP
12 
13 #include <time.h>
14 
15 #include <algorithm>
16 #include <sstream>
17 
18 #include <Xpetra_Matrix.hpp>
19 #include <Xpetra_MultiVectorFactory.hpp>
20 #include <Xpetra_Operator.hpp>
21 #include <Xpetra_IO.hpp>
22 
23 #include "MueLu_Hierarchy_decl.hpp"
24 
25 #include "MueLu_FactoryManager.hpp"
26 #include "MueLu_HierarchyUtils.hpp"
27 #include "MueLu_TopRAPFactory.hpp"
28 #include "MueLu_TopSmootherFactory.hpp"
29 #include "MueLu_Level.hpp"
30 #include "MueLu_Monitor.hpp"
31 #include "MueLu_PerfUtils.hpp"
32 #include "MueLu_PFactory.hpp"
33 #include "MueLu_SmootherFactory.hpp"
34 #include "MueLu_SmootherBase.hpp"
35 
36 #include "Teuchos_TimeMonitor.hpp"
37 
38 namespace MueLu {
39 
40 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
41 Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Hierarchy()
42  : maxCoarseSize_(GetDefaultMaxCoarseSize())
43  , implicitTranspose_(GetDefaultImplicitTranspose())
44  , fuseProlongationAndUpdate_(GetDefaultFuseProlongationAndUpdate())
45  , doPRrebalance_(GetDefaultPRrebalance())
46  , doPRViaCopyrebalance_(false)
47  , isPreconditioner_(true)
48  , Cycle_(GetDefaultCycle())
49  , WCycleStartLevel_(0)
50  , scalingFactor_(Teuchos::ScalarTraits<double>::one())
51  , lib_(Xpetra::UseTpetra)
52  , isDumpingEnabled_(false)
53  , dumpLevel_(-2)
54  , rate_(-1)
55  , sizeOfAllocatedLevelMultiVectors_(0) {
56  AddLevel(rcp(new Level));
57 }
58 
59 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
60 Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Hierarchy(const std::string& label)
61  : Hierarchy() {
62  SetLabel(label);
63  setObjectLabel(label);
64  Levels_[0]->setObjectLabel(label);
65 }
66 
67 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
68 Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Hierarchy(const RCP<Matrix>& A)
69  : maxCoarseSize_(GetDefaultMaxCoarseSize())
70  , implicitTranspose_(GetDefaultImplicitTranspose())
71  , fuseProlongationAndUpdate_(GetDefaultFuseProlongationAndUpdate())
72  , doPRrebalance_(GetDefaultPRrebalance())
73  , doPRViaCopyrebalance_(false)
74  , isPreconditioner_(true)
75  , Cycle_(GetDefaultCycle())
76  , WCycleStartLevel_(0)
77  , scalingFactor_(Teuchos::ScalarTraits<double>::one())
78  , isDumpingEnabled_(false)
79  , dumpLevel_(-2)
80  , rate_(-1)
81  , sizeOfAllocatedLevelMultiVectors_(0) {
82  lib_ = A->getDomainMap()->lib();
83 
84  RCP<Level> Finest = rcp(new Level);
85  AddLevel(Finest);
86 
87  Finest->Set("A", A);
88 }
89 
90 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
91 Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Hierarchy(const RCP<Matrix>& A, const std::string& label)
92  : Hierarchy(A) {
93  setObjectLabel(label);
94  Levels_[0]->setObjectLabel(label);
95 }
96 
97 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
98 Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::~Hierarchy() = default;
99 
100 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
101 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::AddLevel(const RCP<Level>& level) {
102  int levelID = LastLevelID() + 1; // ID of the inserted level
103 
104  if (level->GetLevelID() != -1 && (level->GetLevelID() != levelID))
105  GetOStream(Warnings1) << "Hierarchy::AddLevel(): Level with ID=" << level->GetLevelID() << " have been added at the end of the hierarchy\n but its ID have been redefined"
106  << " because last level ID of the hierarchy was " << LastLevelID() << "." << std::endl;
107 
108  Levels_.push_back(level);
109  level->SetLevelID(levelID);
110  level->setlib(lib_);
111 
112  level->SetPreviousLevel((levelID == 0) ? Teuchos::null : Levels_[LastLevelID() - 1]);
113  level->setObjectLabel(this->getObjectLabel());
114 }
115 
116 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
117 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::AddNewLevel() {
118  RCP<Level> newLevel = Levels_[LastLevelID()]->Build(); // new coarse level, using copy constructor
119  newLevel->setlib(lib_);
120  this->AddLevel(newLevel); // add to hierarchy
121 }
122 
123 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
124 RCP<Level>& Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetLevel(const int levelID) {
125  TEUCHOS_TEST_FOR_EXCEPTION(levelID < 0 || levelID > LastLevelID(), Exceptions::RuntimeError,
126  "MueLu::Hierarchy::GetLevel(): invalid input parameter value: LevelID = " << levelID);
127  return Levels_[levelID];
128 }
129 
130 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
131 int Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetNumLevels() const {
132  return Levels_.size();
133 }
134 
135 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
136 int Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetGlobalNumLevels() const {
137  RCP<Operator> A = Levels_[0]->template Get<RCP<Operator>>("A");
138  RCP<const Teuchos::Comm<int>> comm = A->getDomainMap()->getComm();
139 
140  int numLevels = GetNumLevels();
141  int numGlobalLevels;
142  Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, numLevels, Teuchos::ptr(&numGlobalLevels));
143 
144  return numGlobalLevels;
145 }
146 
147 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
148 double Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetOperatorComplexity() const {
149  double totalNnz = 0, lev0Nnz = 1;
150  for (int i = 0; i < GetNumLevels(); ++i) {
151  TEUCHOS_TEST_FOR_EXCEPTION(!(Levels_[i]->IsAvailable("A")), Exceptions::RuntimeError,
152  "Operator complexity cannot be calculated because A is unavailable on level " << i);
153  RCP<Operator> A = Levels_[i]->template Get<RCP<Operator>>("A");
154  if (A.is_null())
155  break;
156 
157  RCP<Matrix> Am = rcp_dynamic_cast<Matrix>(A);
158  if (Am.is_null()) {
159  GetOStream(Warnings0) << "Some level operators are not matrices, operator complexity calculation aborted" << std::endl;
160  return 0.0;
161  }
162 
163  totalNnz += as<double>(Am->getGlobalNumEntries());
164  if (i == 0)
165  lev0Nnz = totalNnz;
166  }
167  return totalNnz / lev0Nnz;
168 }
169 
170 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
171 double Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetSmootherComplexity() const {
172  double node_sc = 0, global_sc = 0;
173  double a0_nnz = 0;
174  const size_t INVALID = Teuchos::OrdinalTraits<size_t>::invalid();
175  // Get cost of fine matvec
176  if (GetNumLevels() <= 0) return -1.0;
177  if (!Levels_[0]->IsAvailable("A")) return -1.0;
178 
179  RCP<Operator> A = Levels_[0]->template Get<RCP<Operator>>("A");
180  if (A.is_null()) return -1.0;
181  RCP<Matrix> Am = rcp_dynamic_cast<Matrix>(A);
182  if (Am.is_null()) return -1.0;
183  a0_nnz = as<double>(Am->getGlobalNumEntries());
184 
185  // Get smoother complexity at each level
186  for (int i = 0; i < GetNumLevels(); ++i) {
187  size_t level_sc = 0;
188  if (!Levels_[i]->IsAvailable("PreSmoother")) continue;
189  RCP<SmootherBase> S = Levels_[i]->template Get<RCP<SmootherBase>>("PreSmoother");
190  if (S.is_null()) continue;
191  level_sc = S->getNodeSmootherComplexity();
192  if (level_sc == INVALID) {
193  global_sc = -1.0;
194  break;
195  }
196 
197  node_sc += as<double>(level_sc);
198  }
199 
200  double min_sc = 0.0;
201  RCP<const Teuchos::Comm<int>> comm = A->getDomainMap()->getComm();
202  Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, node_sc, Teuchos::ptr(&global_sc));
203  Teuchos::reduceAll(*comm, Teuchos::REDUCE_MIN, node_sc, Teuchos::ptr(&min_sc));
204 
205  if (min_sc < 0.0)
206  return -1.0;
207  else
208  return global_sc / a0_nnz;
209 }
210 
211 // Coherence checks todo in Setup() (using an helper function):
212 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
213 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::CheckLevel(Level& level, int levelID) {
214  TEUCHOS_TEST_FOR_EXCEPTION(level.lib() != lib_, Exceptions::RuntimeError,
215  "MueLu::Hierarchy::CheckLevel(): wrong underlying linear algebra library.");
216  TEUCHOS_TEST_FOR_EXCEPTION(level.GetLevelID() != levelID, Exceptions::RuntimeError,
217  "MueLu::Hierarchy::CheckLevel(): wrong level ID");
218  TEUCHOS_TEST_FOR_EXCEPTION(levelID != 0 && level.GetPreviousLevel() != Levels_[levelID - 1], Exceptions::RuntimeError,
219  "MueLu::Hierarchy::Setup(): wrong level parent");
220 }
221 
222 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
223 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::SetMatvecParams(RCP<ParameterList> matvecParams) {
224  for (int i = 0; i < GetNumLevels(); ++i) {
225  RCP<Level> level = Levels_[i];
226  if (level->IsAvailable("A")) {
227  RCP<Operator> Aop = level->Get<RCP<Operator>>("A");
228  RCP<Matrix> A = rcp_dynamic_cast<Matrix>(Aop);
229  if (!A.is_null()) {
230  RCP<const Import> xpImporter = A->getCrsGraph()->getImporter();
231  if (!xpImporter.is_null())
232  xpImporter->setDistributorParameters(matvecParams);
233  RCP<const Export> xpExporter = A->getCrsGraph()->getExporter();
234  if (!xpExporter.is_null())
235  xpExporter->setDistributorParameters(matvecParams);
236  }
237  }
238  if (level->IsAvailable("P")) {
239  RCP<Matrix> P = level->Get<RCP<Matrix>>("P");
240  RCP<const Import> xpImporter = P->getCrsGraph()->getImporter();
241  if (!xpImporter.is_null())
242  xpImporter->setDistributorParameters(matvecParams);
243  RCP<const Export> xpExporter = P->getCrsGraph()->getExporter();
244  if (!xpExporter.is_null())
245  xpExporter->setDistributorParameters(matvecParams);
246  }
247  if (level->IsAvailable("R")) {
248  RCP<Matrix> R = level->Get<RCP<Matrix>>("R");
249  RCP<const Import> xpImporter = R->getCrsGraph()->getImporter();
250  if (!xpImporter.is_null())
251  xpImporter->setDistributorParameters(matvecParams);
252  RCP<const Export> xpExporter = R->getCrsGraph()->getExporter();
253  if (!xpExporter.is_null())
254  xpExporter->setDistributorParameters(matvecParams);
255  }
256  if (level->IsAvailable("Importer")) {
257  RCP<const Import> xpImporter = level->Get<RCP<const Import>>("Importer");
258  if (!xpImporter.is_null())
259  xpImporter->setDistributorParameters(matvecParams);
260  }
261  }
262 }
263 
264 // The function uses three managers: fine, coarse and next coarse
265 // We construct the data for the coarse level, and do requests for the next coarse
266 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
267 bool Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Setup(int coarseLevelID,
268  const RCP<const FactoryManagerBase> fineLevelManager,
269  const RCP<const FactoryManagerBase> coarseLevelManager,
270  const RCP<const FactoryManagerBase> nextLevelManager) {
271  // Use PrintMonitor/TimerMonitor instead of just a FactoryMonitor to print "Level 0" instead of Hierarchy(0)
272  // Print is done after the requests for next coarse level
273 
274  TEUCHOS_TEST_FOR_EXCEPTION(LastLevelID() < coarseLevelID, Exceptions::RuntimeError,
275  "MueLu::Hierarchy:Setup(): level " << coarseLevelID << " (specified by coarseLevelID argument) "
276  "must be built before calling this function.");
277 
278  Level& level = *Levels_[coarseLevelID];
279 
280  bool useStackedTimer = !Teuchos::TimeMonitor::stackedTimerNameIsDefault();
281 
282  std::string label = FormattingHelper::getColonLabel(level.getObjectLabel());
283  RCP<TimeMonitor> m1;
284  if (!useStackedTimer)
285  m1 = rcp(new TimeMonitor(*this, label + this->ShortClassName() + ": " + "Setup (total)"));
286  TimeMonitor m2(*this, label + this->ShortClassName() + ": " + "Setup" + " (total, level=" + Teuchos::toString(coarseLevelID) + ")");
287 
288  // TODO: pass coarseLevelManager by reference
289  TEUCHOS_TEST_FOR_EXCEPTION(coarseLevelManager == Teuchos::null, Exceptions::RuntimeError,
290  "MueLu::Hierarchy::Setup(): argument coarseLevelManager cannot be null");
291 
292  typedef MueLu::TopRAPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node> TopRAPFactory;
293  typedef MueLu::TopSmootherFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node> TopSmootherFactory;
294 
295  if (levelManagers_.size() < coarseLevelID + 1)
296  levelManagers_.resize(coarseLevelID + 1);
297  levelManagers_[coarseLevelID] = coarseLevelManager;
298 
299  bool isFinestLevel = (fineLevelManager.is_null());
300  bool isLastLevel = (nextLevelManager.is_null());
301 
302  int oldRank = -1;
303  if (isFinestLevel) {
304  RCP<Operator> A = level.Get<RCP<Operator>>("A");
305  RCP<const Map> domainMap = A->getDomainMap();
306  RCP<const Teuchos::Comm<int>> comm = domainMap->getComm();
307 
308  // Initialize random seed for reproducibility
309  Utilities::SetRandomSeed(*comm);
310 
311  // Record the communicator on the level (used for timers sync)
312  level.SetComm(comm);
313  oldRank = SetProcRankVerbose(comm->getRank());
314 
315  // Set the Hierarchy library to match that of the finest level matrix,
316  // even if it was already set
317  lib_ = domainMap->lib();
318  level.setlib(lib_);
319 
320  } else {
321  // Permeate library to a coarser level
322  level.setlib(lib_);
323 
324  Level& prevLevel = *Levels_[coarseLevelID - 1];
325  oldRank = SetProcRankVerbose(prevLevel.GetComm()->getRank());
326  }
327 
328  CheckLevel(level, coarseLevelID);
329 
330  // Attach FactoryManager to the fine level
331  RCP<SetFactoryManager> SFMFine;
332  if (!isFinestLevel)
333  SFMFine = rcp(new SetFactoryManager(Levels_[coarseLevelID - 1], fineLevelManager));
334 
335  if (isFinestLevel && Levels_[coarseLevelID]->IsAvailable("Coordinates"))
336  ReplaceCoordinateMap(*Levels_[coarseLevelID]);
337 
338  // Attach FactoryManager to the coarse level
339  SetFactoryManager SFMCoarse(Levels_[coarseLevelID], coarseLevelManager);
340 
341  if (isDumpingEnabled_ && (dumpLevel_ == 0 || dumpLevel_ == -1) && coarseLevelID == 1)
342  DumpCurrentGraph(0);
343 
344  RCP<TopSmootherFactory> coarseFact;
345  RCP<TopSmootherFactory> smootherFact = rcp(new TopSmootherFactory(coarseLevelManager, "Smoother"));
346 
347  int nextLevelID = coarseLevelID + 1;
348 
349  RCP<SetFactoryManager> SFMNext;
350  if (isLastLevel == false) {
351  // We are not at the coarsest level, so there is going to be another level ("next coarse") after this one ("coarse")
352  if (nextLevelID > LastLevelID())
353  AddNewLevel();
354  CheckLevel(*Levels_[nextLevelID], nextLevelID);
355 
356  // Attach FactoryManager to the next level (level after coarse)
357  SFMNext = rcp(new SetFactoryManager(Levels_[nextLevelID], nextLevelManager));
358  Levels_[nextLevelID]->Request(TopRAPFactory(coarseLevelManager, nextLevelManager));
359 
360  // Do smoother requests here. We don't know whether this is going to be
361  // the coarsest level or not, but we need to DeclareInput before we call
362  // coarseRAPFactory.Build(), otherwise some stuff may be erased after
363  // level releases
364  level.Request(*smootherFact);
365 
366  } else {
367  // Similar to smoother above, do the coarse solver request here. We don't
368  // know whether this is going to be the coarsest level or not, but we
369  // need to DeclareInput before we call coarseRAPFactory.Build(),
370  // otherwise some stuff may be erased after level releases. This is
371  // actually evident on ProjectorSmoother. It requires both "A" and
372  // "Nullspace". However, "Nullspace" is erased after all releases, so if
373  // we call the coarse factory request after RAP build we would not have
374  // any data, and cannot get it as we don't have previous managers. The
375  // typical trace looks like this:
376  //
377  // MueLu::Level(0)::GetFactory(Aggregates, 0): No FactoryManager
378  // during request for data " Aggregates" on level 0 by factory TentativePFactory
379  // during request for data " P" on level 1 by factory EminPFactory
380  // during request for data " P" on level 1 by factory TransPFactory
381  // during request for data " R" on level 1 by factory RAPFactory
382  // during request for data " A" on level 1 by factory TentativePFactory
383  // during request for data " Nullspace" on level 2 by factory NullspaceFactory
384  // during request for data " Nullspace" on level 2 by factory NullspacePresmoothFactory
385  // during request for data " Nullspace" on level 2 by factory ProjectorSmoother
386  // during request for data " PreSmoother" on level 2 by factory NoFactory
387  if (coarseFact.is_null())
388  coarseFact = rcp(new TopSmootherFactory(coarseLevelManager, "CoarseSolver"));
389  level.Request(*coarseFact);
390  }
391 
392  GetOStream(Runtime0) << std::endl;
393  PrintMonitor m0(*this, "Level " + Teuchos::toString(coarseLevelID), static_cast<MsgType>(Runtime0 | Test));
394 
395  // Build coarse level hierarchy
396  RCP<Operator> Ac = Teuchos::null;
397  TopRAPFactory coarseRAPFactory(fineLevelManager, coarseLevelManager);
398 
399  if (level.IsAvailable("A")) {
400  Ac = level.Get<RCP<Operator>>("A");
401  } else if (!isFinestLevel) {
402  // We only build here, the release is done later
403  coarseRAPFactory.Build(*level.GetPreviousLevel(), level);
404  }
405 
406  bool setLastLevelviaMaxCoarseSize = false;
407  if (level.IsAvailable("A"))
408  Ac = level.Get<RCP<Operator>>("A");
409  RCP<Matrix> Acm = rcp_dynamic_cast<Matrix>(Ac);
410 
411  // Record the communicator on the level
412  if (!Ac.is_null())
413  level.SetComm(Ac->getDomainMap()->getComm());
414 
415  // Test if we reach the end of the hierarchy
416  bool isOrigLastLevel = isLastLevel;
417  if (isLastLevel) {
418  // Last level as we have achieved the max limit
419  isLastLevel = true;
420 
421  } else if (Ac.is_null()) {
422  // Last level for this processor, as it does not belong to the next
423  // subcommunicator. Other processors may continue working on the
424  // hierarchy
425  isLastLevel = true;
426 
427  } else {
428  if (!Acm.is_null() && Acm->getGlobalNumRows() <= maxCoarseSize_) {
429  // Last level as the size of the coarse matrix became too small
430  GetOStream(Runtime0) << "Max coarse size (<= " << maxCoarseSize_ << ") achieved" << std::endl;
431  isLastLevel = true;
432  if (Acm->getGlobalNumRows() != 0) setLastLevelviaMaxCoarseSize = true;
433  }
434  }
435 
436  if (!Ac.is_null() && !isFinestLevel) {
437  RCP<Operator> A = Levels_[coarseLevelID - 1]->template Get<RCP<Operator>>("A");
438  RCP<Matrix> Am = rcp_dynamic_cast<Matrix>(A);
439 
440  const double maxCoarse2FineRatio = 0.8;
441  if (!Acm.is_null() && !Am.is_null() && Acm->getGlobalNumRows() > maxCoarse2FineRatio * Am->getGlobalNumRows()) {
442  // We could abort here, but for now we simply notify user.
443  // Couple of additional points:
444  // - if repartitioning is delayed until level K, but the aggregation
445  // procedure stagnates between levels K-1 and K. In this case,
446  // repartitioning could enable faster coarsening once again, but the
447  // hierarchy construction will abort due to the stagnation check.
448  // - if the matrix is small enough, we could move it to one processor.
449  GetOStream(Warnings0) << "Aggregation stagnated. Please check your matrix and/or adjust your configuration file."
450  << "Possible fixes:\n"
451  << " - reduce the maximum number of levels\n"
452  << " - enable repartitioning\n"
453  << " - increase the minimum coarse size." << std::endl;
454  }
455  }
456 
457  if (isLastLevel) {
458  if (!isOrigLastLevel) {
459  // We did not expect to finish this early so we did request a smoother.
460  // We need a coarse solver instead. Do the magic.
461  level.Release(*smootherFact);
462  if (coarseFact.is_null())
463  coarseFact = rcp(new TopSmootherFactory(coarseLevelManager, "CoarseSolver"));
464  level.Request(*coarseFact);
465  }
466 
467  // Do the actual build, if we have any data.
468  // NOTE: this is not a great check, we may want to call Build() regardless.
469  if (!Ac.is_null())
470  coarseFact->Build(level);
471 
472  // Once the dirty deed is done, release stuff. The smoother has already
473  // been released.
474  level.Release(*coarseFact);
475 
476  } else {
477  // isLastLevel = false => isOrigLastLevel = false, meaning that we have
478  // requested the smoother. Now we need to build it and to release it.
479  // We don't need to worry about the coarse solver, as we didn't request it.
480  if (!Ac.is_null())
481  smootherFact->Build(level);
482 
483  level.Release(*smootherFact);
484  }
485 
486  if (isLastLevel == true) {
487  int actualNumLevels = nextLevelID;
488  if (isOrigLastLevel == false) {
489  // Earlier in the function, we constructed the next coarse level, and requested data for the that level,
490  // assuming that we are not at the coarsest level. Now, we changed our mind, so we have to release those.
491  Levels_[nextLevelID]->Release(TopRAPFactory(coarseLevelManager, nextLevelManager));
492 
493  // We truncate/resize the hierarchy and possibly remove the last created level if there is
494  // something wrong with it as indicated by its P not being valid. This might happen
495  // if the global number of aggregates turns out to be zero
496 
497  if (!setLastLevelviaMaxCoarseSize) {
498  if (Levels_[nextLevelID - 1]->IsAvailable("P")) {
499  if (Levels_[nextLevelID - 1]->template Get<RCP<Matrix>>("P") == Teuchos::null) actualNumLevels = nextLevelID - 1;
500  } else
501  actualNumLevels = nextLevelID - 1;
502  }
503  }
504  if (actualNumLevels == nextLevelID - 1) {
505  // Didn't expect to finish early so we requested smoother but need coarse solver instead.
506  Levels_[nextLevelID - 2]->Release(*smootherFact);
507 
508  if (Levels_[nextLevelID - 2]->IsAvailable("PreSmoother")) Levels_[nextLevelID - 2]->RemoveKeepFlag("PreSmoother", NoFactory::get());
509  if (Levels_[nextLevelID - 2]->IsAvailable("PostSmoother")) Levels_[nextLevelID - 2]->RemoveKeepFlag("PostSmoother", NoFactory::get());
510  if (coarseFact.is_null())
511  coarseFact = rcp(new TopSmootherFactory(coarseLevelManager, "CoarseSolver"));
512  Levels_[nextLevelID - 2]->Request(*coarseFact);
513  if (!(Levels_[nextLevelID - 2]->template Get<RCP<Matrix>>("A").is_null()))
514  coarseFact->Build(*(Levels_[nextLevelID - 2]));
515  Levels_[nextLevelID - 2]->Release(*coarseFact);
516  }
517  Levels_.resize(actualNumLevels);
518  }
519 
520  // I think this is the proper place for graph so that it shows every dependence
521  if (isDumpingEnabled_ && ((dumpLevel_ > 0 && coarseLevelID == dumpLevel_) || dumpLevel_ == -1))
522  DumpCurrentGraph(coarseLevelID);
523 
524  if (!isFinestLevel) {
525  // Release the hierarchy data
526  // We release so late to help blocked solvers, as the smoothers for them need A blocks
527  // which we construct in RAPFactory
528  level.Release(coarseRAPFactory);
529  }
530 
531  if (oldRank != -1)
532  SetProcRankVerbose(oldRank);
533 
534  return isLastLevel;
535 }
536 
537 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
538 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::SetupRe() {
539  int numLevels = Levels_.size();
540  TEUCHOS_TEST_FOR_EXCEPTION(levelManagers_.size() != numLevels, Exceptions::RuntimeError,
541  "Hierarchy::SetupRe: " << Levels_.size() << " levels, but " << levelManagers_.size() << " level factory managers");
542 
543  const int startLevel = 0;
544  Clear(startLevel);
545 
546 #ifdef HAVE_MUELU_DEBUG
547  // Reset factories' data used for debugging
548  for (int i = 0; i < numLevels; i++)
549  levelManagers_[i]->ResetDebugData();
550 
551 #endif
552 
553  int levelID;
554  for (levelID = startLevel; levelID < numLevels;) {
555  bool r = Setup(levelID,
556  (levelID != 0 ? levelManagers_[levelID - 1] : Teuchos::null),
557  levelManagers_[levelID],
558  (levelID + 1 != numLevels ? levelManagers_[levelID + 1] : Teuchos::null));
559  levelID++;
560  if (r) break;
561  }
562  // We may construct fewer levels for some reason, make sure we continue
563  // doing that in the future
564  Levels_.resize(levelID);
565  levelManagers_.resize(levelID);
566 
567  int sizeOfVecs = sizeOfAllocatedLevelMultiVectors_;
568 
569  AllocateLevelMultiVectors(sizeOfVecs, true);
570 
571  // since the # of levels, etc. may have changed, force re-determination of description during next call to description()
572  ResetDescription();
573 
574  describe(GetOStream(Statistics0), GetVerbLevel());
575 
576  CheckForEmptySmoothersAndCoarseSolve();
577 }
578 
579 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
580 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Setup(const FactoryManagerBase& manager, int startLevel, int numDesiredLevels) {
581  // Use MueLu::BaseClass::description() to avoid printing "{numLevels = 1}" (numLevels is increasing...)
582  PrintMonitor m0(*this, "Setup (" + this->MueLu::BaseClass::description() + ")", Runtime0);
583 
584  Clear(startLevel);
585 
586  // Check Levels_[startLevel] exists.
587  TEUCHOS_TEST_FOR_EXCEPTION(Levels_.size() <= startLevel, Exceptions::RuntimeError,
588  "MueLu::Hierarchy::Setup(): fine level (" << startLevel << ") does not exist");
589 
590  TEUCHOS_TEST_FOR_EXCEPTION(numDesiredLevels <= 0, Exceptions::RuntimeError,
591  "Constructing non-positive (" << numDesiredLevels << ") number of levels does not make sense.");
592 
593  // Check for fine level matrix A
594  TEUCHOS_TEST_FOR_EXCEPTION(!Levels_[startLevel]->IsAvailable("A"), Exceptions::RuntimeError,
595  "MueLu::Hierarchy::Setup(): fine level (" << startLevel << ") has no matrix A! "
596  "Set fine level matrix A using Level.Set()");
597 
598  RCP<Operator> A = Levels_[startLevel]->template Get<RCP<Operator>>("A");
599  lib_ = A->getDomainMap()->lib();
600 
601  if (IsPrint(Statistics2)) {
602  RCP<Matrix> Amat = rcp_dynamic_cast<Matrix>(A);
603 
604  if (!Amat.is_null()) {
605  RCP<ParameterList> params = rcp(new ParameterList());
606  params->set("printLoadBalancingInfo", true);
607  params->set("printCommInfo", true);
608 
609  GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*Amat, "A0", params);
610  } else {
611  GetOStream(Warnings1) << "Fine level operator is not a matrix, statistics are not available" << std::endl;
612  }
613  }
614 
615  RCP<const FactoryManagerBase> rcpmanager = rcpFromRef(manager);
616 
617  const int lastLevel = startLevel + numDesiredLevels - 1;
618  GetOStream(Runtime0) << "Setup loop: startLevel = " << startLevel << ", lastLevel = " << lastLevel
619  << " (stop if numLevels = " << numDesiredLevels << " or Ac.size() < " << maxCoarseSize_ << ")" << std::endl;
620 
621  // Setup multigrid levels
622  int iLevel = 0;
623  if (numDesiredLevels == 1) {
624  iLevel = 0;
625  Setup(startLevel, Teuchos::null, rcpmanager, Teuchos::null); // setup finest==coarsest level (first and last managers are Teuchos::null)
626 
627  } else {
628  bool bIsLastLevel = Setup(startLevel, Teuchos::null, rcpmanager, rcpmanager); // setup finest level (level 0) (first manager is Teuchos::null)
629  if (bIsLastLevel == false) {
630  for (iLevel = startLevel + 1; iLevel < lastLevel; iLevel++) {
631  bIsLastLevel = Setup(iLevel, rcpmanager, rcpmanager, rcpmanager); // setup intermediate levels
632  if (bIsLastLevel == true)
633  break;
634  }
635  if (bIsLastLevel == false)
636  Setup(lastLevel, rcpmanager, rcpmanager, Teuchos::null); // setup coarsest level (last manager is Teuchos::null)
637  }
638  }
639 
640  // TODO: some check like this should be done at the beginning of the routine
641  TEUCHOS_TEST_FOR_EXCEPTION(iLevel != Levels_.size() - 1, Exceptions::RuntimeError,
642  "MueLu::Hierarchy::Setup(): number of level");
643 
644  // TODO: this is not exception safe: manager will still hold default
645  // factories if you exit this function with an exception
646  manager.Clean();
647 
648  describe(GetOStream(Statistics0), GetVerbLevel());
649 
650  CheckForEmptySmoothersAndCoarseSolve();
651 }
652 
653 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
654 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::CheckForEmptySmoothersAndCoarseSolve() {
655  for (LO levelNo = 0; levelNo < GetNumLevels(); ++levelNo) {
656  auto level = Levels_[levelNo];
657  if ((level->IsAvailable("A") && !level->template Get<RCP<Operator>>("A").is_null()) && (!level->IsAvailable("PreSmoother")) && (!level->IsAvailable("PostSmoother"))) {
658  GetOStream(Warnings1) << "No " << (levelNo == as<LO>(Levels_.size()) - 1 ? "coarse grid solver" : "smoother") << " on level " << level->GetLevelID() << std::endl;
659  }
660  }
661 }
662 
663 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
664 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Clear(int startLevel) {
665  if (startLevel < GetNumLevels())
666  GetOStream(Runtime0) << "Clearing old data (if any)" << std::endl;
667 
668  for (int iLevel = startLevel; iLevel < GetNumLevels(); iLevel++)
669  Levels_[iLevel]->Clear();
670 }
671 
672 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
673 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ExpertClear() {
674  GetOStream(Runtime0) << "Clearing old data (expert)" << std::endl;
675  for (int iLevel = 0; iLevel < GetNumLevels(); iLevel++)
676  Levels_[iLevel]->ExpertClear();
677 }
678 
679 #if defined(HAVE_MUELU_EXPERIMENTAL) && defined(HAVE_MUELU_ADDITIVE_VARIANT)
680 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
681 ConvergenceStatus Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Iterate(const MultiVector& B, MultiVector& X, ConvData conv,
682  bool InitialGuessIsZero, LO startLevel) {
683  LO nIts = conv.maxIts_;
684  MagnitudeType tol = conv.tol_;
685 
686  std::string prefix = this->ShortClassName() + ": ";
687  std::string levelSuffix = " (level=" + toString(startLevel) + ")";
688  std::string levelSuffix1 = " (level=" + toString(startLevel + 1) + ")";
689 
690  using namespace Teuchos;
691  RCP<Time> CompTime = Teuchos::TimeMonitor::getNewCounter(prefix + "Computational Time (total)");
692  RCP<Time> Concurrent = Teuchos::TimeMonitor::getNewCounter(prefix + "Concurrent portion");
693  RCP<Time> ApplyR = Teuchos::TimeMonitor::getNewCounter(prefix + "R: Computational Time");
694  RCP<Time> ApplyPbar = Teuchos::TimeMonitor::getNewCounter(prefix + "Pbar: Computational Time");
695  RCP<Time> CompFine = Teuchos::TimeMonitor::getNewCounter(prefix + "Fine: Computational Time");
696  RCP<Time> CompCoarse = Teuchos::TimeMonitor::getNewCounter(prefix + "Coarse: Computational Time");
697  RCP<Time> ApplySum = Teuchos::TimeMonitor::getNewCounter(prefix + "Sum: Computational Time");
698  RCP<Time> Synchronize_beginning = Teuchos::TimeMonitor::getNewCounter(prefix + "Synchronize_beginning");
699  RCP<Time> Synchronize_center = Teuchos::TimeMonitor::getNewCounter(prefix + "Synchronize_center");
700  RCP<Time> Synchronize_end = Teuchos::TimeMonitor::getNewCounter(prefix + "Synchronize_end");
701 
702  RCP<Level> Fine = Levels_[0];
703  RCP<Level> Coarse;
704 
705  RCP<Operator> A = Fine->Get<RCP<Operator>>("A");
706  Teuchos::RCP<const Teuchos::Comm<int>> communicator = A->getDomainMap()->getComm();
707 
708  // Synchronize_beginning->start();
709  // communicator->barrier();
710  // Synchronize_beginning->stop();
711 
712  CompTime->start();
713 
714  SC one = STS::one(), zero = STS::zero();
715 
716  bool zeroGuess = InitialGuessIsZero;
717 
718  // ======= UPFRONT DEFINITION OF COARSE VARIABLES ===========
719 
720  // RCP<const Map> origMap;
721  RCP<Operator> P;
722  RCP<Operator> Pbar;
723  RCP<Operator> R;
724  RCP<MultiVector> coarseRhs, coarseX;
725  RCP<Operator> Ac;
726  RCP<SmootherBase> preSmoo_coarse, postSmoo_coarse;
727  bool emptyCoarseSolve = true;
728  RCP<MultiVector> coarseX_prolonged = MultiVectorFactory::Build(X.getMap(), X.getNumVectors(), true);
729 
730  RCP<const Import> importer;
731 
732  if (Levels_.size() > 1) {
733  Coarse = Levels_[1];
734  if (Coarse->IsAvailable("Importer"))
735  importer = Coarse->Get<RCP<const Import>>("Importer");
736 
737  R = Coarse->Get<RCP<Operator>>("R");
738  P = Coarse->Get<RCP<Operator>>("P");
739 
740  // if(Coarse->IsAvailable("Pbar"))
741  Pbar = Coarse->Get<RCP<Operator>>("Pbar");
742 
743  coarseRhs = MultiVectorFactory::Build(R->getRangeMap(), B.getNumVectors(), true);
744 
745  Ac = Coarse->Get<RCP<Operator>>("A");
746 
747  ApplyR->start();
748  R->apply(B, *coarseRhs, Teuchos::NO_TRANS, one, zero);
749  // P->apply(B, *coarseRhs, Teuchos::TRANS, one, zero);
750  ApplyR->stop();
751 
752  if (doPRrebalance_ || importer.is_null()) {
753  coarseX = MultiVectorFactory::Build(coarseRhs->getMap(), X.getNumVectors(), true);
754 
755  } else {
756  RCP<TimeMonitor> ITime = rcp(new TimeMonitor(*this, prefix + "Solve : import (total)", Timings0));
757  RCP<TimeMonitor> ILevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : import" + levelSuffix1, Timings0));
758 
759  // Import: range map of R --> domain map of rebalanced Ac (before subcomm replacement)
760  RCP<MultiVector> coarseTmp = MultiVectorFactory::Build(importer->getTargetMap(), coarseRhs->getNumVectors());
761  coarseTmp->doImport(*coarseRhs, *importer, Xpetra::INSERT);
762  coarseRhs.swap(coarseTmp);
763 
764  coarseX = MultiVectorFactory::Build(importer->getTargetMap(), X.getNumVectors(), true);
765  }
766 
767  if (Coarse->IsAvailable("PreSmoother"))
768  preSmoo_coarse = Coarse->Get<RCP<SmootherBase>>("PreSmoother");
769  if (Coarse->IsAvailable("PostSmoother"))
770  postSmoo_coarse = Coarse->Get<RCP<SmootherBase>>("PostSmoother");
771  }
772 
773  // ==========================================================
774 
775  MagnitudeType prevNorm = STS::magnitude(STS::one()), curNorm = STS::magnitude(STS::one());
776  rate_ = 1.0;
777 
778  for (LO i = 1; i <= nIts; i++) {
779 #ifdef HAVE_MUELU_DEBUG
780  if (A->getDomainMap()->isCompatible(*(X.getMap())) == false) {
781  std::ostringstream ss;
782  ss << "Level " << startLevel << ": level A's domain map is not compatible with X";
783  throw Exceptions::Incompatible(ss.str());
784  }
785 
786  if (A->getRangeMap()->isCompatible(*(B.getMap())) == false) {
787  std::ostringstream ss;
788  ss << "Level " << startLevel << ": level A's range map is not compatible with B";
789  throw Exceptions::Incompatible(ss.str());
790  }
791 #endif
792  }
793 
794  bool emptyFineSolve = true;
795 
796  RCP<MultiVector> fineX;
797  fineX = MultiVectorFactory::Build(X.getMap(), X.getNumVectors(), true);
798 
799  // Synchronize_center->start();
800  // communicator->barrier();
801  // Synchronize_center->stop();
802 
803  Concurrent->start();
804 
805  // NOTE: we need to check using IsAvailable before Get here to avoid building default smoother
806  if (Fine->IsAvailable("PreSmoother")) {
807  RCP<SmootherBase> preSmoo = Fine->Get<RCP<SmootherBase>>("PreSmoother");
808  CompFine->start();
809  preSmoo->Apply(*fineX, B, zeroGuess);
810  CompFine->stop();
811  emptyFineSolve = false;
812  }
813  if (Fine->IsAvailable("PostSmoother")) {
814  RCP<SmootherBase> postSmoo = Fine->Get<RCP<SmootherBase>>("PostSmoother");
815  CompFine->start();
816  postSmoo->Apply(*fineX, B, zeroGuess);
817  CompFine->stop();
818 
819  emptyFineSolve = false;
820  }
821  if (emptyFineSolve == true) {
822  // Fine grid smoother is identity
823  fineX->update(one, B, zero);
824  }
825 
826  if (Levels_.size() > 1) {
827  // NOTE: we need to check using IsAvailable before Get here to avoid building default smoother
828  if (Coarse->IsAvailable("PreSmoother")) {
829  CompCoarse->start();
830  preSmoo_coarse->Apply(*coarseX, *coarseRhs, zeroGuess);
831  CompCoarse->stop();
832  emptyCoarseSolve = false;
833  }
834  if (Coarse->IsAvailable("PostSmoother")) {
835  CompCoarse->start();
836  postSmoo_coarse->Apply(*coarseX, *coarseRhs, zeroGuess);
837  CompCoarse->stop();
838  emptyCoarseSolve = false;
839  }
840  if (emptyCoarseSolve == true) {
841  // Coarse operator is identity
842  coarseX->update(one, *coarseRhs, zero);
843  }
844  Concurrent->stop();
845  // Synchronize_end->start();
846  // communicator->barrier();
847  // Synchronize_end->stop();
848 
849  if (!doPRrebalance_ && !importer.is_null()) {
850  RCP<TimeMonitor> ITime = rcp(new TimeMonitor(*this, prefix + "Solve : export (total)", Timings0));
851  RCP<TimeMonitor> ILevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : export" + levelSuffix1, Timings0));
852 
853  // Import: range map of rebalanced Ac (before subcomm replacement) --> domain map of P
854  RCP<MultiVector> coarseTmp = MultiVectorFactory::Build(importer->getSourceMap(), coarseX->getNumVectors());
855  coarseTmp->doExport(*coarseX, *importer, Xpetra::INSERT);
856  coarseX.swap(coarseTmp);
857  }
858 
859  ApplyPbar->start();
860  Pbar->apply(*coarseX, *coarseX_prolonged, Teuchos::NO_TRANS, one, zero);
861  ApplyPbar->stop();
862  }
863 
864  ApplySum->start();
865  X.update(1.0, *fineX, 1.0, *coarseX_prolonged, 0.0);
866  ApplySum->stop();
867 
868  CompTime->stop();
869 
870  // communicator->barrier();
871 
872  return (tol > 0 ? ConvergenceStatus::Unconverged : ConvergenceStatus::Undefined);
873 }
874 #else
875 // ---------------------------------------- Iterate -------------------------------------------------------
876 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
877 ConvergenceStatus Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Iterate(const MultiVector& B, MultiVector& X, ConvData conv,
878  bool InitialGuessIsZero, LO startLevel) {
879  LO nIts = conv.maxIts_;
880  MagnitudeType tol = conv.tol_;
881 
882  // These timers work as follows. "iterateTime" records total time spent in
883  // iterate. "levelTime" records time on a per level basis. The label is
884  // crafted to mimic the per-level messages used in Monitors. Note that a
885  // given level is timed with a TimeMonitor instead of a Monitor or
886  // SubMonitor. This is mainly because I want to time each level
887  // separately, and Monitors/SubMonitors print "(total) xx yy zz" ,
888  // "(sub,total) xx yy zz", respectively, which is subject to
889  // misinterpretation. The per-level TimeMonitors are stopped/started
890  // manually before/after a recursive call to Iterate. A side artifact to
891  // this approach is that the counts for intermediate level timers are twice
892  // the counts for the finest and coarsest levels.
893 
894  RCP<Level> Fine = Levels_[startLevel];
895 
896  std::string label = FormattingHelper::getColonLabel(Fine->getObjectLabel());
897  std::string prefix = label + this->ShortClassName() + ": ";
898  std::string levelSuffix = " (level=" + toString(startLevel) + ")";
899  std::string levelSuffix1 = " (level=" + toString(startLevel + 1) + ")";
900 
901  bool useStackedTimer = !Teuchos::TimeMonitor::stackedTimerNameIsDefault();
902 
903  RCP<Monitor> iterateTime;
904  RCP<TimeMonitor> iterateTime1;
905  if (startLevel == 0)
906  iterateTime = rcp(new Monitor(*this, "Solve", label, (nIts == 1) ? None : Runtime0, Timings0));
907  else if (!useStackedTimer)
908  iterateTime1 = rcp(new TimeMonitor(*this, prefix + "Solve (total, level=" + toString(startLevel) + ")", Timings0));
909 
910  std::string iterateLevelTimeLabel = prefix + "Solve" + levelSuffix;
911  RCP<TimeMonitor> iterateLevelTime = rcp(new TimeMonitor(*this, iterateLevelTimeLabel, Timings0));
912 
913  bool zeroGuess = InitialGuessIsZero;
914 
915  RCP<Operator> A = Fine->Get<RCP<Operator>>("A");
916  using namespace Teuchos;
917  RCP<Time> CompCoarse = Teuchos::TimeMonitor::getNewCounter(prefix + "Coarse: Computational Time");
918 
919  if (A.is_null()) {
920  // This processor does not have any data for this process on coarser
921  // levels. This can only happen when there are multiple processors and
922  // we use repartitioning.
923  return ConvergenceStatus::Undefined;
924  }
925 
926  // If we switched the number of vectors, we'd need to reallocate here.
927  // If the number of vectors is unchanged, this is a noop.
928  // NOTE: We need to check against B because the tests in AllocateLevelMultiVectors
929  // will fail on Stokhos Scalar types (due to the so-called 'hidden dimension')
930  const BlockedMultiVector* Bblocked = dynamic_cast<const BlockedMultiVector*>(&B);
931  if (residual_.size() > startLevel &&
932  ((Bblocked && !Bblocked->isSameSize(*residual_[startLevel])) ||
933  (!Bblocked && !residual_[startLevel]->isSameSize(B))))
934  DeleteLevelMultiVectors();
935  AllocateLevelMultiVectors(X.getNumVectors());
936 
937  // Print residual information before iterating
938  typedef Teuchos::ScalarTraits<typename STS::magnitudeType> STM;
939  MagnitudeType prevNorm = STM::one();
940  rate_ = 1.0;
941  if (IsCalculationOfResidualRequired(startLevel, conv)) {
942  ConvergenceStatus convergenceStatus = ComputeResidualAndPrintHistory(*A, X, B, Teuchos::ScalarTraits<LO>::zero(), startLevel, conv, prevNorm);
943  if (convergenceStatus == MueLu::ConvergenceStatus::Converged)
944  return convergenceStatus;
945  }
946 
947  SC one = STS::one(), zero = STS::zero();
948  for (LO iteration = 1; iteration <= nIts; iteration++) {
949 #ifdef HAVE_MUELU_DEBUG
950 #if 0 // TODO fix me
951  if (A->getDomainMap()->isCompatible(*(X.getMap())) == false) {
952  std::ostringstream ss;
953  ss << "Level " << startLevel << ": level A's domain map is not compatible with X";
954  throw Exceptions::Incompatible(ss.str());
955  }
956 
957  if (A->getRangeMap()->isCompatible(*(B.getMap())) == false) {
958  std::ostringstream ss;
959  ss << "Level " << startLevel << ": level A's range map is not compatible with B";
960  throw Exceptions::Incompatible(ss.str());
961  }
962 #endif
963 #endif
964 
965  if (startLevel == as<LO>(Levels_.size()) - 1) {
966  // On the coarsest level, we do either smoothing (if defined) or a direct solve.
967  RCP<TimeMonitor> CLevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : coarse" + levelSuffix, Timings0));
968 
969  bool emptySolve = true;
970 
971  // NOTE: we need to check using IsAvailable before Get here to avoid building default smoother
972  if (Fine->IsAvailable("PreSmoother")) {
973  RCP<SmootherBase> preSmoo = Fine->Get<RCP<SmootherBase>>("PreSmoother");
974  CompCoarse->start();
975  preSmoo->Apply(X, B, zeroGuess);
976  CompCoarse->stop();
977  zeroGuess = false;
978  emptySolve = false;
979  }
980  if (Fine->IsAvailable("PostSmoother")) {
981  RCP<SmootherBase> postSmoo = Fine->Get<RCP<SmootherBase>>("PostSmoother");
982  CompCoarse->start();
983  postSmoo->Apply(X, B, zeroGuess);
984  CompCoarse->stop();
985  emptySolve = false;
986  zeroGuess = false;
987  }
988  if (emptySolve == true) {
989  // Coarse operator is identity
990  X.update(one, B, zero);
991  }
992 
993  } else {
994  // On intermediate levels, we do cycles
995  RCP<Level> Coarse = Levels_[startLevel + 1];
996  {
997  // ============== PRESMOOTHING ==============
998  RCP<TimeMonitor> STime;
999  if (!useStackedTimer)
1000  STime = rcp(new TimeMonitor(*this, prefix + "Solve : smoothing (total)", Timings0));
1001  RCP<TimeMonitor> SLevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : smoothing" + levelSuffix, Timings0));
1002 
1003  if (Fine->IsAvailable("PreSmoother")) {
1004  RCP<SmootherBase> preSmoo = Fine->Get<RCP<SmootherBase>>("PreSmoother");
1005  preSmoo->Apply(X, B, zeroGuess);
1006  zeroGuess = false;
1007  }
1008  }
1009 
1010  RCP<MultiVector> residual;
1011  {
1012  RCP<TimeMonitor> ATime;
1013  if (!useStackedTimer)
1014  ATime = rcp(new TimeMonitor(*this, prefix + "Solve : residual calculation (total)", Timings0));
1015  RCP<TimeMonitor> ALevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : residual calculation" + levelSuffix, Timings0));
1016  if (zeroGuess) {
1017  // If there's a pre-smoother, then zeroGuess is false. If there isn't and people still have zeroGuess set,
1018  // then then X still has who knows what, so we need to zero it before we go to the coarse grid.
1019  X.putScalar(zero);
1020  }
1021 
1022  Utilities::Residual(*A, X, B, *residual_[startLevel]);
1023  residual = residual_[startLevel];
1024  }
1025 
1026  RCP<Operator> P = Coarse->Get<RCP<Operator>>("P");
1027  if (Coarse->IsAvailable("Pbar"))
1028  P = Coarse->Get<RCP<Operator>>("Pbar");
1029 
1030  RCP<MultiVector> coarseRhs, coarseX;
1031  // const bool initializeWithZeros = true;
1032  {
1033  // ============== RESTRICTION ==============
1034  RCP<TimeMonitor> RTime;
1035  if (!useStackedTimer)
1036  RTime = rcp(new TimeMonitor(*this, prefix + "Solve : restriction (total)", Timings0));
1037  RCP<TimeMonitor> RLevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : restriction" + levelSuffix, Timings0));
1038  coarseRhs = coarseRhs_[startLevel];
1039 
1040  if (implicitTranspose_) {
1041  P->apply(*residual, *coarseRhs, Teuchos::TRANS, one, zero);
1042 
1043  } else {
1044  RCP<Operator> R = Coarse->Get<RCP<Operator>>("R");
1045  R->apply(*residual, *coarseRhs, Teuchos::NO_TRANS, one, zero);
1046  }
1047  }
1048 
1049  RCP<const Import> importer;
1050  if (Coarse->IsAvailable("Importer"))
1051  importer = Coarse->Get<RCP<const Import>>("Importer");
1052 
1053  coarseX = coarseX_[startLevel];
1054  if (!doPRrebalance_ && !importer.is_null()) {
1055  RCP<TimeMonitor> ITime;
1056  if (!useStackedTimer)
1057  ITime = rcp(new TimeMonitor(*this, prefix + "Solve : import (total)", Timings0));
1058  RCP<TimeMonitor> ILevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : import" + levelSuffix1, Timings0));
1059 
1060  // Import: range map of R --> domain map of rebalanced Ac (before subcomm replacement)
1061  RCP<MultiVector> coarseTmp = coarseImport_[startLevel];
1062  coarseTmp->doImport(*coarseRhs, *importer, Xpetra::INSERT);
1063  coarseRhs.swap(coarseTmp);
1064  }
1065 
1066  RCP<Operator> Ac = Coarse->Get<RCP<Operator>>("A");
1067  if (!Ac.is_null()) {
1068  RCP<const Map> origXMap = coarseX->getMap();
1069  RCP<const Map> origRhsMap = coarseRhs->getMap();
1070 
1071  // Replace maps with maps with a subcommunicator
1072  coarseRhs->replaceMap(Ac->getRangeMap());
1073  coarseX->replaceMap(Ac->getDomainMap());
1074 
1075  {
1076  iterateLevelTime = Teuchos::null; // stop timing this level
1077 
1078  Iterate(*coarseRhs, *coarseX, 1, true, startLevel + 1);
1079  // ^^ zero initial guess
1080  if (Cycle_ == WCYCLE && WCycleStartLevel_ >= startLevel)
1081  Iterate(*coarseRhs, *coarseX, 1, false, startLevel + 1);
1082  // ^^ nonzero initial guess
1083 
1084  iterateLevelTime = rcp(new TimeMonitor(*this, iterateLevelTimeLabel)); // restart timing this level
1085  }
1086  coarseX->replaceMap(origXMap);
1087  coarseRhs->replaceMap(origRhsMap);
1088  }
1089 
1090  if (!doPRrebalance_ && !importer.is_null()) {
1091  RCP<TimeMonitor> ITime;
1092  if (!useStackedTimer)
1093  ITime = rcp(new TimeMonitor(*this, prefix + "Solve : export (total)", Timings0));
1094  RCP<TimeMonitor> ILevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : export" + levelSuffix1, Timings0));
1095 
1096  // Import: range map of rebalanced Ac (before subcomm replacement) --> domain map of P
1097  RCP<MultiVector> coarseTmp = coarseExport_[startLevel];
1098  coarseTmp->doExport(*coarseX, *importer, Xpetra::INSERT);
1099  coarseX.swap(coarseTmp);
1100  }
1101 
1102  {
1103  // ============== PROLONGATION ==============
1104  RCP<TimeMonitor> PTime;
1105  if (!useStackedTimer)
1106  PTime = rcp(new TimeMonitor(*this, prefix + "Solve : prolongation (total)", Timings0));
1107  RCP<TimeMonitor> PLevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : prolongation" + levelSuffix, Timings0));
1108  // Update X += P * coarseX
1109  // Note that due to what may be round-off error accumulation, use of the fused kernel
1110  // P->apply(*coarseX, X, Teuchos::NO_TRANS, one, one);
1111  // can in some cases result in slightly higher iteration counts.
1112  if (fuseProlongationAndUpdate_) {
1113  P->apply(*coarseX, X, Teuchos::NO_TRANS, scalingFactor_, one);
1114  } else {
1115  RCP<MultiVector> correction = correction_[startLevel];
1116  P->apply(*coarseX, *correction, Teuchos::NO_TRANS, one, zero);
1117  X.update(scalingFactor_, *correction, one);
1118  }
1119  }
1120 
1121  {
1122  // ============== POSTSMOOTHING ==============
1123  RCP<TimeMonitor> STime;
1124  if (!useStackedTimer)
1125  STime = rcp(new TimeMonitor(*this, prefix + "Solve : smoothing (total)", Timings0));
1126  RCP<TimeMonitor> SLevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : smoothing" + levelSuffix, Timings0));
1127 
1128  if (Fine->IsAvailable("PostSmoother")) {
1129  RCP<SmootherBase> postSmoo = Fine->Get<RCP<SmootherBase>>("PostSmoother");
1130  postSmoo->Apply(X, B, false);
1131  }
1132  }
1133  }
1134  zeroGuess = false;
1135 
1136  if (IsCalculationOfResidualRequired(startLevel, conv)) {
1137  ConvergenceStatus convergenceStatus = ComputeResidualAndPrintHistory(*A, X, B, iteration, startLevel, conv, prevNorm);
1138  if (convergenceStatus == MueLu::ConvergenceStatus::Converged)
1139  return convergenceStatus;
1140  }
1141  }
1142  return (tol > 0 ? ConvergenceStatus::Unconverged : ConvergenceStatus::Undefined);
1143 }
1144 #endif
1145 
1146 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1147 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(const LO& start, const LO& end, const std::string& suffix) {
1148  LO startLevel = (start != -1 ? start : 0);
1149  LO endLevel = (end != -1 ? end : Levels_.size() - 1);
1150 
1151  TEUCHOS_TEST_FOR_EXCEPTION(startLevel > endLevel, Exceptions::RuntimeError,
1152  "MueLu::Hierarchy::Write : startLevel must be <= endLevel");
1153 
1154  TEUCHOS_TEST_FOR_EXCEPTION(startLevel < 0 || endLevel >= Levels_.size(), Exceptions::RuntimeError,
1155  "MueLu::Hierarchy::Write bad start or end level");
1156 
1157  for (LO i = startLevel; i < endLevel + 1; i++) {
1158  RCP<Matrix> A = rcp_dynamic_cast<Matrix>(Levels_[i]->template Get<RCP<Operator>>("A")), P, R;
1159  if (i > 0) {
1160  P = rcp_dynamic_cast<Matrix>(Levels_[i]->template Get<RCP<Operator>>("P"));
1161  if (!implicitTranspose_)
1162  R = rcp_dynamic_cast<Matrix>(Levels_[i]->template Get<RCP<Operator>>("R"));
1163  }
1164 
1165  if (!A.is_null()) Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("A_" + toString(i) + suffix + ".m", *A);
1166  if (!P.is_null()) {
1167  Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("P_" + toString(i) + suffix + ".m", *P);
1168  }
1169  if (!R.is_null()) {
1170  Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("R_" + toString(i) + suffix + ".m", *R);
1171  }
1172  }
1173 }
1174 
1175 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1176 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Keep(const std::string& ename, const FactoryBase* factory) {
1177  for (Array<RCP<Level>>::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1178  (*it)->Keep(ename, factory);
1179 }
1180 
1181 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1182 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Delete(const std::string& ename, const FactoryBase* factory) {
1183  for (Array<RCP<Level>>::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1184  (*it)->Delete(ename, factory);
1185 }
1186 
1187 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1188 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::AddKeepFlag(const std::string& ename, const FactoryBase* factory, KeepType keep) {
1189  for (Array<RCP<Level>>::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1190  (*it)->AddKeepFlag(ename, factory, keep);
1191 }
1192 
1193 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1194 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::RemoveKeepFlag(const std::string& ename, const FactoryBase* factory, KeepType keep) {
1195  for (Array<RCP<Level>>::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1196  (*it)->RemoveKeepFlag(ename, factory, keep);
1197 }
1198 
1199 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1200 std::string Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::description() const {
1201  if (description_.empty()) {
1202  std::ostringstream out;
1203  out << BaseClass::description();
1204  out << "{#levels = " << GetGlobalNumLevels() << ", complexity = " << GetOperatorComplexity() << "}";
1205  description_ = out.str();
1206  }
1207  return description_;
1208 }
1209 
1210 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1211 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::describe(Teuchos::FancyOStream& out, const Teuchos::EVerbosityLevel tVerbLevel) const {
1212  describe(out, toMueLuVerbLevel(tVerbLevel));
1213 }
1214 
1215 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1216 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::describe(Teuchos::FancyOStream& out, const VerbLevel verbLevel) const {
1217  RCP<Operator> A0 = Levels_[0]->template Get<RCP<Operator>>("A");
1218  RCP<const Teuchos::Comm<int>> comm = A0->getDomainMap()->getComm();
1219 
1220  int numLevels = GetNumLevels();
1221  RCP<Operator> Ac = Levels_[numLevels - 1]->template Get<RCP<Operator>>("A");
1222  if (Ac.is_null()) {
1223  // It may happen that we do repartition on the last level, but the matrix
1224  // is small enough to satisfy "max coarse size" requirement. Then, even
1225  // though we have the level, the matrix would be null on all but one processors
1226  numLevels--;
1227  }
1228  int root = comm->getRank();
1229 
1230 #ifdef HAVE_MPI
1231  int smartData = numLevels * comm->getSize() + comm->getRank(), maxSmartData;
1232  reduceAll(*comm, Teuchos::REDUCE_MAX, smartData, Teuchos::ptr(&maxSmartData));
1233  root = maxSmartData % comm->getSize();
1234 #endif
1235 
1236  // Compute smoother complexity, if needed
1237  double smoother_comp = -1.0;
1238  if (verbLevel & (Statistics0 | Test))
1239  smoother_comp = GetSmootherComplexity();
1240 
1241  std::string outstr;
1242  if (comm->getRank() == root && verbLevel & (Statistics0 | Test)) {
1243  std::vector<Xpetra::global_size_t> nnzPerLevel;
1244  std::vector<Xpetra::global_size_t> rowsPerLevel;
1245  std::vector<int> numProcsPerLevel;
1246  bool someOpsNotMatrices = false;
1247  const Xpetra::global_size_t INVALID = Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid();
1248  for (int i = 0; i < numLevels; i++) {
1249  TEUCHOS_TEST_FOR_EXCEPTION(!(Levels_[i]->IsAvailable("A")), Exceptions::RuntimeError,
1250  "Operator A is not available on level " << i);
1251 
1252  RCP<Operator> A = Levels_[i]->template Get<RCP<Operator>>("A");
1253  TEUCHOS_TEST_FOR_EXCEPTION(A.is_null(), Exceptions::RuntimeError,
1254  "Operator A on level " << i << " is null.");
1255 
1256  RCP<Matrix> Am = rcp_dynamic_cast<Matrix>(A);
1257  if (Am.is_null()) {
1258  someOpsNotMatrices = true;
1259  nnzPerLevel.push_back(INVALID);
1260  rowsPerLevel.push_back(A->getDomainMap()->getGlobalNumElements());
1261  numProcsPerLevel.push_back(A->getDomainMap()->getComm()->getSize());
1262  } else {
1263  LO storageblocksize = Am->GetStorageBlockSize();
1264  Xpetra::global_size_t nnz = Am->getGlobalNumEntries() * storageblocksize * storageblocksize;
1265  nnzPerLevel.push_back(nnz);
1266  rowsPerLevel.push_back(Am->getGlobalNumRows() * storageblocksize);
1267  numProcsPerLevel.push_back(Am->getRowMap()->getComm()->getSize());
1268  }
1269  }
1270  if (someOpsNotMatrices)
1271  GetOStream(Warnings0) << "Some level operators are not matrices, statistics calculation are incomplete" << std::endl;
1272 
1273  {
1274  std::string label = Levels_[0]->getObjectLabel();
1275  std::ostringstream oss;
1276  oss << std::setfill(' ');
1277  oss << "\n--------------------------------------------------------------------------------\n";
1278  oss << "--- Multigrid Summary " << std::setw(32) << "---\n";
1279  oss << "--------------------------------------------------------------------------------" << std::endl;
1280  if (hierarchyLabel_ != "") oss << "Label = " << hierarchyLabel_ << std::endl;
1281  if (verbLevel & Parameters1)
1282  oss << "Scalar = " << Teuchos::ScalarTraits<Scalar>::name() << std::endl;
1283  oss << "Number of levels = " << numLevels << std::endl;
1284  oss << "Operator complexity = " << std::setprecision(2) << std::setiosflags(std::ios::fixed);
1285  if (!someOpsNotMatrices)
1286  oss << GetOperatorComplexity() << std::endl;
1287  else
1288  oss << "not available (Some operators in hierarchy are not matrices.)" << std::endl;
1289 
1290  if (smoother_comp != -1.0) {
1291  oss << "Smoother complexity = " << std::setprecision(2) << std::setiosflags(std::ios::fixed)
1292  << smoother_comp << std::endl;
1293  }
1294 
1295  switch (Cycle_) {
1296  case VCYCLE:
1297  oss << "Cycle type = V" << std::endl;
1298  break;
1299  case WCYCLE:
1300  oss << "Cycle type = W" << std::endl;
1301  if (WCycleStartLevel_ > 0)
1302  oss << "Cycle start level = " << WCycleStartLevel_ << std::endl;
1303  break;
1304  default:
1305  break;
1306  };
1307  oss << std::endl;
1308 
1309  Xpetra::global_size_t tt = rowsPerLevel[0];
1310  int rowspacer = 2;
1311  while (tt != 0) {
1312  tt /= 10;
1313  rowspacer++;
1314  }
1315  for (size_t i = 0; i < nnzPerLevel.size(); ++i) {
1316  tt = nnzPerLevel[i];
1317  if (tt != INVALID)
1318  break;
1319  tt = 100; // This will get used if all levels are operators.
1320  }
1321  int nnzspacer = 2;
1322  while (tt != 0) {
1323  tt /= 10;
1324  nnzspacer++;
1325  }
1326  tt = numProcsPerLevel[0];
1327  int npspacer = 2;
1328  while (tt != 0) {
1329  tt /= 10;
1330  npspacer++;
1331  }
1332  oss << "level " << std::setw(rowspacer) << " rows " << std::setw(nnzspacer) << " nnz "
1333  << " nnz/row" << std::setw(npspacer) << " c ratio"
1334  << " procs" << std::endl;
1335  for (size_t i = 0; i < nnzPerLevel.size(); ++i) {
1336  oss << " " << i << " ";
1337  oss << std::setw(rowspacer) << rowsPerLevel[i];
1338  if (nnzPerLevel[i] != INVALID) {
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  } else {
1343  oss << std::setw(nnzspacer) << "Operator";
1344  oss << std::setprecision(2) << std::setiosflags(std::ios::fixed);
1345  oss << std::setw(9) << " ";
1346  }
1347  if (i)
1348  oss << std::setw(9) << as<double>(rowsPerLevel[i - 1]) / rowsPerLevel[i];
1349  else
1350  oss << std::setw(9) << " ";
1351  oss << " " << std::setw(npspacer) << numProcsPerLevel[i] << std::endl;
1352  }
1353  oss << std::endl;
1354  for (int i = 0; i < GetNumLevels(); ++i) {
1355  RCP<SmootherBase> preSmoo, postSmoo;
1356  if (Levels_[i]->IsAvailable("PreSmoother"))
1357  preSmoo = Levels_[i]->template Get<RCP<SmootherBase>>("PreSmoother");
1358  if (Levels_[i]->IsAvailable("PostSmoother"))
1359  postSmoo = Levels_[i]->template Get<RCP<SmootherBase>>("PostSmoother");
1360 
1361  if (preSmoo != null && preSmoo == postSmoo)
1362  oss << "Smoother (level " << i << ") both : " << preSmoo->description() << std::endl;
1363  else {
1364  oss << "Smoother (level " << i << ") pre : "
1365  << (preSmoo != null ? preSmoo->description() : "no smoother") << std::endl;
1366  oss << "Smoother (level " << i << ") post : "
1367  << (postSmoo != null ? postSmoo->description() : "no smoother") << std::endl;
1368  }
1369 
1370  oss << std::endl;
1371  }
1372 
1373  outstr = oss.str();
1374  }
1375  }
1376 
1377 #ifdef HAVE_MPI
1378  RCP<const Teuchos::MpiComm<int>> mpiComm = rcp_dynamic_cast<const Teuchos::MpiComm<int>>(comm);
1379  MPI_Comm rawComm = (*mpiComm->getRawMpiComm())();
1380 
1381  int strLength = outstr.size();
1382  MPI_Bcast(&strLength, 1, MPI_INT, root, rawComm);
1383  if (comm->getRank() != root)
1384  outstr.resize(strLength);
1385  MPI_Bcast(&outstr[0], strLength, MPI_CHAR, root, rawComm);
1386 #endif
1387 
1388  out << outstr;
1389 }
1390 
1391 // NOTE: at some point this should be replaced by a friend operator <<
1392 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1393 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::print(std::ostream& out, const VerbLevel verbLevel) const {
1394  Teuchos::OSTab tab2(out);
1395  for (int i = 0; i < GetNumLevels(); ++i)
1396  Levels_[i]->print(out, verbLevel);
1397 }
1398 
1399 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1400 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::IsPreconditioner(const bool flag) {
1401  isPreconditioner_ = flag;
1402 }
1403 
1404 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1405 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DumpCurrentGraph(int currLevel) const {
1406  if (GetProcRankVerbose() != 0)
1407  return;
1408 #if defined(HAVE_MUELU_BOOST) && defined(HAVE_MUELU_BOOST_FOR_REAL) && defined(BOOST_VERSION) && (BOOST_VERSION >= 104400)
1409 
1410  BoostGraph graph;
1411 
1412  BoostProperties dp;
1413  dp.property("label", boost::get(boost::vertex_name, graph));
1414  dp.property("id", boost::get(boost::vertex_index, graph));
1415  dp.property("label", boost::get(boost::edge_name, graph));
1416  dp.property("color", boost::get(boost::edge_color, graph));
1417 
1418  // create local maps
1419  std::map<const FactoryBase*, BoostVertex> vindices;
1420  typedef std::map<std::pair<BoostVertex, BoostVertex>, std::string> emap;
1421  emap edges;
1422 
1423  static int call_id = 0;
1424 
1425  RCP<Operator> A = Levels_[0]->template Get<RCP<Operator>>("A");
1426  int rank = A->getDomainMap()->getComm()->getRank();
1427 
1428  // printf("[%d] CMS: ----------------------\n",rank);
1429  for (int i = currLevel; i <= currLevel + 1 && i < GetNumLevels(); i++) {
1430  edges.clear();
1431  Levels_[i]->UpdateGraph(vindices, edges, dp, graph);
1432 
1433  for (emap::const_iterator eit = edges.begin(); eit != edges.end(); eit++) {
1434  std::pair<BoostEdge, bool> boost_edge = boost::add_edge(eit->first.first, eit->first.second, graph);
1435  // printf("[%d] CMS: Hierarchy, adding edge (%d->%d) %d\n",rank,(int)eit->first.first,(int)eit->first.second,(int)boost_edge.second);
1436  // Because xdot.py views 'Graph' as a keyword
1437  if (eit->second == std::string("Graph"))
1438  boost::put("label", dp, boost_edge.first, std::string("Graph_"));
1439  else
1440  boost::put("label", dp, boost_edge.first, eit->second);
1441  if (i == currLevel)
1442  boost::put("color", dp, boost_edge.first, std::string("red"));
1443  else
1444  boost::put("color", dp, boost_edge.first, std::string("blue"));
1445  }
1446  }
1447 
1448  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"));
1449  boost::write_graphviz_dp(out, graph, dp, std::string("id"));
1450  out.close();
1451  call_id++;
1452 #else
1453  GetOStream(Errors) << "Dependency graph output requires boost and MueLu_ENABLE_Boost_for_real" << std::endl;
1454 #endif
1455 }
1456 
1457 // Enforce that coordinate vector's map is consistent with that of A
1458 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1459 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ReplaceCoordinateMap(Level& level) {
1460  RCP<Operator> Ao = level.Get<RCP<Operator>>("A");
1461  RCP<Matrix> A = rcp_dynamic_cast<Matrix>(Ao);
1462  if (A.is_null()) {
1463  GetOStream(Runtime1) << "Hierarchy::ReplaceCoordinateMap: operator is not a matrix, skipping..." << std::endl;
1464  return;
1465  }
1466  if (Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A) != Teuchos::null) {
1467  GetOStream(Runtime1) << "Hierarchy::ReplaceCoordinateMap: operator is a BlockedCrsMatrix, skipping..." << std::endl;
1468  return;
1469  }
1470 
1471  typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LO, GO, NO> xdMV;
1472 
1473  RCP<xdMV> coords = level.Get<RCP<xdMV>>("Coordinates");
1474 
1475  if (A->getRowMap()->isSameAs(*(coords->getMap()))) {
1476  GetOStream(Runtime1) << "Hierarchy::ReplaceCoordinateMap: matrix and coordinates maps are same, skipping..." << std::endl;
1477  return;
1478  }
1479 
1480  if (A->IsView("stridedMaps") && rcp_dynamic_cast<const StridedMap>(A->getRowMap("stridedMaps")) != Teuchos::null) {
1481  RCP<const StridedMap> stridedRowMap = rcp_dynamic_cast<const StridedMap>(A->getRowMap("stridedMaps"));
1482 
1483  // It is better to through an exceptions if maps may be inconsistent, than to ignore it and experience unfathomable breakdowns
1484  TEUCHOS_TEST_FOR_EXCEPTION(stridedRowMap->getStridedBlockId() != -1 || stridedRowMap->getOffset() != 0,
1485  Exceptions::RuntimeError, "Hierarchy::ReplaceCoordinateMap: nontrivial maps (block id = " << stridedRowMap->getStridedBlockId() << ", offset = " << stridedRowMap->getOffset() << ")");
1486  }
1487 
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");
1490 
1491  size_t blkSize = A->GetFixedBlockSize() / A->GetStorageBlockSize();
1492 
1493  RCP<const Map> nodeMap = A->getRowMap();
1494  if (blkSize > 1) {
1495  // Create a nodal map, as coordinates have not been expanded to a DOF map yet.
1496  RCP<const Map> dofMap = A->getRowMap();
1497  GO indexBase = dofMap->getIndexBase();
1498  size_t numLocalDOFs = dofMap->getLocalNumElements();
1499  TEUCHOS_TEST_FOR_EXCEPTION(numLocalDOFs % blkSize, Exceptions::RuntimeError,
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();
1502 
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;
1506 
1507  Xpetra::global_size_t INVALID = Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid();
1508  nodeMap = MapFactory::Build(dofMap->lib(), INVALID, nodeGIDs(), indexBase, dofMap->getComm());
1509  } else {
1510  // blkSize == 1
1511  // Check whether the length of vectors fits to the size of A
1512  // If yes, make sure that the maps are matching
1513  // If no, throw a warning but do not touch the Coordinates
1514  if (coords->getLocalLength() != A->getRowMap()->getLocalNumElements()) {
1515  GetOStream(Warnings) << "Coordinate vector does not match row map of matrix A!" << std::endl;
1516  return;
1517  }
1518  }
1519 
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]());
1525  }
1526 
1527  RCP<xdMV> newCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LO, GO, NO>::Build(nodeMap, coordDataView(), coords->getNumVectors());
1528  level.Set("Coordinates", newCoords);
1529 }
1530 
1531 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1532 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::AllocateLevelMultiVectors(int numvecs, bool forceMapCheck) {
1533  int N = Levels_.size();
1534  if (((sizeOfAllocatedLevelMultiVectors_ == numvecs && residual_.size() == N) || numvecs <= 0) && !forceMapCheck) return;
1535 
1536  // If, somehow, we changed the number of levels, delete everything first
1537  if (residual_.size() != N) {
1538  DeleteLevelMultiVectors();
1539 
1540  residual_.resize(N);
1541  coarseRhs_.resize(N);
1542  coarseX_.resize(N);
1543  coarseImport_.resize(N);
1544  coarseExport_.resize(N);
1545  correction_.resize(N);
1546  }
1547 
1548  for (int i = 0; i < N; i++) {
1549  RCP<Operator> A = Levels_[i]->template Get<RCP<Operator>>("A");
1550  if (!A.is_null()) {
1551  // This dance is because we allow A to have a BlockedMap and X/B to have (compatible) non-blocked map
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();
1557  }
1558 
1559  if (residual_[i].is_null() || !residual_[i]->getMap()->isSameAs(*Arm))
1560  // This is zero'd by default since it is filled via an operator apply
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);
1564  }
1565 
1566  if (i + 1 < N) {
1567  // This is zero'd by default since it is filled via an operator apply
1568  if (implicitTranspose_) {
1569  RCP<Operator> P = Levels_[i + 1]->template Get<RCP<Operator>>("P");
1570  if (!P.is_null()) {
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);
1574  }
1575  } else {
1576  RCP<Operator> R = Levels_[i + 1]->template Get<RCP<Operator>>("R");
1577  if (!R.is_null()) {
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);
1581  }
1582  }
1583 
1584  RCP<const Import> importer;
1585  if (Levels_[i + 1]->IsAvailable("Importer"))
1586  importer = Levels_[i + 1]->template Get<RCP<const Import>>("Importer");
1587  if (doPRrebalance_ || importer.is_null()) {
1588  RCP<const Map> map = coarseRhs_[i]->getMap();
1589  if (coarseX_[i].is_null() || !coarseX_[i]->getMap()->isSameAs(*map))
1590  coarseX_[i] = MultiVectorFactory::Build(map, numvecs, true);
1591  } else {
1592  RCP<const Map> map;
1593  map = importer->getTargetMap();
1594  if (coarseImport_[i].is_null() || !coarseImport_[i]->getMap()->isSameAs(*map)) {
1595  coarseImport_[i] = MultiVectorFactory::Build(map, numvecs, false);
1596  coarseX_[i] = MultiVectorFactory::Build(map, numvecs, false);
1597  }
1598  map = importer->getSourceMap();
1599  if (coarseExport_[i].is_null() || !coarseExport_[i]->getMap()->isSameAs(*map))
1600  coarseExport_[i] = MultiVectorFactory::Build(map, numvecs, false);
1601  }
1602  }
1603  }
1604  sizeOfAllocatedLevelMultiVectors_ = numvecs;
1605 }
1606 
1607 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1608 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeleteLevelMultiVectors() {
1609  if (sizeOfAllocatedLevelMultiVectors_ == 0) return;
1610  residual_.resize(0);
1611  coarseRhs_.resize(0);
1612  coarseX_.resize(0);
1613  coarseImport_.resize(0);
1614  coarseExport_.resize(0);
1615  correction_.resize(0);
1616  sizeOfAllocatedLevelMultiVectors_ = 0;
1617 }
1618 
1619 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1620 bool Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::IsCalculationOfResidualRequired(
1621  const LO startLevel, const ConvData& conv) const {
1622  return (startLevel == 0 && !isPreconditioner_ && (IsPrint(Statistics1) || conv.tol_ > 0));
1623 }
1624 
1625 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1626 ConvergenceStatus Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::IsConverged(
1627  const Teuchos::Array<MagnitudeType>& residualNorm, const MagnitudeType convergenceTolerance) const {
1628  ConvergenceStatus convergenceStatus = ConvergenceStatus::Undefined;
1629 
1630  if (convergenceTolerance > Teuchos::ScalarTraits<MagnitudeType>::zero()) {
1631  bool passed = true;
1632  for (LO k = 0; k < residualNorm.size(); k++)
1633  if (residualNorm[k] >= convergenceTolerance)
1634  passed = false;
1635 
1636  if (passed)
1637  convergenceStatus = ConvergenceStatus::Converged;
1638  else
1639  convergenceStatus = ConvergenceStatus::Unconverged;
1640  }
1641 
1642  return convergenceStatus;
1643 }
1644 
1645 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1646 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::PrintResidualHistory(
1647  const LO iteration, const Teuchos::Array<MagnitudeType>& residualNorm) const {
1648  GetOStream(Statistics1) << "iter: "
1649  << std::setiosflags(std::ios::left)
1650  << std::setprecision(3) << std::setw(4) << iteration
1651  << " residual = "
1652  << std::setprecision(10) << residualNorm
1653  << std::endl;
1654 }
1655 
1656 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1657 ConvergenceStatus Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ComputeResidualAndPrintHistory(
1658  const Operator& A, const MultiVector& X, const MultiVector& B, const LO iteration,
1659  const LO startLevel, const ConvData& conv, MagnitudeType& previousResidualNorm) {
1660  Teuchos::Array<MagnitudeType> residualNorm;
1661  residualNorm = Utilities::ResidualNorm(A, X, B, *residual_[startLevel]);
1662 
1663  const MagnitudeType currentResidualNorm = residualNorm[0];
1664  rate_ = currentResidualNorm / previousResidualNorm;
1665  previousResidualNorm = currentResidualNorm;
1666 
1667  if (IsPrint(Statistics1))
1668  PrintResidualHistory(iteration, residualNorm);
1669 
1670  return IsConverged(residualNorm, conv.tol_);
1671 }
1672 
1673 } // namespace MueLu
1674 
1675 #endif // MUELU_HIERARCHY_DEF_HPP
MueLu::Warnings0
Important warning messages (one line)
Definition: MueLu_VerbosityLevel.hpp:20
MueLu::TopSmootherFactory::Build
void Build(Level &level) const
Build an object with this factory.
Definition: MueLu_TopSmootherFactory_def.hpp:78
MueLu::Hierarchy::IsPreconditioner
void IsPreconditioner(const bool flag)
Definition: MueLu_Hierarchy_def.hpp:1400
MueLu::Timings0
High level timing information (use Teuchos::TimeMonitor::summarize() to print)
Definition: MueLu_VerbosityLevel.hpp:37
Teuchos::Describable::getObjectLabel
virtual std::string getObjectLabel() const
MueLu_SmootherBase.hpp
MueLu::Level::Get
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-&gt;Get&lt; RCP&lt;Matrix&gt; &gt;(&quot;A&quot;, factory) if factory == NULL =&gt; use default factory.
Definition: MueLu_Level.hpp:154
is_null
bool is_null(const boost::shared_ptr< T > &p)
MueLu::Hierarchy::GetLevel
RCP< Level > & GetLevel(const int levelID=0)
Retrieve a certain level from hierarchy.
Definition: MueLu_Hierarchy_def.hpp:124
MueLu::toString
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Definition: MueLu_Utilities_decl.hpp:599
MueLu::Hierarchy::Hierarchy
Hierarchy()
Default constructor.
Definition: MueLu_Hierarchy_def.hpp:41
MueLu::SmootherBase::getNodeSmootherComplexity
virtual size_t getNodeSmootherComplexity() const =0
Compute a rough estimate of the cost to apply this smoother on this MPI rank. Return Teuchos::Ordinal...
MueLu::Level::Release
void Release(const FactoryBase &factory)
Decrement the storage counter for all the inputs of a factory.
Definition: MueLu_Level.cpp:149
GO
GlobalOrdinal GO
Teuchos::TimeMonitor::getNewCounter
static RCP< Time > getNewCounter(const std::string &name)
MueLu::Hierarchy::AddLevel
void AddLevel(const RCP< Level > &level)
Add a level at the end of the hierarchy.
Definition: MueLu_Hierarchy_def.hpp:101
MueLu::Hierarchy::CheckForEmptySmoothersAndCoarseSolve
void CheckForEmptySmoothersAndCoarseSolve()
Definition: MueLu_Hierarchy_def.hpp:654
MueLu::Hierarchy::PrintResidualHistory
void PrintResidualHistory(const LO iteration, const Teuchos::Array< MagnitudeType > &residualNorm) const
Print residualNorm for this iteration to the screen.
Definition: MueLu_Hierarchy_def.hpp:1646
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
MueLu::KeepType
short KeepType
Definition: MueLu_KeepType.hpp:27
Xpetra::IO::Write
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
MueLu::Hierarchy::DumpCurrentGraph
void DumpCurrentGraph(int level) const
Definition: MueLu_Hierarchy_def.hpp:1405
MueLu::Statistics1
Print more statistics.
Definition: MueLu_VerbosityLevel.hpp:34
MueLu::Hierarchy::AddNewLevel
void AddNewLevel()
Add a new level at the end of the hierarchy.
Definition: MueLu_Hierarchy_def.hpp:117
MueLu::toMueLuVerbLevel
VerbLevel toMueLuVerbLevel(const Teuchos::EVerbosityLevel verbLevel)
Translate Teuchos verbosity level to MueLu verbosity level.
Definition: MueLu_VerbosityLevel.cpp:18
Teuchos::RCP::swap
void swap(RCP< T > &r_ptr)
LO
LocalOrdinal LO
MueLu::Runtime0
One-liner description of what is happening.
Definition: MueLu_VerbosityLevel.hpp:25
MueLu_Level.hpp
MueLu::Hierarchy::Clear
void Clear(int startLevel=0)
Clear impermanent data from previous setup.
Definition: MueLu_Hierarchy_def.hpp:664
Teuchos::EVerbosityLevel
EVerbosityLevel
MueLu::Level::SetPreviousLevel
void SetPreviousLevel(const RCP< Level > &previousLevel)
Definition: MueLu_Level.cpp:62
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Teuchos::basic_OSTab
Xpetra_Matrix.hpp
MueLu::Exceptions::Incompatible
Exception throws to report incompatible objects (like maps).
Definition: MueLu_Exceptions.hpp:50
MueLu::UtilitiesBase::SetRandomSeed
static void SetRandomSeed(const Teuchos::Comm< int > &comm)
Set seed for random number generator.
Definition: MueLu_UtilitiesBase_def.hpp:1756
MueLu::TimeMonitor
Integrates Teuchos::TimeMonitor with MueLu verbosity system.
Definition: MueLu_TimeMonitor.hpp:31
Teuchos::ScalarTraits::name
static std::string name()
MueLu::NoFactory::get
static const NoFactory * get()
Definition: MueLu_NoFactory.cpp:34
MueLu::Hierarchy::Iterate
ConvergenceStatus Iterate(const MultiVector &B, MultiVector &X, ConvData conv=ConvData(), bool InitialGuessIsZero=false, LO startLevel=0)
Apply the multigrid preconditioner.
Definition: MueLu_Hierarchy_def.hpp:877
MueLu::Hierarchy::ReplaceCoordinateMap
void ReplaceCoordinateMap(Level &level)
Definition: MueLu_Hierarchy_def.hpp:1459
MueLu::FormattingHelper::getColonLabel
static std::string getColonLabel(const std::string &label)
Helper function for object label.
Definition: MueLu_Monitor.hpp:30
MueLu::Hierarchy::AddKeepFlag
void AddKeepFlag(const std::string &ename, const FactoryBase *factory=NoFactory::get(), KeepType keep=MueLu::Keep)
Call Level::AddKeepFlag for each level of the Hierarchy.
Definition: MueLu_Hierarchy_def.hpp:1188
MueLu::Statistics2
Print even more statistics.
Definition: MueLu_VerbosityLevel.hpp:35
MueLu::Level::setlib
void setlib(Xpetra::UnderlyingLib lib2)
Definition: MueLu_Level.hpp:326
MueLu::Warnings1
Additional warnings.
Definition: MueLu_VerbosityLevel.hpp:22
Xpetra_IO.hpp
MueLu::Hierarchy::print
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...
Definition: MueLu_Hierarchy_def.hpp:1393
MueLu::Hierarchy::GetSmootherComplexity
double GetSmootherComplexity() const
Definition: MueLu_Hierarchy_def.hpp:171
MueLu::FactoryBase
Base class for factories (e.g., R, P, and A_coarse).
Definition: MueLu_FactoryBase.hpp:24
MueLu::Statistics0
Print statistics that do not involve significant additional computation.
Definition: MueLu_VerbosityLevel.hpp:33
MueLu::Hierarchy::lib_
Xpetra::UnderlyingLib lib_
Epetra/Tpetra mode.
Definition: MueLu_Hierarchy_decl.hpp:397
MueLu::UtilitiesBase::Residual
static RCP< MultiVector > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS)
Definition: MueLu_UtilitiesBase_def.hpp:925
rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::Hierarchy::Write
void Write(const LO &start=-1, const LO &end=-1, const std::string &suffix="")
Print matrices in the multigrid hierarchy to file.
Definition: MueLu_Hierarchy_def.hpp:1147
MueLu::Level
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:63
MueLu_Hierarchy_decl.hpp
MueLu::Hierarchy::description
std::string description() const
Return a simple one-line description of this object.
Definition: MueLu_Hierarchy_def.hpp:1200
MueLu::FactoryManagerBase
Class that provides default factories within Needs class.
Definition: MueLu_FactoryManagerBase.hpp:27
Teuchos::basic_FancyOStream
MueLu::Hierarchy::RemoveKeepFlag
void RemoveKeepFlag(const std::string &ename, const FactoryBase *factory, KeepType keep=MueLu::All)
Call Level::RemoveKeepFlag for each level of the Hierarchy.
Definition: MueLu_Hierarchy_def.hpp:1194
MueLu::Hierarchy::CheckLevel
void CheckLevel(Level &level, int levelID)
Helper function.
Definition: MueLu_Hierarchy_def.hpp:213
MueLu::Level::lib
Xpetra::UnderlyingLib lib()
Definition: MueLu_Level.hpp:327
Teuchos::Describable::setObjectLabel
virtual void setObjectLabel(const std::string &objectLabel)
MueLu::Hierarchy::AllocateLevelMultiVectors
void AllocateLevelMultiVectors(int numvecs, bool forceMapCheck=false)
Definition: MueLu_Hierarchy_def.hpp:1532
MueLu::PerfUtils::PrintMatrixInfo
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
Definition: MueLu_PerfUtils_def.hpp:92
Xpetra::global_size_t
size_t global_size_t
MueLu::Level::Set
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
Definition: MueLu_Level.hpp:120
MueLu::Hierarchy::SetLabel
void SetLabel(const std::string &hierarchyLabel)
Definition: MueLu_Hierarchy_decl.hpp:258
MueLu::Hierarchy::SetMatvecParams
void SetMatvecParams(RCP< ParameterList > matvecParams)
Definition: MueLu_Hierarchy_def.hpp:223
MueLu::TopRAPFactory::Build
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Definition: MueLu_TopRAPFactory_def.hpp:51
MueLu::Hierarchy::GetGlobalNumLevels
int GetGlobalNumLevels() const
Definition: MueLu_Hierarchy_def.hpp:136
Teuchos::Array::push_back
void push_back(const value_type &x)
MueLu::Level::GetPreviousLevel
RCP< Level > & GetPreviousLevel()
Previous level.
Definition: MueLu_Level.cpp:60
MueLu::UtilitiesBase::ResidualNorm
static Teuchos::Array< Magnitude > ResidualNorm(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS)
Definition: MueLu_UtilitiesBase_def.hpp:901
MueLu::Hierarchy::ConvData
Data struct for defining stopping criteria of multigrid iteration.
Definition: MueLu_Hierarchy_decl.hpp:71
MueLu::Level::SetLevelID
void SetLevelID(int levelID)
Set level number.
Definition: MueLu_Level.cpp:53
MueLu::Hierarchy::MagnitudeType
STS::magnitudeType MagnitudeType
Definition: MueLu_Hierarchy_decl.hpp:68
MueLu::Hierarchy::Delete
void Delete(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Call Level::Delete(ename, factory) for each level of the Hierarchy.
Definition: MueLu_Hierarchy_def.hpp:1182
MueLu::Monitor
Timer to be used in non-factories.
Definition: MueLu_Monitor.hpp:68
Teuchos::Array::size
size_type size() const
MueLu::Hierarchy::Levels_
Array< RCP< Level > > Levels_
Container for Level objects.
Definition: MueLu_Hierarchy_decl.hpp:359
SC
Scalar SC
MueLu_Monitor.hpp
MueLu::Hierarchy::IsCalculationOfResidualRequired
bool IsCalculationOfResidualRequired(const LO startLevel, const ConvData &conv) const
Decide if the residual needs to be computed.
Definition: MueLu_Hierarchy_def.hpp:1620
NO
Node NO
Teuchos_TimeMonitor.hpp
MueLu::Level::GetLevelID
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:51
MueLu::Parameters1
Print class parameters (more parameters, more verbose)
Definition: MueLu_VerbosityLevel.hpp:31
MueLu::Hierarchy::~Hierarchy
virtual ~Hierarchy()
Destructor.
residual
void residual(const Operator< SC, LO, GO, NO > &Aop, const MultiVector< SC, LO, GO, NO > &X_in, const MultiVector< SC, LO, GO, NO > &B_in, MultiVector< SC, LO, GO, NO > &R_in)
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: MueLu_Exceptions.hpp:36
MueLu::SetFactoryManager
An exception safe way to call the method &#39;Level::SetFactoryManager()&#39;.
Definition: MueLu_HierarchyUtils_decl.hpp:29
MueLu::Level::SetComm
void SetComm(RCP< const Teuchos::Comm< int > > const &comm)
Definition: MueLu_Level.hpp:329
MueLu::Hierarchy::describe
void describe(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the Hierarchy with some verbosity level to a FancyOStream object.
Definition: MueLu_Hierarchy_def.hpp:1216
MueLu::Warnings
Print all warning messages.
Definition: MueLu_VerbosityLevel.hpp:52
MueLu::Runtime1
Description of what is happening (more verbose)
Definition: MueLu_VerbosityLevel.hpp:26
MueLu::Hierarchy::Keep
void Keep(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Call Level::Keep(ename, factory) for each level of the Hierarchy.
Definition: MueLu_Hierarchy_def.hpp:1176
MueLu::Hierarchy::IsConverged
ConvergenceStatus IsConverged(const Teuchos::Array< MagnitudeType > &residualNorm, const MagnitudeType convergenceTolerance) const
Decide if the multigrid iteration is converged.
Definition: MueLu_Hierarchy_def.hpp:1626
MueLu::Hierarchy::GetOperatorComplexity
double GetOperatorComplexity() const
Definition: MueLu_Hierarchy_def.hpp:148
MueLu::Hierarchy::Setup
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.
Definition: MueLu_Hierarchy_def.hpp:267
MueLu::Level::GetComm
RCP< const Teuchos::Comm< int > > GetComm() const
Definition: MueLu_Level.hpp:330
Xpetra_Operator.hpp
MueLu::Hierarchy::ComputeResidualAndPrintHistory
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.
Definition: MueLu_Hierarchy_def.hpp:1657
MueLu::Describable::description
virtual std::string description() const
Return a simple one-line description of this object.
Definition: MueLu_Describable.cpp:25
MueLu::Hierarchy
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
Definition: MueLu_Hierarchy_decl.hpp:63
Teuchos::toString
std::string toString(const T &t)
MueLu::Level::IsAvailable
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need&#39;s value has been saved.
Definition: MueLu_Level.cpp:281
MueLu::Level::Request
void Request(const FactoryBase &factory)
Increment the storage counter for all the inputs of a factory.
Definition: MueLu_Level.cpp:142
Teuchos::RCP::is_null
bool is_null() const