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);