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