46 #ifndef MUELU_RAPSHIFTFACTORY_DEF_HPP
47 #define MUELU_RAPSHIFTFACTORY_DEF_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>
77 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
85 #undef SET_VALID_ENTRY
87 validParamList->
set<
RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A used during the prolongator smoothing process");
88 validParamList->
set<
RCP<const FactoryBase> >(
"M", Teuchos::null,
"Generating factory of the matrix M used during the non-Galerkin RAP");
89 validParamList->
set<
RCP<const FactoryBase> >(
"Mdiag", Teuchos::null,
"Generating factory of the matrix Mdiag used during the non-Galerkin RAP");
90 validParamList->
set<
RCP<const FactoryBase> >(
"K", Teuchos::null,
"Generating factory of the matrix K used during the non-Galerkin RAP");
94 validParamList->
set<
bool > (
"CheckMainDiagonal",
false,
"Check main diagonal for zeros");
95 validParamList->
set<
bool > (
"RepairMainDiagonal",
false,
"Repair zeros on main diagonal");
99 validParamList->
set<
RCP<const FactoryBase> > (
"cfl-based shift array", Teuchos::null,
"MueLu-generated shift array for CFL-based shifting");
103 norecurse.disableRecursiveValidation();
104 validParamList->
set<
ParameterList> (
"matrixmatrix: kernel params", norecurse,
"MatrixMatrix kernel parameters");
106 return validParamList;
110 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
114 bool use_mdiag =
false;
116 use_mdiag = pL.
get<
bool>(
"rap: shift diagonal M");
119 bool use_low_storage =
false;
121 use_low_storage = pL.
get<
bool>(
"rap: shift low storage");
122 use_mdiag = use_low_storage ?
true : use_mdiag;
125 if (implicitTranspose_ ==
false) {
126 Input(coarseLevel,
"R");
129 if(!use_low_storage) Input(fineLevel,
"K");
130 else Input(fineLevel,
"A");
131 Input(coarseLevel,
"P");
133 if(!use_mdiag) Input(fineLevel,
"M");
134 else Input(fineLevel,
"Mdiag");
144 "deltaT was not provided by the user on level0!");
152 "cfl was not provided by the user on level0!");
156 Input(fineLevel,
"cfl-based shift array");
161 for(std::vector<
RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it!=transferFacts_.end(); ++it) {
162 (*it)->CallDeclareInput(coarseLevel);
166 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
172 bool M_is_diagonal =
false;
174 M_is_diagonal = pL.
get<
bool>(
"rap: shift diagonal M");
177 bool use_low_storage =
false;
179 use_low_storage = pL.
get<
bool>(
"rap: shift low storage");
180 M_is_diagonal = use_low_storage ?
true : M_is_diagonal;
193 double dt = Get<double>(fineLevel,
"deltaT");
194 double cfl = Get<double>(fineLevel,
"cfl");
195 double ts_at_cfl1 = dt / cfl;
201 for(
int i=1; i<(int)CFLs.
size(); i++)
202 myCFLs[i] = (CFLs[i]> cfl) ? cfl : CFLs[i];
205 std::ostringstream ofs;
206 ofs<<
"RAPShiftFactory: CFL schedule = ";
207 for(
int i=0; i<(int)CFLs.
size(); i++)
211 GetOStream(
Statistics0)<<
"RAPShiftFactory: Timestep at CFL=1 is "<< ts_at_cfl1 <<
" " <<std::endl;
214 for(
int i=0; i<(int)myshifts.
size(); i++)
215 myshifts[i] = 1.0 / (ts_at_cfl1*myCFLs[i]);
216 doubleShifts = myshifts();
219 std::ostringstream ofs;
220 ofs<<
"RAPShiftFactory: shift schedule = ";
221 for(
int i=0; i<(int)doubleShifts.
size(); i++)
222 ofs<<
" "<<doubleShifts[i];
225 Set(coarseLevel,
"cfl-based shift array",myshifts);
229 doubleShifts = myshifts();
230 Set(coarseLevel,
"cfl-based shift array",myshifts);
241 if(use_low_storage) K = Get< RCP<Matrix> >(fineLevel,
"A");
242 else K = Get< RCP<Matrix> >(fineLevel,
"K");
243 if(!M_is_diagonal) M = Get< RCP<Matrix> >(fineLevel,
"M");
244 else Mdiag = Get< RCP<Vector> >(fineLevel,
"Mdiag");
246 RCP<Matrix> P = Get< RCP<Matrix> >(coarseLevel,
"P");
266 MP->leftScale(*Mdiag);
269 Set(coarseLevel,
"AP Pattern", KP);
272 bool doOptimizedStorage =
true;
280 bool doFillComplete=
true;
281 if (implicitTranspose_) {
283 Kc =
Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*P,
true, *KP,
false, Kc, GetOStream(
Statistics2), doFillComplete, doOptimizedStorage);
284 Mc =
Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*P,
true, *MP,
false, Mc, GetOStream(
Statistics2), doFillComplete, doOptimizedStorage);
287 RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel,
"R");
289 Kc =
Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*R,
false, *KP,
false, Kc, GetOStream(
Statistics2), doFillComplete, doOptimizedStorage);
290 Mc =
Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*R,
false, *MP,
false, Mc, GetOStream(
Statistics2), doFillComplete, doOptimizedStorage);
297 int level = coarseLevel.GetLevelID();
299 if(!use_low_storage) {
301 if(level < (
int)shifts_.size()) shift = shifts_[level];
302 else shift = Teuchos::as<Scalar>(pL.
get<
double>(
"rap: shift"));
306 if(level < (
int)shifts_.size()) {
307 if(level==1) shift = shifts_[level];
310 for(
int i=1; i < level-1; i++) {
313 shift = (prod1 * shifts_[level] - prod1);
316 else if(doubleShifts.
size() != 0) {
317 double d_shift = 0.0;
318 if(level < doubleShifts.
size())
319 d_shift = doubleShifts[level] - doubleShifts[level-1];
322 GetOStream(
Warnings1) <<
"WARNING: RAPShiftFactory has detected a negative shift... This implies a less stable coarse grid."<<std::endl;
323 shift = Teuchos::as<Scalar>(d_shift);
326 double base_shift = pL.
get<
double>(
"rap: shift");
327 if(level == 1) shift = Teuchos::as<Scalar>(base_shift);
328 else shift = Teuchos::as<Scalar>(pow(base_shift,level) - pow(base_shift,level-1));
331 GetOStream(
Runtime0) <<
"RAPShiftFactory: Using shift " << shift << std::endl;
337 Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixAdd(*Kc,
false,
Teuchos::ScalarTraits<Scalar>::one(), *Mc,
false, shift, Ac, GetOStream(
Statistics2));
341 bool repairZeroDiagonals = pL.
get<
bool>(
"RepairMainDiagonal") || pL.
get<
bool>(
"rap: fix zero diagonals");
342 bool checkAc = pL.
get<
bool>(
"CheckMainDiagonal")|| pL.
get<
bool>(
"rap: fix zero diagonals"); ;
343 if (checkAc || repairZeroDiagonals)
347 params->
set(
"printLoadBalancingInfo",
true);
350 Set(coarseLevel,
"A", Ac);
353 Set(coarseLevel,
"K", Kc);
356 Set(coarseLevel,
"M", Mc);
362 Mc->getLocalDiagCopy(*Mcv);
363 Set(coarseLevel,
"Mdiag", Mcv);
369 if (transferFacts_.begin() != transferFacts_.end()) {
373 for (std::vector<
RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it) {
375 GetOStream(
Runtime0) <<
"RAPShiftFactory: call transfer factory: " << fac->
description() << std::endl;
384 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
387 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)");
388 transferFacts_.push_back(factory);
393 #define MUELU_RAPSHIFTFACTORY_SHORT
394 #endif // MUELU_RAPSHIFTFACTORY_DEF_HPP
Exception indicating invalid cast attempted.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
virtual void CallBuild(Level &requestedLevel) const =0
void Release(const FactoryBase &factory)
Decrement the storage counter for all the inputs of a factory.
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
One-liner description of what is happening.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
static void CheckRepairMainDiagonal(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> &Ac, bool const &repairZeroDiagonals, Teuchos::FancyOStream &fos, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType threshold=Teuchos::ScalarTraits< typename Teuchos::ScalarTraits< Scalar >::magnitudeType >::zero())
static const NoFactory * get()
Print even more statistics.
bool isParameter(const std::string &name) const
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Print statistics that do not involve significant additional computation.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > BuildCopy(const RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > A)
void resize(size_type new_size, const value_type &x=value_type())
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
static RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
static void TwoMatrixAdd(const Matrix &A, bool transposeA, SC alpha, Matrix &B, SC beta)
void AddTransferFactory(const RCP< const FactoryBase > &factory)
Add transfer factory in the end of list of transfer factories in RepartitionAcFactory.
static void Multiply(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, Matrix &C, bool call_FillComplete_on_result=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > ¶ms=null)
int GetLevelID() const
Return level number.
Exception throws to report errors in the internal logical of the program.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
virtual std::string description() const
Return a simple one-line description of this object.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
#define SET_VALID_ENTRY(name)