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_Vector.hpp> 54 #include <Xpetra_VectorFactory.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>
73 RCP<ParameterList> validParamList = rcp(
new ParameterList());
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");
82 validParamList->set< RCP<const FactoryBase> >(
"P", Teuchos::null,
"Prolongator factory");
83 validParamList->set< RCP<const FactoryBase> >(
"R", Teuchos::null,
"Restrictor factory");
85 validParamList->set<
bool > (
"CheckMainDiagonal",
false,
"Check main diagonal for zeros");
86 validParamList->set<
bool > (
"RepairMainDiagonal",
false,
"Repair zeros on main diagonal");
89 ParameterList norecurse;
90 norecurse.disableRecursiveValidation();
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>
117 const Teuchos::ParameterList& pL = GetParameterList();
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");
136 KP = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*K,
false, *P,
false, KP, GetOStream(
Statistics2));
137 MP = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*M,
false, *P,
false, MP, GetOStream(
Statistics2));
138 Set(coarseLevel,
"AP Pattern", KP);
143 bool doOptimizedStorage = !checkAc_;
145 RCP<Matrix> Ac, Kc, Mc;
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();
167 Scalar shift = Teuchos::ScalarTraits<Scalar>::zero();
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);
183 RCP<ParameterList> params = rcp(
new ParameterList());;
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) {
198 RCP<const FactoryBase> fac = *it;
199 GetOStream(
Runtime0) <<
"RAPShiftFactory: call transfer factory: " << fac->description() << std::endl;
200 fac->CallBuild(coarseLevel);
208 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
212 RCP<Vector> diagVec = VectorFactory::Build(Ac->getRowMap());
213 Ac->getLocalDiagCopy(*diagVec);
214 Teuchos::ArrayRCP< Scalar > diagVal = diagVec->getDataNonConst(0);
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);
221 Teuchos::ArrayRCP<LocalOrdinal> indout(1,lcid);
222 Teuchos::ArrayRCP<Scalar> valout(1,Teuchos::ScalarTraits<Scalar>::one());
223 Ac->insertLocalValues(r, indout.view(0,indout.size()), valout.view(0,valout.size()));
229 const RCP<const Teuchos::Comm<int> > & comm = Ac->getRowMap()->getComm();
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)
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.
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.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
#define SET_VALID_ENTRY(name)