MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_EminPFactory_def.hpp
Go to the documentation of this file.
1 #ifndef MUELU_EMINPFACTORY_DEF_HPP
2 #define MUELU_EMINPFACTORY_DEF_HPP
3 
4 #include <Xpetra_Matrix.hpp>
5 #include <Xpetra_StridedMapFactory.hpp>
6 
8 
9 #include "MueLu_CGSolver.hpp"
10 #include "MueLu_Constraint.hpp"
11 #include "MueLu_GMRESSolver.hpp"
12 #include "MueLu_MasterList.hpp"
13 #include "MueLu_Monitor.hpp"
14 #include "MueLu_PerfUtils.hpp"
15 #include "MueLu_SolverBase.hpp"
16 #include "MueLu_SteepestDescentSolver.hpp"
17 
18 namespace MueLu {
19 
20 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
22  RCP<ParameterList> validParamList = rcp(new ParameterList());
23 
24 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
25  SET_VALID_ENTRY("emin: num iterations");
26  SET_VALID_ENTRY("emin: num reuse iterations");
27  SET_VALID_ENTRY("emin: iterative method");
28  {
30  validParamList->getEntry("emin: iterative method").setValidator(rcp(new validatorType(Teuchos::tuple<std::string>("cg", "sd", "gmres"), "emin: iterative method")));
31  }
32 #undef SET_VALID_ENTRY
33 
34  validParamList->set<RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory for the matrix A used during internal iterations");
35  validParamList->set<RCP<const FactoryBase> >("P", Teuchos::null, "Generating factory for the initial guess");
36  validParamList->set<RCP<const FactoryBase> >("Constraint", Teuchos::null, "Generating factory for constraints");
37 
38  validParamList->set<RCP<Matrix> >("P0", Teuchos::null, "Initial guess at P");
39  validParamList->set<bool>("Keep P0", false, "Keep an initial P0 (for reuse)");
40 
41  validParamList->set<RCP<Constraint> >("Constraint0", Teuchos::null, "Initial Constraint");
42  validParamList->set<bool>("Keep Constraint0", false, "Keep an initial Constraint (for reuse)");
43 
44  return validParamList;
45 }
46 
47 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
49  Input(fineLevel, "A");
50 
51  static bool isAvailableP0 = false;
52  static bool isAvailableConstraint0 = false;
53 
54  // Here is a tricky little piece of code
55  // We don't want to request (aka call Input) when we reuse and P0 is available
56  // However, we cannot run something simple like this:
57  // if (!coarseLevel.IsAvailable("P0", this))
58  // Input(coarseLevel, "P");
59  // The reason is that it works fine during the request stage, but fails
60  // in the release stage as we _construct_ P0 during Build process. Therefore,
61  // we need to understand whether we are in Request or Release mode
62  // NOTE: This is a very unique situation, please try not to propagate the
63  // mode check any further
64 
65  if (coarseLevel.GetRequestMode() == Level::REQUEST) {
66  isAvailableP0 = coarseLevel.IsAvailable("P0", this);
67  isAvailableConstraint0 = coarseLevel.IsAvailable("Constraint0", this);
68  }
69 
70  if (isAvailableP0 == false)
71  Input(coarseLevel, "P");
72 
73  if (isAvailableConstraint0 == false)
74  Input(coarseLevel, "Constraint");
75 }
76 
77 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
79  BuildP(fineLevel, coarseLevel);
80 }
81 
82 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
84  FactoryMonitor m(*this, "Prolongator minimization", coarseLevel);
85 
86  const ParameterList& pL = GetParameterList();
87 
88  // Get the matrix
89  RCP<Matrix> A = Get<RCP<Matrix> >(fineLevel, "A");
90 
91  if (restrictionMode_) {
92  SubFactoryMonitor m2(*this, "Transpose A", coarseLevel);
93 
94  A = Utilities::Transpose(*A, true);
95  }
96 
97  // Get/make initial guess
98  RCP<Matrix> P0;
99  int numIts;
100  if (coarseLevel.IsAvailable("P0", this)) {
101  // Reuse data
102  P0 = coarseLevel.Get<RCP<Matrix> >("P0", this);
103  numIts = pL.get<int>("emin: num reuse iterations");
104  GetOStream(Runtime0) << "Reusing P0" << std::endl;
105 
106  } else {
107  // Construct data
108  P0 = Get<RCP<Matrix> >(coarseLevel, "P");
109  numIts = pL.get<int>("emin: num iterations");
110  }
111  // NOTE: the main assumption here that P0 satisfies both constraints:
112  // - nonzero pattern
113  // - nullspace preservation
114 
115  // Get/make constraint operator
116  RCP<Constraint> X;
117  if (coarseLevel.IsAvailable("Constraint0", this)) {
118  // Reuse data
119  X = coarseLevel.Get<RCP<Constraint> >("Constraint0", this);
120  GetOStream(Runtime0) << "Reusing Constraint0" << std::endl;
121 
122  } else {
123  // Construct data
124  X = Get<RCP<Constraint> >(coarseLevel, "Constraint");
125  }
126  GetOStream(Runtime0) << "Number of emin iterations = " << numIts << std::endl;
127 
128  std::string solverType = pL.get<std::string>("emin: iterative method");
129  RCP<SolverBase> solver;
130  if (solverType == "cg")
131  solver = rcp(new CGSolver(numIts));
132  else if (solverType == "sd")
133  solver = rcp(new SteepestDescentSolver(numIts));
134  else if (solverType == "gmres")
135  solver = rcp(new GMRESSolver(numIts));
136 
137  RCP<Matrix> P;
138  solver->Iterate(*A, *X, *P0, P);
139 
140  // NOTE: EXPERIMENTAL and FRAGILE
141  if (!P->IsView("stridedMaps")) {
142  if (A->IsView("stridedMaps") == true) {
143  GetOStream(Runtime1) << "Using A to fillComplete P" << std::endl;
144 
145  // FIXME: X->GetPattern() actually returns a CrsGraph.
146  // CrsGraph has no knowledge of Xpetra's sup/Matrix views. As such,
147  // it has no idea about strided maps. We create one, which is
148  // most likely incorrect for many use cases.
149  std::vector<size_t> stridingInfo(1, 1);
150  RCP<const StridedMap> dMap = StridedMapFactory::Build(X->GetPattern()->getDomainMap(), stridingInfo);
151 
152  P->CreateView("stridedMaps", A->getRowMap("stridedMaps"), dMap);
153 
154  } else {
155  P->CreateView("stridedMaps", P->getRangeMap(), P->getDomainMap());
156  }
157  }
158 
159  // Level Set
160  if (!restrictionMode_) {
161  // The factory is in prolongation mode
162  Set(coarseLevel, "P", P);
163 
164  if (pL.get<bool>("Keep P0")) {
165  // NOTE: we must do Keep _before_ set as the Needs class only sets if
166  // a) data has been requested (which is not the case here), or
167  // b) data has some keep flag
168  coarseLevel.Keep("P0", this);
169  Set(coarseLevel, "P0", P);
170  }
171  if (pL.get<bool>("Keep Constraint0")) {
172  // NOTE: we must do Keep _before_ set as the Needs class only sets if
173  // a) data has been requested (which is not the case here), or
174  // b) data has some keep flag
175  coarseLevel.Keep("Constraint0", this);
176  Set(coarseLevel, "Constraint0", X);
177  }
178 
179  if (IsPrint(Statistics2)) {
180  RCP<ParameterList> params = rcp(new ParameterList());
181  params->set("printLoadBalancingInfo", true);
182  params->set("printCommInfo", true);
183  GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*P, "P", params);
184  }
185 
186  } else {
187  // The factory is in restriction mode
188  RCP<Matrix> R;
189  {
190  SubFactoryMonitor m2(*this, "Transpose P", coarseLevel);
191 
192  R = Utilities::Transpose(*P, true);
193  }
194 
195  Set(coarseLevel, "R", R);
196 
197  if (IsPrint(Statistics2)) {
198  RCP<ParameterList> params = rcp(new ParameterList());
199  params->set("printLoadBalancingInfo", true);
200  params->set("printCommInfo", true);
201  GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*R, "R", params);
202  }
203  }
204 }
205 
206 } // namespace MueLu
207 
208 #endif // MUELU_EMINPFACTORY_DEF_HPP
void Keep(const std::string &ename, const FactoryBase *factory)
Request to keep variable &#39;ename&#39; generated by &#39;factory&#39; after the setup phase.
#define SET_VALID_ENTRY(name)
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-&gt;Get&lt; RCP&lt;Matrix&gt; &gt;(&quot;A&quot;, factory) if factory == NULL =&gt; use default factory.
void Build(Level &fineLevel, Level &coarseLevel) const
Build method.
void setValidator(RCP< const ParameterEntryValidator > const &validator)
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Implements conjugate gradient algorithm for energy-minimization.
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
Timer to be used in factories. Similar to Monitor but with additional timers.
One-liner description of what is happening.
RCP< const CrsGraph > GetPattern() const
Print even more statistics.
RequestMode GetRequestMode() const
Implements conjugate gradient algorithm for energy-minimization.
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
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Transpose(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, bool optimizeTranspose=false, const std::string &label=std::string(), const Teuchos::RCP< Teuchos::ParameterList > &params=Teuchos::null)
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Implements steepest descent algorithm for energy-minimization.
Description of what is happening (more verbose)
ParameterEntry & getEntry(const std::string &name)
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need&#39;s value has been saved.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.