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