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 #include "MueLu_PerfUtils.hpp"
65 
66 #ifdef HAVE_MUELU_INTREPID2
67 #include "Kokkos_DynRankView.hpp"
68 #endif
69 
70 namespace MueLu {
71 
72  // This class stores the configuration of a Hierarchy.
73  // The class also provides an algorithm to build a Hierarchy from the configuration.
74  //
75  // See also: FactoryManager
76  //
77  template <class Scalar = double, class LocalOrdinal = int, class GlobalOrdinal = LocalOrdinal, class Node = KokkosClassic::DefaultNode::DefaultNodeType>
78  class HierarchyManager : public HierarchyFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node> {
79 #undef MUELU_HIERARCHYMANAGER_SHORT
80 #include "MueLu_UseShortNames.hpp"
81  typedef std::pair<std::string, const FactoryBase*> keep_pair;
82 
83  public:
84 
86  HierarchyManager(int numDesiredLevel = MasterList::getDefault<int>("max levels")) :
87  numDesiredLevel_ (numDesiredLevel),
88  maxCoarseSize_ (MasterList::getDefault<int>("coarse: max size")),
90  doPRrebalance_ (MasterList::getDefault<bool>("repartition: rebalance P and R")),
91  implicitTranspose_ (MasterList::getDefault<bool>("transpose: use implicit")),
92  graphOutputLevel_(-1) { }
93 
95  virtual ~HierarchyManager() { }
96 
98  void AddFactoryManager(int startLevel, int numDesiredLevel, RCP<FactoryManagerBase> manager) {
99  const int lastLevel = startLevel + numDesiredLevel - 1;
100  if (levelManagers_.size() < lastLevel + 1)
101  levelManagers_.resize(lastLevel + 1);
102 
103  for (int iLevel = startLevel; iLevel <= lastLevel; iLevel++)
104  levelManagers_[iLevel] = manager;
105  }
106 
109  // NOTE: last levelManager is used for all the remaining levels
110  return (levelID >= levelManagers_.size() ? levelManagers_[levelManagers_.size()-1] : levelManagers_[levelID]);
111  }
112 
114  size_t getNumFactoryManagers() const {
115  return levelManagers_.size();
116  }
117 
119  void CheckConfig() {
120  for (int i = 0; i < levelManagers_.size(); i++)
121  TEUCHOS_TEST_FOR_EXCEPTION(levelManagers_[i] == Teuchos::null, Exceptions::RuntimeError, "MueLu:HierarchyConfig::CheckConfig(): Undefined configuration for level:");
122  }
123 
125 
126  virtual RCP<Hierarchy> CreateHierarchy() const {
127  return rcp(new Hierarchy());
128  }
129 
130  virtual RCP<Hierarchy> CreateHierarchy(const std::string& label) const {
131  return rcp(new Hierarchy(label));
132  }
133 
135  virtual void SetupHierarchy(Hierarchy& H) const {
136  TEUCHOS_TEST_FOR_EXCEPTION(!H.GetLevel(0)->IsAvailable("A"), Exceptions::RuntimeError, "No fine level operator");
137 
138 #ifdef HAVE_MUELU_DEBUG
139  // Reset factories' data used for debugging
140  for (int i = 0; i < levelManagers_.size(); i++)
141  levelManagers_[i]->ResetDebugData();
142 
143 #endif
144 
145  // Setup Matrix
146  // TODO: I should certainly undo this somewhere...
147  RCP<Level> l0 = H.GetLevel(0);
148  RCP<Operator> Op = l0->Get<RCP<Operator> >("A");
149 
150  Xpetra::UnderlyingLib lib = Op->getDomainMap()->lib();
151  H.setlib(lib);
152 
153  SetupOperator(*Op);
154  SetupExtra(H);
155 
156  // Setup Hierarchy
159  if (graphOutputLevel_ >= 0)
160  H.EnableGraphDumping("dep_graph.dot", graphOutputLevel_);
161 
163  RCP<Matrix> Amat = rcp_dynamic_cast<Matrix>(Op);
164 
165  if (!Amat.is_null()) {
166  RCP<ParameterList> params = rcp(new ParameterList());
167  params->set("printLoadBalancingInfo", true);
168  params->set("printCommInfo", true);
169 
171  } else {
172  VerboseObject::GetOStream(Warnings1) << "Fine level operator is not a matrix, statistics are not available" << std::endl;
173  }
174  }
175 
178 
179  H.Clear();
180 
181  // There are few issues with using Keep in the interpreter:
182  // 1. Hierarchy::Keep interface takes a name and a factory. If
183  // factories are different on different levels, the AddNewLevel() call
184  // in Hierarchy does not work properly, as it assume that factories are
185  // the same.
186  // 2. FactoryManager does not have a Keep option, only Hierarchy and
187  // Level have it
188  // 3. Interpreter constructs factory managers, but not levels. So we
189  // cannot set up Keep flags there.
190  //
191  // The solution implemented here does the following:
192  // 1. Construct hierarchy with dummy levels. This avoids
193  // Hierarchy::AddNewLevel() calls which will propagate wrong
194  // inheritance.
195  // 2. Interpreter constructs keep_ array with names and factories for
196  // that level
197  // 3. For each level, we call Keep(name, factory) for each keep_
198  for (int i = 0; i < numDesiredLevel_; i++) {
199  std::map<int, std::vector<keep_pair> >::const_iterator it = keep_.find(i);
200  if (it != keep_.end()) {
201  RCP<Level> l = H.GetLevel(i);
202  const std::vector<keep_pair>& keeps = it->second;
203  for (size_t j = 0; j < keeps.size(); j++)
204  l->Keep(keeps[j].first, keeps[j].second);
205  }
206  if (i < numDesiredLevel_-1) {
207  RCP<Level> newLevel = rcp(new Level());
208  H.AddLevel(newLevel);
209  }
210  }
214  ExportDataSetKeepFlags(H, nullspaceToPrint_, "Nullspace");
215  ExportDataSetKeepFlags(H, coordinatesToPrint_, "Coordinates");
216 #ifdef HAVE_MUELU_INTREPID2
217  ExportDataSetKeepFlags(H,elementToNodeMapsToPrint_, "pcoarsen: element to node map");
218 #endif
219 
220  int levelID = 0;
221  int lastLevelID = numDesiredLevel_ - 1;
222  bool isLastLevel = false;
223 
224  while (!isLastLevel) {
225  bool r = H.Setup(levelID,
226  LvlMngr(levelID-1, lastLevelID),
227  LvlMngr(levelID, lastLevelID),
228  LvlMngr(levelID+1, lastLevelID));
229 
230  isLastLevel = r || (levelID == lastLevelID);
231  levelID++;
232  }
233  // FIXME: Should allow specification of NumVectors on parameterlist
235 
237  fos->setOutputToRootOnly(0);
238  H.describe(*fos, verbosity_);
239 
240  // When we reuse hierarchy, it is necessary that we don't
241  // change the number of levels. We also cannot make requests
242  // for coarser levels, because we don't construct all the
243  // data on previous levels. For instance, let's say our first
244  // run constructed three levels. If we try to do requests during
245  // next setup for the fourth level, it would need Aggregates
246  // which we didn't construct for level 3 because we reused P.
247  // To fix this situation, we change the number of desired levels
248  // here.
249  numDesiredLevel_ = levelID;
250 
251  WriteData<Matrix>(H, matricesToPrint_, "A");
252  WriteData<Matrix>(H, prolongatorsToPrint_, "P");
253  WriteData<Matrix>(H, restrictorsToPrint_, "R");
254  WriteData<MultiVector>(H, nullspaceToPrint_, "Nullspace");
255  WriteData<MultiVector>(H, coordinatesToPrint_, "Coordinates");
256 #ifdef HAVE_MUELU_INTREPID2
257  typedef Kokkos::DynRankView<LocalOrdinal,typename Node::device_type> FCi;
258  WriteDataFC<FCi>(H,elementToNodeMapsToPrint_, "pcoarsen: element to node map","el2node");
259 #endif
260 
261 
262  } //SetupHierarchy
263 
265 
266  typedef std::map<std::string, RCP<const FactoryBase> > FactoryMap;
267 
268  protected: //TODO: access function
269 
271  virtual void SetupOperator(Operator& Op) const { }
272 
274  // TODO: merge with SetupMatrix ?
275  virtual void SetupExtra(Hierarchy& H) const { }
276 
277  // TODO this was private
278  // Used in SetupHierarchy() to access levelManagers_
279  // Inputs i=-1 and i=size() are allowed to simplify calls to hierarchy->Setup()
280  Teuchos::RCP<FactoryManagerBase> LvlMngr(int levelID, int lastLevelID) const {
281  // NOTE: the order of 'if' statements is important
282  if (levelID == -1) // levelID = -1 corresponds to the finest level
283  return Teuchos::null;
284 
285  if (levelID == lastLevelID+1) // levelID = 'lastLevelID+1' corresponds to the last level (i.e., no nextLevel)
286  return Teuchos::null;
287 
288  if (levelManagers_.size() == 0) { // default factory manager.
289  // The default manager is shared across levels, initialized only if needed and deleted with the HierarchyManager
290  static RCP<FactoryManagerBase> defaultMngr = rcp(new FactoryManager());
291  return defaultMngr;
292  }
293 
294  return GetFactoryManager(levelID);
295  }
296 
297  // Hierarchy parameters
298  mutable int numDesiredLevel_;
310 
311  std::map<int, std::vector<keep_pair> > keep_;
312 
313  private:
314  // Set the keep flags for Export Data
315  void ExportDataSetKeepFlags(Hierarchy& H, const Teuchos::Array<int>& data, const std::string& name) const {
316  for (int i = 0; i < data.size(); ++i) {
317  if (data[i] < H.GetNumLevels()) {
318  RCP<Level> L = H.GetLevel(data[i]);
319  if(!L.is_null() && data[i] < levelManagers_.size())
320  L->AddKeepFlag(name, &*levelManagers_[data[i]]->GetFactory(name));
321  }
322  }
323  }
324 
325 
326  template<class T>
327  void WriteData(Hierarchy& H, const Teuchos::Array<int>& data, const std::string& name) const {
328  for (int i = 0; i < data.size(); ++i) {
329  std::string fileName = name + "_" + Teuchos::toString(data[i]) + ".m";
330 
331  if (data[i] < H.GetNumLevels()) {
332  RCP<Level> L = H.GetLevel(data[i]);
333  if (data[i] < levelManagers_.size() && L->IsAvailable(name,&*levelManagers_[i]->GetFactory(name))) {
334  // Try generating factory
335  RCP<T> M = L->template Get< RCP<T> >(name,&*levelManagers_[i]->GetFactory(name));
336  if (!M.is_null()) {
338  }
339  }
340  else if (L->IsAvailable(name)) {
341  // Try nofactory
342  RCP<T> M = L->template Get< RCP<T> >(name);
343  if (!M.is_null()) {
345  }
346  }
347 
348  }
349  }
350  }
351 
352 
353  template<class T>
354  void WriteDataFC(Hierarchy& H, const Teuchos::Array<int>& data, const std::string& name, const std::string & ofname) const {
355  for (int i = 0; i < data.size(); ++i) {
356  const std::string fileName = ofname + "_" + Teuchos::toString(data[i]) + ".m";
357 
358  if (data[i] < H.GetNumLevels()) {
359  RCP<Level> L = H.GetLevel(data[i]);
360 
361  if (L->IsAvailable(name)) {
362  RCP<T> M = L->template Get< RCP<T> >(name);
363  if (!M.is_null()) {
364  RCP<Matrix> A = L->template Get<RCP<Matrix> >("A");
365  RCP<const CrsGraph> AG = A->getCrsGraph();
366  WriteFieldContainer<T>(fileName,*M,*AG->getColMap());
367  }
368  }
369  }
370  }
371  }
372 
373  // For dumping an IntrepidPCoarsening element-to-node map to disk
374  template<class T>
375  void WriteFieldContainer(const std::string& fileName, T & fcont,const Map &colMap) const {
376  typedef LocalOrdinal LO;
377  typedef GlobalOrdinal GO;
378  typedef Node NO;
379  typedef Xpetra::MultiVector<GO,LO,GO,NO> GOMultiVector;
380 
381  size_t num_els = (size_t) fcont.extent(0);
382  size_t num_vecs =(size_t) fcont.extent(1);
383 
384  // Generate rowMap
385  Teuchos::RCP<const Map> rowMap = Xpetra::MapFactory<LO,GO,NO>::Build(colMap.lib(),Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),fcont.extent(0),colMap.getIndexBase(),colMap.getComm());
386 
387  // Fill multivector to use *petra dump routines
389 
390  for(size_t j=0; j<num_vecs; j++) {
391  Teuchos::ArrayRCP<GO> v = vec->getDataNonConst(j);
392  for(size_t i=0; i<num_els; i++)
393  v[i] = colMap.getGlobalElement(fcont(i,j));
394  }
395 
396  Xpetra::IO<GO,LO,GO,NO>::Write(fileName,*vec);
397  }
398 
399 
400 
401  // Levels
402  Array<RCP<FactoryManagerBase> > levelManagers_; // one FactoryManager per level (the last levelManager is used for all the remaining levels)
403 
404  }; // class HierarchyManager
405 
406 } // namespace MueLu
407 
408 #define MUELU_HIERARCHYMANAGER_SHORT
409 #endif // MUELU_HIERARCHYMANAGER_HPP
410 
411 //TODO: split into _decl/_def
412 // TODO: default value for first param (FactoryManager()) should not be duplicated (code maintainability)
virtual void SetupOperator(Operator &Op) const
Setup Matrix object.
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.
GlobalOrdinal GO
void AddLevel(const RCP< Level > &level)
Add a level at the end of the hierarchy.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::Array< int > elementToNodeMapsToPrint_
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
void SetMaxCoarseSize(Xpetra::global_size_t maxCoarseSize)
void AddFactoryManager(int startLevel, int numDesiredLevel, RCP< FactoryManagerBase > manager)
LocalOrdinal LO
virtual RCP< FancyOStream > getOStream() const
void WriteData(Hierarchy &H, const Teuchos::Array< int > &data, const std::string &name) const
void Clear(int startLevel=0)
Clear impermanent data from previous setup.
Namespace for MueLu classes and methods.
RCP< FactoryManagerBase > GetFactoryManager(int levelID) const
std::map< std::string, RCP< const FactoryBase > > FactoryMap
Static class that holds the complete list of valid MueLu parameters.
bool is_null() const
virtual void SetupHierarchy(Hierarchy &H) const
Setup Hierarchy object.
void AllocateLevelMultiVectors(int numvecs)
Print even more statistics.
Additional warnings.
Teuchos::FancyOStream & GetOStream(MsgType type, int thisProcRankOnly=0) const
Get an output stream for outputting the input message type.
Teuchos::RCP< FactoryManagerBase > LvlMngr(int levelID, int lastLevelID) const
static void SetDefaultVerbLevel(const VerbLevel defaultVerbLevel)
Set the default (global) verbosity level.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
Teuchos::Array< int > nullspaceToPrint_
bool IsPrint(MsgType type, int thisProcRankOnly=-1) const
Find out whether we need to print out information for a specific message type.
size_type size() const
Array< RCP< FactoryManagerBase > > levelManagers_
void ExportDataSetKeepFlags(Hierarchy &H, const Teuchos::Array< int > &data, const std::string &name) const
void describe(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the Hierarchy with some verbosity level to a FancyOStream object.
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
void SetPRrebalance(bool doPRrebalance)
size_t global_size_t
void SetImplicitTranspose(const bool &implicit)
void WriteDataFC(Hierarchy &H, const Teuchos::Array< int > &data, const std::string &name, const std::string &ofname) const
virtual RCP< Hierarchy > CreateHierarchy() const
Create an empty Hierarchy object.
Teuchos::Array< int > prolongatorsToPrint_
virtual void SetupExtra(Hierarchy &H) const
Setup extra data.
virtual RCP< Hierarchy > CreateHierarchy(const std::string &label) const
Create a labeled empty Hierarchy object.
void WriteFieldContainer(const std::string &fileName, T &fcont, const Map &colMap) const
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &map, size_t NumVectors, bool zeroOut=true)
Teuchos::Array< int > coordinatesToPrint_
Node NO
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
size_t getNumFactoryManagers() const
returns number of factory managers stored in levelManagers_ vector.
void setlib(Xpetra::UnderlyingLib inlib)
Xpetra::global_size_t maxCoarseSize_
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.
static Teuchos::RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, LocalGlobal lg=Xpetra::GloballyDistributed, const Teuchos::RCP< Node > &=Teuchos::null)
Teuchos::Array< int > restrictorsToPrint_
void EnableGraphDumping(const std::string &filename, int levelID=1)
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
std::string toString(const T &t)