46 #ifndef MUELU_RAPSHIFTFACTORY_DEF_HPP 47 #define MUELU_RAPSHIFTFACTORY_DEF_HPP 59 #include "MueLu_PerfUtils.hpp" 60 #include "MueLu_Utilities.hpp" 64 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
66 : implicitTranspose_(false), checkAc_(false), repairZeroDiagonals_(false) { }
68 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
70 if (implicitTranspose_ ==
false) {
71 Input(coarseLevel,
"R");
74 Input(fineLevel,
"K");
75 Input(fineLevel,
"M");
76 Input(coarseLevel,
"P");
79 for(std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it!=transferFacts_.end(); ++it) {
80 (*it)->CallDeclareInput(coarseLevel);
84 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
90 RCP<Matrix> K = Get< RCP<Matrix> >(fineLevel,
"K");
91 RCP<Matrix> M = Get< RCP<Matrix> >(fineLevel,
"M");
92 RCP<Matrix> P = Get< RCP<Matrix> >(coarseLevel,
"P");
98 if (IsAvailable(coarseLevel,
"AP Pattern")) {
99 KP = Get< RCP<Matrix> >(coarseLevel,
"AP Pattern");
100 MP = Get< RCP<Matrix> >(coarseLevel,
"AP Pattern");
107 Set(coarseLevel,
"AP Pattern", KP);
112 bool doOptimizedStorage = !checkAc_;
114 RCP<Matrix> Ac, Kc, Mc;
120 bool doFillComplete=
true;
121 if (implicitTranspose_) {
123 Kc =
Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*P,
true, *KP,
false, Kc, GetOStream(
Statistics2), doFillComplete, doOptimizedStorage);
124 Mc =
Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*P,
true, *MP,
false, Mc, GetOStream(
Statistics2), doFillComplete, doOptimizedStorage);
127 RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel,
"R");
129 Kc =
Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*R,
false, *KP,
false, Kc, GetOStream(
Statistics2), doFillComplete, doOptimizedStorage);
130 Mc =
Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*R,
false, *MP,
false, Mc, GetOStream(
Statistics2), doFillComplete, doOptimizedStorage);
134 int level = coarseLevel.GetLevelID();
135 Scalar shift = shifts_[level];
136 Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixAdd(*Kc,
false, (Scalar) 1.0, *Mc,
false, shift, Ac, GetOStream(
Statistics2));
140 CheckMainDiagonal(Ac);
142 RCP<ParameterList> params = rcp(
new ParameterList());;
143 params->set(
"printLoadBalancingInfo",
true);
146 Set(coarseLevel,
"A", Ac);
147 Set(coarseLevel,
"K", Kc);
148 Set(coarseLevel,
"M", Mc);
149 Set(coarseLevel,
"RAP Pattern", Ac);
152 if (transferFacts_.begin() != transferFacts_.end()) {
156 for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it) {
157 RCP<const FactoryBase> fac = *it;
158 GetOStream(
Runtime0) <<
"RAPShiftFactory: call transfer factory: " << fac->description() << std::endl;
159 fac->CallBuild(coarseLevel);
167 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
171 RCP<Vector> diagVec = VectorFactory::Build(Ac->getRowMap());
172 Ac->getLocalDiagCopy(*diagVec);
173 Teuchos::ArrayRCP< Scalar > diagVal = diagVec->getDataNonConst(0);
174 for (
size_t r=0; r<Ac->getRowMap()->getNodeNumElements(); r++) {
175 if(diagVal[r]==0.0) {
177 if(repairZeroDiagonals_) {
178 GlobalOrdinal grid = Ac->getRowMap()->getGlobalElement(r);
179 LocalOrdinal lcid = Ac->getColMap()->getLocalElement(grid);
180 Teuchos::ArrayRCP<LocalOrdinal> indout(1,lcid);
181 Teuchos::ArrayRCP<Scalar> valout(1,Teuchos::ScalarTraits<Scalar>::one());
182 Ac->insertLocalValues(r, indout.view(0,indout.size()), valout.view(0,valout.size()));
188 const RCP<const Teuchos::Comm<int> > & comm = Ac->getRowMap()->getComm();
189 GO lZeroDiagsGO = lZeroDiags;
192 if(repairZeroDiagonals_) GetOStream(
Warnings0) <<
"RAPShiftFactory (WARNING): repaired " << gZeroDiags <<
" zeros on main diagonal of Ac." << std::endl;
193 else GetOStream(
Warnings0) <<
"RAPShiftFactory (WARNING): found " << gZeroDiags <<
" zeros on main diagonal of Ac." << std::endl;
197 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
200 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)");
201 transferFacts_.push_back(factory);
206 #define MUELU_RAPSHIFTFACTORY_SHORT 207 #endif // MUELU_RAPSHIFTFACTORY_DEF_HPP Important warning messages (one line)
Exception indicating invalid cast attempted.
#define MueLu_sumAll(rcpComm, in, out)
void CheckMainDiagonal(RCP< Matrix > &Ac) const
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.
Print even more statistics.
Print statistics that do not involve significant additional computation.
Class that holds all level-specific information.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
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)
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.