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