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