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