MueLu  Version of the Day
MueLu_Hierarchy_decl.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_HIERARCHY_DECL_HPP
47 #define MUELU_HIERARCHY_DECL_HPP
48 
49 #include <Teuchos_ParameterList.hpp>
50 #include <Teuchos_Ptr.hpp>
51 
52 #include <Xpetra_ConfigDefs.hpp> // global_size_t
53 #include <Xpetra_Matrix_fwd.hpp>
54 #include <Xpetra_MultiVector_fwd.hpp>
55 #include <Xpetra_MultiVectorFactory_fwd.hpp>
56 #include <Xpetra_Operator_fwd.hpp>
57 
58 #include <Xpetra_Cloner.hpp>
59 #include <MueLu_SmootherCloner.hpp>
60 #include "MueLu_ConfigDefs.hpp"
61 #include "MueLu_BaseClass.hpp"
62 #include "MueLu_Hierarchy_fwd.hpp"
63 
64 #include "MueLu_Types.hpp"
65 
67 #include "MueLu_FactoryManager.hpp" // no fwd declaration because constructor of FactoryManager is used as a default parameter of Setup()
69 #include "MueLu_KeepType.hpp"
70 #include "MueLu_Level_fwd.hpp"
71 #include "MueLu_MasterList.hpp"
72 #include "MueLu_NoFactory.hpp"
73 #include "MueLu_PerfUtils_fwd.hpp"
74 #include "MueLu_PFactory_fwd.hpp"
75 #include "MueLu_RFactory_fwd.hpp"
80 #include "MueLu_Utilities_fwd.hpp"
81 
82 namespace MueLu {
83 
84  enum ReturnType {
88  };
89 
100  template <class Scalar = Xpetra::Operator<>::scalar_type,
101  class LocalOrdinal = typename Xpetra::Operator<Scalar>::local_ordinal_type,
102  class GlobalOrdinal = typename Xpetra::Operator<Scalar, LocalOrdinal>::global_ordinal_type,
103  class Node = typename Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal>::node_type>
104  class Hierarchy : public BaseClass {
105 #undef MUELU_HIERARCHY_SHORT
106 #include "MueLu_UseShortNames.hpp"
107 
108  typedef Teuchos::ScalarTraits<SC> STS;
109  typedef typename STS::magnitudeType MagnitudeType;
110 
112  struct ConvData {
113  ConvData() : maxIts_(1), tol_(-STS::magnitude(STS::one())) { }
114  ConvData(LO maxIts) : maxIts_(maxIts), tol_(-STS::magnitude(STS::one())) { }
115  ConvData(MagnitudeType tol) : maxIts_(10000), tol_(tol) { }
116  ConvData(std::pair<LO,MagnitudeType> p) : maxIts_(p.first), tol_(p.second) { }
117 
120  };
121 
122  public:
123 
125 
126 
128  Hierarchy();
130  Hierarchy(const std::string& label);
131 
133  Hierarchy(const RCP<Matrix> & A);
134 
136  Hierarchy(const RCP<Matrix> & A, const std::string& label);
137 
139  virtual ~Hierarchy() { }
140 
142 
144 
145 
147  static CycleType GetDefaultCycle() { return MasterList::getDefault<std::string>("cycle type") == "V" ? VCYCLE : WCYCLE; }
148  static bool GetDefaultImplicitTranspose() { return MasterList::getDefault<bool>("transpose: use implicit"); }
149  static Xpetra::global_size_t GetDefaultMaxCoarseSize() { return MasterList::getDefault<int>("coarse: max size"); }
150  static int GetDefaultMaxLevels() { return MasterList::getDefault<int>("max levels"); }
151  static bool GetDefaultPRrebalance() { return MasterList::getDefault<bool>("repartition: rebalance P and R"); }
152 
153  Xpetra::global_size_t GetMaxCoarseSize() const { return maxCoarseSize_; }
154  bool GetImplicitTranspose() const { return implicitTranspose_; }
155 
156  void SetMaxCoarseSize(Xpetra::global_size_t maxCoarseSize) { maxCoarseSize_ = maxCoarseSize; }
157  void SetPRrebalance(bool doPRrebalance) { doPRrebalance_ = doPRrebalance; }
158  void SetImplicitTranspose(const bool& implicit) { implicitTranspose_ = implicit; }
159 
161 
163 
164  template<class S2, class LO2, class GO2, class N2>
165  friend class Hierarchy;
166 
167  private:
168  int LastLevelID() const { return Levels_.size() - 1; }
169  void DumpCurrentGraph() const;
170 
171  public:
172 
174  void AddLevel(const RCP<Level> & level);
175 
177  void AddNewLevel();
178 
180  RCP<Level> & GetLevel(const int levelID = 0);
181 
182  int GetNumLevels() const;
183  int GetGlobalNumLevels() const;
184 
185  MagnitudeType GetRate() const { return rate_; }
186 
187  // This function is global
188  double GetOperatorComplexity() const;
189 
190  // This function is global
191  double GetSmootherComplexity() const;
192 
194  void CheckLevel(Level& level, int levelID);
195 
197 
234  bool Setup(int coarseLevelID, const RCP<const FactoryManagerBase> fineLevelManager /* = Teuchos::null */, const RCP<const FactoryManagerBase> coarseLevelManager,
235  const RCP<const FactoryManagerBase> nextLevelManager = Teuchos::null);
236 
238  void Setup(const FactoryManagerBase& manager = FactoryManager(), int startLevel = 0, int numDesiredLevels = GetDefaultMaxLevels());
239 
240  void SetupRe();
241 
243  void Clear(int startLevel = 0);
244  void ExpertClear();
245 
247  CycleType GetCycle() const { return Cycle_; }
248 
250  void SetCycle(CycleType Cycle) { Cycle_ = Cycle; }
251 
253  void SetProlongatorScalingFactor(double scalingFactor) { scalingFactor_ = scalingFactor; }
254 
267  ReturnType Iterate(const MultiVector& B, MultiVector& X, ConvData conv = ConvData(),
268  bool InitialGuessIsZero = false, LO startLevel = 0);
269 
279  void Write(const LO &start=-1, const LO &end=-1, const std::string &suffix="");
280 
282 
284 
285 
287  void Keep(const std::string & ename, const FactoryBase* factory = NoFactory::get());
288 
290  void Delete(const std::string& ename, const FactoryBase* factory = NoFactory::get());
291 
293  void AddKeepFlag(const std::string & ename, const FactoryBase* factory = NoFactory::get(), KeepType keep = MueLu::Keep);
294 
296  void RemoveKeepFlag(const std::string & ename, const FactoryBase* factory, KeepType keep = MueLu::All);
297 
299 
301 
302 
304  std::string description() const;
305 
311  void describe(Teuchos::FancyOStream& out, const VerbLevel verbLevel = Default) const;
312  void describe(Teuchos::FancyOStream& out, const Teuchos::EVerbosityLevel verbLevel = Teuchos::VERB_HIGH) const;
313 
314  // Hierarchy::print is local hierarchy function, thus the statistics can be different from global ones
315  void print(std::ostream& out = std::cout, const VerbLevel verbLevel = (MueLu::Parameters | MueLu::Statistics0)) const;
316 
321  void IsPreconditioner(const bool flag);
322 
324 
325  void EnableGraphDumping(const std::string& filename, int levelID = 1) {
326  isDumpingEnabled_ = true;
327  dumpLevel_ = levelID;
328  dumpFile_ = filename;
329  }
330 
331  template<class Node2>
332  Teuchos::RCP< Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node2> >
333  clone (const RCP<Node2> &node2) const;
334 
335  void setlib(Xpetra::UnderlyingLib inlib) { lib_ = inlib; }
336  Xpetra::UnderlyingLib lib() { return lib_; }
337 
338  // force recreation of cached description_ next time description() is called:
340  description_ = "";
341  }
342 
343  void AllocateLevelMultiVectors(int numvecs);
345 
346  protected:
347  const RCP<const FactoryManagerBase>& GetLevelManager(const int levelID) const {
348  return levelManagers_[levelID];
349  }
350 
351  private:
353  Hierarchy(const Hierarchy &h);
354 
356  Array<RCP<Level> > Levels_;
357 
358  // We replace coordinates GIDs to make them consistent with matrix GIDs,
359  // even if user does not do that. Ideally, though, we should completely
360  // remove any notion of coordinate GIDs, and deal only with LIDs, assuming
361  // that they are consistent with matrix block IDs
362  void ReplaceCoordinateMap(Level& level);
363 
364  // Minimum size of a matrix on any level. If we fall below that, we stop
365  // the coarsening
366  Xpetra::global_size_t maxCoarseSize_;
367 
368  // Potential speed up of the setup by skipping R construction, and using
369  // transpose matrix-matrix product for RAP
371 
372  // Potential speed up of the setup by skipping rebalancing of P and R, and
373  // doing extra import during solve
375 
376  // Hierarchy may be used in a standalone mode, or as a preconditioner
378 
379  // V- or W-cycle
381 
382  // Scaling factor to be applied to coarse grid correction.
384 
385  // Epetra/Tpetra mode
386  Xpetra::UnderlyingLib lib_;
387 
388  // cache description to avoid recreating in each call to description() - use ResetDescription() to force recreation in Setup, SetupRe, etc.
389  mutable std::string description_ = ""; // mutable so that we can lazily initialize in description(), which is declared const
390 
392  // If enabled, we dump the graph on a specified level into a specified file
395  std::string dumpFile_;
396 
399 
400  // Level managers used during the Setup
401  Array<RCP<const FactoryManagerBase> > levelManagers_;
402 
403  // Caching (Multi)Vectors used in Hierarchy::Iterate()
406 
407 
408  }; //class Hierarchy
409 
410  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
411  template<typename Node2>
412  Teuchos::RCP<Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node2> >
414  clone (const Teuchos::RCP<Node2> &node2) const {
416  typedef Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node2> CloneOperator;
418 
419  Teuchos::RCP<New_H_Type> new_h = Teuchos::rcp(new New_H_Type());
420  new_h->Levels_.resize(this->GetNumLevels());
421  new_h->maxCoarseSize_ = maxCoarseSize_;
422  new_h->implicitTranspose_ = implicitTranspose_;
423  new_h->isPreconditioner_ = isPreconditioner_;
424  new_h->isDumpingEnabled_ = isDumpingEnabled_;
425  new_h->dumpLevel_ = dumpLevel_;
426  new_h->dumpFile_ = dumpFile_;
427 
428  RCP<SmootherBase> Pre, Post;
429  RCP<CloneSmoother> clonePre, clonePost;
430  RCP<CloneOperator> cloneA, cloneR, cloneP;
431  RCP<Operator> A, R, P;
432  for (int i = 0; i < GetNumLevels(); i++) {
433  RCP<Level> level = this->Levels_[i];
434  RCP<Level> clonelevel = rcp(new Level());
435 
436  if (level->IsAvailable("A")) {
437  A = level->template Get<RCP<Operator> >("A");
438  cloneA = Xpetra::clone<Scalar, LocalOrdinal, GlobalOrdinal, Node, Node2>(*A, node2);
439  clonelevel->template Set<RCP<CloneOperator> >("A", cloneA);
440  }
441  if (level->IsAvailable("R")){
442  R = level->template Get<RCP<Operator> >("R");
443  cloneR = Xpetra::clone<Scalar, LocalOrdinal, GlobalOrdinal, Node, Node2>(*R, node2);
444  clonelevel->template Set<RCP<CloneOperator> >("R", cloneR);
445  }
446  if (level->IsAvailable("P")){
447  P = level->template Get<RCP<Operator> >("P");
448  cloneP = Xpetra::clone<Scalar, LocalOrdinal, GlobalOrdinal, Node, Node2>(*P, node2);
449  clonelevel->template Set<RCP<CloneOperator> >("P", cloneP);
450  }
451  if (level->IsAvailable("PreSmoother")){
452  Pre = level->template Get<RCP<SmootherBase> >("PreSmoother");
453  clonePre = MueLu::clone<Scalar, LocalOrdinal, GlobalOrdinal, Node, Node2> (Pre, cloneA, node2);
454  clonelevel->template Set<RCP<CloneSmoother> >("PreSmoother", clonePre);
455  }
456  if (level->IsAvailable("PostSmoother")){
457  Post = level->template Get<RCP<SmootherBase> >("PostSmoother");
458  clonePost = MueLu::clone<Scalar, LocalOrdinal, GlobalOrdinal, Node, Node2> (Post, cloneA, node2);
459  clonelevel-> template Set<RCP<CloneSmoother> >("PostSmoother", clonePost);
460  }
461  new_h->Levels_[i] = clonelevel;
462  }
463 
464  return new_h;
465  }
466 
467 } //namespace MueLu
468 
469 #define MUELU_HIERARCHY_SHORT
470 #endif // MUELU_HIERARCHY_DECL_HPP
void IsPreconditioner(const bool flag)
Array< RCP< MultiVector > > coarseExport_
This class specifies the default factory that should generate some data on a Level if the data does n...
double GetSmootherComplexity() const
RCP< Level > & GetLevel(const int levelID=0)
Retrieve a certain level from hierarchy.
Teuchos::RCP< Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node2 > > clone(const RCP< Node2 > &node2) const
Hierarchy()
Default constructor.
Array< RCP< MultiVector > > coarseImport_
Print parameters.
void AddLevel(const RCP< Level > &level)
Add a level at the end of the hierarchy.
short KeepType
void SetMaxCoarseSize(Xpetra::global_size_t maxCoarseSize)
void AddNewLevel()
Add a new level at the end of the hierarchy.
static Xpetra::global_size_t GetDefaultMaxCoarseSize()
void SetCycle(CycleType Cycle)
Supports VCYCLE and WCYCLE types.
void Clear(int startLevel=0)
Clear impermanent data from previous setup.
Namespace for MueLu classes and methods.
Teuchos::ScalarTraits< SC > STS
MagnitudeType rate_
Convergece rate.
static const NoFactory * get()
void ReplaceCoordinateMap(Level &level)
void AllocateLevelMultiVectors(int numvecs)
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 print(std::ostream &out=std::cout, const VerbLevel verbLevel=(MueLu::Parameters|MueLu::Statistics0)) const
Xpetra::UnderlyingLib lib()
Base class for factories (e.g., R, P, and A_coarse).
Print statistics that do not involve significant additional computation.
ConvData(std::pair< LO, MagnitudeType > p)
Xpetra::UnderlyingLib lib_
CycleType GetCycle() const
Returns multigrid cycle type (supports VCYCLE and WCYCLE)
static CycleType GetDefaultCycle()
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.
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.
const RCP< const FactoryManagerBase > & GetLevelManager(const int levelID) const
Always keep data, even accross run. This flag is set by Level::Keep(). This flag is propagated to coa...
static bool GetDefaultImplicitTranspose()
void describe(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the Hierarchy with some verbosity level to a FancyOStream object.
bool GetImplicitTranspose() const
void SetPRrebalance(bool doPRrebalance)
Base class for smoothers.
double GetOperatorComplexity() const
static int GetDefaultMaxLevels()
Base class for MueLu classes.
void SetImplicitTranspose(const bool &implicit)
virtual ~Hierarchy()
Destructor.
Data struct for defining stopping criteria of multigrid iteration.
Array< RCP< MultiVector > > coarseRhs_
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.
MagnitudeType GetRate() const
Array< RCP< Level > > Levels_
Container for Level objects.
static bool GetDefaultPRrebalance()
Xpetra::global_size_t GetMaxCoarseSize() const
void SetProlongatorScalingFactor(double scalingFactor)
Specify damping factor alpha such that x = x + alpha*P*c, where c is the coarse grid correction...
ReturnType Iterate(const MultiVector &B, MultiVector &X, ConvData conv=ConvData(), bool InitialGuessIsZero=false, LO startLevel=0)
Apply the multigrid preconditioner.
Array< RCP< MultiVector > > residual_
void Keep(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Call Level::Keep(ename, factory) for each level of the Hierarchy.
bool isDumpingEnabled_
Graph dumping.
void setlib(Xpetra::UnderlyingLib inlib)
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.
Array< RCP< MultiVector > > coarseX_
Array< RCP< MultiVector > > correction_
void EnableGraphDumping(const std::string &filename, int levelID=1)
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
Array< RCP< const FactoryManagerBase > > levelManagers_
Xpetra::global_size_t maxCoarseSize_