46#ifndef MUELU_RAPSHIFTFACTORY_DEF_HPP
47#define MUELU_RAPSHIFTFACTORY_DEF_HPP
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>
61#include "MueLu_PerfUtils.hpp"
62#include "MueLu_Utilities.hpp"
67 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
69 : implicitTranspose_(false) { }
73 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
75 RCP<ParameterList> validParamList = rcp(
new ParameterList());
77#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
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");
95 validParamList->set<
bool > (
"CheckMainDiagonal",
false,
"Check main diagonal for zeros");
96 validParamList->set<
bool > (
"RepairMainDiagonal",
false,
"Repair zeros on main diagonal");
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");
103 ParameterList norecurse;
104 norecurse.disableRecursiveValidation();
105 validParamList->set<ParameterList> (
"matrixmatrix: kernel params", norecurse,
"MatrixMatrix kernel parameters");
107 return validParamList;
111 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
113 const Teuchos::ParameterList& pL = GetParameterList();
115 bool use_mdiag =
false;
116 if(pL.isParameter(
"rap: shift diagonal M"))
117 use_mdiag = pL.get<
bool>(
"rap: shift diagonal M");
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;
126 if (implicitTranspose_ ==
false) {
127 Input(coarseLevel,
"R");
130 if(!use_low_storage) Input(fineLevel,
"K");
131 else Input(fineLevel,
"A");
132 Input(coarseLevel,
"P");
134 if(!use_mdiag) Input(fineLevel,
"M");
135 else Input(fineLevel,
"Mdiag");
138 if(pL.isParameter(
"rap: cfl array") && pL.get<Teuchos::Array<double> >(
"rap: cfl array").size() > 0) {
145 "deltaT was not provided by the user on level0!");
153 "cfl was not provided by the user on level0!");
157 Input(fineLevel,
"cfl-based shift array");
162 for(std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it!=transferFacts_.end(); ++it) {
163 (*it)->CallDeclareInput(coarseLevel);
167 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
171 const Teuchos::ParameterList& pL = GetParameterList();
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");
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;
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 ) {
188 doubleShifts = pL.get<Teuchos::Array<double> >(
"rap: shift array")();
190 if(pL.isParameter(
"rap: cfl array") && pL.get<Teuchos::Array<double> >(
"rap: cfl array").size() > 0) {
192 Teuchos::ArrayView<const double> CFLs = pL.get<Teuchos::Array<double> >(
"rap: cfl array")();
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());
202 for(
int i=1; i<(int)CFLs.size(); i++)
203 myCFLs[i] = (CFLs[i]> cfl) ? cfl : CFLs[i];
206 std::ostringstream ofs;
207 ofs<<
"RAPShiftFactory: CFL schedule = ";
208 for(
int i=0; i<(int)CFLs.size(); i++)
212 GetOStream(
Statistics0)<<
"RAPShiftFactory: Timestep at CFL=1 is "<< ts_at_cfl1 <<
" " <<std::endl;
215 for(
int i=0; i<(int)myshifts.size(); i++)
216 myshifts[i] = 1.0 / (ts_at_cfl1*myCFLs[i]);
217 doubleShifts = myshifts();
220 std::ostringstream ofs;
221 ofs<<
"RAPShiftFactory: shift schedule = ";
222 for(
int i=0; i<(int)doubleShifts.size(); i++)
223 ofs<<
" "<<doubleShifts[i];
226 Set(coarseLevel,
"cfl-based shift array",myshifts);
229 myshifts = Get<Teuchos::ArrayRCP<double> > (fineLevel,
"cfl-based shift array");
230 doubleShifts = myshifts();
231 Set(coarseLevel,
"cfl-based shift array",myshifts);
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");
247 RCP<Matrix> P = Get< RCP<Matrix> >(coarseLevel,
"P");
261 KP = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*K,
false, *P,
false, KP, GetOStream(
Statistics2));
263 MP = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*M,
false, *P,
false, MP, GetOStream(
Statistics2));
266 MP = Xpetra::MatrixFactory2<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildCopy(P);
267 MP->leftScale(*Mdiag);
270 Set(coarseLevel,
"AP Pattern", KP);
273 bool doOptimizedStorage =
true;
275 RCP<Matrix> Ac, Kc, Mc;
281 bool doFillComplete=
true;
282 if (implicitTranspose_) {
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);
288 RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel,
"R");
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);
298 int level = coarseLevel.GetLevelID();
299 Scalar shift = Teuchos::ScalarTraits<Scalar>::zero();
300 if(!use_low_storage) {
302 if(level < (
int)shifts_.size()) shift = shifts_[level];
303 else shift = Teuchos::as<Scalar>(pL.get<
double>(
"rap: shift"));
307 if(level < (
int)shifts_.size()) {
308 if(level==1) shift = shifts_[level];
310 Scalar prod1 = Teuchos::ScalarTraits<Scalar>::one();
311 for(
int i=1; i < level-1; i++) {
314 shift = (prod1 * shifts_[level] - prod1);
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];
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);
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));
332 GetOStream(
Runtime0) <<
"RAPShiftFactory: Using shift " << shift << std::endl;
338 Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixAdd(*Kc,
false, Teuchos::ScalarTraits<Scalar>::one(), *Mc,
false, shift, Ac, GetOStream(
Statistics2));
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));
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));
352 RCP<ParameterList> params = rcp(
new ParameterList());;
353 params->set(
"printLoadBalancingInfo",
true);
356 Set(coarseLevel,
"A", Ac);
359 Set(coarseLevel,
"K", Kc);
362 Set(coarseLevel,
"M", Mc);
367 RCP<Vector> Mcv = Xpetra::VectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Mc->getRowMap(),
false);
368 Mc->getLocalDiagCopy(*Mcv);
369 Set(coarseLevel,
"Mdiag", Mcv);
375 if (transferFacts_.begin() != transferFacts_.end()) {
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);
390 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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);
399#define MUELU_RAPSHIFTFACTORY_SHORT
#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.
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.
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.