MueLu  Version of the Day
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()),
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()),
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  // NOTE: All reuse cases leave all of the maps the same, meaning that we do not
556  // need to reallocated the cached multivectors for Iterate(). If this were to change,
557  // we'd want to do a DeleteLevelMultiVectors() and AllocateLevelMultiVectors() here.
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  }
989  if (emptySolve == true) {
990  GetOStream(Warnings1) << "No coarse grid solver" << std::endl;
991  // Coarse operator is identity
992  X.update(one, B, zero);
993  }
994 
995  } else {
996  // On intermediate levels, we do cycles
997  RCP<Level> Coarse = Levels_[startLevel+1];
998  {
999  // ============== PRESMOOTHING ==============
1000  RCP<TimeMonitor> STime;
1001  if (!useStackedTimer)
1002  STime = rcp(new TimeMonitor(*this, prefix + "Solve : smoothing (total)" , Timings0));
1003  RCP<TimeMonitor> SLevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : smoothing" + levelSuffix, Timings0));
1004 
1005  if (Fine->IsAvailable("PreSmoother")) {
1006  RCP<SmootherBase> preSmoo = Fine->Get< RCP<SmootherBase> >("PreSmoother");
1007  preSmoo->Apply(X, B, zeroGuess);
1008  }
1009  }
1010 
1011  RCP<MultiVector> residual;
1012  {
1013  RCP<TimeMonitor> ATime;
1014  if (!useStackedTimer)
1015  ATime = rcp(new TimeMonitor(*this, prefix + "Solve : residual calculation (total)" , Timings0));
1016  RCP<TimeMonitor> ALevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : residual calculation" + levelSuffix, Timings0));
1017  Utilities::Residual(*A, X, B,*residual_[startLevel]);
1018  residual = residual_[startLevel];
1019  }
1020 
1021  RCP<Operator> P = Coarse->Get< RCP<Operator> >("P");
1022  if (Coarse->IsAvailable("Pbar"))
1023  P = Coarse->Get< RCP<Operator> >("Pbar");
1024 
1025  RCP<MultiVector> coarseRhs, coarseX;
1026  // const bool initializeWithZeros = true;
1027  {
1028  // ============== RESTRICTION ==============
1029  RCP<TimeMonitor> RTime;
1030  if (!useStackedTimer)
1031  RTime = rcp(new TimeMonitor(*this, prefix + "Solve : restriction (total)" , Timings0));
1032  RCP<TimeMonitor> RLevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : restriction" + levelSuffix, Timings0));
1033  coarseRhs = coarseRhs_[startLevel];
1034 
1035  if (implicitTranspose_) {
1036  P->apply(*residual, *coarseRhs, Teuchos::TRANS, one, zero);
1037 
1038  } else {
1039  RCP<Operator> R = Coarse->Get< RCP<Operator> >("R");
1040  R->apply(*residual, *coarseRhs, Teuchos::NO_TRANS, one, zero);
1041  }
1042  }
1043 
1044  RCP<const Import> importer;
1045  if (Coarse->IsAvailable("Importer"))
1046  importer = Coarse->Get< RCP<const Import> >("Importer");
1047 
1048  coarseX = coarseX_[startLevel];
1049  if (!doPRrebalance_ && !importer.is_null()) {
1050  RCP<TimeMonitor> ITime;
1051  if (!useStackedTimer)
1052  ITime = rcp(new TimeMonitor(*this, prefix + "Solve : import (total)" , Timings0));
1053  RCP<TimeMonitor> ILevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : import" + levelSuffix1, Timings0));
1054 
1055  // Import: range map of R --> domain map of rebalanced Ac (before subcomm replacement)
1056  RCP<MultiVector> coarseTmp = coarseImport_[startLevel];
1057  coarseTmp->doImport(*coarseRhs, *importer, Xpetra::INSERT);
1058  coarseRhs.swap(coarseTmp);
1059  }
1060 
1061  RCP<Operator> Ac = Coarse->Get< RCP<Operator> >("A");
1062  if (!Ac.is_null()) {
1063  RCP<const Map> origXMap = coarseX->getMap();
1064  RCP<const Map> origRhsMap = coarseRhs->getMap();
1065 
1066  // Replace maps with maps with a subcommunicator
1067  coarseRhs->replaceMap(Ac->getRangeMap());
1068  coarseX ->replaceMap(Ac->getDomainMap());
1069 
1070  {
1071  iterateLevelTime = Teuchos::null; // stop timing this level
1072 
1073  Iterate(*coarseRhs, *coarseX, 1, true, startLevel+1);
1074  // ^^ zero initial guess
1075  if (Cycle_ == WCYCLE)
1076  Iterate(*coarseRhs, *coarseX, 1, false, startLevel+1);
1077  // ^^ nonzero initial guess
1078 
1079  iterateLevelTime = rcp(new TimeMonitor(*this, iterateLevelTimeLabel)); // restart timing this level
1080  }
1081  coarseX->replaceMap(origXMap);
1082  coarseRhs->replaceMap(origRhsMap);
1083  }
1084 
1085  if (!doPRrebalance_ && !importer.is_null()) {
1086  RCP<TimeMonitor> ITime;
1087  if (!useStackedTimer)
1088  ITime = rcp(new TimeMonitor(*this, prefix + "Solve : export (total)" , Timings0));
1089  RCP<TimeMonitor> ILevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : export" + levelSuffix1, Timings0));
1090 
1091  // Import: range map of rebalanced Ac (before subcomm replacement) --> domain map of P
1092  RCP<MultiVector> coarseTmp = coarseExport_[startLevel];
1093  coarseTmp->doExport(*coarseX, *importer, Xpetra::INSERT);
1094  coarseX.swap(coarseTmp);
1095  }
1096 
1097  {
1098  // ============== PROLONGATION ==============
1099  RCP<TimeMonitor> PTime;
1100  if (!useStackedTimer)
1101  PTime = rcp(new TimeMonitor(*this, prefix + "Solve : prolongation (total)" , Timings0));
1102  RCP<TimeMonitor> PLevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : prolongation" + levelSuffix, Timings0));
1103  // Update X += P * coarseX
1104  // Note that due to what may be round-off error accumulation, use of the fused kernel
1105  // P->apply(*coarseX, X, Teuchos::NO_TRANS, one, one);
1106  // can in some cases result in slightly higher iteration counts.
1107  if (fuseProlongationAndUpdate_) {
1108  P->apply(*coarseX, X, Teuchos::NO_TRANS, scalingFactor_, one);
1109  } else {
1110  RCP<MultiVector> correction = correction_[startLevel];
1111  P->apply(*coarseX, *correction, Teuchos::NO_TRANS, one, zero);
1112  X.update(scalingFactor_, *correction, one);
1113  }
1114  }
1115 
1116  {
1117  // ============== POSTSMOOTHING ==============
1118  RCP<TimeMonitor> STime;
1119  if (!useStackedTimer)
1120  STime = rcp(new TimeMonitor(*this, prefix + "Solve : smoothing (total)" , Timings0));
1121  RCP<TimeMonitor> SLevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : smoothing" + levelSuffix, Timings0));
1122 
1123  if (Fine->IsAvailable("PostSmoother")) {
1124  RCP<SmootherBase> postSmoo = Fine->Get< RCP<SmootherBase> >("PostSmoother");
1125  postSmoo->Apply(X, B, false);
1126  }
1127  }
1128  }
1129  zeroGuess = false;
1130 
1131  if (startLevel == 0 && !isPreconditioner_ &&
1132  (IsPrint(Statistics1) || tol > 0)) {
1133  // We calculate the residual only if we want to print it out, or if we
1134  // want to stop once we achive the tolerance
1135  Teuchos::Array<MagnitudeType> rn;
1136  rn = Utilities::ResidualNorm(*A, X, B,*residual_[startLevel]);
1137 
1138  prevNorm = curNorm;
1139  curNorm = rn[0];
1140  rate_ = as<MagnitudeType>(curNorm / prevNorm);
1141 
1142  if (IsPrint(Statistics1))
1143  GetOStream(Statistics1) << "iter: "
1144  << std::setiosflags(std::ios::left)
1145  << std::setprecision(3) << i
1146  << " residual = "
1147  << std::setprecision(10) << rn
1148  << std::endl;
1149 
1150  if (tol > 0) {
1151  bool passed = true;
1152  for (LO k = 0; k < rn.size(); k++)
1153  if (rn[k] >= tol)
1154  passed = false;
1155 
1156  if (passed)
1157  return Converged;
1158  }
1159  }
1160  }
1161  return (tol > 0 ? Unconverged : Undefined);
1162  }
1163 #endif
1164 
1165 
1166  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1167  void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(const LO& start, const LO& end, const std::string &suffix) {
1168  LO startLevel = (start != -1 ? start : 0);
1169  LO endLevel = (end != -1 ? end : Levels_.size()-1);
1170 
1171  TEUCHOS_TEST_FOR_EXCEPTION(startLevel > endLevel, Exceptions::RuntimeError,
1172  "MueLu::Hierarchy::Write : startLevel must be <= endLevel");
1173 
1174  TEUCHOS_TEST_FOR_EXCEPTION(startLevel < 0 || endLevel >= Levels_.size(), Exceptions::RuntimeError,
1175  "MueLu::Hierarchy::Write bad start or end level");
1176 
1177  for (LO i = startLevel; i < endLevel + 1; i++) {
1178  RCP<Matrix> A = rcp_dynamic_cast<Matrix>(Levels_[i]-> template Get< RCP< Operator> >("A")), P, R;
1179  if (i > 0) {
1180  P = rcp_dynamic_cast<Matrix>(Levels_[i]-> template Get< RCP< Operator> >("P"));
1181  if (!implicitTranspose_)
1182  R = rcp_dynamic_cast<Matrix>(Levels_[i]-> template Get< RCP< Operator> >("R"));
1183  }
1184 
1185  if (!A.is_null()) Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("A_" + toString(i) + suffix + ".m", *A);
1186  if (!P.is_null()) {
1187  Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("P_" + toString(i) + suffix + ".m", *P);
1188  }
1189  if (!R.is_null()) {
1190  Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("R_" + toString(i) + suffix + ".m", *R);
1191  }
1192  }
1193  }
1194 
1195  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1196  void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Keep(const std::string & ename, const FactoryBase* factory) {
1197  for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1198  (*it)->Keep(ename, factory);
1199  }
1200 
1201  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1202  void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Delete(const std::string& ename, const FactoryBase* factory) {
1203  for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1204  (*it)->Delete(ename, factory);
1205  }
1206 
1207  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1208  void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::AddKeepFlag(const std::string & ename, const FactoryBase* factory, KeepType keep) {
1209  for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1210  (*it)->AddKeepFlag(ename, factory, keep);
1211  }
1212 
1213  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1214  void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::RemoveKeepFlag(const std::string & ename, const FactoryBase* factory, KeepType keep) {
1215  for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1216  (*it)->RemoveKeepFlag(ename, factory, keep);
1217  }
1218 
1219  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1220  std::string Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::description() const {
1221  if ( description_.empty() )
1222  {
1223  std::ostringstream out;
1224  out << BaseClass::description();
1225  out << "{#levels = " << GetGlobalNumLevels() << ", complexity = " << GetOperatorComplexity() << "}";
1226  description_ = out.str();
1227  }
1228  return description_;
1229  }
1230 
1231  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1232  void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::describe(Teuchos::FancyOStream& out, const Teuchos::EVerbosityLevel tVerbLevel) const {
1233  describe(out, toMueLuVerbLevel(tVerbLevel));
1234  }
1235 
1236  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1237  void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::describe(Teuchos::FancyOStream& out, const VerbLevel verbLevel) const {
1238  RCP<Operator> A0 = Levels_[0]->template Get<RCP<Operator> >("A");
1239  RCP<const Teuchos::Comm<int> > comm = A0->getDomainMap()->getComm();
1240 
1241  int numLevels = GetNumLevels();
1242  RCP<Operator> Ac = Levels_[numLevels-1]->template Get<RCP<Operator> >("A");
1243  if (Ac.is_null()) {
1244  // It may happen that we do repartition on the last level, but the matrix
1245  // is small enough to satisfy "max coarse size" requirement. Then, even
1246  // though we have the level, the matrix would be null on all but one processors
1247  numLevels--;
1248  }
1249  int root = comm->getRank();
1250 
1251 #ifdef HAVE_MPI
1252  int smartData = numLevels*comm->getSize() + comm->getRank(), maxSmartData;
1253  reduceAll(*comm, Teuchos::REDUCE_MAX, smartData, Teuchos::ptr(&maxSmartData));
1254  root = maxSmartData % comm->getSize();
1255 #endif
1256 
1257  // Compute smoother complexity, if needed
1258  double smoother_comp =-1.0;
1259  if (verbLevel & (Statistics0 | Test))
1260  smoother_comp = GetSmootherComplexity();
1261 
1262  std::string outstr;
1263  if (comm->getRank() == root && verbLevel & (Statistics0 | Test)) {
1264  std::vector<Xpetra::global_size_t> nnzPerLevel;
1265  std::vector<Xpetra::global_size_t> rowsPerLevel;
1266  std::vector<int> numProcsPerLevel;
1267  bool aborted = false;
1268  for (int i = 0; i < numLevels; i++) {
1269  TEUCHOS_TEST_FOR_EXCEPTION(!(Levels_[i]->IsAvailable("A")) , Exceptions::RuntimeError,
1270  "Operator A is not available on level " << i);
1271 
1272  RCP<Operator> A = Levels_[i]->template Get<RCP<Operator> >("A");
1273  TEUCHOS_TEST_FOR_EXCEPTION(A.is_null(), Exceptions::RuntimeError,
1274  "Operator A on level " << i << " is null.");
1275 
1276  RCP<Matrix> Am = rcp_dynamic_cast<Matrix>(A);
1277  if (Am.is_null()) {
1278  GetOStream(Warnings0) << "Some level operators are not matrices, statistics calculation aborted" << std::endl;
1279  aborted = true;
1280  break;
1281  }
1282 
1283  Xpetra::global_size_t nnz = Am->getGlobalNumEntries();
1284  nnzPerLevel .push_back(nnz);
1285  rowsPerLevel .push_back(Am->getGlobalNumRows());
1286  numProcsPerLevel.push_back(Am->getRowMap()->getComm()->getSize());
1287  }
1288 
1289  if (!aborted) {
1290  std::string label = Levels_[0]->getObjectLabel();
1291  std::ostringstream oss;
1292  oss << std::setfill(' ');
1293  oss << "\n--------------------------------------------------------------------------------\n";
1294  oss << "--- Multigrid Summary " << std::setw(28) << std::left << label << "---\n";
1295  oss << "--------------------------------------------------------------------------------" << std::endl;
1296  oss << "Number of levels = " << numLevels << std::endl;
1297  oss << "Operator complexity = " << std::setprecision(2) << std::setiosflags(std::ios::fixed)
1298  << GetOperatorComplexity() << std::endl;
1299 
1300  if(smoother_comp!=-1.0) {
1301  oss << "Smoother complexity = " << std::setprecision(2) << std::setiosflags(std::ios::fixed)
1302  << smoother_comp << std::endl;
1303  }
1304 
1305  switch (Cycle_) {
1306  case VCYCLE:
1307  oss << "Cycle type = V" << std::endl;
1308  break;
1309  case WCYCLE:
1310  oss << "Cycle type = W" << std::endl;
1311  break;
1312  default:
1313  break;
1314  };
1315  oss << std::endl;
1316 
1317  Xpetra::global_size_t tt = rowsPerLevel[0];
1318  int rowspacer = 2; while (tt != 0) { tt /= 10; rowspacer++; }
1319  tt = nnzPerLevel[0];
1320  int nnzspacer = 2; while (tt != 0) { tt /= 10; nnzspacer++; }
1321  tt = numProcsPerLevel[0];
1322  int npspacer = 2; while (tt != 0) { tt /= 10; npspacer++; }
1323  oss << "level " << std::setw(rowspacer) << " rows " << std::setw(nnzspacer) << " nnz " << " nnz/row" << std::setw(npspacer) << " c ratio" << " procs" << std::endl;
1324  for (size_t i = 0; i < nnzPerLevel.size(); ++i) {
1325  oss << " " << i << " ";
1326  oss << std::setw(rowspacer) << rowsPerLevel[i];
1327  oss << std::setw(nnzspacer) << nnzPerLevel[i];
1328  oss << std::setprecision(2) << std::setiosflags(std::ios::fixed);
1329  oss << std::setw(9) << as<double>(nnzPerLevel[i]) / rowsPerLevel[i];
1330  if (i) oss << std::setw(9) << as<double>(rowsPerLevel[i-1])/rowsPerLevel[i];
1331  else oss << std::setw(9) << " ";
1332  oss << " " << std::setw(npspacer) << numProcsPerLevel[i] << std::endl;
1333  }
1334  oss << std::endl;
1335  for (int i = 0; i < GetNumLevels(); ++i) {
1336  RCP<SmootherBase> preSmoo, postSmoo;
1337  if (Levels_[i]->IsAvailable("PreSmoother"))
1338  preSmoo = Levels_[i]->template Get< RCP<SmootherBase> >("PreSmoother");
1339  if (Levels_[i]->IsAvailable("PostSmoother"))
1340  postSmoo = Levels_[i]->template Get< RCP<SmootherBase> >("PostSmoother");
1341 
1342  if (preSmoo != null && preSmoo == postSmoo)
1343  oss << "Smoother (level " << i << ") both : " << preSmoo->description() << std::endl;
1344  else {
1345  oss << "Smoother (level " << i << ") pre : "
1346  << (preSmoo != null ? preSmoo->description() : "no smoother") << std::endl;
1347  oss << "Smoother (level " << i << ") post : "
1348  << (postSmoo != null ? postSmoo->description() : "no smoother") << std::endl;
1349  }
1350 
1351  oss << std::endl;
1352  }
1353 
1354  outstr = oss.str();
1355  }
1356  }
1357 
1358 #ifdef HAVE_MPI
1359  RCP<const Teuchos::MpiComm<int> > mpiComm = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(comm);
1360  MPI_Comm rawComm = (*mpiComm->getRawMpiComm())();
1361 
1362  int strLength = outstr.size();
1363  MPI_Bcast(&strLength, 1, MPI_INT, root, rawComm);
1364  if (comm->getRank() != root)
1365  outstr.resize(strLength);
1366  MPI_Bcast(&outstr[0], strLength, MPI_CHAR, root, rawComm);
1367 #endif
1368 
1369  out << outstr;
1370  }
1371 
1372  // NOTE: at some point this should be replaced by a friend operator <<
1373  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1374  void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::print(std::ostream& out, const VerbLevel verbLevel) const {
1375  Teuchos::OSTab tab2(out);
1376  for (int i = 0; i < GetNumLevels(); ++i)
1377  Levels_[i]->print(out, verbLevel);
1378  }
1379 
1380  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1381  void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::IsPreconditioner(const bool flag) {
1382  isPreconditioner_ = flag;
1383  }
1384 
1385  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1386  void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DumpCurrentGraph() const {
1387  if (GetProcRankVerbose() != 0)
1388  return;
1389 #if defined(HAVE_MUELU_BOOST) && defined(HAVE_MUELU_BOOST_FOR_REAL) && defined(BOOST_VERSION) && (BOOST_VERSION >= 104400)
1390  BoostGraph graph;
1391 
1392  BoostProperties dp;
1393  dp.property("label", boost::get(boost::vertex_name, graph));
1394  dp.property("id", boost::get(boost::vertex_index, graph));
1395  dp.property("label", boost::get(boost::edge_name, graph));
1396  dp.property("color", boost::get(boost::edge_color, graph));
1397 
1398  // create local maps
1399  std::map<const FactoryBase*, BoostVertex> vindices;
1400  typedef std::map<std::pair<BoostVertex,BoostVertex>, std::string> emap; emap edges;
1401 
1402  static int call_id=0;
1403 
1404  RCP<Operator> A = Levels_[0]->template Get<RCP<Operator> >("A");
1405  int rank = A->getDomainMap()->getComm()->getRank();
1406 
1407  // printf("[%d] CMS: ----------------------\n",rank);
1408  for (int i = dumpLevel_; i <= dumpLevel_+1 && i < GetNumLevels(); i++) {
1409  edges.clear();
1410  Levels_[i]->UpdateGraph(vindices, edges, dp, graph);
1411 
1412  for (emap::const_iterator eit = edges.begin(); eit != edges.end(); eit++) {
1413  std::pair<BoostEdge, bool> boost_edge = boost::add_edge(eit->first.first, eit->first.second, graph);
1414  // printf("[%d] CMS: Hierarchy, adding edge (%d->%d) %d\n",rank,(int)eit->first.first,(int)eit->first.second,(int)boost_edge.second);
1415  // Because xdot.py views 'Graph' as a keyword
1416  if(eit->second==std::string("Graph")) boost::put("label", dp, boost_edge.first, std::string("Graph_"));
1417  else boost::put("label", dp, boost_edge.first, eit->second);
1418  if (i == dumpLevel_)
1419  boost::put("color", dp, boost_edge.first, std::string("red"));
1420  else
1421  boost::put("color", dp, boost_edge.first, std::string("blue"));
1422  }
1423  }
1424 
1425 #if 0
1426  std::ostringstream legend;
1427  legend << "< <TABLE BORDER=\"0\" CELLBORDER=\"1\" CELLSPACING=\"0\" CELLPADDING=\"4\"> \
1428  <TR><TD COLSPAN=\"2\">Legend</TD></TR> \
1429  <TR><TD><FONT color=\"red\">Level " << dumpLevel_ << "</FONT></TD><TD><FONT color=\"blue\">Level " << dumpLevel_+1 << "</FONT></TD></TR> \
1430  </TABLE> >";
1431  BoostVertex boost_vertex = boost::add_vertex(graph);
1432  boost::put("label", dp, boost_vertex, legend.str());
1433 #endif
1434 
1435  std::ofstream out(dumpFile_.c_str() +std::to_string(call_id)+std::string("_")+ std::to_string(rank) + std::string(".dot"));
1436  boost::write_graphviz_dp(out, graph, dp, std::string("id"));
1437  out.close();
1438  call_id++;
1439 #else
1440  GetOStream(Errors) << "Dependency graph output requires boost and MueLu_ENABLE_Boost_for_real" << std::endl;
1441 #endif
1442  }
1443 
1444  // Enforce that coordinate vector's map is consistent with that of A
1445  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1446  void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ReplaceCoordinateMap(Level& level) {
1447  RCP<Operator> Ao = level.Get<RCP<Operator> >("A");
1448  RCP<Matrix> A = rcp_dynamic_cast<Matrix>(Ao);
1449  if (A.is_null()) {
1450  GetOStream(Warnings1) << "Hierarchy::ReplaceCoordinateMap: operator is not a matrix, skipping..." << std::endl;
1451  return;
1452  }
1453  if(Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A) != Teuchos::null) {
1454  GetOStream(Warnings1) << "Hierarchy::ReplaceCoordinateMap: operator is a BlockedCrsMatrix, skipping..." << std::endl;
1455  return;
1456  }
1457 
1458  typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO> xdMV;
1459 
1460  RCP<xdMV> coords = level.Get<RCP<xdMV> >("Coordinates");
1461 
1462  if (A->getRowMap()->isSameAs(*(coords->getMap()))) {
1463  GetOStream(Warnings1) << "Hierarchy::ReplaceCoordinateMap: matrix and coordinates maps are same, skipping..." << std::endl;
1464  return;
1465  }
1466 
1467  if (A->IsView("stridedMaps") && rcp_dynamic_cast<const StridedMap>(A->getRowMap("stridedMaps")) != Teuchos::null) {
1468  RCP<const StridedMap> stridedRowMap = rcp_dynamic_cast<const StridedMap>(A->getRowMap("stridedMaps"));
1469 
1470  // It is better to through an exceptions if maps may be inconsistent, than to ignore it and experience unfathomable breakdowns
1471  TEUCHOS_TEST_FOR_EXCEPTION(stridedRowMap->getStridedBlockId() != -1 || stridedRowMap->getOffset() != 0,
1472  Exceptions::RuntimeError, "Hierarchy::ReplaceCoordinateMap: nontrivial maps (block id = " << stridedRowMap->getStridedBlockId()
1473  << ", offset = " << stridedRowMap->getOffset() << ")");
1474  }
1475 
1476  GetOStream(Runtime1) << "Replacing coordinate map" << std::endl;
1477 
1478  size_t blkSize = A->GetFixedBlockSize();
1479 
1480  RCP<const Map> nodeMap = A->getRowMap();
1481  if (blkSize > 1) {
1482  // Create a nodal map, as coordinates have not been expanded to a DOF map yet.
1483  RCP<const Map> dofMap = A->getRowMap();
1484  GO indexBase = dofMap->getIndexBase();
1485  size_t numLocalDOFs = dofMap->getNodeNumElements();
1486  TEUCHOS_TEST_FOR_EXCEPTION(numLocalDOFs % blkSize, Exceptions::RuntimeError,
1487  "Hierarchy::ReplaceCoordinateMap: block size (" << blkSize << ") is incompatible with the number of local dofs in a row map (" << numLocalDOFs);
1488  ArrayView<const GO> GIDs = dofMap->getNodeElementList();
1489 
1490  Array<GO> nodeGIDs(numLocalDOFs/blkSize);
1491  for (size_t i = 0; i < numLocalDOFs; i += blkSize)
1492  nodeGIDs[i/blkSize] = (GIDs[i] - indexBase)/blkSize + indexBase;
1493 
1494  Xpetra::global_size_t INVALID = Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid();
1495  nodeMap = MapFactory::Build(dofMap->lib(), INVALID, nodeGIDs(), indexBase, dofMap->getComm());
1496  } else {
1497  // blkSize == 1
1498  // Check whether the length of vectors fits to the size of A
1499  // If yes, make sure that the maps are matching
1500  // If no, throw a warning but do not touch the Coordinates
1501  if(coords->getLocalLength() != A->getRowMap()->getNodeNumElements()) {
1502  GetOStream(Warnings) << "Coordinate vector does not match row map of matrix A!" << std::endl;
1503  return;
1504  }
1505  }
1506 
1507  Array<ArrayView<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> > coordDataView;
1508  std::vector<ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> > coordData;
1509  for (size_t i = 0; i < coords->getNumVectors(); i++) {
1510  coordData.push_back(coords->getData(i));
1511  coordDataView.push_back(coordData[i]());
1512  }
1513 
1514  RCP<xdMV> newCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO>::Build(nodeMap, coordDataView(), coords->getNumVectors());
1515  level.Set("Coordinates", newCoords);
1516  }
1517 
1518  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1519  void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::AllocateLevelMultiVectors(int numvecs) {
1520  int N = Levels_.size();
1521  if( (sizeOfAllocatedLevelMultiVectors_ == numvecs && residual_.size() == N) || numvecs<=0 ) return;
1522 
1523  // If, somehow, we changed the number of levels, delete everything first
1524  if(residual_.size() != N) {
1525  DeleteLevelMultiVectors();
1526 
1527  residual_.resize(N);
1528  coarseRhs_.resize(N);
1529  coarseX_.resize(N);
1530  coarseImport_.resize(N);
1531  coarseExport_.resize(N);
1532  correction_.resize(N);
1533  }
1534 
1535  for(int i=0; i<N; i++) {
1536  RCP<Operator> A = Levels_[i]->template Get< RCP<Operator> >("A");
1537  if(!A.is_null()) {
1538  // This dance is because we allow A to have a BlockedMap and X/B to have (compatible) non-blocked map
1539  RCP<const BlockedCrsMatrix> A_as_blocked = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(A);
1540  RCP<const Map> Arm = A->getRangeMap();
1541  RCP<const Map> Adm = A->getDomainMap();
1542  if(!A_as_blocked.is_null()) {
1543  Adm = A_as_blocked->getFullDomainMap();
1544  }
1545 
1546  // This is zero'd by default since it is filled via an operator apply
1547  residual_[i] = MultiVectorFactory::Build(Arm, numvecs, true);
1548  correction_[i] = MultiVectorFactory::Build(Adm, numvecs, false);
1549  }
1550 
1551  if(i+1<N) {
1552  // This is zero'd by default since it is filled via an operator apply
1553  if(implicitTranspose_) {
1554  RCP<Operator> P = Levels_[i+1]->template Get< RCP<Operator> >("P");
1555  if(!P.is_null()) coarseRhs_[i] = MultiVectorFactory::Build(P->getDomainMap(),numvecs,true);
1556  } else {
1557  RCP<Operator> R = Levels_[i+1]->template Get< RCP<Operator> >("R");
1558  if(!R.is_null()) coarseRhs_[i] = MultiVectorFactory::Build(R->getRangeMap(),numvecs,true);
1559  }
1560 
1561 
1562  RCP<const Import> importer;
1563  if(Levels_[i+1]->IsAvailable("Importer"))
1564  importer = Levels_[i+1]->template Get< RCP<const Import> >("Importer");
1565  if (doPRrebalance_ || importer.is_null())
1566  coarseX_[i] = MultiVectorFactory::Build(coarseRhs_[i]->getMap(),numvecs,false);
1567  else {
1568  coarseImport_[i] = MultiVectorFactory::Build(importer->getTargetMap(), numvecs,false);
1569  coarseExport_[i] = MultiVectorFactory::Build(importer->getSourceMap(), numvecs,false);
1570  coarseX_[i] = MultiVectorFactory::Build(importer->getTargetMap(),numvecs,false);
1571  }
1572  }
1573  }
1574  sizeOfAllocatedLevelMultiVectors_ = numvecs;
1575  }
1576 
1577 
1578 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1579 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeleteLevelMultiVectors() {
1580  if(sizeOfAllocatedLevelMultiVectors_==0) return;
1581  residual_.resize(0);
1582  coarseRhs_.resize(0);
1583  coarseX_.resize(0);
1584  coarseImport_.resize(0);
1585  coarseExport_.resize(0);
1586  correction_.resize(0);
1587  sizeOfAllocatedLevelMultiVectors_ = 0;
1588 }
1589 
1590 
1591 } //namespace MueLu
1592 
1593 #endif // MUELU_HIERARCHY_DEF_HPP
MueLu_BoostGraphviz.hpp
MueLu_Hierarchy_decl.hpp
MueLu_Level.hpp
MueLu_Monitor.hpp
MueLu_SmootherBase.hpp
MueLu_SmootherFactoryBase.hpp
MueLu::Describable::description
virtual std::string description() const
Return a simple one-line description of this object.
Definition: MueLu_Describable.cpp:61
MueLu::Exceptions::Incompatible
Exception throws to report incompatible objects (like maps).
Definition: MueLu_Exceptions.hpp:82
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: MueLu_Exceptions.hpp:70
MueLu::FactoryBase
Base class for factories (e.g., R, P, and A_coarse).
Definition: MueLu_FactoryBase.hpp:60
MueLu::FactoryManagerBase
Class that provides default factories within Needs class.
Definition: MueLu_FactoryManagerBase.hpp:63
MueLu::Hierarchy
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
Definition: MueLu_Hierarchy_decl.hpp:104
MueLu::Hierarchy::AddLevel
void AddLevel(const RCP< Level > &level)
Add a level at the end of the hierarchy.
Definition: MueLu_Hierarchy_def.hpp:125
MueLu::Hierarchy::GetSmootherComplexity
double GetSmootherComplexity() const
Definition: MueLu_Hierarchy_def.hpp:196
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:1167
MueLu::Hierarchy::GetLevel
RCP< Level > & GetLevel(const int levelID=0)
Retrieve a certain level from hierarchy.
Definition: MueLu_Hierarchy_def.hpp:149
MueLu::Hierarchy::CheckLevel
void CheckLevel(Level &level, int levelID)
Helper function.
Definition: MueLu_Hierarchy_def.hpp:236
MueLu::Hierarchy::description
std::string description() const
Return a simple one-line description of this object.
Definition: MueLu_Hierarchy_def.hpp:1220
MueLu::Hierarchy::IsPreconditioner
void IsPreconditioner(const bool flag)
Definition: MueLu_Hierarchy_def.hpp:1381
MueLu::Hierarchy::Levels_
Array< RCP< Level > > Levels_
Container for Level objects.
Definition: MueLu_Hierarchy_decl.hpp:364
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::Hierarchy::MagnitudeType
STS::magnitudeType MagnitudeType
Definition: MueLu_Hierarchy_decl.hpp:109
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:1237
MueLu::Hierarchy::SetMatvecParams
void SetMatvecParams(RCP< ParameterList > matvecParams)
Definition: MueLu_Hierarchy_def.hpp:246
MueLu::Hierarchy::lib_
Xpetra::UnderlyingLib lib_
Definition: MueLu_Hierarchy_decl.hpp:398
MueLu::Hierarchy::Clear
void Clear(int startLevel=0)
Clear impermanent data from previous setup.
Definition: MueLu_Hierarchy_def.hpp:638
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::Hierarchy::GetOperatorComplexity
double GetOperatorComplexity() const
Definition: MueLu_Hierarchy_def.hpp:173
MueLu::Hierarchy::print
void print(std::ostream &out=std::cout, const VerbLevel verbLevel=(MueLu::Parameters|MueLu::Statistics0)) const
Definition: MueLu_Hierarchy_def.hpp:1374
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:1202
MueLu::Hierarchy::GetGlobalNumLevels
int GetGlobalNumLevels() const
Definition: MueLu_Hierarchy_def.hpp:161
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:1196
MueLu::Hierarchy::AllocateLevelMultiVectors
void AllocateLevelMultiVectors(int numvecs)
Definition: MueLu_Hierarchy_def.hpp:1519
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:1208
MueLu::Hierarchy::AddNewLevel
void AddNewLevel()
Add a new level at the end of the hierarchy.
Definition: MueLu_Hierarchy_def.hpp:142
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:1214
MueLu::Hierarchy::ReplaceCoordinateMap
void ReplaceCoordinateMap(Level &level)
Definition: MueLu_Hierarchy_def.hpp:1446
MueLu::Level
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
MueLu::Level::IsAvailable
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
Definition: MueLu_Level.hpp:333
MueLu::Level::SetComm
void SetComm(RCP< const Teuchos::Comm< int > > const &comm)
Definition: MueLu_Level.hpp:381
MueLu::Level::setlib
void setlib(Xpetra::UnderlyingLib lib2)
Definition: MueLu_Level.hpp:378
MueLu::Level::Release
void Release(const FactoryBase &factory)
Decrement the storage counter for all the inputs of a factory.
Definition: MueLu_Level.cpp:143
MueLu::Level::GetLevelID
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
MueLu::Level::GetComm
RCP< const Teuchos::Comm< int > > GetComm() const
Definition: MueLu_Level.hpp:382
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)....
Definition: MueLu_Level.hpp:192
MueLu::Level::Set
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
Definition: MueLu_Level.hpp:157
MueLu::Level::Request
void Request(const FactoryBase &factory)
Increment the storage counter for all the inputs of a factory.
Definition: MueLu_Level.cpp:136
MueLu::Level::lib
Xpetra::UnderlyingLib lib()
Definition: MueLu_Level.hpp:379
MueLu::Level::GetPreviousLevel
RCP< Level > & GetPreviousLevel()
Previous level.
Definition: MueLu_Level.hpp:132
MueLu::Monitor
Timer to be used in non-factories.
Definition: MueLu_Monitor.hpp:107
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:102
MueLu::SetFactoryManager
An exception safe way to call the method 'Level::SetFactoryManager()'.
Definition: MueLu_HierarchyUtils_decl.hpp:68
MueLu::TimeMonitor
Integrates Teuchos::TimeMonitor with MueLu verbosity system.
Definition: MueLu_TimeMonitor.hpp:68
MueLu::TopRAPFactory::Build
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Definition: MueLu_TopRAPFactory_def.hpp:55
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:212
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:214
MueLu::Utilities::SetRandomSeed
static void SetRandomSeed(const Teuchos::Comm< int > &comm)
Definition: MueLu_Utilities_decl.hpp:222
MueLu
Namespace for MueLu classes and methods.
Definition: MueLu_BrickAggregationFactory_decl.hpp:76
MueLu::Warnings0
@ Warnings0
Important warning messages (one line)
Definition: MueLu_VerbosityLevel.hpp:57
MueLu::Statistics2
@ Statistics2
Print even more statistics.
Definition: MueLu_VerbosityLevel.hpp:72
MueLu::Warnings
@ Warnings
Print all warning messages.
Definition: MueLu_VerbosityLevel.hpp:85
MueLu::Errors
@ Errors
Errors.
Definition: MueLu_VerbosityLevel.hpp:55
MueLu::Statistics1
@ Statistics1
Print more statistics.
Definition: MueLu_VerbosityLevel.hpp:71
MueLu::Timings0
@ Timings0
High level timing information (use Teuchos::TimeMonitor::summarize() to print)
Definition: MueLu_VerbosityLevel.hpp:74
MueLu::Runtime0
@ Runtime0
One-liner description of what is happening.
Definition: MueLu_VerbosityLevel.hpp:62
MueLu::Runtime1
@ Runtime1
Description of what is happening (more verbose)
Definition: MueLu_VerbosityLevel.hpp:63
MueLu::Warnings1
@ Warnings1
Additional warnings.
Definition: MueLu_VerbosityLevel.hpp:59
MueLu::Statistics0
@ Statistics0
Print statistics that do not involve significant additional computation.
Definition: MueLu_VerbosityLevel.hpp:70
MueLu::KeepType
short KeepType
Definition: MueLu_KeepType.hpp:63
MueLu::toString
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Definition: MueLu_Utilities_decl.hpp:1000
MueLu::toMueLuVerbLevel
VerbLevel toMueLuVerbLevel(const Teuchos::EVerbosityLevel verbLevel)
Definition: MueLu_VerbosityLevel.cpp:51
MueLu::FormattingHelper::getColonLabel
static std::string getColonLabel(const std::string &label)
Helper function for object label.
Definition: MueLu_Monitor.hpp:66
MueLu::Hierarchy::ConvData
Data struct for defining stopping criteria of multigrid iteration.
Definition: MueLu_Hierarchy_decl.hpp:112