MueLu  Version of the Day
MueLu_HierarchyManager.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_HIERARCHYMANAGER_DECL_HPP
47 #define MUELU_HIERARCHYMANAGER_DECL_HPP
48 
49 #include <string>
50 #include <map>
51 
52 #include <Teuchos_Array.hpp>
53 
54 #include <Xpetra_Operator.hpp>
55 #include <Xpetra_IO.hpp>
56 
57 #include "MueLu_ConfigDefs.hpp"
58 
59 #include "MueLu_Exceptions.hpp"
60 #include "MueLu_Hierarchy.hpp"
62 #include "MueLu_Level.hpp"
63 #include "MueLu_MasterList.hpp"
64 
65 namespace MueLu {
66 
67  // This class stores the configuration of a Hierarchy.
68  // The class also provides an algorithm to build a Hierarchy from the configuration.
69  //
70  // See also: FactoryManager
71  //
72  template <class Scalar = double, class LocalOrdinal = int, class GlobalOrdinal = LocalOrdinal, class Node = KokkosClassic::DefaultNode::DefaultNodeType>
73  class HierarchyManager : public HierarchyFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node> {
74 #undef MUELU_HIERARCHYMANAGER_SHORT
75 #include "MueLu_UseShortNames.hpp"
76  typedef std::pair<std::string, const FactoryBase*> keep_pair;
77 
78  public:
79 
81  HierarchyManager(int numDesiredLevel = MasterList::getDefault<int>("max levels")) :
82  numDesiredLevel_ (numDesiredLevel),
83  maxCoarseSize_ (MasterList::getDefault<int>("coarse: max size")),
85  doPRrebalance_ (MasterList::getDefault<bool>("repartition: rebalance P and R")),
86  implicitTranspose_ (MasterList::getDefault<bool>("transpose: use implicit")),
87  graphOutputLevel_(-1) { }
88 
90  virtual ~HierarchyManager() { }
91 
93  void AddFactoryManager(int startLevel, int numDesiredLevel, RCP<FactoryManagerBase> manager) {
94  const int lastLevel = startLevel + numDesiredLevel - 1;
95  if (levelManagers_.size() < lastLevel + 1)
96  levelManagers_.resize(lastLevel + 1);
97 
98  for (int iLevel = startLevel; iLevel <= lastLevel; iLevel++)
99  levelManagers_[iLevel] = manager;
100  }
101 
103  RCP<FactoryManagerBase> GetFactoryManager(int levelID) const {
104  // NOTE: last levelManager is used for all the remaining levels
105  return (levelID >= levelManagers_.size() ? levelManagers_[levelManagers_.size()-1] : levelManagers_[levelID]);
106  }
107 
109  size_t getNumFactoryManagers() const {
110  return levelManagers_.size();
111  }
112 
114  void CheckConfig() {
115  for (int i = 0; i < levelManagers_.size(); i++)
116  TEUCHOS_TEST_FOR_EXCEPTION(levelManagers_[i] == Teuchos::null, Exceptions::RuntimeError, "MueLu:HierarchyConfig::CheckConfig(): Undefined configuration for level:");
117  }
118 
120 
121  virtual RCP<Hierarchy> CreateHierarchy() const {
122  return rcp(new Hierarchy());
123  }
124 
126  virtual void SetupHierarchy(Hierarchy& H) const {
127  TEUCHOS_TEST_FOR_EXCEPTION(!H.GetLevel(0)->IsAvailable("A"), Exceptions::RuntimeError, "No fine level operator");
128 
129 #ifdef HAVE_MUELU_DEBUG
130  // Reset factories' data used for debugging
131  for (int i = 0; i < levelManagers_.size(); i++)
132  levelManagers_[i]->ResetDebugData();
133 
134 #endif
135 
136  // Setup Matrix
137  // TODO: I should certainly undo this somewhere...
138  RCP<Level> l0 = H.GetLevel(0);
139  RCP<Operator> Op = l0->Get<RCP<Operator> >("A");
140 
141  Xpetra::UnderlyingLib lib = Op->getDomainMap()->lib();
142  H.setlib(lib);
143 
144  SetupOperator(*Op);
145  SetupExtra(H);
146 
147  // Setup Hierarchy
150  if (graphOutputLevel_ >= 0)
151  H.EnableGraphDumping("dep_graph.dot", graphOutputLevel_);
152 
155 
156  H.Clear();
157 
158  // There are few issues with using Keep in the interpreter:
159  // 1. Hierarchy::Keep interface takes a name and a factory. If
160  // factories are different on different levels, the AddNewLevel() call
161  // in Hierarchy does not work properly, as it assume that factories are
162  // the same.
163  // 2. FactoryManager does not have a Keep option, only Hierarchy and
164  // Level have it
165  // 3. Interpreter constructs factory managers, but not levels. So we
166  // cannot set up Keep flags there.
167  //
168  // The solution implemented here does the following:
169  // 1. Construct hierarchy with dummy levels. This avoids
170  // Hierarchy::AddNewLevel() calls which will propagate wrong
171  // inheritance.
172  // 2. Interpreter constructs keep_ array with names and factories for
173  // that level
174  // 3. For each level, we call Keep(name, factory) for each keep_
175  for (int i = 0; i < numDesiredLevel_; i++) {
176  std::map<int, std::vector<keep_pair> >::const_iterator it = keep_.find(i);
177  if (it != keep_.end()) {
178  RCP<Level> l = H.GetLevel(i);
179  const std::vector<keep_pair>& keeps = it->second;
180  for (size_t j = 0; j < keeps.size(); j++)
181  l->Keep(keeps[j].first, keeps[j].second);
182  }
183  if (i < numDesiredLevel_-1) {
184  RCP<Level> newLevel = rcp(new Level());
185  H.AddLevel(newLevel);
186  }
187  }
188 
189  int levelID = 0;
190  int lastLevelID = numDesiredLevel_ - 1;
191  bool isLastLevel = false;
192 
193  while (!isLastLevel) {
194  bool r = H.Setup(levelID,
195  LvlMngr(levelID-1, lastLevelID),
196  LvlMngr(levelID, lastLevelID),
197  LvlMngr(levelID+1, lastLevelID));
198 
199  isLastLevel = r || (levelID == lastLevelID);
200  levelID++;
201  }
202 
203  RCP<Teuchos::FancyOStream> fos = this->getOStream();
204  fos->setOutputToRootOnly(0);
205  H.describe(*fos, verbosity_);
206 
207  // When we reuse hierarchy, it is necessary that we don't
208  // change the number of levels. We also cannot make requests
209  // for coarser levels, because we don't construct all the
210  // data on previous levels. For instance, let's say our first
211  // run constructed three levels. If we try to do requests during
212  // next setup for the fourth level, it would need Aggregates
213  // which we didn't construct for level 3 because we reused P.
214  // To fix this situation, we change the number of desired levels
215  // here.
216  numDesiredLevel_ = levelID;
217 
218  WriteData<Matrix>(H, matricesToPrint_, "A");
219  WriteData<Matrix>(H, prolongatorsToPrint_, "P");
220  WriteData<Matrix>(H, restrictorsToPrint_, "R");
221 
222  } //SetupHierarchy
223 
225 
226  typedef std::map<std::string, RCP<const FactoryBase> > FactoryMap;
227 
228  protected: //TODO: access function
229 
231  virtual void SetupOperator(Operator& Op) const { }
232 
234  // TODO: merge with SetupMatrix ?
235  virtual void SetupExtra(Hierarchy& H) const { }
236 
237  // TODO this was private
238  // Used in SetupHierarchy() to access levelManagers_
239  // Inputs i=-1 and i=size() are allowed to simplify calls to hierarchy->Setup()
240  Teuchos::RCP<FactoryManagerBase> LvlMngr(int levelID, int lastLevelID) const {
241  // NOTE: the order of 'if' statements is important
242  if (levelID == -1) // levelID = -1 corresponds to the finest level
243  return Teuchos::null;
244 
245  if (levelID == lastLevelID+1) // levelID = 'lastLevelID+1' corresponds to the last level (i.e., no nextLevel)
246  return Teuchos::null;
247 
248  if (levelManagers_.size() == 0) { // default factory manager.
249  // The default manager is shared across levels, initialized only if needed and deleted with the HierarchyManager
250  static RCP<FactoryManagerBase> defaultMngr = rcp(new FactoryManager());
251  return defaultMngr;
252  }
253 
254  return GetFactoryManager(levelID);
255  }
256 
257  // Hierarchy parameters
258  mutable int numDesiredLevel_;
259  Xpetra::global_size_t maxCoarseSize_;
264  Teuchos::Array<int> matricesToPrint_;
265  Teuchos::Array<int> prolongatorsToPrint_;
266  Teuchos::Array<int> restrictorsToPrint_;
267 
268  std::map<int, std::vector<keep_pair> > keep_;
269 
270  private:
271 
272  template<class T>
273  void WriteData(Hierarchy& H, const Teuchos::Array<int>& data, const std::string& name) const {
274  for (int i = 0; i < data.size(); ++i) {
275  std::string fileName = name + "_" + Teuchos::toString(data[i]) + ".m";
276 
277  if (data[i] < H.GetNumLevels()) {
278  RCP<Level> L = H.GetLevel(data[i]);
279 
280  if (L->IsAvailable(name)) {
281  RCP<T> M = L->template Get< RCP<T> >(name);
282  if (!M.is_null())
283  Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(fileName,* M);
284  }
285  }
286  }
287  }
288 
289  // Levels
290  Array<RCP<FactoryManagerBase> > levelManagers_; // one FactoryManager per level (the last levelManager is used for all the remaining levels)
291 
292  }; // class HierarchyManager
293 
294 } // namespace MueLu
295 
296 #define MUELU_HIERARCHYMANAGER_SHORT
297 #endif // MUELU_HIERARCHYMANAGER_HPP
298 
299 //TODO: split into _decl/_def
300 // TODO: default value for first param (FactoryManager()) should not be duplicated (code maintainability)
virtual void SetupExtra(Hierarchy &H) const
Setup extra data.
This class specifies the default factory that should generate some data on a Level if the data does n...
Teuchos::Array< int > matricesToPrint_
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.
static void SetDefaultVerbLevel(const VerbLevel defaultVerbLevel)
Set the default (global) verbosity level.
void AddLevel(const RCP< Level > &level)
Add a level at the end of the hierarchy.
void SetMaxCoarseSize(Xpetra::global_size_t maxCoarseSize)
void AddFactoryManager(int startLevel, int numDesiredLevel, RCP< FactoryManagerBase > manager)
void Clear(int startLevel=0)
Clear impermanent data from previous setup.
Namespace for MueLu classes and methods.
std::map< std::string, RCP< const FactoryBase > > FactoryMap
virtual void SetupOperator(Operator &Op) const
Setup Matrix object.
Static class that holds the complete list of valid MueLu parameters.
Teuchos::RCP< FactoryManagerBase > LvlMngr(int levelID, int lastLevelID) const
void WriteData(Hierarchy &H, const Teuchos::Array< int > &data, const std::string &name) const
size_t getNumFactoryManagers() const
returns number of factory managers stored in levelManagers_ vector.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
Array< RCP< FactoryManagerBase > > levelManagers_
RCP< FactoryManagerBase > GetFactoryManager(int levelID) const
void SetPRrebalance(bool doPRrebalance)
void SetImplicitTranspose(const bool &implicit)
Teuchos::Array< int > prolongatorsToPrint_
std::map< int, std::vector< keep_pair > > keep_
HierarchyManager(int numDesiredLevel=MasterList::getDefault< int >("max levels"))
Exception throws to report errors in the internal logical of the program.
std::pair< std::string, const FactoryBase * > keep_pair
void describe(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the Hierarchy with some verbosity level to a FancyOStream object.
void setlib(Xpetra::UnderlyingLib inlib)
Xpetra::global_size_t maxCoarseSize_
virtual RCP< Hierarchy > CreateHierarchy() const
Create an empty Hierarchy 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.
Teuchos::Array< int > restrictorsToPrint_
void EnableGraphDumping(const std::string &filename, int levelID=1)
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
virtual void SetupHierarchy(Hierarchy &H) const
Setup Hierarchy object.