46 #ifndef MUELU_RAPFACTORY_DEF_HPP 47 #define MUELU_RAPFACTORY_DEF_HPP 52 #include <Xpetra_Matrix.hpp> 53 #include <Xpetra_MatrixFactory.hpp> 54 #include <Xpetra_MatrixMatrix.hpp> 55 #include <Xpetra_MatrixUtils.hpp> 56 #include <Xpetra_TripleMatrixMultiply.hpp> 57 #include <Xpetra_Vector.hpp> 58 #include <Xpetra_VectorFactory.hpp> 64 #include "MueLu_PerfUtils.hpp" 70 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
72 : hasDeclaredInput_(false) { }
74 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
76 RCP<ParameterList> validParamList = rcp(
new ParameterList());
78 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 85 #undef SET_VALID_ENTRY 86 validParamList->set< RCP<const FactoryBase> >(
"A", null,
"Generating factory of the matrix A used during the prolongator smoothing process");
87 validParamList->set< RCP<const FactoryBase> >(
"P", null,
"Prolongator factory");
88 validParamList->set< RCP<const FactoryBase> >(
"R", null,
"Restrictor factory");
90 validParamList->set<
bool > (
"CheckMainDiagonal",
false,
"Check main diagonal for zeros");
91 validParamList->set<
bool > (
"RepairMainDiagonal",
false,
"Repair zeros on main diagonal");
94 ParameterList norecurse;
95 norecurse.disableRecursiveValidation();
96 validParamList->set<ParameterList> (
"matrixmatrix: kernel params", norecurse,
"MatrixMatrix kernel parameters");
98 return validParamList;
101 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
103 const Teuchos::ParameterList& pL = GetParameterList();
104 if (pL.get<
bool>(
"transpose: use implicit") ==
false)
105 Input(coarseLevel,
"R");
107 Input(fineLevel,
"A");
108 Input(coarseLevel,
"P");
111 for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it)
112 (*it)->CallDeclareInput(coarseLevel);
114 hasDeclaredInput_ =
true;
117 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
119 const bool doTranspose =
true;
120 const bool doFillComplete =
true;
121 const bool doOptimizeStorage =
true;
125 std::ostringstream levelstr;
130 "MueLu::RAPFactory::Build(): CallDeclareInput has not been called before Build!");
132 const Teuchos::ParameterList& pL = GetParameterList();
133 RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel,
"A");
134 RCP<Matrix> P = Get< RCP<Matrix> >(coarseLevel,
"P"), AP;
136 bool isEpetra = A->getRowMap()->lib() == Xpetra::UseEpetra;
138 #ifdef KOKKOS_ENABLE_CUDA 139 (
typeid(
Node).name() ==
typeid(Kokkos::Compat::KokkosCudaWrapperNode).name()) ||
141 #ifdef KOKKOS_ENABLE_HIP
142 (
typeid(
Node).name() ==
typeid(Kokkos::Compat::KokkosHIPWrapperNode).name()) ||
146 if (pL.get<
bool>(
"rap: triple product") ==
false || isEpetra || isGPU) {
147 if (pL.get<
bool>(
"rap: triple product") && isEpetra)
148 GetOStream(
Warnings1) <<
"Switching from triple product to R x (A x P) since triple product has not been implemented for Epetra.\n";
149 #if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) 150 if (pL.get<
bool>(
"rap: triple product") && isGPU)
151 GetOStream(
Warnings1) <<
"Switching from triple product to R x (A x P) since triple product has not been implemented for " 152 << Node::execution_space::name() << std::endl;
156 RCP<ParameterList> APparams = rcp(
new ParameterList);
157 if(pL.isSublist(
"matrixmatrix: kernel params"))
158 APparams->sublist(
"matrixmatrix: kernel params") = pL.sublist(
"matrixmatrix: kernel params");
161 APparams->set(
"compute global constants: temporaries",APparams->get(
"compute global constants: temporaries",
false));
162 APparams->set(
"compute global constants",APparams->get(
"compute global constants",
false));
164 if (coarseLevel.IsAvailable(
"AP reuse data",
this)) {
165 GetOStream(static_cast<MsgType>(
Runtime0 |
Test)) <<
"Reusing previous AP data" << std::endl;
167 APparams = coarseLevel.Get< RCP<ParameterList> >(
"AP reuse data",
this);
169 if (APparams->isParameter(
"graph"))
170 AP = APparams->get< RCP<Matrix> >(
"graph");
176 AP = MatrixMatrix::Multiply(*A, !doTranspose, *P, !doTranspose, AP, GetOStream(
Statistics2),
177 doFillComplete, doOptimizeStorage, labelstr+std::string(
"MueLu::A*P-")+levelstr.str(), APparams);
181 RCP<ParameterList> RAPparams = rcp(
new ParameterList);
182 if(pL.isSublist(
"matrixmatrix: kernel params"))
183 RAPparams->sublist(
"matrixmatrix: kernel params") = pL.sublist(
"matrixmatrix: kernel params");
185 if (coarseLevel.IsAvailable(
"RAP reuse data",
this)) {
186 GetOStream(static_cast<MsgType>(
Runtime0 |
Test)) <<
"Reusing previous RAP data" << std::endl;
188 RAPparams = coarseLevel.Get< RCP<ParameterList> >(
"RAP reuse data",
this);
190 if (RAPparams->isParameter(
"graph"))
191 Ac = RAPparams->get< RCP<Matrix> >(
"graph");
195 Ac->SetMaxEigenvalueEstimate(-Teuchos::ScalarTraits<SC>::one());
199 RAPparams->set(
"compute global constants: temporaries",RAPparams->get(
"compute global constants: temporaries",
false));
200 RAPparams->set(
"compute global constants",
true);
206 if (pL.get<
bool>(
"transpose: use implicit") ==
true) {
209 Ac = MatrixMatrix::Multiply(*P, doTranspose, *AP, !doTranspose, Ac, GetOStream(
Statistics2),
210 doFillComplete, doOptimizeStorage, labelstr+std::string(
"MueLu::R*(AP)-implicit-")+levelstr.str(), RAPparams);
213 RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel,
"R");
217 Ac = MatrixMatrix::Multiply(*R, !doTranspose, *AP, !doTranspose, Ac, GetOStream(
Statistics2),
218 doFillComplete, doOptimizeStorage, labelstr+std::string(
"MueLu::R*(AP)-explicit-")+levelstr.str(), RAPparams);
221 Teuchos::ArrayView<const double> relativeFloor = pL.get<Teuchos::Array<double> >(
"rap: relative diagonal floor")();
222 if(relativeFloor.size() > 0) {
223 Xpetra::MatrixUtils<SC,LO,GO,NO>::RelativeDiagonalBoost(Ac, relativeFloor,GetOStream(
Statistics2));
226 bool repairZeroDiagonals = pL.get<
bool>(
"RepairMainDiagonal") || pL.get<
bool>(
"rap: fix zero diagonals");
227 bool checkAc = pL.get<
bool>(
"CheckMainDiagonal")|| pL.get<
bool>(
"rap: fix zero diagonals"); ;
228 if (checkAc || repairZeroDiagonals) {
229 using magnitudeType =
typename Teuchos::ScalarTraits<Scalar>::magnitudeType;
230 magnitudeType threshold;
231 if (pL.isType<magnitudeType>(
"rap: fix zero diagonals threshold"))
232 threshold = pL.get<magnitudeType>(
"rap: fix zero diagonals threshold");
234 threshold = Teuchos::as<magnitudeType>(pL.get<
double>(
"rap: fix zero diagonals threshold"));
235 Scalar replacement = Teuchos::as<Scalar>(pL.get<
double>(
"rap: fix zero diagonals replacement"));
236 Xpetra::MatrixUtils<SC,LO,GO,NO>::CheckRepairMainDiagonal(Ac, repairZeroDiagonals, GetOStream(
Warnings1), threshold, replacement);
240 RCP<ParameterList> params = rcp(
new ParameterList());;
241 params->set(
"printLoadBalancingInfo",
true);
242 params->set(
"printCommInfo",
true);
246 if(!Ac.is_null()) {std::ostringstream oss; oss <<
"A_" << coarseLevel.GetLevelID(); Ac->setObjectLabel(oss.str());}
247 Set(coarseLevel,
"A", Ac);
249 APparams->set(
"graph", AP);
250 Set(coarseLevel,
"AP reuse data", APparams);
251 RAPparams->set(
"graph", Ac);
252 Set(coarseLevel,
"RAP reuse data", RAPparams);
254 RCP<ParameterList> RAPparams = rcp(
new ParameterList);
255 if(pL.isSublist(
"matrixmatrix: kernel params"))
256 RAPparams->sublist(
"matrixmatrix: kernel params") = pL.sublist(
"matrixmatrix: kernel params");
258 if (coarseLevel.IsAvailable(
"RAP reuse data",
this)) {
259 GetOStream(static_cast<MsgType>(
Runtime0 |
Test)) <<
"Reusing previous RAP data" << std::endl;
261 RAPparams = coarseLevel.Get< RCP<ParameterList> >(
"RAP reuse data",
this);
263 if (RAPparams->isParameter(
"graph"))
264 Ac = RAPparams->get< RCP<Matrix> >(
"graph");
268 Ac->SetMaxEigenvalueEstimate(-Teuchos::ScalarTraits<SC>::one());
272 RAPparams->set(
"compute global constants: temporaries",RAPparams->get(
"compute global constants: temporaries",
false));
273 RAPparams->set(
"compute global constants",
true);
275 if (pL.get<
bool>(
"transpose: use implicit") ==
true) {
277 Ac = MatrixFactory::Build(P->getDomainMap(), Teuchos::as<LO>(0));
281 Xpetra::TripleMatrixMultiply<SC,LO,GO,NO>::
282 MultiplyRAP(*P, doTranspose, *A, !doTranspose, *P, !doTranspose, *Ac, doFillComplete,
283 doOptimizeStorage, labelstr+std::string(
"MueLu::R*A*P-implicit-")+levelstr.str(),
287 RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel,
"R");
288 Ac = MatrixFactory::Build(R->getRowMap(), Teuchos::as<LO>(0));
292 Xpetra::TripleMatrixMultiply<SC,LO,GO,NO>::
293 MultiplyRAP(*R, !doTranspose, *A, !doTranspose, *P, !doTranspose, *Ac, doFillComplete,
294 doOptimizeStorage, labelstr+std::string(
"MueLu::R*A*P-explicit-")+levelstr.str(),
298 Teuchos::ArrayView<const double> relativeFloor = pL.get<Teuchos::Array<double> >(
"rap: relative diagonal floor")();
299 if(relativeFloor.size() > 0) {
300 Xpetra::MatrixUtils<SC,LO,GO,NO>::RelativeDiagonalBoost(Ac, relativeFloor,GetOStream(
Statistics2));
303 bool repairZeroDiagonals = pL.get<
bool>(
"RepairMainDiagonal") || pL.get<
bool>(
"rap: fix zero diagonals");
304 bool checkAc = pL.get<
bool>(
"CheckMainDiagonal")|| pL.get<
bool>(
"rap: fix zero diagonals"); ;
305 if (checkAc || repairZeroDiagonals) {
306 using magnitudeType =
typename Teuchos::ScalarTraits<Scalar>::magnitudeType;
307 magnitudeType threshold;
308 if (pL.isType<magnitudeType>(
"rap: fix zero diagonals threshold"))
309 threshold = pL.get<magnitudeType>(
"rap: fix zero diagonals threshold");
311 threshold = Teuchos::as<magnitudeType>(pL.get<
double>(
"rap: fix zero diagonals threshold"));
312 Scalar replacement = Teuchos::as<Scalar>(pL.get<
double>(
"rap: fix zero diagonals replacement"));
313 Xpetra::MatrixUtils<SC,LO,GO,NO>::CheckRepairMainDiagonal(Ac, repairZeroDiagonals, GetOStream(
Warnings1), threshold, replacement);
318 RCP<ParameterList> params = rcp(
new ParameterList());;
319 params->set(
"printLoadBalancingInfo",
true);
320 params->set(
"printCommInfo",
true);
324 if(!Ac.is_null()) {std::ostringstream oss; oss <<
"A_" << coarseLevel.GetLevelID(); Ac->setObjectLabel(oss.str());}
325 Set(coarseLevel,
"A", Ac);
327 RAPparams->set(
"graph", Ac);
328 Set(coarseLevel,
"RAP reuse data", RAPparams);
334 #ifdef HAVE_MUELU_DEBUG 335 MatrixUtils::checkLocalRowMapMatchesColMap(*Ac);
336 #endif // HAVE_MUELU_DEBUG 338 if (transferFacts_.begin() != transferFacts_.end()) {
342 for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it) {
343 RCP<const FactoryBase> fac = *it;
344 GetOStream(
Runtime0) <<
"RAPFactory: call transfer factory: " << fac->description() << std::endl;
345 fac->CallBuild(coarseLevel);
380 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
383 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const TwoLevelFactoryBase>(factory) == Teuchos::null,
Exceptions::BadCast,
384 "MueLu::RAPFactory::AddTransferFactory: Transfer factory is not derived from TwoLevelFactoryBase. " 385 "This is very strange. (Note: you can remove this exception if there's a good reason for)");
386 TEUCHOS_TEST_FOR_EXCEPTION(hasDeclaredInput_,
Exceptions::RuntimeError,
"MueLu::RAPFactory::AddTransferFactory: Factory is being added after we have already declared input");
387 transferFacts_.push_back(factory);
392 #define MUELU_RAPFACTORY_SHORT 393 #endif // MUELU_RAPFACTORY_DEF_HPP #define SET_VALID_ENTRY(name)
Exception indicating invalid cast attempted.
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.
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.
MueLu::DefaultScalar Scalar
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.