MueLu Version of the Day
Loading...
Searching...
No Matches
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_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
21namespace 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 {
32 typedef Teuchos::StringToIntegralParameterEntryValidator<int> validatorType;
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
#define SET_VALID_ENTRY(name)
Implements conjugate gradient algorithm for energy-minimization.
void Build(Level &fineLevel, Level &coarseLevel) const
Build method.
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Timer to be used in factories. Similar to Monitor but with additional timers.
Implements conjugate gradient algorithm for energy-minimization.
Class that holds all level-specific information.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
RequestMode GetRequestMode() const
void Keep(const std::string &ename, const FactoryBase *factory)
Request to keep variable 'ename' generated by 'factory' after the setup phase.
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
Implements steepest descent algorithm for energy-minimization.
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)
Namespace for MueLu classes and methods.
@ Statistics2
Print even more statistics.
@ Runtime0
One-liner description of what is happening.
@ Runtime1
Description of what is happening (more verbose)