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>
59#include <Xpetra_IO.hpp>
65#include "MueLu_PerfUtils.hpp"
71 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
73 : hasDeclaredInput_(false) { }
75 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
77 RCP<ParameterList> validParamList = rcp(
new ParameterList());
79#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
87 validParamList->set< RCP<const FactoryBase> >(
"A", null,
"Generating factory of the matrix A used during the prolongator smoothing process");
88 validParamList->set< RCP<const FactoryBase> >(
"P", null,
"Prolongator factory");
89 validParamList->set< RCP<const FactoryBase> >(
"R", null,
"Restrictor factory");
91 validParamList->set<
bool > (
"CheckMainDiagonal",
false,
"Check main diagonal for zeros");
92 validParamList->set<
bool > (
"RepairMainDiagonal",
false,
"Repair zeros on main diagonal");
95 ParameterList norecurse;
96 norecurse.disableRecursiveValidation();
97 validParamList->set<ParameterList> (
"matrixmatrix: kernel params", norecurse,
"MatrixMatrix kernel parameters");
99 return validParamList;
102 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
104 const Teuchos::ParameterList& pL = GetParameterList();
105 if (pL.get<
bool>(
"transpose: use implicit") ==
false)
106 Input(coarseLevel,
"R");
108 Input(fineLevel,
"A");
109 Input(coarseLevel,
"P");
112 for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it)
113 (*it)->CallDeclareInput(coarseLevel);
115 hasDeclaredInput_ =
true;
118 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
120 const bool doTranspose =
true;
121 const bool doFillComplete =
true;
122 const bool doOptimizeStorage =
true;
126 std::ostringstream levelstr;
131 "MueLu::RAPFactory::Build(): CallDeclareInput has not been called before Build!");
133 const Teuchos::ParameterList& pL = GetParameterList();
134 RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel,
"A");
135 RCP<Matrix> P = Get< RCP<Matrix> >(coarseLevel,
"P"), AP;
138 if (P == Teuchos::null) {
140 Set(coarseLevel,
"A", Ac);
144 bool isEpetra = A->getRowMap()->lib() == Xpetra::UseEpetra;
146#ifdef KOKKOS_ENABLE_CUDA
147 (
typeid(
Node).name() ==
typeid(Kokkos::Compat::KokkosCudaWrapperNode).name()) ||
149#ifdef KOKKOS_ENABLE_HIP
150 (
typeid(
Node).name() ==
typeid(Kokkos::Compat::KokkosHIPWrapperNode).name()) ||
154 if (pL.get<
bool>(
"rap: triple product") ==
false || isEpetra || isGPU) {
155 if (pL.get<
bool>(
"rap: triple product") && isEpetra)
156 GetOStream(
Warnings1) <<
"Switching from triple product to R x (A x P) since triple product has not been implemented for Epetra.\n";
157#if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP)
158 if (pL.get<
bool>(
"rap: triple product") && isGPU)
159 GetOStream(
Warnings1) <<
"Switching from triple product to R x (A x P) since triple product has not been implemented for "
160 << Node::execution_space::name() << std::endl;
164 RCP<ParameterList> APparams = rcp(
new ParameterList);
165 if(pL.isSublist(
"matrixmatrix: kernel params"))
166 APparams->sublist(
"matrixmatrix: kernel params") = pL.sublist(
"matrixmatrix: kernel params");
169 APparams->set(
"compute global constants: temporaries",APparams->get(
"compute global constants: temporaries",
false));
170 APparams->set(
"compute global constants",APparams->get(
"compute global constants",
false));
172 if (coarseLevel.IsAvailable(
"AP reuse data",
this)) {
173 GetOStream(
static_cast<MsgType>(
Runtime0 |
Test)) <<
"Reusing previous AP data" << std::endl;
175 APparams = coarseLevel.Get< RCP<ParameterList> >(
"AP reuse data",
this);
177 if (APparams->isParameter(
"graph"))
178 AP = APparams->get< RCP<Matrix> >(
"graph");
184 AP = MatrixMatrix::Multiply(*A, !doTranspose, *P, !doTranspose, AP, GetOStream(
Statistics2),
185 doFillComplete, doOptimizeStorage, labelstr+std::string(
"MueLu::A*P-")+levelstr.str(), APparams);
189 RCP<ParameterList> RAPparams = rcp(
new ParameterList);
190 if(pL.isSublist(
"matrixmatrix: kernel params"))
191 RAPparams->sublist(
"matrixmatrix: kernel params") = pL.sublist(
"matrixmatrix: kernel params");
193 if (coarseLevel.IsAvailable(
"RAP reuse data",
this)) {
194 GetOStream(
static_cast<MsgType>(
Runtime0 |
Test)) <<
"Reusing previous RAP data" << std::endl;
196 RAPparams = coarseLevel.Get< RCP<ParameterList> >(
"RAP reuse data",
this);
198 if (RAPparams->isParameter(
"graph"))
199 Ac = RAPparams->get< RCP<Matrix> >(
"graph");
203 Ac->SetMaxEigenvalueEstimate(-Teuchos::ScalarTraits<SC>::one());
207 RAPparams->set(
"compute global constants: temporaries",RAPparams->get(
"compute global constants: temporaries",
false));
208 RAPparams->set(
"compute global constants",
true);
214 if (pL.get<
bool>(
"transpose: use implicit") ==
true) {
217 Ac = MatrixMatrix::Multiply(*P, doTranspose, *AP, !doTranspose, Ac, GetOStream(
Statistics2),
218 doFillComplete, doOptimizeStorage, labelstr+std::string(
"MueLu::R*(AP)-implicit-")+levelstr.str(), RAPparams);
221 RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel,
"R");
225 Ac = MatrixMatrix::Multiply(*R, !doTranspose, *AP, !doTranspose, Ac, GetOStream(
Statistics2),
226 doFillComplete, doOptimizeStorage, labelstr+std::string(
"MueLu::R*(AP)-explicit-")+levelstr.str(), RAPparams);
229 Teuchos::ArrayView<const double> relativeFloor = pL.get<Teuchos::Array<double> >(
"rap: relative diagonal floor")();
230 if(relativeFloor.size() > 0) {
231 Xpetra::MatrixUtils<SC,LO,GO,NO>::RelativeDiagonalBoost(Ac, relativeFloor,GetOStream(
Statistics2));
234 bool repairZeroDiagonals = pL.get<
bool>(
"RepairMainDiagonal") || pL.get<
bool>(
"rap: fix zero diagonals");
235 bool checkAc = pL.get<
bool>(
"CheckMainDiagonal")|| pL.get<
bool>(
"rap: fix zero diagonals"); ;
236 if (checkAc || repairZeroDiagonals) {
237 using magnitudeType =
typename Teuchos::ScalarTraits<Scalar>::magnitudeType;
238 magnitudeType threshold;
239 if (pL.isType<magnitudeType>(
"rap: fix zero diagonals threshold"))
240 threshold = pL.get<magnitudeType>(
"rap: fix zero diagonals threshold");
242 threshold = Teuchos::as<magnitudeType>(pL.get<
double>(
"rap: fix zero diagonals threshold"));
243 Scalar replacement = Teuchos::as<Scalar>(pL.get<
double>(
"rap: fix zero diagonals replacement"));
244 Xpetra::MatrixUtils<SC,LO,GO,NO>::CheckRepairMainDiagonal(Ac, repairZeroDiagonals, GetOStream(
Warnings1), threshold, replacement);
248 RCP<ParameterList> params = rcp(
new ParameterList());;
249 params->set(
"printLoadBalancingInfo",
true);
250 params->set(
"printCommInfo",
true);
254 if(!Ac.is_null()) {std::ostringstream oss; oss <<
"A_" << coarseLevel.GetLevelID(); Ac->setObjectLabel(oss.str());}
255 Set(coarseLevel,
"A", Ac);
258 APparams->set(
"graph", AP);
259 Set(coarseLevel,
"AP reuse data", APparams);
262 RAPparams->set(
"graph", Ac);
263 Set(coarseLevel,
"RAP reuse data", RAPparams);
266 RCP<ParameterList> RAPparams = rcp(
new ParameterList);
267 if(pL.isSublist(
"matrixmatrix: kernel params"))
268 RAPparams->sublist(
"matrixmatrix: kernel params") = pL.sublist(
"matrixmatrix: kernel params");
270 if (coarseLevel.IsAvailable(
"RAP reuse data",
this)) {
271 GetOStream(
static_cast<MsgType>(
Runtime0 |
Test)) <<
"Reusing previous RAP data" << std::endl;
273 RAPparams = coarseLevel.Get< RCP<ParameterList> >(
"RAP reuse data",
this);
275 if (RAPparams->isParameter(
"graph"))
276 Ac = RAPparams->get< RCP<Matrix> >(
"graph");
280 Ac->SetMaxEigenvalueEstimate(-Teuchos::ScalarTraits<SC>::one());
284 RAPparams->set(
"compute global constants: temporaries",RAPparams->get(
"compute global constants: temporaries",
false));
285 RAPparams->set(
"compute global constants",
true);
287 if (pL.get<
bool>(
"transpose: use implicit") ==
true) {
289 Ac = MatrixFactory::Build(P->getDomainMap(), Teuchos::as<LO>(0));
293 Xpetra::TripleMatrixMultiply<SC,LO,GO,NO>::
294 MultiplyRAP(*P, doTranspose, *A, !doTranspose, *P, !doTranspose, *Ac, doFillComplete,
295 doOptimizeStorage, labelstr+std::string(
"MueLu::R*A*P-implicit-")+levelstr.str(),
298 RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel,
"R");
299 Ac = MatrixFactory::Build(R->getRowMap(), Teuchos::as<LO>(0));
303 Xpetra::TripleMatrixMultiply<SC,LO,GO,NO>::
304 MultiplyRAP(*R, !doTranspose, *A, !doTranspose, *P, !doTranspose, *Ac, doFillComplete,
305 doOptimizeStorage, labelstr+std::string(
"MueLu::R*A*P-explicit-")+levelstr.str(),
309 Teuchos::ArrayView<const double> relativeFloor = pL.get<Teuchos::Array<double> >(
"rap: relative diagonal floor")();
310 if(relativeFloor.size() > 0) {
311 Xpetra::MatrixUtils<SC,LO,GO,NO>::RelativeDiagonalBoost(Ac, relativeFloor,GetOStream(
Statistics2));
314 bool repairZeroDiagonals = pL.get<
bool>(
"RepairMainDiagonal") || pL.get<
bool>(
"rap: fix zero diagonals");
315 bool checkAc = pL.get<
bool>(
"CheckMainDiagonal")|| pL.get<
bool>(
"rap: fix zero diagonals"); ;
316 if (checkAc || repairZeroDiagonals) {
317 using magnitudeType =
typename Teuchos::ScalarTraits<Scalar>::magnitudeType;
318 magnitudeType threshold;
319 if (pL.isType<magnitudeType>(
"rap: fix zero diagonals threshold"))
320 threshold = pL.get<magnitudeType>(
"rap: fix zero diagonals threshold");
322 threshold = Teuchos::as<magnitudeType>(pL.get<
double>(
"rap: fix zero diagonals threshold"));
323 Scalar replacement = Teuchos::as<Scalar>(pL.get<
double>(
"rap: fix zero diagonals replacement"));
324 Xpetra::MatrixUtils<SC,LO,GO,NO>::CheckRepairMainDiagonal(Ac, repairZeroDiagonals, GetOStream(
Warnings1), threshold, replacement);
329 RCP<ParameterList> params = rcp(
new ParameterList());;
330 params->set(
"printLoadBalancingInfo",
true);
331 params->set(
"printCommInfo",
true);
335 if(!Ac.is_null()) {std::ostringstream oss; oss <<
"A_" << coarseLevel.GetLevelID(); Ac->setObjectLabel(oss.str());}
336 Set(coarseLevel,
"A", Ac);
339 RAPparams->set(
"graph", Ac);
340 Set(coarseLevel,
"RAP reuse data", RAPparams);
347#ifdef HAVE_MUELU_DEBUG
348 MatrixUtils::checkLocalRowMapMatchesColMap(*Ac);
351 if (transferFacts_.begin() != transferFacts_.end()) {
355 for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it) {
356 RCP<const FactoryBase> fac = *it;
357 GetOStream(
Runtime0) <<
"RAPFactory: call transfer factory: " << fac->description() << std::endl;
358 fac->CallBuild(coarseLevel);
393 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
396 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const TwoLevelFactoryBase>(factory) == Teuchos::null,
Exceptions::BadCast,
397 "MueLu::RAPFactory::AddTransferFactory: Transfer factory is not derived from TwoLevelFactoryBase. "
398 "This is very strange. (Note: you can remove this exception if there's a good reason for)");
399 TEUCHOS_TEST_FOR_EXCEPTION(hasDeclaredInput_,
Exceptions::RuntimeError,
"MueLu::RAPFactory::AddTransferFactory: Factory is being added after we have already declared input");
400 transferFacts_.push_back(factory);
405#define MUELU_RAPFACTORY_SHORT
#define SET_VALID_ENTRY(name)
MueLu::DefaultScalar Scalar
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
void Release(const FactoryBase &factory)
Decrement the storage counter for all the inputs of a factory.
int GetLevelID() const
Return level number.
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
void AddTransferFactory(const RCP< const FactoryBase > &factory)
Add transfer factory in the end of list of transfer factories in RepartitionAcFactory.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Namespace for MueLu classes and methods.
@ Statistics2
Print even more statistics.
@ Runtime0
One-liner description of what is happening.
@ Warnings1
Additional warnings.