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