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