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 #ifndef MUELU_HIERARCHY_DEF_HPP
47 #define MUELU_HIERARCHY_DEF_HPP
48 
49 #include <algorithm>
50 #include <sstream>
51 
52 #include <Xpetra_Matrix.hpp>
53 #include <Xpetra_MultiVectorFactory.hpp>
54 #include <Xpetra_Operator.hpp>
55 #include <Xpetra_IO.hpp>
56 
57 #include "MueLu_Hierarchy_decl.hpp"
58 
59 #include "MueLu_BoostGraphviz.hpp"
60 #include "MueLu_FactoryManager.hpp"
61 #include "MueLu_HierarchyHelpers.hpp"
62 #include "MueLu_Level.hpp"
63 #include "MueLu_Monitor.hpp"
64 #include "MueLu_PFactory.hpp"
66 #include "MueLu_SmootherFactory.hpp"
67 #include "MueLu_SmootherBase.hpp"
68 
69 namespace MueLu {
70 
71  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
73  : maxCoarseSize_(GetDefaultMaxCoarseSize()), implicitTranspose_(GetDefaultImplicitTranspose()),
74  doPRrebalance_(GetDefaultPRrebalance()), isPreconditioner_(true), Cycle_(GetDefaultCycle()),
75  lib_(Xpetra::UseTpetra), isDumpingEnabled_(false), dumpLevel_(-1), rate_(-1)
76  {
77  AddLevel(rcp(new Level));
78  }
79 
80  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
84  isDumpingEnabled_(false), dumpLevel_(-1), rate_(-1)
85  {
86  lib_ = A->getDomainMap()->lib();
87 
88  RCP<Level> Finest = rcp(new Level);
89  AddLevel(Finest);
90 
91  Finest->Set("A", A);
92  }
93 
94  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
96  int levelID = LastLevelID() + 1; // ID of the inserted level
97 
98  if (level->GetLevelID() != -1 && (level->GetLevelID() != levelID))
99  GetOStream(Warnings1) << "Hierarchy::AddLevel(): Level with ID=" << level->GetLevelID() <<
100  " have been added at the end of the hierarchy\n but its ID have been redefined" <<
101  " because last level ID of the hierarchy was " << LastLevelID() << "." << std::endl;
102 
103  Levels_.push_back(level);
104  level->SetLevelID(levelID);
105  level->setlib(lib_);
106 
107  level->SetPreviousLevel( (levelID == 0) ? Teuchos::null : Levels_[LastLevelID() - 1] );
108  }
109 
110  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
112  RCP<Level> newLevel = Levels_[LastLevelID()]->Build(); // new coarse level, using copy constructor
113  newLevel->setlib(lib_);
114  this->AddLevel(newLevel); // add to hierarchy
115  }
116 
117  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
119  TEUCHOS_TEST_FOR_EXCEPTION(levelID < 0 || levelID > LastLevelID(), Exceptions::RuntimeError,
120  "MueLu::Hierarchy::GetLevel(): invalid input parameter value: LevelID = " << levelID);
121  return Levels_[levelID];
122  }
123 
124  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
126  return Levels_.size();
127  }
128 
129  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
131  RCP<Operator> A = Levels_[0]->template Get<RCP<Operator> >("A");
132  RCP<const Teuchos::Comm<int> > comm = A->getDomainMap()->getComm();
133 
134  int numLevels = GetNumLevels();
135  int numGlobalLevels;
136  Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, numLevels, Teuchos::ptr(&numGlobalLevels));
137 
138  return numGlobalLevels;
139  }
140 
141  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
143  double totalNnz = 0, lev0Nnz = 1;
144  for (int i = 0; i < GetNumLevels(); ++i) {
145  TEUCHOS_TEST_FOR_EXCEPTION(!(Levels_[i]->IsAvailable("A")) , Exceptions::RuntimeError,
146  "Operator complexity cannot be calculated because A is unavailable on level " << i);
147  RCP<Operator> A = Levels_[i]->template Get<RCP<Operator> >("A");
148  if (A.is_null())
149  break;
150 
151  RCP<Matrix> Am = rcp_dynamic_cast<Matrix>(A);
152  if (Am.is_null()) {
153  GetOStream(Warnings0) << "Some level operators are not matrices, operator complexity calculation aborted" << std::endl;
154  return 0.0;
155  }
156 
157  totalNnz += as<double>(Am->getGlobalNumEntries());
158  if (i == 0)
159  lev0Nnz = totalNnz;
160  }
161  return totalNnz / lev0Nnz;
162  }
163 
164  // Coherence checks todo in Setup() (using an helper function):
165  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
167  TEUCHOS_TEST_FOR_EXCEPTION(level.lib() != lib_, Exceptions::RuntimeError,
168  "MueLu::Hierarchy::CheckLevel(): wrong underlying linear algebra library.");
169  TEUCHOS_TEST_FOR_EXCEPTION(level.GetLevelID() != levelID, Exceptions::RuntimeError,
170  "MueLu::Hierarchy::CheckLevel(): wrong level ID");
171  TEUCHOS_TEST_FOR_EXCEPTION(levelID != 0 && level.GetPreviousLevel() != Levels_[levelID-1], Exceptions::RuntimeError,
172  "MueLu::Hierarchy::Setup(): wrong level parent");
173  }
174 
175  // The function uses three managers: fine, coarse and next coarse
176  // We construct the data for the coarse level, and do requests for the next coarse
177  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
179  const RCP<const FactoryManagerBase> fineLevelManager,
180  const RCP<const FactoryManagerBase> coarseLevelManager,
181  const RCP<const FactoryManagerBase> nextLevelManager) {
182  // Use PrintMonitor/TimerMonitor instead of just a FactoryMonitor to print "Level 0" instead of Hierarchy(0)
183  // Print is done after the requests for next coarse level
184  TimeMonitor m1(*this, this->ShortClassName() + ": " + "Setup (total)");
185  TimeMonitor m2(*this, this->ShortClassName() + ": " + "Setup" + " (total, level=" + Teuchos::toString(coarseLevelID) + ")");
186 
187  // TODO: pass coarseLevelManager by reference
188  TEUCHOS_TEST_FOR_EXCEPTION(coarseLevelManager == Teuchos::null, Exceptions::RuntimeError,
189  "MueLu::Hierarchy::Setup(): argument coarseLevelManager cannot be null");
190 
193 
194  TEUCHOS_TEST_FOR_EXCEPTION(LastLevelID() < coarseLevelID, Exceptions::RuntimeError,
195  "MueLu::Hierarchy:Setup(): level " << coarseLevelID << " (specified by coarseLevelID argument) "
196  "must be built before calling this function.");
197 
198  Level& level = *Levels_[coarseLevelID];
199 
200  if (levelManagers_.size() < coarseLevelID+1)
201  levelManagers_.resize(coarseLevelID+1);
202  levelManagers_[coarseLevelID] = coarseLevelManager;
203 
204  bool isFinestLevel = (fineLevelManager.is_null());
205  bool isLastLevel = (nextLevelManager.is_null());
206 
207  int oldRank = -1;
208  if (isFinestLevel) {
209  RCP<Operator> A = level.Get< RCP<Operator> >("A");
210  RCP<const Map> domainMap = A->getDomainMap();
211  RCP<const Teuchos::Comm<int> > comm = domainMap->getComm();
212 
213  // Initialize random seed for reproducibility
215 
216  // Record the communicator on the level (used for timers sync)
217  level.SetComm(comm);
218  oldRank = SetProcRankVerbose(comm->getRank());
219 
220  // Set the Hierarchy library to match that of the finest level matrix,
221  // even if it was already set
222  lib_ = domainMap->lib();
223  level.setlib(lib_);
224 
225  } else {
226  // Permeate library to a coarser level
227  level.setlib(lib_);
228 
229  Level& prevLevel = *Levels_[coarseLevelID-1];
230  oldRank = SetProcRankVerbose(prevLevel.GetComm()->getRank());
231  }
232 
233  CheckLevel(level, coarseLevelID);
234 
235  // Attach FactoryManager to the fine level
236  RCP<SetFactoryManager> SFMFine;
237  if (!isFinestLevel)
238  SFMFine = rcp(new SetFactoryManager(Levels_[coarseLevelID-1], fineLevelManager));
239 
240  if (isFinestLevel && Levels_[coarseLevelID]->IsAvailable("Coordinates"))
241  ReplaceCoordinateMap(*Levels_[coarseLevelID]);
242 
243  // Attach FactoryManager to the coarse level
244  SetFactoryManager SFMCoarse(Levels_[coarseLevelID], coarseLevelManager);
245 
246  if (isDumpingEnabled_ && dumpLevel_ == 0 && coarseLevelID == 1)
248 
249  RCP<TopSmootherFactory> coarseFact = rcp(new TopSmootherFactory(coarseLevelManager, "CoarseSolver"));
250  RCP<TopSmootherFactory> smootherFact = rcp(new TopSmootherFactory(coarseLevelManager, "Smoother"));
251 
252  int nextLevelID = coarseLevelID + 1;
253 
254  RCP<SetFactoryManager> SFMNext;
255  if (isLastLevel == false) {
256  // We are not at the coarsest level, so there is going to be another level ("next coarse") after this one ("coarse")
257  if (nextLevelID > LastLevelID())
258  AddNewLevel();
259  CheckLevel(*Levels_[nextLevelID], nextLevelID);
260 
261  // Attach FactoryManager to the next level (level after coarse)
262  SFMNext = rcp(new SetFactoryManager(Levels_[nextLevelID], nextLevelManager));
263  Levels_[nextLevelID]->Request(TopRAPFactory(coarseLevelManager, nextLevelManager));
264 
265  // Do smoother requests here. We don't know whether this is going to be
266  // the coarsest level or not, but we need to DeclareInput before we call
267  // coarseRAPFactory.Build(), otherwise some stuff may be erased after
268  // level releases
269  level.Request(*smootherFact);
270 
271  } else {
272  // Similar to smoother above, do the coarse solver request here. We don't
273  // know whether this is going to be the coarsest level or not, but we
274  // need to DeclareInput before we call coarseRAPFactory.Build(),
275  // otherwise some stuff may be erased after level releases. This is
276  // actually evident on ProjectorSmoother. It requires both "A" and
277  // "Nullspace". However, "Nullspace" is erased after all releases, so if
278  // we call the coarse factory request after RAP build we would not have
279  // any data, and cannot get it as we don't have previous managers. The
280  // typical trace looks like this:
281  //
282  // MueLu::Level(0)::GetFactory(Aggregates, 0): No FactoryManager
283  // during request for data " Aggregates" on level 0 by factory TentativePFactory
284  // during request for data " P" on level 1 by factory EminPFactory
285  // during request for data " P" on level 1 by factory TransPFactory
286  // during request for data " R" on level 1 by factory RAPFactory
287  // during request for data " A" on level 1 by factory TentativePFactory
288  // during request for data " Nullspace" on level 2 by factory NullspaceFactory
289  // during request for data " Nullspace" on level 2 by factory NullspacePresmoothFactory
290  // during request for data " Nullspace" on level 2 by factory ProjectorSmoother
291  // during request for data " PreSmoother" on level 2 by factory NoFactory
292  level.Request(*coarseFact);
293  }
294 
295  PrintMonitor m0(*this, "Level " + Teuchos::toString(coarseLevelID), static_cast<MsgType>(GetVerbLevel()));
296 
297  // Build coarse level hierarchy
298  RCP<Operator> Ac = Teuchos::null;
299  TopRAPFactory coarseRAPFactory(fineLevelManager, coarseLevelManager);
300 
301  if (level.IsAvailable("A")) {
302  Ac = level.Get<RCP<Operator> >("A");
303  } else if (!isFinestLevel) {
304  // We only build here, the release is done later
305  coarseRAPFactory.Build(*level.GetPreviousLevel(), level);
306  }
307 
308  if (level.IsAvailable("A"))
309  Ac = level.Get<RCP<Operator> >("A");
310  RCP<Matrix> Acm = rcp_dynamic_cast<Matrix>(Ac);
311 
312  // Record the communicator on the level
313  if (!Ac.is_null())
314  level.SetComm(Ac->getDomainMap()->getComm());
315 
316  // Test if we reach the end of the hierarchy
317  bool isOrigLastLevel = isLastLevel;
318  if (isLastLevel) {
319  // Last level as we have achieved the max limit
320  isLastLevel = true;
321 
322  } else if (Ac.is_null()) {
323  // Last level for this processor, as it does not belong to the next
324  // subcommunicator. Other processors may continue working on the
325  // hierarchy
326  isLastLevel = true;
327 
328  } else {
329  if (!Acm.is_null() && Acm->getGlobalNumRows() <= maxCoarseSize_) {
330  // Last level as the size of the coarse matrix became too small
331  GetOStream(Runtime0) << "Max coarse size (<= " << maxCoarseSize_ << ") achieved" << std::endl;
332  isLastLevel = true;
333  }
334  }
335 
336  if (!Ac.is_null() && !isFinestLevel) {
337  RCP<Operator> A = Levels_[coarseLevelID-1]->template Get< RCP<Operator> >("A");
338  RCP<Matrix> Am = rcp_dynamic_cast<Matrix>(A);
339 
340  const double maxCoarse2FineRatio = 0.8;
341  if (!Acm.is_null() && !Am.is_null() && Acm->getGlobalNumRows() > maxCoarse2FineRatio * Am->getGlobalNumRows()) {
342  // We could abort here, but for now we simply notify user.
343  // Couple of additional points:
344  // - if repartitioning is delayed until level K, but the aggregation
345  // procedure stagnates between levels K-1 and K. In this case,
346  // repartitioning could enable faster coarsening once again, but the
347  // hierarchy construction will abort due to the stagnation check.
348  // - if the matrix is small enough, we could move it to one processor.
349  GetOStream(Warnings0) << "Aggregation stagnated. Please check your matrix and/or adjust your configuration file."
350  << "Possible fixes:\n"
351  << " - reduce the maximum number of levels\n"
352  << " - enable repartitioning\n"
353  << " - increase the minimum coarse size." << std::endl;
354 
355  }
356  }
357 
358  if (isLastLevel) {
359  if (!isOrigLastLevel) {
360  // We did not expect to finish this early so we did request a smoother.
361  // We need a coarse solver instead. Do the magic.
362  level.Release(*smootherFact);
363  level.Request(*coarseFact);
364  }
365 
366  // Do the actual build, if we have any data.
367  // NOTE: this is not a great check, we may want to call Build() regardless.
368  if (!Ac.is_null())
369  coarseFact->Build(level);
370 
371  // Once the dirty deed is done, release stuff. The smoother has already
372  // been released.
373  level.Release(*coarseFact);
374 
375  } else {
376  // isLastLevel = false => isOrigLastLevel = false, meaning that we have
377  // requested the smoother. Now we need to build it and to release it.
378  // We don't need to worry about the coarse solver, as we didn't request it.
379  if (!Ac.is_null())
380  smootherFact->Build(level);
381 
382  level.Release(*smootherFact);
383  }
384 
385  if (isLastLevel == true) {
386  if (isOrigLastLevel == false) {
387  // Earlier in the function, we constructed the next coarse level, and requested data for the that level,
388  // assuming that we are not at the coarsest level. Now, we changed our mind, so we have to release those.
389  Levels_[nextLevelID]->Release(TopRAPFactory(coarseLevelManager, nextLevelManager));
390  }
391  Levels_.resize(nextLevelID);
392  }
393 
394  // I think this is the proper place for graph so that it shows every dependence
395  if (isDumpingEnabled_ && dumpLevel_ > 0 && coarseLevelID == dumpLevel_)
397 
398  if (!isFinestLevel) {
399  // Release the hierarchy data
400  // We release so late to help blocked solvers, as the smoothers for them need A blocks
401  // which we construct in RAPFactory
402  level.Release(coarseRAPFactory);
403  }
404 
405  if (oldRank != -1)
406  SetProcRankVerbose(oldRank);
407 
408  return isLastLevel;
409  }
410 
411  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
413  int numLevels = Levels_.size();
414  TEUCHOS_TEST_FOR_EXCEPTION(levelManagers_.size() != numLevels, Exceptions::RuntimeError,
415  "Hierarchy::SetupRe: " << Levels_.size() << " levels, but " << levelManagers_.size() << " level factory managers");
416 
417 #ifdef HAVE_MUELU_DEBUG
418  // Reset factories' data used for debugging
419  for (int i = 0; i < numLevels; i++)
420  levelManagers_[i]->ResetDebugData();
421 
422 #endif
423 
424  int levelID;
425  for (levelID = 0; levelID < numLevels;) {
426  bool r = Setup(levelID,
427  (levelID != 0 ? levelManagers_[levelID-1] : Teuchos::null),
428  levelManagers_[levelID],
429  (levelID+1 != numLevels ? levelManagers_[levelID+1] : Teuchos::null));
430  levelID++;
431  if (r) break;
432  }
433  // We may construct fewer levels for some reason, make sure we continue
434  // doing that in the future
435  Levels_ .resize(levelID);
436  levelManagers_.resize(levelID);
437 
439  }
440 
441  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
442  void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Setup(const FactoryManagerBase& manager, int startLevel, int numDesiredLevels) {
443  // Use MueLu::BaseClass::description() to avoid printing "{numLevels = 1}" (numLevels is increasing...)
444  PrintMonitor m0(*this, "Setup (" + this->MueLu::BaseClass::description() + ")");
445 
446  Clear(startLevel);
447 
448  // Check Levels_[startLevel] exists.
449  TEUCHOS_TEST_FOR_EXCEPTION(Levels_.size() <= startLevel, Exceptions::RuntimeError,
450  "MueLu::Hierarchy::Setup(): fine level (" << startLevel << ") does not exist");
451 
452  TEUCHOS_TEST_FOR_EXCEPTION(numDesiredLevels <= 0, Exceptions::RuntimeError,
453  "Constructing non-positive (" << numDesiredLevels << ") number of levels does not make sense.");
454 
455  // Check for fine level matrix A
456  TEUCHOS_TEST_FOR_EXCEPTION(!Levels_[startLevel]->IsAvailable("A"), Exceptions::RuntimeError,
457  "MueLu::Hierarchy::Setup(): fine level (" << startLevel << ") has no matrix A! "
458  "Set fine level matrix A using Level.Set()");
459 
460  RCP<Operator> A = Levels_[startLevel]->template Get<RCP<Operator> >("A");
461  lib_ = A->getDomainMap()->lib();
462 
463  RCP<const FactoryManagerBase> rcpmanager = rcpFromRef(manager);
464 
465  const int lastLevel = startLevel + numDesiredLevels - 1;
466  GetOStream(Runtime0) << "Setup loop: startLevel = " << startLevel << ", lastLevel = " << lastLevel
467  << " (stop if numLevels = " << numDesiredLevels << " or Ac.size() < " << maxCoarseSize_ << ")" << std::endl;
468 
469  // Setup multigrid levels
470  int iLevel = 0;
471  if (numDesiredLevels == 1) {
472  iLevel = 0;
473  Setup(startLevel, Teuchos::null, rcpmanager, Teuchos::null); // setup finest==coarsest level (first and last managers are Teuchos::null)
474 
475  } else {
476  bool bIsLastLevel = Setup(startLevel, Teuchos::null, rcpmanager, rcpmanager); // setup finest level (level 0) (first manager is Teuchos::null)
477  if (bIsLastLevel == false) {
478  for (iLevel = startLevel + 1; iLevel < lastLevel; iLevel++) {
479  bIsLastLevel = Setup(iLevel, rcpmanager, rcpmanager, rcpmanager); // setup intermediate levels
480  if (bIsLastLevel == true)
481  break;
482  }
483  if (bIsLastLevel == false)
484  Setup(lastLevel, rcpmanager, rcpmanager, Teuchos::null); // setup coarsest level (last manager is Teuchos::null)
485  }
486  }
487 
488  // TODO: some check like this should be done at the beginning of the routine
489  TEUCHOS_TEST_FOR_EXCEPTION(iLevel != Levels_.size() - 1, Exceptions::RuntimeError,
490  "MueLu::Hierarchy::Setup(): number of level");
491 
492  // TODO: this is not exception safe: manager will still hold default
493  // factories if you exit this function with an exception
494  manager.Clean();
495 
497  }
498 
499  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
501  if (startLevel < GetNumLevels())
502  GetOStream(Runtime0) << "Clearing old data (if any)" << std::endl;
503 
504  for (int iLevel = startLevel; iLevel < GetNumLevels(); iLevel++)
505  Levels_[iLevel]->Clear();
506  }
507 
508  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
510  GetOStream(Runtime0) << "Clearing old data (expert)" << std::endl;
511  for (int iLevel = 0; iLevel < GetNumLevels(); iLevel++)
512  Levels_[iLevel]->ExpertClear();
513  }
514 
515  // ---------------------------------------- Iterate -------------------------------------------------------
516 
517  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
519  bool InitialGuessIsZero, LO startLevel) {
520  LO nIts = conv.maxIts_;
521  MagnitudeType tol = conv.tol_;
522 
523  // These timers work as follows. "iterateTime" records total time spent in
524  // iterate. "levelTime" records time on a per level basis. The label is
525  // crafted to mimic the per-level messages used in Monitors. Note that a
526  // given level is timed with a TimeMonitor instead of a Monitor or
527  // SubMonitor. This is mainly because I want to time each level
528  // separately, and Monitors/SubMonitors print "(total) xx yy zz" ,
529  // "(sub,total) xx yy zz", respectively, which is subject to
530  // misinterpretation. The per-level TimeMonitors are stopped/started
531  // manually before/after a recursive call to Iterate. A side artifact to
532  // this approach is that the counts for intermediate level timers are twice
533  // the counts for the finest and coarsest levels.
534  std::string prefix = this->ShortClassName() + ": ";
535  std::string levelSuffix = " (level=" + toString(startLevel) + ")";
536  RCP<Monitor> iterateTime;
537  RCP<TimeMonitor> iterateTime1;
538  if (startLevel == 0)
539  iterateTime = rcp(new Monitor(*this, "Solve", (nIts == 1) ? None : Runtime0, Timings0));
540  else
541  iterateTime1 = rcp(new TimeMonitor(*this, prefix + "Solve (total, level=" + toString(startLevel) + ")", Timings0));
542 
543  std::string iterateLevelTimeLabel = prefix + "Solve" + levelSuffix;
544  RCP<TimeMonitor> iterateLevelTime = rcp(new TimeMonitor(*this, iterateLevelTimeLabel, Timings0));
545 
546  bool zeroGuess = InitialGuessIsZero;
547 
548  RCP<Level> Fine = Levels_[startLevel];
549  RCP<Operator> A = Fine->Get< RCP<Operator> >("A");
550 
551  if (A.is_null()) {
552  // This processor does not have any data for this process on coarser
553  // levels. This can only happen when there are multiple processors and
554  // we use repartitioning.
555  return Undefined;
556  }
557 
558  // Print residual information before iterating
559  MagnitudeType prevNorm = STS::magnitude(STS::one()), curNorm = STS::magnitude(STS::one());
560  rate_ = 1.0;
561  if (startLevel == 0 && !isPreconditioner_ &&
562  (IsPrint(Statistics1) || tol > 0)) {
563  // We calculate the residual only if we want to print it out, or if we
564  // want to stop once we achive the tolerance
565  Teuchos::Array<MagnitudeType> rn;
566  rn = Utilities::ResidualNorm(*A, X, B);
567 
568  if (tol > 0) {
569  bool passed = true;
570  for (LO k = 0; k < rn.size(); k++)
571  if (rn[k] >= tol)
572  passed = false;
573 
574  if (passed)
575  return Converged;
576  }
577 
578  if (IsPrint(Statistics1))
579  GetOStream(Statistics1) << "iter: "
580  << std::setiosflags(std::ios::left)
581  << std::setprecision(3) << 0 // iter 0
582  << " residual = "
583  << std::setprecision(10) << rn
584  << std::endl;
585  }
586 
587  SC one = STS::one(), zero = STS::zero();
588  for (LO i = 1; i <= nIts; i++) {
589 #ifdef HAVE_MUELU_DEBUG
590  if (A->getDomainMap()->isCompatible(*(X.getMap())) == false) {
591  std::ostringstream ss;
592  ss << "Level " << startLevel << ": level A's domain map is not compatible with X";
593  throw Exceptions::Incompatible(ss.str());
594  }
595 
596  if (A->getRangeMap()->isCompatible(*(B.getMap())) == false) {
597  std::ostringstream ss;
598  ss << "Level " << startLevel << ": level A's range map is not compatible with B";
599  throw Exceptions::Incompatible(ss.str());
600  }
601 #endif
602 
603  if (startLevel == as<LO>(Levels_.size())-1) {
604  // On the coarsest level, we do either smoothing (if defined) or a direct solve.
605  RCP<TimeMonitor> CLevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : coarse" + levelSuffix, Timings0));
606 
607  bool emptySolve = true;
608 
609  // NOTE: we need to check using IsAvailable before Get here to avoid building default smoother
610  if (Fine->IsAvailable("PreSmoother")) {
611  RCP<SmootherBase> preSmoo = Fine->Get< RCP<SmootherBase> >("PreSmoother");
612  preSmoo->Apply(X, B, zeroGuess);
613  zeroGuess = false;
614  emptySolve = false;
615  }
616  if (Fine->IsAvailable("PostSmoother")) {
617  RCP<SmootherBase> postSmoo = Fine->Get< RCP<SmootherBase> >("PostSmoother");
618  postSmoo->Apply(X, B, zeroGuess);
619  emptySolve = false;
620  }
621  if (emptySolve == true)
622  GetOStream(Warnings0) << "No coarse grid solver" << std::endl;
623 
624  } else {
625  // On intermediate levels, we do cycles
626  RCP<Level> Coarse = Levels_[startLevel+1];
627 
628  {
629  // ============== PRESMOOTHING ==============
630  RCP<TimeMonitor> STime = rcp(new TimeMonitor(*this, prefix + "Solve : smoothing (total)" , Timings0));
631  RCP<TimeMonitor> SLevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : smoothing" + levelSuffix, Timings0));
632 
633  if (Fine->IsAvailable("PreSmoother")) {
634  RCP<SmootherBase> preSmoo = Fine->Get< RCP<SmootherBase> >("PreSmoother");
635  preSmoo->Apply(X, B, zeroGuess);
636  } else {
637  GetOStream(Warnings1) << "Level " << startLevel << ": No PreSmoother!" << std::endl;
638  }
639  }
640 
641  RCP<MultiVector> residual;
642  {
643  RCP<TimeMonitor> ATime = rcp(new TimeMonitor(*this, prefix + "Solve : residual calculation (total)" , Timings0));
644  RCP<TimeMonitor> ALevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : residual calculation" + levelSuffix, Timings0));
645  residual = Utilities::Residual(*A, X, B);
646  }
647 
648  RCP<Operator> P = Coarse->Get< RCP<Operator> >("P");
649  RCP<MultiVector> coarseRhs, coarseX;
650  const bool initializeWithZeros = true;
651  {
652  // ============== RESTRICTION ==============
653  RCP<TimeMonitor> RTime = rcp(new TimeMonitor(*this, prefix + "Solve : restriction (total)" , Timings0));
654  RCP<TimeMonitor> RLevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : restriction" + levelSuffix, Timings0));
655 
656  if (implicitTranspose_) {
657  coarseRhs = MultiVectorFactory::Build(P->getDomainMap(), X.getNumVectors(), !initializeWithZeros);
658  P->apply(*residual, *coarseRhs, Teuchos::TRANS, one, zero);
659 
660  } else {
661  RCP<Operator> R = Coarse->Get< RCP<Operator> >("R");
662  coarseRhs = MultiVectorFactory::Build(R->getRangeMap(), X.getNumVectors(), !initializeWithZeros);
663  R->apply(*residual, *coarseRhs, Teuchos::NO_TRANS, one, zero);
664  }
665  }
666 
667  RCP<const Import> importer;
668  if (Coarse->IsAvailable("Importer"))
669  importer = Coarse->Get< RCP<const Import> >("Importer");
670 
671  if (doPRrebalance_ || importer.is_null()) {
672  coarseX = MultiVectorFactory::Build(coarseRhs->getMap(), X.getNumVectors(), initializeWithZeros);
673 
674  } else {
675  RCP<TimeMonitor> ITime = rcp(new TimeMonitor(*this, prefix + "Solve : import (total)" , Timings0));
676  RCP<TimeMonitor> ILevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : import" + levelSuffix, Timings0));
677 
678  // Import: range map of R --> domain map of rebalanced Ac (before subcomm replacement)
679  RCP<MultiVector> coarseTmp = MultiVectorFactory::Build(importer->getTargetMap(), coarseRhs->getNumVectors(), false);
680  coarseTmp->doImport(*coarseRhs, *importer, Xpetra::INSERT);
681  coarseRhs.swap(coarseTmp);
682 
683  coarseX = MultiVectorFactory::Build(importer->getTargetMap(), X.getNumVectors(), initializeWithZeros);
684  }
685 
686  RCP<Operator> Ac = Coarse->Get< RCP<Operator> >("A");
687  if (!Ac.is_null()) {
688  RCP<const Map> origXMap = coarseX->getMap();
689 
690  // Replace maps with maps with a subcommunicator
691  coarseRhs->replaceMap(Ac->getRangeMap());
692  coarseX ->replaceMap(Ac->getDomainMap());
693 
694  {
695  iterateLevelTime = Teuchos::null; // stop timing this level
696 
697  Iterate(*coarseRhs, *coarseX, 1, true, startLevel+1);
698  // ^^ zero initial guess
699  if (Cycle_ == WCYCLE)
700  Iterate(*coarseRhs, *coarseX, 1, false, startLevel+1);
701  // ^^ nonzero initial guess
702 
703  iterateLevelTime = rcp(new TimeMonitor(*this, iterateLevelTimeLabel)); // restart timing this level
704  }
705  coarseX->replaceMap(origXMap);
706  }
707 
708  if (!doPRrebalance_ && !importer.is_null()) {
709  RCP<TimeMonitor> ITime = rcp(new TimeMonitor(*this, prefix + "Solve : export (total)" , Timings0));
710  RCP<TimeMonitor> ILevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : export" + levelSuffix, Timings0));
711 
712  // Import: range map of rebalanced Ac (before subcomm replacement) --> domain map of P
713  RCP<MultiVector> coarseTmp = MultiVectorFactory::Build(importer->getSourceMap(), coarseX->getNumVectors(), false);
714  coarseTmp->doExport(*coarseX, *importer, Xpetra::INSERT);
715  coarseX.swap(coarseTmp);
716  }
717 
718  // Update X += P * coarseX
719  // Note that due to what may be round-off error accumulation, use of the fused kernel
720  // P->apply(*coarseX, X, Teuchos::NO_TRANS, one, one);
721  // can in some cases result in slightly higher iteration counts.
722  RCP<MultiVector> correction = MultiVectorFactory::Build(P->getRangeMap(), X.getNumVectors(),false);
723  {
724  // ============== PROLONGATION ==============
725  RCP<TimeMonitor> PTime = rcp(new TimeMonitor(*this, prefix + "Solve : prolongation (total)" , Timings0));
726  RCP<TimeMonitor> PLevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : prolongation" + levelSuffix, Timings0));
727  P->apply(*coarseX, *correction, Teuchos::NO_TRANS, one, zero);
728  }
729  X.update(one, *correction, one);
730 
731  {
732  // ============== POSTSMOOTHING ==============
733  RCP<TimeMonitor> STime = rcp(new TimeMonitor(*this, prefix + "Solve : smoothing (total)" , Timings0));
734  RCP<TimeMonitor> SLevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : smoothing" + levelSuffix, Timings0));
735 
736  if (Fine->IsAvailable("PostSmoother")) {
737  RCP<SmootherBase> postSmoo = Fine->Get< RCP<SmootherBase> >("PostSmoother");
738  postSmoo->Apply(X, B, false);
739 
740  } else {
741  GetOStream(Warnings1) << "Level " << startLevel << ": No PostSmoother!" << std::endl;
742  }
743  }
744  }
745  zeroGuess = false;
746 
747  if (startLevel == 0 && !isPreconditioner_ &&
748  (IsPrint(Statistics1) || tol > 0)) {
749  // We calculate the residual only if we want to print it out, or if we
750  // want to stop once we achive the tolerance
751  Teuchos::Array<MagnitudeType> rn;
752  rn = Utilities::ResidualNorm(*A, X, B);
753 
754  prevNorm = curNorm;
755  curNorm = rn[0];
756  rate_ = as<MagnitudeType>(curNorm / prevNorm);
757 
758  if (tol > 0) {
759  bool passed = true;
760  for (LO k = 0; k < rn.size(); k++)
761  if (rn[k] >= tol)
762  passed = false;
763 
764  if (passed)
765  return Converged;
766  }
767 
768  if (IsPrint(Statistics1))
769  GetOStream(Statistics1) << "iter: "
770  << std::setiosflags(std::ios::left)
771  << std::setprecision(3) << i
772  << " residual = "
773  << std::setprecision(10) << rn
774  << std::endl;
775  }
776  }
777  return (tol > 0 ? Unconverged : Undefined);
778  }
779 
780  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
781  void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(const LO& start, const LO& end) {
782  LO startLevel = (start != -1 ? start : 0);
783  LO endLevel = (end != -1 ? end : Levels_.size()-1);
784 
785  TEUCHOS_TEST_FOR_EXCEPTION(startLevel > endLevel, Exceptions::RuntimeError,
786  "MueLu::Hierarchy::Write : startLevel must be <= endLevel");
787 
788  TEUCHOS_TEST_FOR_EXCEPTION(startLevel < 0 || endLevel >= Levels_.size(), Exceptions::RuntimeError,
789  "MueLu::Hierarchy::Write bad start or end level");
790 
791  for (LO i = startLevel; i < endLevel + 1; i++) {
792  RCP<Matrix> A = rcp_dynamic_cast<Matrix>(Levels_[i]-> template Get< RCP< Operator> >("A")), P, R;
793  if (i > 0) {
794  P = rcp_dynamic_cast<Matrix>(Levels_[i]-> template Get< RCP< Operator> >("P"));
795  if (!implicitTranspose_)
796  R = rcp_dynamic_cast<Matrix>(Levels_[i]-> template Get< RCP< Operator> >("R"));
797  }
798 
799  if (!A.is_null()) Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("A_" + toString(i) + ".m", *A);
800  if (!P.is_null()) Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("P_" + toString(i) + ".m", *P);
801  if (!R.is_null()) Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("R_" + toString(i) + ".m", *R);
802  }
803  }
804 
805  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
806  void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Keep(const std::string & ename, const FactoryBase* factory) {
807  for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
808  (*it)->Keep(ename, factory);
809  }
810 
811  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
812  void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Delete(const std::string& ename, const FactoryBase* factory) {
813  for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
814  (*it)->Delete(ename, factory);
815  }
816 
817  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
818  void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::AddKeepFlag(const std::string & ename, const FactoryBase* factory, KeepType keep) {
819  for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
820  (*it)->AddKeepFlag(ename, factory, keep);
821  }
822 
823  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
824  void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::RemoveKeepFlag(const std::string & ename, const FactoryBase* factory, KeepType keep) {
825  for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
826  (*it)->RemoveKeepFlag(ename, factory, keep);
827  }
828 
829  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
831  std::ostringstream out;
832  out << BaseClass::description();
833  out << "{#levels = " << GetGlobalNumLevels() << ", complexity = " << GetOperatorComplexity() << "}";
834  return out.str();
835  }
836 
837  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
838  void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::describe(Teuchos::FancyOStream& out, const Teuchos::EVerbosityLevel tVerbLevel) const {
839  describe(out, toMueLuVerbLevel(tVerbLevel));
840  }
841 
842  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
843  void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::describe(Teuchos::FancyOStream& out, const VerbLevel verbLevel) const {
844  RCP<Operator> A0 = Levels_[0]->template Get<RCP<Operator> >("A");
845  RCP<const Teuchos::Comm<int> > comm = A0->getDomainMap()->getComm();
846 
847  int numLevels = GetNumLevels();
848  RCP<Operator> Ac = Levels_[numLevels-1]->template Get<RCP<Operator> >("A");
849  if (Ac.is_null()) {
850  // It may happen that we do repartition on the last level, but the matrix
851  // is small enough to satisfy "max coarse size" requirement. Then, even
852  // though we have the level, the matrix would be null on all but one processors
853  numLevels--;
854  }
855  int root = comm->getRank();
856 
857 #ifdef HAVE_MPI
858  int smartData = numLevels*comm->getSize() + comm->getRank(), maxSmartData;
859  reduceAll(*comm, Teuchos::REDUCE_MAX, smartData, Teuchos::ptr(&maxSmartData));
860  root = maxSmartData % comm->getSize();
861 #endif
862 
863  std::string outstr;
864  if (comm->getRank() == root && verbLevel & (Statistics0 | Test)) {
865  std::vector<Xpetra::global_size_t> nnzPerLevel;
866  std::vector<Xpetra::global_size_t> rowsPerLevel;
867  std::vector<int> numProcsPerLevel;
868  bool aborted = false;
869  for (int i = 0; i < numLevels; i++) {
870  TEUCHOS_TEST_FOR_EXCEPTION(!(Levels_[i]->IsAvailable("A")) , Exceptions::RuntimeError,
871  "Operator A is not available on level " << i);
872 
873  RCP<Operator> A = Levels_[i]->template Get<RCP<Operator> >("A");
874  TEUCHOS_TEST_FOR_EXCEPTION(A.is_null(), Exceptions::RuntimeError,
875  "Operator A on level " << i << " is null.");
876 
877  RCP<Matrix> Am = rcp_dynamic_cast<Matrix>(A);
878  if (Am.is_null()) {
879  GetOStream(Warnings0) << "Some level operators are not matrices, statistics calculation aborted" << std::endl;
880  aborted = true;
881  break;
882  }
883 
884  Xpetra::global_size_t nnz = Am->getGlobalNumEntries();
885  nnzPerLevel .push_back(nnz);
886  rowsPerLevel .push_back(Am->getGlobalNumRows());
887  numProcsPerLevel.push_back(Am->getRowMap()->getComm()->getSize());
888  }
889 
890  if (!aborted) {
891  std::ostringstream oss;
892  oss << "\n--------------------------------------------------------------------------------\n" <<
893  "--- Multigrid Summary ---\n"
894  "--------------------------------------------------------------------------------" << std::endl;
895  oss << "Number of levels = " << numLevels << std::endl;
896  oss << "Operator complexity = " << std::setprecision(2) << std::setiosflags(std::ios::fixed)
897  << GetOperatorComplexity() << std::endl;
898  oss << std::endl;
899 
900  Xpetra::global_size_t tt = rowsPerLevel[0];
901  int rowspacer = 2; while (tt != 0) { tt /= 10; rowspacer++; }
902  tt = nnzPerLevel[0];
903  int nnzspacer = 2; while (tt != 0) { tt /= 10; nnzspacer++; }
904  tt = numProcsPerLevel[0];
905  int npspacer = 2; while (tt != 0) { tt /= 10; npspacer++; }
906  oss << "level " << std::setw(rowspacer) << " rows " << std::setw(nnzspacer) << " nnz " << " nnz/row" << std::setw(npspacer) << " c ratio" << " procs" << std::endl;
907  for (size_t i = 0; i < nnzPerLevel.size(); ++i) {
908  oss << " " << i << " ";
909  oss << std::setw(rowspacer) << rowsPerLevel[i];
910  oss << std::setw(nnzspacer) << nnzPerLevel[i];
911  oss << std::setprecision(2) << std::setiosflags(std::ios::fixed);
912  oss << std::setw(9) << as<double>(nnzPerLevel[i]) / rowsPerLevel[i];
913  if (i) oss << std::setw(9) << as<double>(rowsPerLevel[i-1])/rowsPerLevel[i];
914  else oss << std::setw(9) << " ";
915  oss << " " << std::setw(npspacer) << numProcsPerLevel[i] << std::endl;
916  }
917  oss << std::endl;
918  for (int i = 0; i < GetNumLevels(); ++i) {
919  RCP<SmootherBase> preSmoo, postSmoo;
920  if (Levels_[i]->IsAvailable("PreSmoother"))
921  preSmoo = Levels_[i]->template Get< RCP<SmootherBase> >("PreSmoother");
922  if (Levels_[i]->IsAvailable("PostSmoother"))
923  postSmoo = Levels_[i]->template Get< RCP<SmootherBase> >("PostSmoother");
924 
925  if (preSmoo != null && preSmoo == postSmoo)
926  oss << "Smoother (level " << i << ") both : " << preSmoo->description() << std::endl;
927  else {
928  oss << "Smoother (level " << i << ") pre : "
929  << (preSmoo != null ? preSmoo->description() : "no smoother") << std::endl;
930  oss << "Smoother (level " << i << ") post : "
931  << (postSmoo != null ? postSmoo->description() : "no smoother") << std::endl;
932  }
933 
934  oss << std::endl;
935  }
936 
937  outstr = oss.str();
938  }
939  }
940 
941 #ifdef HAVE_MPI
942  RCP<const Teuchos::MpiComm<int> > mpiComm = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(comm);
943  MPI_Comm rawComm = (*mpiComm->getRawMpiComm())();
944 
945  int strLength = outstr.size();
946  MPI_Bcast(&strLength, 1, MPI_INT, root, rawComm);
947  if (comm->getRank() != root)
948  outstr.resize(strLength);
949  MPI_Bcast(&outstr[0], strLength, MPI_CHAR, root, rawComm);
950 #endif
951 
952  out << outstr;
953  }
954 
955  // NOTE: at some point this should be replaced by a friend operator <<
956  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
957  void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::print(std::ostream& out, const VerbLevel verbLevel) const {
958  Teuchos::OSTab tab2(out);
959  for (int i = 0; i < GetNumLevels(); ++i)
960  Levels_[i]->print(out, verbLevel);
961  }
962 
963  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
965  isPreconditioner_ = flag;
966  }
967 
968  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
970  if (GetProcRankVerbose() != 0)
971  return;
972 #if defined(HAVE_MUELU_BOOST) && defined(HAVE_MUELU_BOOST_FOR_REAL) && defined(BOOST_VERSION) && (BOOST_VERSION >= 104400)
973  BoostGraph graph;
974 
975  BoostProperties dp;
976  dp.property("label", boost::get(boost::vertex_name, graph));
977  dp.property("id", boost::get(boost::vertex_index, graph));
978  dp.property("label", boost::get(boost::edge_name, graph));
979  dp.property("color", boost::get(boost::edge_color, graph));
980 
981  // create local maps
982  std::map<const FactoryBase*, BoostVertex> vindices;
983  typedef std::map<std::pair<BoostVertex,BoostVertex>, std::string> emap; emap edges;
984 
985  for (int i = dumpLevel_; i <= dumpLevel_+1 && i < GetNumLevels(); i++) {
986  edges.clear();
987  Levels_[i]->UpdateGraph(vindices, edges, dp, graph);
988 
989  for (emap::const_iterator eit = edges.begin(); eit != edges.end(); eit++) {
990  std::pair<BoostEdge, bool> boost_edge = boost::add_edge(eit->first.first, eit->first.second, graph);
991  boost::put("label", dp, boost_edge.first, eit->second);
992  if (i == dumpLevel_)
993  boost::put("color", dp, boost_edge.first, std::string("red"));
994  else
995  boost::put("color", dp, boost_edge.first, std::string("blue"));
996  }
997  }
998 
999  // add legend
1000  std::ostringstream legend;
1001  legend << "< <TABLE BORDER=\"0\" CELLBORDER=\"1\" CELLSPACING=\"0\" CELLPADDING=\"4\"> \
1002  <TR><TD COLSPAN=\"2\">Legend</TD></TR> \
1003  <TR><TD><FONT color=\"red\">Level " << dumpLevel_ << "</FONT></TD><TD><FONT color=\"blue\">Level " << dumpLevel_+1 << "</FONT></TD></TR> \
1004  </TABLE> >";
1005  BoostVertex boost_vertex = boost::add_vertex(graph);
1006  boost::put("label", dp, boost_vertex, legend.str());
1007 
1008  std::ofstream out(dumpFile_.c_str());
1009  boost::write_graphviz_dp(out, graph, dp, std::string("id"));
1010 #else
1011  GetOStream(Errors) << "Dependency graph output requires boost and MueLu_ENABLE_Boost_for_real" << std::endl;
1012 #endif
1013  }
1014 
1015  // Enforce that coordinate vector's map is consistent with that of A
1016  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1018  RCP<Operator> Ao = level.Get<RCP<Operator> >("A");
1019  RCP<Matrix> A = rcp_dynamic_cast<Matrix>(Ao);
1020  if (A.is_null()) {
1021  GetOStream(Warnings0) << "Hierarchy::ReplaceCoordinateMap: operator is not a matrix, skipping..." << std::endl;
1022  return;
1023  }
1024 
1025  typedef Xpetra::MultiVector<double,LO,GO,NO> xdMV;
1026 
1027  RCP<xdMV> coords = level.Get<RCP<xdMV> >("Coordinates");
1028 
1029  if (A->getRowMap()->isSameAs(*(coords->getMap()))) {
1030  GetOStream(Warnings0) << "Hierarchy::ReplaceCoordinateMap: matrix and coordinates maps are same, skipping..." << std::endl;
1031  return;
1032  }
1033 
1034  if (A->IsView("stridedMaps") && rcp_dynamic_cast<const StridedMap>(A->getRowMap("stridedMaps")) != Teuchos::null) {
1035  RCP<const StridedMap> stridedRowMap = rcp_dynamic_cast<const StridedMap>(A->getRowMap("stridedMaps"));
1036 
1037  // It is better to through an exceptions if maps may be inconsistent, than to ignore it and experience unfathomable breakdowns
1038  TEUCHOS_TEST_FOR_EXCEPTION(stridedRowMap->getStridedBlockId() != -1 || stridedRowMap->getOffset() != 0,
1039  Exceptions::RuntimeError, "Hierarchy::ReplaceCoordinateMap: nontrivial maps (block id = " << stridedRowMap->getStridedBlockId()
1040  << ", offset = " << stridedRowMap->getOffset() << ")");
1041  }
1042 
1043  GetOStream(Runtime1) << "Replacing coordinate map" << std::endl;
1044 
1045  size_t blkSize = A->GetFixedBlockSize();
1046 
1047  RCP<const Map> nodeMap = A->getRowMap();
1048  if (blkSize > 1) {
1049  // Create a nodal map, as coordinates have not been expanded to a DOF map yet.
1050  RCP<const Map> dofMap = A->getRowMap();
1051  GO indexBase = dofMap->getIndexBase();
1052  size_t numLocalDOFs = dofMap->getNodeNumElements();
1053  TEUCHOS_TEST_FOR_EXCEPTION(numLocalDOFs % blkSize, Exceptions::RuntimeError,
1054  "Hierarchy::ReplaceCoordinateMap: block size (" << blkSize << ") is incompatible with the number of local dofs in a row map (" << numLocalDOFs);
1055  ArrayView<const GO> GIDs = dofMap->getNodeElementList();
1056 
1057  Array<GO> nodeGIDs(numLocalDOFs/blkSize);
1058  for (size_t i = 0; i < numLocalDOFs; i += blkSize)
1059  nodeGIDs[i/blkSize] = (GIDs[i] - indexBase)/blkSize + indexBase;
1060 
1061  Xpetra::global_size_t INVALID = Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid();
1062  nodeMap = MapFactory::Build(dofMap->lib(), INVALID, nodeGIDs(), indexBase, dofMap->getComm());
1063  }
1064 
1065  Array<ArrayView<const double> > coordDataView;
1066  std::vector<ArrayRCP<const double> > coordData;
1067  for (size_t i = 0; i < coords->getNumVectors(); i++) {
1068  coordData.push_back(coords->getData(i));
1069  coordDataView.push_back(coordData[i]());
1070  }
1071 
1072  RCP<xdMV> newCoords = Xpetra::MultiVectorFactory<double,LO,GO,NO>::Build(nodeMap, coordDataView(), coords->getNumVectors());
1073  level.Set("Coordinates", newCoords);
1074  }
1075 
1076 } //namespace MueLu
1077 
1078 #endif // MUELU_HIERARCHY_DEF_HPP
Important warning messages (one line)
void IsPreconditioner(const bool flag)
RCP< Level > & GetPreviousLevel()
Previous level.
static void SetRandomSeed(const Teuchos::Comm< int > &comm)
High level timing information (use Teuchos::TimeMonitor::summarize() to print)
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->Get< RCP<Matrix> >("A", factory) if factory == NULL => use default factory.
RCP< Level > & GetLevel(const int levelID=0)
Retrieve a certain level from hierarchy.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Hierarchy()
Default constructor.
void Release(const FactoryBase &factory)
Decrement the storage counter for all the inputs of a factory.
void AddLevel(const RCP< Level > &level)
Add a level at the end of the hierarchy.
short KeepType
Print more statistics.
void AddNewLevel()
Add a new level at the end of the hierarchy.
VerbLevel toMueLuVerbLevel(const Teuchos::EVerbosityLevel verbLevel)
One-liner description of what is happening.
static Xpetra::global_size_t GetDefaultMaxCoarseSize()
void DumpCurrentGraph() const
void Clear(int startLevel=0)
Clear impermanent data from previous setup.
Namespace for MueLu classes and methods.
void Write(const LO &start=-1, const LO &end=-1)
Print matrices in the multigrid hierarchy to file.
Exception throws to report incompatible objects (like maps).
Integrates Teuchos::TimeMonitor with MueLu verbosity system.
Print skeleton for the run, i.e. factory calls and used parameters.
MagnitudeType rate_
Convergece rate.
virtual std::string ShortClassName() const
Return the class name of the object, without template parameters and without namespace.
void ReplaceCoordinateMap(Level &level)
void AddKeepFlag(const std::string &ename, const FactoryBase *factory=NoFactory::get(), KeepType keep=MueLu::Keep)
Call Level::AddKeepFlag for each level of the Hierarchy.
void setlib(Xpetra::UnderlyingLib lib2)
int SetProcRankVerbose(int procRank) const
Set proc rank used for printing.
Additional warnings.
void print(std::ostream &out=std::cout, const VerbLevel verbLevel=(MueLu::Parameters|MueLu::Statistics0)) const
VerbLevel GetVerbLevel() const
Get the verbosity level.
Base class for factories (e.g., R, P, and A_coarse).
Print statistics that do not involve significant additional computation.
Xpetra::UnderlyingLib lib_
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)
static CycleType GetDefaultCycle()
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
std::string description() const
Return a simple one-line description of this object.
Class that provides default factories within Needs class.
void RemoveKeepFlag(const std::string &ename, const FactoryBase *factory, KeepType keep=MueLu::All)
Call Level::RemoveKeepFlag for each level of the Hierarchy.
void CheckLevel(Level &level, int levelID)
Helper function.
Xpetra::UnderlyingLib lib()
static bool GetDefaultImplicitTranspose()
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
int GetGlobalNumLevels() const
Teuchos::FancyOStream & GetOStream(MsgType type, int thisProcRankOnly=0) const
Get an output stream for outputting the input message type.
bool IsPrint(MsgType type, int thisProcRankOnly=-1) const
Find out whether we need to print out information for a specific message type.
Data struct for defining stopping criteria of multigrid iteration.
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)
STS::magnitudeType MagnitudeType
void Delete(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Call Level::Delete(ename, factory) for each level of the Hierarchy.
Timer to be used in non-factories.
Array< RCP< Level > > Levels_
Container for Level objects.
static bool GetDefaultPRrebalance()
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
Exception throws to report errors in the internal logical of the program.
An exception safe way to call the method &#39;Level::SetFactoryManager()&#39;.
void SetComm(RCP< const Teuchos::Comm< int > > const &comm)
void describe(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the Hierarchy with some verbosity level to a FancyOStream object.
ReturnType Iterate(const MultiVector &B, MultiVector &X, ConvData conv=ConvData(), bool InitialGuessIsZero=false, LO startLevel=0)
Apply the multigrid preconditioner.
Description of what is happening (more verbose)
void Keep(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Call Level::Keep(ename, factory) for each level of the Hierarchy.
bool isDumpingEnabled_
Graph dumping.
double GetOperatorComplexity() const
bool Setup(int coarseLevelID, const RCP< const FactoryManagerBase > fineLevelManager, const RCP< const FactoryManagerBase > coarseLevelManager, const RCP< const FactoryManagerBase > nextLevelManager=Teuchos::null)
Multi-level setup phase: build a new level of the hierarchy.
RCP< const Teuchos::Comm< int > > GetComm() const
virtual std::string description() const
Return a simple one-line description of this object.
Array< RCP< const FactoryManagerBase > > levelManagers_
Xpetra::global_size_t maxCoarseSize_
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need&#39;s value has been saved.
int GetProcRankVerbose() const
Get proc rank used for printing. Do not use this information for any other purpose.
void Request(const FactoryBase &factory)
Increment the storage counter for all the inputs of a factory.