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" 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>
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");
105 KP = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*K,
false, *P,
false, KP, GetOStream(
Statistics2));
106 MP = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*M,
false, *P,
false, MP, GetOStream(
Statistics2));
107 Set(coarseLevel,
"AP Pattern", KP);
112 bool doOptimizedStorage = !checkAc_;
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);
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) {
158 GetOStream(
Runtime0) <<
"RAPShiftFactory: call transfer factory: " << fac->
description() << std::endl;
167 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
171 RCP<Vector> diagVec = VectorFactory::Build(Ac->getRowMap());
172 Ac->getLocalDiagCopy(*diagVec);
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);
182 Ac->insertLocalValues(r, indout.
view(0,indout.
size()), valout.
view(0,valout.
size()));
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)
virtual void CallBuild(Level &requestedLevel) const =0
void CheckMainDiagonal(RCP< Matrix > &Ac) const
void Release(const FactoryBase &factory)
Decrement the storage counter for all the inputs of a factory.
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.
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.
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.