MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_RAPShiftFactory_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_RAPSHIFTFACTORY_DEF_HPP
47#define MUELU_RAPSHIFTFACTORY_DEF_HPP
48
49#include <sstream>
50
51#include <Xpetra_Matrix.hpp>
52#include <Xpetra_MatrixMatrix.hpp>
53#include <Xpetra_MatrixUtils.hpp>
54#include <Xpetra_Vector.hpp>
55#include <Xpetra_VectorFactory.hpp>
56
57
59#include "MueLu_MasterList.hpp"
60#include "MueLu_Monitor.hpp"
61#include "MueLu_PerfUtils.hpp"
62#include "MueLu_Utilities.hpp"
63
64namespace MueLu {
65
66 /*********************************************************************************************************/
67 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
69 : implicitTranspose_(false) { }
70
71
72 /*********************************************************************************************************/
73 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
75 RCP<ParameterList> validParamList = rcp(new ParameterList());
76
77#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
78 SET_VALID_ENTRY("transpose: use implicit");
79 SET_VALID_ENTRY("rap: fix zero diagonals");
80 SET_VALID_ENTRY("rap: shift");
81 SET_VALID_ENTRY("rap: shift array");
82 SET_VALID_ENTRY("rap: cfl array");
83 SET_VALID_ENTRY("rap: shift diagonal M");
84 SET_VALID_ENTRY("rap: shift low storage");
85 SET_VALID_ENTRY("rap: relative diagonal floor");
86#undef SET_VALID_ENTRY
87
88 validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A used during the prolongator smoothing process");
89 validParamList->set< RCP<const FactoryBase> >("M", Teuchos::null, "Generating factory of the matrix M used during the non-Galerkin RAP");
90 validParamList->set< RCP<const FactoryBase> >("Mdiag", Teuchos::null, "Generating factory of the matrix Mdiag used during the non-Galerkin RAP");
91 validParamList->set< RCP<const FactoryBase> >("K", Teuchos::null, "Generating factory of the matrix K used during the non-Galerkin RAP");
92 validParamList->set< RCP<const FactoryBase> >("P", Teuchos::null, "Prolongator factory");
93 validParamList->set< RCP<const FactoryBase> >("R", Teuchos::null, "Restrictor factory");
94
95 validParamList->set< bool > ("CheckMainDiagonal", false, "Check main diagonal for zeros");
96 validParamList->set< bool > ("RepairMainDiagonal", false, "Repair zeros on main diagonal");
97
98 validParamList->set<RCP<const FactoryBase> > ("deltaT", Teuchos::null, "user deltaT");
99 validParamList->set<RCP<const FactoryBase> > ("cfl", Teuchos::null, "user cfl");
100 validParamList->set<RCP<const FactoryBase> > ("cfl-based shift array", Teuchos::null, "MueLu-generated shift array for CFL-based shifting");
101
102 // Make sure we don't recursively validate options for the matrixmatrix kernels
103 ParameterList norecurse;
104 norecurse.disableRecursiveValidation();
105 validParamList->set<ParameterList> ("matrixmatrix: kernel params", norecurse, "MatrixMatrix kernel parameters");
106
107 return validParamList;
108 }
109
110 /*********************************************************************************************************/
111 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
113 const Teuchos::ParameterList& pL = GetParameterList();
114
115 bool use_mdiag = false;
116 if(pL.isParameter("rap: shift diagonal M"))
117 use_mdiag = pL.get<bool>("rap: shift diagonal M");
118
119 // The low storage version requires mdiag
120 bool use_low_storage = false;
121 if(pL.isParameter("rap: shift low storage")) {
122 use_low_storage = pL.get<bool>("rap: shift low storage");
123 use_mdiag = use_low_storage ? true : use_mdiag;
124 }
125
126 if (implicitTranspose_ == false) {
127 Input(coarseLevel, "R");
128 }
129
130 if(!use_low_storage) Input(fineLevel, "K");
131 else Input(fineLevel, "A");
132 Input(coarseLevel, "P");
133
134 if(!use_mdiag) Input(fineLevel, "M");
135 else Input(fineLevel, "Mdiag");
136
137 // CFL array stuff
138 if(pL.isParameter("rap: cfl array") && pL.get<Teuchos::Array<double> >("rap: cfl array").size() > 0) {
139 if(fineLevel.GetLevelID() == 0) {
140 if(fineLevel.IsAvailable("deltaT", NoFactory::get())) {
141 fineLevel.DeclareInput("deltaT", NoFactory::get(), this);
142 } else {
143 TEUCHOS_TEST_FOR_EXCEPTION(fineLevel.IsAvailable("fine deltaT", NoFactory::get()),
145 "deltaT was not provided by the user on level0!");
146 }
147
148 if(fineLevel.IsAvailable("cfl", NoFactory::get())) {
149 fineLevel.DeclareInput("cfl", NoFactory::get(), this);
150 } else {
151 TEUCHOS_TEST_FOR_EXCEPTION(fineLevel.IsAvailable("fine cfl", NoFactory::get()),
153 "cfl was not provided by the user on level0!");
154 }
155 }
156 else {
157 Input(fineLevel,"cfl-based shift array");
158 }
159 }
160
161 // call DeclareInput of all user-given transfer factories
162 for(std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it!=transferFacts_.end(); ++it) {
163 (*it)->CallDeclareInput(coarseLevel);
164 }
165 }
166
167 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
168 void RAPShiftFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level &fineLevel, Level &coarseLevel) const { // FIXME make fineLevel const
169 {
170 FactoryMonitor m(*this, "Computing Ac", coarseLevel);
171 const Teuchos::ParameterList& pL = GetParameterList();
172
173 bool M_is_diagonal = false;
174 if(pL.isParameter("rap: shift diagonal M"))
175 M_is_diagonal = pL.get<bool>("rap: shift diagonal M");
176
177 // The low storage version requires mdiag
178 bool use_low_storage = false;
179 if(pL.isParameter("rap: shift low storage")) {
180 use_low_storage = pL.get<bool>("rap: shift low storage");
181 M_is_diagonal = use_low_storage ? true : M_is_diagonal;
182 }
183
184 Teuchos::ArrayView<const double> doubleShifts;
185 Teuchos::ArrayRCP<double> myshifts;
186 if(pL.isParameter("rap: shift array") && pL.get<Teuchos::Array<double> >("rap: shift array").size() > 0 ) {
187 // Do we have an array of shifts? If so, we set doubleShifts_
188 doubleShifts = pL.get<Teuchos::Array<double> >("rap: shift array")();
189 }
190 if(pL.isParameter("rap: cfl array") && pL.get<Teuchos::Array<double> >("rap: cfl array").size() > 0) {
191 // Do we have an array of CFLs? If so, we calculated the shifts from them.
192 Teuchos::ArrayView<const double> CFLs = pL.get<Teuchos::Array<double> >("rap: cfl array")();
193 if(fineLevel.GetLevelID() == 0) {
194 double dt = Get<double>(fineLevel,"deltaT");
195 double cfl = Get<double>(fineLevel,"cfl");
196 double ts_at_cfl1 = dt / cfl;
197 myshifts.resize(CFLs.size());
198 Teuchos::Array<double> myCFLs(CFLs.size());
199 myCFLs[0] = cfl;
200
201 // Never make the CFL bigger
202 for(int i=1; i<(int)CFLs.size(); i++)
203 myCFLs[i] = (CFLs[i]> cfl) ? cfl : CFLs[i];
204
205 {
206 std::ostringstream ofs;
207 ofs<<"RAPShiftFactory: CFL schedule = ";
208 for(int i=0; i<(int)CFLs.size(); i++)
209 ofs<<" "<<myCFLs[i];
210 GetOStream(Statistics0) <<ofs.str() << std::endl;
211 }
212 GetOStream(Statistics0)<< "RAPShiftFactory: Timestep at CFL=1 is "<< ts_at_cfl1 << " " <<std::endl;
213
214 // The shift array needs to be 1/dt
215 for(int i=0; i<(int)myshifts.size(); i++)
216 myshifts[i] = 1.0 / (ts_at_cfl1*myCFLs[i]);
217 doubleShifts = myshifts();
218
219 {
220 std::ostringstream ofs;
221 ofs<<"RAPShiftFactory: shift schedule = ";
222 for(int i=0; i<(int)doubleShifts.size(); i++)
223 ofs<<" "<<doubleShifts[i];
224 GetOStream(Statistics0) <<ofs.str() << std::endl;
225 }
226 Set(coarseLevel,"cfl-based shift array",myshifts);
227 }
228 else {
229 myshifts = Get<Teuchos::ArrayRCP<double> > (fineLevel,"cfl-based shift array");
230 doubleShifts = myshifts();
231 Set(coarseLevel,"cfl-based shift array",myshifts);
232 // NOTE: If we're not on level zero, then we should have a shift array
233 }
234 }
235
236 // Inputs: K, M, P
237 // Note: In the low-storage case we do not keep a separate "K", we just use A
238 RCP<Matrix> K;
239 RCP<Matrix> M;
240 RCP<Vector> Mdiag;
241
242 if(use_low_storage) K = Get< RCP<Matrix> >(fineLevel, "A");
243 else K = Get< RCP<Matrix> >(fineLevel, "K");
244 if(!M_is_diagonal) M = Get< RCP<Matrix> >(fineLevel, "M");
245 else Mdiag = Get< RCP<Vector> >(fineLevel, "Mdiag");
246
247 RCP<Matrix> P = Get< RCP<Matrix> >(coarseLevel, "P");
248
249 // Build Kc = RKP, Mc = RMP
250 RCP<Matrix> KP, MP;
251
252 // Reuse pattern if available (multiple solve)
253 // FIXME: Old style reuse doesn't work any more
254 // if (IsAvailable(coarseLevel, "AP Pattern")) {
255 // KP = Get< RCP<Matrix> >(coarseLevel, "AP Pattern");
256 // MP = Get< RCP<Matrix> >(coarseLevel, "AP Pattern");
257 // }
258
259 {
260 SubFactoryMonitor subM(*this, "MxM: K x P", coarseLevel);
261 KP = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*K, false, *P, false, KP, GetOStream(Statistics2));
262 if(!M_is_diagonal) {
263 MP = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*M, false, *P, false, MP, GetOStream(Statistics2));
264 }
265 else {
266 MP = Xpetra::MatrixFactory2<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildCopy(P);
267 MP->leftScale(*Mdiag);
268 }
269
270 Set(coarseLevel, "AP Pattern", KP);
271 }
272
273 bool doOptimizedStorage = true;
274
275 RCP<Matrix> Ac, Kc, Mc;
276
277 // Reuse pattern if available (multiple solve)
278 // if (IsAvailable(coarseLevel, "RAP Pattern"))
279 // Ac = Get< RCP<Matrix> >(coarseLevel, "RAP Pattern");
280
281 bool doFillComplete=true;
282 if (implicitTranspose_) {
283 SubFactoryMonitor m2(*this, "MxM: P' x (KP) (implicit)", coarseLevel);
284 Kc = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*P, true, *KP, false, Kc, GetOStream(Statistics2), doFillComplete, doOptimizedStorage);
285 Mc = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*P, true, *MP, false, Mc, GetOStream(Statistics2), doFillComplete, doOptimizedStorage);
286 }
287 else {
288 RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel, "R");
289 SubFactoryMonitor m2(*this, "MxM: R x (KP) (explicit)", coarseLevel);
290 Kc = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*R, false, *KP, false, Kc, GetOStream(Statistics2), doFillComplete, doOptimizedStorage);
291 Mc = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*R, false, *MP, false, Mc, GetOStream(Statistics2), doFillComplete, doOptimizedStorage);
292 }
293
294 // Get the shift
295 // FIXME - We should really get rid of the shifts array and drive this the same way everything else works
296 // If we're using the recursive "low storage" version, we need to shift by ( \prod_{i=1}^k shift[i] - \prod_{i=1}^{k-1} shift[i]) to
297 // get the recursive relationships correct
298 int level = coarseLevel.GetLevelID();
299 Scalar shift = Teuchos::ScalarTraits<Scalar>::zero();
300 if(!use_low_storage) {
301 // High Storage version
302 if(level < (int)shifts_.size()) shift = shifts_[level];
303 else shift = Teuchos::as<Scalar>(pL.get<double>("rap: shift"));
304 }
305 else {
306 // Low Storage Version
307 if(level < (int)shifts_.size()) {
308 if(level==1) shift = shifts_[level];
309 else {
310 Scalar prod1 = Teuchos::ScalarTraits<Scalar>::one();
311 for(int i=1; i < level-1; i++) {
312 prod1 *= shifts_[i];
313 }
314 shift = (prod1 * shifts_[level] - prod1);
315 }
316 }
317 else if(doubleShifts.size() != 0) {
318 double d_shift = 0.0;
319 if(level < doubleShifts.size())
320 d_shift = doubleShifts[level] - doubleShifts[level-1];
321
322 if(d_shift < 0.0)
323 GetOStream(Warnings1) << "WARNING: RAPShiftFactory has detected a negative shift... This implies a less stable coarse grid."<<std::endl;
324 shift = Teuchos::as<Scalar>(d_shift);
325 }
326 else {
327 double base_shift = pL.get<double>("rap: shift");
328 if(level == 1) shift = Teuchos::as<Scalar>(base_shift);
329 else shift = Teuchos::as<Scalar>(pow(base_shift,level) - pow(base_shift,level-1));
330 }
331 }
332 GetOStream(Runtime0) << "RAPShiftFactory: Using shift " << shift << std::endl;
333
334
335 // recombine to get K+shift*M
336 {
337 SubFactoryMonitor m2(*this, "Add: RKP + s*RMP", coarseLevel);
338 Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixAdd(*Kc, false, Teuchos::ScalarTraits<Scalar>::one(), *Mc, false, shift, Ac, GetOStream(Statistics2));
339 Ac->fillComplete();
340 }
341
342 Teuchos::ArrayView<const double> relativeFloor = pL.get<Teuchos::Array<double> >("rap: relative diagonal floor")();
343 if(relativeFloor.size() > 0)
344 Xpetra::MatrixUtils<SC,LO,GO,NO>::RelativeDiagonalBoost(Ac, relativeFloor,GetOStream(Statistics2));
345
346
347 bool repairZeroDiagonals = pL.get<bool>("RepairMainDiagonal") || pL.get<bool>("rap: fix zero diagonals");
348 bool checkAc = pL.get<bool>("CheckMainDiagonal")|| pL.get<bool>("rap: fix zero diagonals"); ;
349 if (checkAc || repairZeroDiagonals)
350 Xpetra::MatrixUtils<SC,LO,GO,NO>::CheckRepairMainDiagonal(Ac, repairZeroDiagonals, GetOStream(Warnings1));
351
352 RCP<ParameterList> params = rcp(new ParameterList());;
353 params->set("printLoadBalancingInfo", true);
354 GetOStream(Statistics0) << PerfUtils::PrintMatrixInfo(*Ac, "Ac", params);
355
356 Set(coarseLevel, "A", Ac);
357 // We only need K in the 'high storage' mode
358 if(!use_low_storage)
359 Set(coarseLevel, "K", Kc);
360
361 if(!M_is_diagonal) {
362 Set(coarseLevel, "M", Mc);
363 }
364 else {
365 // If M is diagonal, then we only pass that part down the hierarchy
366 // NOTE: Should we be doing some kind of rowsum instead?
367 RCP<Vector> Mcv = Xpetra::VectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Mc->getRowMap(),false);
368 Mc->getLocalDiagCopy(*Mcv);
369 Set(coarseLevel, "Mdiag", Mcv);
370 }
371
372 // Set(coarseLevel, "RAP Pattern", Ac);
373 }
374
375 if (transferFacts_.begin() != transferFacts_.end()) {
376 SubFactoryMonitor m(*this, "Projections", coarseLevel);
377
378 // call Build of all user-given transfer factories
379 for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it) {
380 RCP<const FactoryBase> fac = *it;
381 GetOStream(Runtime0) << "RAPShiftFactory: call transfer factory: " << fac->description() << std::endl;
382 fac->CallBuild(coarseLevel);
383 // AP (11/11/13): I am not sure exactly why we need to call Release, but we do need it to get rid
384 // of dangling data for CoordinatesTransferFactory
385 coarseLevel.Release(*fac);
386 }
387 }
388 }
389
390 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
392 // check if it's a TwoLevelFactoryBase based transfer factory
393 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const TwoLevelFactoryBase>(factory) == Teuchos::null, Exceptions::BadCast, "MueLu::RAPShiftFactory::AddTransferFactory: Transfer factory is not derived from TwoLevelFactoryBase. This is very strange. (Note: you can remove this exception if there's a good reason for)");
394 transferFacts_.push_back(factory);
395 }
396
397} //namespace MueLu
398
399#define MUELU_RAPSHIFTFACTORY_SHORT
400#endif // MUELU_RAPSHIFTFACTORY_DEF_HPP
#define SET_VALID_ENTRY(name)
MueLu::DefaultScalar Scalar
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
void Release(const FactoryBase &factory)
Decrement the storage counter for all the inputs of a factory.
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
static const NoFactory * get()
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.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void AddTransferFactory(const RCP< const FactoryBase > &factory)
Add transfer factory in the end of list of transfer factories in RepartitionAcFactory.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Namespace for MueLu classes and methods.
@ Statistics2
Print even more statistics.
@ Runtime0
One-liner description of what is happening.
@ Warnings1
Additional warnings.
@ Statistics0
Print statistics that do not involve significant additional computation.