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>
74 RCP<ParameterList> validParamList = rcp(
new ParameterList());
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");
80 validParamList->set< RCP<const FactoryBase> >(
"P", null,
"Prolongator factory");
81 validParamList->set< RCP<const FactoryBase> >(
"R", null,
"Restrictor factory");
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>
91 const Teuchos::ParameterList& pL = GetParameterList();
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!");
118 const Teuchos::ParameterList& pL = GetParameterList();
119 RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel,
"A");
120 RCP<Matrix> P = Get< RCP<Matrix> >(coarseLevel,
"P"), AP, Ac;
123 RCP<ParameterList> APparams = rcp(
new ParameterList);
124 if (coarseLevel.IsAvailable(
"AP reuse data",
this)) {
125 GetOStream(static_cast<MsgType>(
Runtime0 |
Test)) <<
"Reusing previous AP data" << std::endl;
127 APparams = coarseLevel.Get< RCP<ParameterList> >(
"AP reuse data",
this);
129 if (APparams->isParameter(
"graph"))
130 AP = APparams->get< RCP<Matrix> >(
"graph");
136 AP = MatrixMatrix::Multiply(*A, !doTranspose, *P, !doTranspose, AP, GetOStream(
Statistics2),
137 doFillComplete, doOptimizeStorage, std::string(
"MueLu::A*P-")+levelstr.str(), APparams);
141 RCP<ParameterList> RAPparams = rcp(
new ParameterList);
142 if (coarseLevel.IsAvailable(
"RAP reuse data",
this)) {
143 GetOStream(static_cast<MsgType>(
Runtime0 |
Test)) <<
"Reusing previous RAP data" << std::endl;
145 RAPparams = coarseLevel.Get< RCP<ParameterList> >(
"RAP reuse data",
this);
147 if (RAPparams->isParameter(
"graph"))
148 Ac = RAPparams->get< RCP<Matrix> >(
"graph");
152 Ac->SetMaxEigenvalueEstimate(-Teuchos::ScalarTraits<SC>::one());
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);
177 RCP<ParameterList> params = rcp(
new ParameterList());;
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) {
196 RCP<const FactoryBase> fac = *it;
197 GetOStream(
Runtime0) <<
"RAPFactory: call transfer factory: " << fac->description() << std::endl;
198 fac->CallBuild(coarseLevel);
234 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
236 const Teuchos::ParameterList& pL = GetParameterList();
237 bool repairZeroDiagonals = pL.get<
bool>(
"RepairMainDiagonal");
238 bool checkAc = pL.get<
bool>(
"CheckMainDiagonal");
240 if (!checkAc && !repairZeroDiagonals)
243 SC zero = Teuchos::ScalarTraits<SC>::zero(), one = Teuchos::ScalarTraits<SC>::one();
245 Teuchos::RCP<Teuchos::ParameterList> p = Teuchos::rcp(
new Teuchos::ParameterList());
246 p->set(
"DoOptimizeStorage",
true);
248 RCP<const Map> rowMap = Ac->getRowMap();
249 RCP<Vector> diagVec = VectorFactory::Build(rowMap);
250 Ac->getLocalDiagCopy(*diagVec);
253 Teuchos::ArrayRCP< Scalar > diagVal = diagVec->getDataNonConst(0);
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);
282 Teuchos::ArrayRCP<GO> indout(1,grid);
283 Teuchos::ArrayRCP<SC> valout(1, one);
284 fixDiagMatrix->insertGlobalValues(grid,indout.view(0, 1), valout.view(0, 1));
288 Teuchos::TimeMonitor m1(*Teuchos::TimeMonitor::getNewTimer(
"CheckRepairMainDiagonal: fillComplete1"));
292 MatrixMatrix::TwoMatrixAdd(*Ac,
false, 1.0, *fixDiagMatrix, 1.0);
293 if (Ac->IsView(
"stridedMaps"))
294 fixDiagMatrix->CreateView(
"stridedMaps", Ac);
302 Teuchos::TimeMonitor m1(*Teuchos::TimeMonitor::getNewTimer(
"CheckRepairMainDiagonal: fillComplete2"));
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);
318 Teuchos::ArrayRCP< Scalar > diagVal2 = diagVec->getDataNonConst(0);
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>
331 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const TwoLevelFactoryBase>(factory) == Teuchos::null,
Exceptions::BadCast,
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)");
334 TEUCHOS_TEST_FOR_EXCEPTION(hasDeclaredInput_,
Exceptions::RuntimeError,
"MueLu::RAPFactory::AddTransferFactory: Factory is being added after we have already declared input");
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)
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 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.
int GetLevelID() const
Return level number.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
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 DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Exception throws to report errors in the internal logical of the program.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.