46 #ifndef MUELU_RAPSHIFTFACTORY_DEF_HPP 47 #define MUELU_RAPSHIFTFACTORY_DEF_HPP 59 #include "MueLu_PerfUtils.hpp" 60 #include "MueLu_Utilities.hpp" 65 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
67 : implicitTranspose_(false), checkAc_(false), repairZeroDiagonals_(false) { }
71 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
75 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 78 #undef SET_VALID_ENTRY 80 validParamList->
set<
RCP<const FactoryBase> >(
"M", Teuchos::null,
"Generating factory of the matrix M used during the non-Galerkin RAP");
81 validParamList->
set<
RCP<const FactoryBase> >(
"K", Teuchos::null,
"Generating factory of the matrix K used during the non-Galerkin RAP");
85 validParamList->
set<
bool > (
"CheckMainDiagonal",
false,
"Check main diagonal for zeros");
86 validParamList->
set<
bool > (
"RepairMainDiagonal",
false,
"Repair zeros on main diagonal");
91 validParamList->
set<
ParameterList> (
"matrixmatrix: kernel params", norecurse,
"MatrixMatrix kernel parameters");
93 return validParamList;
97 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
99 if (implicitTranspose_ ==
false) {
100 Input(coarseLevel,
"R");
103 Input(fineLevel,
"K");
104 Input(fineLevel,
"M");
105 Input(coarseLevel,
"P");
108 for(std::vector<
RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it!=transferFacts_.end(); ++it) {
109 (*it)->CallDeclareInput(coarseLevel);
113 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
120 RCP<Matrix> K = Get< RCP<Matrix> >(fineLevel,
"K");
121 RCP<Matrix> M = Get< RCP<Matrix> >(fineLevel,
"M");
122 RCP<Matrix> P = Get< RCP<Matrix> >(coarseLevel,
"P");
138 Set(coarseLevel,
"AP Pattern", KP);
143 bool doOptimizedStorage = !checkAc_;
151 bool doFillComplete=
true;
152 if (implicitTranspose_) {
154 Kc =
Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*P,
true, *KP,
false, Kc, GetOStream(
Statistics2), doFillComplete, doOptimizedStorage);
155 Mc =
Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*P,
true, *MP,
false, Mc, GetOStream(
Statistics2), doFillComplete, doOptimizedStorage);
158 RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel,
"R");
160 Kc =
Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*R,
false, *KP,
false, Kc, GetOStream(
Statistics2), doFillComplete, doOptimizedStorage);
161 Mc =
Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*R,
false, *MP,
false, Mc, GetOStream(
Statistics2), doFillComplete, doOptimizedStorage);
166 int level = coarseLevel.GetLevelID();
168 if(level < (
int)shifts_.size())
169 shift = shifts_[level];
171 shift = Teuchos::as<Scalar>(pL.
get<
double>(
"rap: shift"));
176 Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixAdd(*Kc,
false, (Scalar) 1.0, *Mc,
false, shift, Ac, GetOStream(
Statistics2));
181 CheckMainDiagonal(Ac);
184 params->
set(
"printLoadBalancingInfo",
true);
187 Set(coarseLevel,
"A", Ac);
188 Set(coarseLevel,
"K", Kc);
189 Set(coarseLevel,
"M", Mc);
193 if (transferFacts_.begin() != transferFacts_.end()) {
197 for (std::vector<
RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it) {
199 GetOStream(
Runtime0) <<
"RAPShiftFactory: call transfer factory: " << fac->
description() << std::endl;
208 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
212 RCP<Vector> diagVec = VectorFactory::Build(Ac->getRowMap());
213 Ac->getLocalDiagCopy(*diagVec);
215 for (
size_t r=0; r<Ac->getRowMap()->getNodeNumElements(); r++) {
216 if(diagVal[r]==0.0) {
218 if(repairZeroDiagonals_) {
219 GlobalOrdinal grid = Ac->getRowMap()->getGlobalElement(r);
220 LocalOrdinal lcid = Ac->getColMap()->getLocalElement(grid);
223 Ac->insertLocalValues(r, indout.
view(0,indout.
size()), valout.
view(0,valout.
size()));
230 GO lZeroDiagsGO = lZeroDiags;
233 if(repairZeroDiagonals_) GetOStream(
Warnings0) <<
"RAPShiftFactory (WARNING): repaired " << gZeroDiags <<
" zeros on main diagonal of Ac." << std::endl;
234 else GetOStream(
Warnings0) <<
"RAPShiftFactory (WARNING): found " << gZeroDiags <<
" zeros on main diagonal of Ac." << std::endl;
238 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
241 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)");
242 transferFacts_.push_back(factory);
247 #define MUELU_RAPSHIFTFACTORY_SHORT 248 #endif // MUELU_RAPSHIFTFACTORY_DEF_HPP Important warning messages (one line)
Exception indicating invalid cast attempted.
#define MueLu_sumAll(rcpComm, in, out)
virtual void CallBuild(Level &requestedLevel) const =0
ParameterList & disableRecursiveValidation()
void CheckMainDiagonal(RCP< Matrix > &Ac) const
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.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Namespace for MueLu classes and methods.
Print even more statistics.
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.
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)
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)
ArrayView< T > view(size_type lowerOffset, size_type size) const
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
virtual std::string description() const
Return a simple one-line description of this object.
#define SET_VALID_ENTRY(name)