MueLu  Version of the Day
MueLu_HierarchyHelpers_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 #ifndef MUELU_HIERARCHYHELPERS_DEF_HPP
47 #define MUELU_HIERARCHYHELPERS_DEF_HPP
48 
49 #include <Xpetra_Matrix.hpp>
50 #include <Xpetra_Operator.hpp>
51 
54 #include "MueLu_SmootherBase.hpp"
55 #include "MueLu_SmootherFactory.hpp"
56 #include "MueLu_FactoryManager.hpp"
57 
58 //TODO/FIXME: DeclareInput(, **this**) cannot be used here
59 
60 namespace MueLu {
61 
62  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
63  TopRAPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TopRAPFactory(RCP<const FactoryManagerBase> parentFactoryManager) :
64  PFact_ (parentFactoryManager->GetFactory("P")),
65  RFact_ (parentFactoryManager->GetFactory("R")),
66  AcFact_(parentFactoryManager->GetFactory("A"))
67  { }
68 
69  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
70  TopRAPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TopRAPFactory(RCP<const FactoryManagerBase> parentFactoryManagerFine, RCP<const FactoryManagerBase> parentFactoryManagerCoarse) :
71  PFact_ (parentFactoryManagerCoarse->GetFactory("P")),
72  RFact_ (parentFactoryManagerCoarse->GetFactory("R")),
73  AcFact_(parentFactoryManagerCoarse->GetFactory("A"))
74  { }
75 
76  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
78 
79  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
81  if (PFact_ != Teuchos::null) coarseLevel.DeclareInput("P", PFact_.get());
82  if (RFact_ != Teuchos::null) coarseLevel.DeclareInput("R", RFact_.get());
83  if ((AcFact_ != Teuchos::null) && (AcFact_ != NoFactory::getRCP())) coarseLevel.DeclareInput("A", AcFact_.get());
84  }
85 
86  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
88  if ((PFact_ != Teuchos::null) && (PFact_ != NoFactory::getRCP())) {
89  RCP<Operator> oP = coarseLevel.Get<RCP<Operator> >("P", PFact_.get());
90  RCP<Matrix> P = rcp_dynamic_cast<Matrix>(oP);
91  if (!P.is_null()) coarseLevel.Set("P", P, NoFactory::get());
92  else coarseLevel.Set("P", oP, NoFactory::get());
93  coarseLevel.AddKeepFlag ("P", NoFactory::get(), MueLu::Final); // FIXME2: Order of Remove/Add matter (data removed otherwise). Should do something about this
94  coarseLevel.RemoveKeepFlag("P", NoFactory::get(), MueLu::UserData); // FIXME: This is a hack, I should change behavior of Level::Set() instead. FIXME3: Should not be removed if flag was there already
95 
96  }
97 
98  if ((RFact_ != Teuchos::null) && (RFact_ != NoFactory::getRCP()) ) {
99  RCP<Operator> oR = coarseLevel.Get<RCP<Operator> >("R", RFact_.get());
100  RCP<Matrix> R = rcp_dynamic_cast<Matrix>(oR);
101  if (!R.is_null()) coarseLevel.Set("R", R, NoFactory::get());
102  else coarseLevel.Set("R", oR, NoFactory::get());
103  coarseLevel.AddKeepFlag ("R", NoFactory::get(), MueLu::Final);
104  coarseLevel.RemoveKeepFlag("R", NoFactory::get(), MueLu::UserData); // FIXME: This is a hack
105  }
106 
107  if ((AcFact_ != Teuchos::null) && (AcFact_ != NoFactory::getRCP())) {
108  RCP<Operator> oA = coarseLevel.Get<RCP<Operator> >("A", AcFact_.get());
109  RCP<Matrix> A = rcp_dynamic_cast<Matrix>(oA);
110  if (!A.is_null()) coarseLevel.Set("A", A, NoFactory::get());
111  else coarseLevel.Set("A", oA, NoFactory::get());
112  coarseLevel.AddKeepFlag ("A", NoFactory::get(), MueLu::Final);
113  coarseLevel.RemoveKeepFlag("A", NoFactory::get(), MueLu::UserData); // FIXME: This is a hack
114  }
115  }
116 
117  //
118  //
119  //
120 
121  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
122  TopSmootherFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TopSmootherFactory(RCP<const FactoryManagerBase> parentFactoryManager, const std::string& varName) {
123  TEUCHOS_TEST_FOR_EXCEPTION(varName != "CoarseSolver" && varName != "Smoother", Exceptions::RuntimeError, "varName should be either \"CoarseSolver\" or \"Smoother\"");
124 
125  if (varName == "CoarseSolver") {
126  // For coarsest level, we only need one smoother/solver
127  // If a user wants to do something weird there (like, solve coarsest system by using 2 forward
128  // GS and 1 backward GS), one can use MergedSmoother
129  RCP<const FactoryBase> coarseSolverFactory = parentFactoryManager->GetFactory("CoarseSolver");
130  RCP<const SmootherFactory> coarseSmootherFactory = Teuchos::rcp_dynamic_cast<const SmootherFactory>(coarseSolverFactory);
131  if (coarseSmootherFactory != Teuchos::null) {
132  RCP<SmootherPrototype> preProto;
133  RCP<SmootherPrototype> postProto;
134  coarseSmootherFactory->GetSmootherPrototypes(preProto, postProto);
135 
136  if (preProto == postProto)
137  preSmootherFact_ = parentFactoryManager->GetFactory("CoarseSolver");
138  else {
139  // check whether pre- and/or post-smoothing is desired on coarsest level
140  if(preProto != Teuchos::null)
141  preSmootherFact_ = parentFactoryManager->GetFactory("CoarseSolver");
142  if(postProto != Teuchos::null)
143  postSmootherFact_ = parentFactoryManager->GetFactory("CoarseSolver");
144  }
145  }
146  else // default handling: get default direct solver as presmoother on coarsest level
147  preSmootherFact_ = parentFactoryManager->GetFactory("CoarseSolver");
148 
149  } else {
150  preSmootherFact_ = parentFactoryManager->GetFactory("PreSmoother");
151  postSmootherFact_ = parentFactoryManager->GetFactory("PostSmoother");
152  }
153  }
154 
155  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
157 
158  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
160  if (preSmootherFact_ != Teuchos::null)
161  level.DeclareInput("PreSmoother", preSmootherFact_.get());
162  if (postSmootherFact_ != Teuchos::null)
163  level.DeclareInput("PostSmoother", postSmootherFact_.get());
164  }
165 
166  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
168  if (preSmootherFact_.is_null() && postSmootherFact_.is_null())
169  return;
170 
171  // NOTE 1: We need to set at least some keep flag for the smoothers, otherwise it is going to be removed as soon as all requests are released.
172  // We choose to set the Final flag for the data. In addition, we allow this data to be retrieved by only using the name by the means
173  // of using NoFactory. However, any data set with NoFactory gets UserData flag by default. We don't really want that flag, so we remove it.
174 
175  // NOTE 2: some smoother factories are tricky (see comments in MueLu::SmootherFactory
176  // Sometimes, we don't know whether the factory is able to generate "PreSmoother" or "PostSmoother"
177  // For the SmootherFactory, however, we are able to check that.
178 
179  if (!preSmootherFact_.is_null()) {
180  // Checking for null is not sufficient, as SmootherFactory(null, something) does not generate "PreSmoother"
181  bool isAble = true;
182  RCP<const SmootherFactory> s = rcp_dynamic_cast<const SmootherFactory>(preSmootherFact_);
183  if (!s.is_null()) {
184  RCP<SmootherPrototype> pre, post;
185  s->GetSmootherPrototypes(pre, post);
186  if (pre.is_null())
187  isAble = false;
188  } else {
189  // We assume that if presmoother factory is not SmootherFactory, it *is* able to generate "PreSmoother"
190  }
191 
192  if (isAble) {
193  RCP<SmootherBase> Pre = level.Get<RCP<SmootherBase> >("PreSmoother", preSmootherFact_.get());
194 
195  level.Set ("PreSmoother", Pre, NoFactory::get());
196 
197  level.AddKeepFlag ("PreSmoother", NoFactory::get(), MueLu::Final);
198  level.RemoveKeepFlag("PreSmoother", NoFactory::get(), MueLu::UserData);
199  }
200  }
201 
202  if (!postSmootherFact_.is_null()) {
203  // Checking for null is not sufficient, as SmootherFactory(something, null) does not generate "PostSmoother"
204  bool isAble = true;
205  RCP<const SmootherFactory> s = rcp_dynamic_cast<const SmootherFactory>(postSmootherFact_);
206  if (!s.is_null()) {
207  RCP<SmootherPrototype> pre, post;
208  s->GetSmootherPrototypes(pre, post);
209  if (post.is_null())
210  isAble = false;
211  } else {
212  // We assume that if presmoother factory is not SmootherFactory, it *is* able to generate "PreSmoother"
213  }
214 
215  if (isAble) {
216  RCP<SmootherBase> Post = level.Get<RCP<SmootherBase> >("PostSmoother", postSmootherFact_.get());
217 
218  level.Set ("PostSmoother", Post, NoFactory::get());
219 
220  level.AddKeepFlag ("PostSmoother", NoFactory::get(), MueLu::Final);
221  level.RemoveKeepFlag("PostSmoother", NoFactory::get(), MueLu::UserData);
222  }
223  }
224  }
225 
226 
227  // Adds the following non-serializable data (A,P,R,Nullspace,Coordinates) from level-specific sublist nonSerialList,
228  // calling AddNewLevel as appropriate.
229  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
231  for (ParameterList::ConstIterator it = paramList.begin(); it != paramList.end(); it++) {
232  const std::string& levelName = it->first;
233 
234  // Check for mach of the form "level X" where X is a positive integer
235  if (paramList.isSublist(levelName) && levelName.find("level ") == 0 && levelName.size() > 6) {
236  int levelID = strtol(levelName.substr(6).c_str(), 0, 0);
237  if (levelID > 0)
238  {
239  // Do enough level adding so we can be sure to add the data to the right place
240  for (int i = H.GetNumLevels(); i <= levelID; i++)
241  H.AddNewLevel();
242  }
243  RCP<Level> level = H.GetLevel(levelID);
244 
245  RCP<FactoryManager> M = Teuchos::rcp_dynamic_cast<FactoryManager>(HM.GetFactoryManager(levelID));
246  TEUCHOS_TEST_FOR_EXCEPTION(M.is_null(), Exceptions::InvalidArgument, "MueLu::Utils::AddNonSerializableDataToHierarchy: cannot get FactoryManager");
247 
248  // Grab the level sublist & loop over parameters
249  const ParameterList& levelList = paramList.sublist(levelName);
250  for (ParameterList::ConstIterator it2 = levelList.begin(); it2 != levelList.end(); it2++) {
251  const std::string& name = it2->first;
252  TEUCHOS_TEST_FOR_EXCEPTION(name != "A" && name != "P" && name != "R" &&
253  name != "Nullspace" && name != "Coordinates" &&
255  "MueLu::Utils::AddNonSerializableDataToHierarchy: parameter list contains unknown data type");
256 
257  if (name == "A") {
258  level->Set(name, Teuchos::getValue<RCP<Matrix > > (it2->second),NoFactory::get());
259  M->SetFactory(name, NoFactory::getRCP()); // TAW: not sure about this: be aware that this affects all levels
260  // However, A is accessible through NoFactory anyway, so it should
261  // be fine here.
262  }
263  else if( name == "P" || name == "R") {
264  level->AddKeepFlag(name,NoFactory::get(),MueLu::UserData);
265  level->Set(name, Teuchos::getValue<RCP<Matrix > > (it2->second), M->GetFactory(name).get());
266  }
267  else if (name == "Nullspace")
268  {
269  level->AddKeepFlag(name,NoFactory::get(),MueLu::UserData);
270  level->Set(name, Teuchos::getValue<RCP<MultiVector > >(it2->second), NoFactory::get());
271  //M->SetFactory(name, NoFactory::getRCP()); // TAW: generally it is a bad idea to overwrite the factory manager data here
272  // One should do this only in very special cases
273  }
274  else if(name == "Coordinates") //Scalar of Coordinates MV is always double
275  {
276  level->AddKeepFlag(name,NoFactory::get(),MueLu::UserData);
277  level->Set(name, Teuchos::getValue<RCP<Xpetra::MultiVector<double, LocalOrdinal, GlobalOrdinal, Node> > >(it2->second), NoFactory::get());
278  //M->SetFactory(name, NoFactory::getRCP()); // TAW: generally it is a bad idea to overwrite the factory manager data here
279  }
280  #ifdef HAVE_MUELU_MATLAB
281  else
282  {
283  //Custom variable for Muemex
284  size_t typeNameStart = name.find_first_not_of(' ');
285  size_t typeNameEnd = name.find(' ', typeNameStart);
286  std::string typeName = name.substr(typeNameStart, typeNameEnd - typeNameStart);
287  std::transform(typeName.begin(), typeName.end(), typeName.begin(), ::tolower);
288  level->AddKeepFlag(name, NoFactory::get(), MueLu::UserData);
289  if(typeName == "matrix")
290  level->Set(name, Teuchos::getValue<RCP<Matrix> >(it2->second), NoFactory::get());
291  else if(typeName == "multivector")
292  level->Set(name, Teuchos::getValue<RCP<MultiVector> >(it2->second), NoFactory::get());
293  else if(typeName == "map")
294  level->Set(name, Teuchos::getValue<RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > >(it2->second), NoFactory::get());
295  else if(typeName == "ordinalvector")
296  level->Set(name, Teuchos::getValue<RCP<Xpetra::Vector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node> > >(it2->second), NoFactory::get());
297  else if(typeName == "scalar")
298  level->Set(name, Teuchos::getValue<Scalar>(it2->second), NoFactory::get());
299  else if(typeName == "double")
300  level->Set(name, Teuchos::getValue<double>(it2->second), NoFactory::get());
301  else if(typeName == "complex")
302  level->Set(name, Teuchos::getValue<std::complex<double> >(it2->second), NoFactory::get());
303  else if(typeName == "int")
304  level->Set(name, Teuchos::getValue<int>(it2->second), NoFactory::get());
305  else if(typeName == "string")
306  level->Set(name, Teuchos::getValue<std::string>(it2->second), NoFactory::get());
307  }
308  #endif
309  }
310  }
311  }
312  }
313 
314 } // namespace MueLu
315 
316 #define MUELU_HIERARCHY_HELPERS_SHORT
317 #endif // MUELU_HIERARCHYHELPERS_DEF_HPP
Generic Smoother Factory for generating the smoothers of the MG hierarchy.
TopRAPFactory(RCP< const FactoryManagerBase > parentFactoryManager)
void Build(Level &level) const
Build an object with this factory.
This class specifies the default factory that should generate some data on a Level if the data does n...
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.
static void AddNonSerializableDataToHierarchy(HierarchyManager &HM, Hierarchy &H, const ParameterList &nonSerialList)
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
void AddNewLevel()
Add a new level at the end of the hierarchy.
User data are always kept. This flag is set automatically when Level::Set("data", data) is used...
Namespace for MueLu classes and methods.
bool IsParamMuemexVariable(const std::string &name)
static const NoFactory * get()
RCP< const FactoryBase > PFact_
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
Keep data only for this run. Used to keep data useful for Hierarchy::Iterate(). Data will be deleted ...
void RemoveKeepFlag(const std::string &ename, const FactoryBase *factory, KeepType keep=MueLu::All)
RCP< FactoryManagerBase > GetFactoryManager(int levelID) const
void AddKeepFlag(const std::string &ename, const FactoryBase *factory=NoFactory::get(), KeepType keep=MueLu::Keep)
RCP< const FactoryBase > RFact_
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
RCP< const FactoryBase > AcFact_
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
void GetSmootherPrototypes(RCP< SmootherPrototype > &preSmootherPrototype, RCP< SmootherPrototype > &postSmootherPrototype) const
Get smoother prototypes.
Exception throws to report errors in the internal logical of the program.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
void DeclareInput(Level &level) const
Specifies the data that this class needs, and the factories that generate that data.
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
TopSmootherFactory(RCP< const FactoryManagerBase > parentFactoryManager, const std::string &varName)
Exception throws to report invalid user entry.
static const RCP< const NoFactory > getRCP()
Static Get() functions.
const RCP< const FactoryBase > GetFactory(const std::string &varName) const
Default implementation of FactoryAcceptor::GetFactory()