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)) 86 #undef SET_VALID_ENTRY 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 400 #endif // MUELU_RAPSHIFTFACTORY_DEF_HPP Exception indicating invalid cast attempted.
void Release(const FactoryBase &factory)
Decrement the storage counter for all the inputs of a factory.
Timer to be used in factories. Similar to Monitor but with additional timers.
One-liner description of what is happening.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Namespace for MueLu classes and methods.
static const NoFactory * get()
Print even more statistics.
int GetLevelID() const
Return level number.
Print statistics that do not involve significant additional computation.
MueLu::DefaultScalar Scalar
Class that holds all level-specific information.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
void AddTransferFactory(const RCP< const FactoryBase > &factory)
Add transfer factory in the end of list of transfer factories in RepartitionAcFactory.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
Exception throws to report errors in the internal logical of the program.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
#define SET_VALID_ENTRY(name)