46 #ifndef MUELU_RAPFACTORY_DEF_HPP 47 #define MUELU_RAPFACTORY_DEF_HPP 62 #include "MueLu_PerfUtils.hpp" 68 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
70 : hasDeclaredInput_(false) { }
72 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
76 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 78 #undef SET_VALID_ENTRY 79 validParamList->
set<
RCP<const FactoryBase> >(
"A", null,
"Generating factory of the matrix A used during the prolongator smoothing process");
83 validParamList->
set<
bool > (
"CheckMainDiagonal",
false,
"Check main diagonal for zeros");
84 validParamList->
set<
bool > (
"RepairMainDiagonal",
false,
"Repair zeros on main diagonal");
86 return validParamList;
89 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
92 if (pL.
get<
bool>(
"transpose: use implicit") ==
false)
93 Input(coarseLevel,
"R");
95 Input(fineLevel,
"A");
96 Input(coarseLevel,
"P");
99 for (std::vector<
RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it)
100 (*it)->CallDeclareInput(coarseLevel);
102 hasDeclaredInput_ =
true;
105 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
107 const bool doTranspose =
true;
108 const bool doFillComplete =
true;
109 const bool doOptimizeStorage =
true;
112 std::ostringstream levelstr;
116 "MueLu::RAPFactory::Build(): CallDeclareInput has not been called before Build!");
119 RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel,
"A");
120 RCP<Matrix> P = Get< RCP<Matrix> >(coarseLevel,
"P"), AP, Ac;
124 if (coarseLevel.IsAvailable(
"AP reuse data",
this)) {
125 GetOStream(static_cast<MsgType>(
Runtime0 |
Test)) <<
"Reusing previous AP data" << std::endl;
129 if (APparams->isParameter(
"graph"))
136 AP = MatrixMatrix::Multiply(*A, !doTranspose, *P, !doTranspose, AP, GetOStream(
Statistics2),
137 doFillComplete, doOptimizeStorage, std::string(
"MueLu::A*P-")+levelstr.str(), APparams);
142 if (coarseLevel.IsAvailable(
"RAP reuse data",
this)) {
143 GetOStream(static_cast<MsgType>(
Runtime0 |
Test)) <<
"Reusing previous RAP data" << std::endl;
159 if (pL.
get<
bool>(
"transpose: use implicit") ==
true) {
162 Ac = MatrixMatrix::Multiply(*P, doTranspose, *AP, !doTranspose, Ac, GetOStream(
Statistics2),
163 doFillComplete, doOptimizeStorage, std::string(
"MueLu::R*(AP)-implicit-")+levelstr.str(), RAPparams);
166 RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel,
"R");
170 Ac = MatrixMatrix::Multiply(*R, !doTranspose, *AP, !doTranspose, Ac, GetOStream(
Statistics2),
171 doFillComplete, doOptimizeStorage, std::string(
"MueLu::R*(AP)-explicit-")+levelstr.str(), RAPparams);
174 CheckRepairMainDiagonal(Ac);
178 params->
set(
"printLoadBalancingInfo",
true);
179 params->
set(
"printCommInfo",
true);
183 Set(coarseLevel,
"A", Ac);
185 APparams->set(
"graph", AP);
186 Set(coarseLevel,
"AP reuse data", APparams);
187 RAPparams->
set(
"graph", Ac);
188 Set(coarseLevel,
"RAP reuse data", RAPparams);
191 if (transferFacts_.begin() != transferFacts_.end()) {
195 for (std::vector<
RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it) {
197 GetOStream(
Runtime0) <<
"RAPFactory: call transfer factory: " << fac->
description() << std::endl;
234 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
237 bool repairZeroDiagonals = pL.
get<
bool>(
"RepairMainDiagonal");
238 bool checkAc = pL.
get<
bool>(
"CheckMainDiagonal");
240 if (!checkAc && !repairZeroDiagonals)
246 p->
set(
"DoOptimizeStorage",
true);
249 RCP<Vector> diagVec = VectorFactory::Build(rowMap);
250 Ac->getLocalDiagCopy(*diagVec);
255 for (
size_t i = 0; i < rowMap->getNodeNumElements(); i++) {
256 if (diagVal[i] == zero) {
261 MueLu_sumAll(rowMap->getComm(), Teuchos::as<GO>(lZeroDiags), gZeroDiags);
263 if (repairZeroDiagonals && gZeroDiags > 0) {
278 RCP<Matrix> fixDiagMatrix = MatrixFactory::Build(rowMap, 1);
279 for (
size_t r = 0; r < rowMap->getNodeNumElements(); r++) {
280 if (diagVal[r] == zero) {
281 GO grid = rowMap->getGlobalElement(r);
284 fixDiagMatrix->insertGlobalValues(grid,indout.
view(0, 1), valout.
view(0, 1));
292 MatrixMatrix::TwoMatrixAdd(*Ac,
false, 1.0, *fixDiagMatrix, 1.0);
293 if (Ac->IsView(
"stridedMaps"))
294 fixDiagMatrix->CreateView(
"stridedMaps", Ac);
312 GetOStream(
Warnings0) <<
"RAPFactory (WARNING): " << (repairZeroDiagonals ?
"repaired " :
"found ")
313 << gZeroDiags <<
" zeros on main diagonal of Ac." << std::endl;
315 #ifdef HAVE_MUELU_DEBUG // only for debugging 317 Ac->getLocalDiagCopy(*diagVec);
319 for (
size_t r = 0; r < Ac->getRowMap()->getNodeNumElements(); r++) {
320 if (diagVal2[r] == zero) {
321 GetOStream(
Errors,-1) <<
"Error: there are zeros left on diagonal after repair..." << std::endl;
328 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
332 "MueLu::RAPFactory::AddTransferFactory: Transfer factory is not derived from TwoLevelFactoryBase. " 333 "This is very strange. (Note: you can remove this exception if there's a good reason for)");
335 transferFacts_.push_back(factory);
340 #define MUELU_RAPFACTORY_SHORT 341 #endif // MUELU_RAPFACTORY_DEF_HPP Important warning messages (one line)
#define SET_VALID_ENTRY(name)
Exception indicating invalid cast attempted.
#define MueLu_sumAll(rcpComm, in, out)
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.
void AddTransferFactory(const RCP< const FactoryBase > &factory)
Add transfer factory in the end of list of transfer factories in RepartitionAcFactory.
Namespace for MueLu classes and methods.
void CheckRepairMainDiagonal(RCP< Matrix > &Ac) const
Print skeleton for the run, i.e. factory calls and used parameters.
Print even more statistics.
static RCP< Time > getNewTimer(const std::string &name)
int GetLevelID() const
Return level number.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
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)
ArrayView< T > view(size_type lowerOffset, size_type size) const
bool isParameter(const std::string &name) const
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Exception throws to report errors in the internal logical of the program.
virtual std::string description() const
Return a simple one-line description of this object.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.