10 #ifndef THYRA_MUELU_PRECONDITIONER_FACTORY_DEF_HPP
11 #define THYRA_MUELU_PRECONDITIONER_FACTORY_DEF_HPP
17 #if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA)
20 #if ((defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) && !defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) && !defined(HAVE_TPETRA_INST_COMPLEX_FLOAT)) || \
21 (!defined(HAVE_TPETRA_INST_DOUBLE) && !defined(HAVE_TPETRA_INST_FLOAT) && defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) && defined(HAVE_TPETRA_INST_COMPLEX_FLOAT)) || \
22 (defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) && defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) && defined(HAVE_TPETRA_INST_COMPLEX_FLOAT)))
23 #define MUELU_CAN_USE_MIXED_PRECISION
31 using Teuchos::rcp_const_cast;
32 using Teuchos::rcp_dynamic_cast;
34 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
35 bool Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(ParameterList& paramList, std::string parameterName) {
38 typedef Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpThyUtils;
46 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
47 typedef Thyra::DiagonalLinearOpBase<Scalar> ThyDiagLinOpBase;
54 typedef Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> thyTpV;
57 #if defined(MUELU_CAN_USE_MIXED_PRECISION)
62 if (paramList.isParameter(parameterName)) {
63 if (paramList.isType<RCP<XpMat> >(parameterName))
65 else if (paramList.isType<RCP<const XpMat> >(parameterName)) {
66 RCP<const XpMat> constM = paramList.get<RCP<const XpMat> >(parameterName);
67 paramList.remove(parameterName);
68 RCP<XpMat> M = rcp_const_cast<XpMat>(constM);
69 paramList.set<RCP<XpMat> >(parameterName, M);
71 }
else if (paramList.isType<RCP<XpMultVec> >(parameterName))
73 else if (paramList.isType<RCP<const XpMultVec> >(parameterName)) {
74 RCP<const XpMultVec> constX = paramList.get<RCP<const XpMultVec> >(parameterName);
75 paramList.remove(parameterName);
76 RCP<XpMultVec> X = rcp_const_cast<XpMultVec>(constX);
77 paramList.set<RCP<XpMultVec> >(parameterName, X);
79 }
else if (paramList.isType<RCP<XpMagMultVec> >(parameterName))
81 else if (paramList.isType<RCP<const XpMagMultVec> >(parameterName)) {
82 RCP<const XpMagMultVec> constX = paramList.get<RCP<const XpMagMultVec> >(parameterName);
83 paramList.remove(parameterName);
84 RCP<XpMagMultVec> X = rcp_const_cast<XpMagMultVec>(constX);
85 paramList.set<RCP<XpMagMultVec> >(parameterName, X);
87 }
else if (paramList.isType<RCP<TpCrsMat> >(parameterName)) {
88 RCP<TpCrsMat> tM = paramList.get<RCP<TpCrsMat> >(parameterName);
89 paramList.remove(parameterName);
90 RCP<XpMat> xM = MueLu::TpetraCrs_To_XpetraMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tM);
91 paramList.set<RCP<XpMat> >(parameterName, xM);
93 }
else if (paramList.isType<RCP<tMV> >(parameterName)) {
94 RCP<tMV> tpetra_X = paramList.get<RCP<tMV> >(parameterName);
95 paramList.remove(parameterName);
96 RCP<XpMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpetra_X);
97 paramList.set<RCP<XpMultVec> >(parameterName, X);
100 }
else if (paramList.isType<RCP<tMagMV> >(parameterName)) {
101 RCP<tMagMV> tpetra_X = paramList.get<RCP<tMagMV> >(parameterName);
102 paramList.remove(parameterName);
103 RCP<XpMagMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node>(tpetra_X);
104 paramList.set<RCP<XpMagMultVec> >(parameterName, X);
108 #if defined(MUELU_CAN_USE_MIXED_PRECISION)
109 else if (paramList.isType<RCP<tHalfMagMV> >(parameterName)) {
110 RCP<tHalfMagMV> tpetra_hX = paramList.get<RCP<tHalfMagMV> >(parameterName);
111 paramList.remove(parameterName);
112 RCP<tMagMV> tpetra_X =
rcp(
new tMagMV(tpetra_hX->getMap(), tpetra_hX->getNumVectors()));
114 RCP<XpMagMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node>(tpetra_X);
115 paramList.set<RCP<XpMagMultVec> >(parameterName, X);
120 else if (paramList.isType<RCP<const ThyDiagLinOpBase> >(parameterName) ||
121 (paramList.isType<RCP<const ThyLinOpBase> >(parameterName) && !rcp_dynamic_cast<const ThyDiagLinOpBase>(paramList.get<RCP<const ThyLinOpBase> >(parameterName)).is_null())) {
122 RCP<const ThyDiagLinOpBase> thyM;
123 if (paramList.isType<RCP<const ThyDiagLinOpBase> >(parameterName))
124 thyM = paramList.get<RCP<const ThyDiagLinOpBase> >(parameterName);
126 thyM = rcp_dynamic_cast<
const ThyDiagLinOpBase>(paramList.get<RCP<const ThyLinOpBase> >(parameterName),
true);
127 paramList.remove(parameterName);
128 RCP<const Thyra::VectorBase<Scalar> > diag = thyM->getDiag();
130 RCP<const XpVec> xpDiag;
131 if (!rcp_dynamic_cast<const thyTpV>(diag).
is_null()) {
132 RCP<const tV> tDiag = Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getConstTpetraVector(diag);
133 if (!tDiag.is_null())
138 paramList.set<RCP<XpMat> >(parameterName, M);
140 }
else if (paramList.isType<RCP<const ThyLinOpBase> >(parameterName)) {
141 RCP<const ThyLinOpBase> thyM = paramList.get<RCP<const ThyLinOpBase> >(parameterName);
142 paramList.remove(parameterName);
144 RCP<XpMat> M = XpThyUtils::toXpetra(Teuchos::rcp_const_cast<ThyLinOpBase>(thyM));
145 paramList.set<RCP<XpMat> >(parameterName, M);
146 }
catch (std::exception& e) {
147 RCP<XpOp> M = XpThyUtils::toXpetraOperator(Teuchos::rcp_const_cast<ThyLinOpBase>(thyM));
151 if (tO->hasDiagonal()) {
152 diag =
rcp(
new tV(tO->getRangeMap()));
153 tO->getLocalDiagCopy(*diag);
157 auto op = rcp_dynamic_cast<XpOp>(tpFOp);
158 paramList.set<RCP<XpOp> >(parameterName, op);
169 #ifdef HAVE_MUELU_EPETRA
170 template <
class GlobalOrdinal>
171 bool Converters<double, int, GlobalOrdinal, Tpetra::KokkosCompat::KokkosSerialWrapperNode>::replaceWithXpetra(ParameterList& paramList, std::string parameterName) {
174 typedef Tpetra::KokkosCompat::KokkosSerialWrapperNode
Node;
177 typedef Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpThyUtils;
185 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
186 typedef Thyra::DiagonalLinearOpBase<Scalar> ThyDiagLinOpBase;
187 typedef Thyra::SpmdVectorSpaceBase<Scalar> ThyVSBase;
192 typedef Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> thyTpV;
195 #if defined(MUELU_CAN_USE_MIXED_PRECISION)
199 #if defined(HAVE_MUELU_EPETRA)
203 if (paramList.isParameter(parameterName)) {
204 if (paramList.isType<RCP<XpMat> >(parameterName))
206 else if (paramList.isType<RCP<const XpMat> >(parameterName)) {
207 RCP<const XpMat> constM = paramList.get<RCP<const XpMat> >(parameterName);
208 paramList.remove(parameterName);
209 RCP<XpMat> M = rcp_const_cast<XpMat>(constM);
210 paramList.set<RCP<XpMat> >(parameterName, M);
212 }
else if (paramList.isType<RCP<XpMultVec> >(parameterName))
214 else if (paramList.isType<RCP<const XpMultVec> >(parameterName)) {
215 RCP<const XpMultVec> constX = paramList.get<RCP<const XpMultVec> >(parameterName);
216 paramList.remove(parameterName);
217 RCP<XpMultVec> X = rcp_const_cast<XpMultVec>(constX);
218 paramList.set<RCP<XpMultVec> >(parameterName, X);
220 }
else if (paramList.isType<RCP<XpMagMultVec> >(parameterName))
222 else if (paramList.isType<RCP<const XpMagMultVec> >(parameterName)) {
223 RCP<const XpMagMultVec> constX = paramList.get<RCP<const XpMagMultVec> >(parameterName);
224 paramList.remove(parameterName);
225 RCP<XpMagMultVec> X = rcp_const_cast<XpMagMultVec>(constX);
226 paramList.set<RCP<XpMagMultVec> >(parameterName, X);
228 }
else if (paramList.isType<RCP<TpCrsMat> >(parameterName)) {
229 RCP<TpCrsMat> tM = paramList.get<RCP<TpCrsMat> >(parameterName);
230 paramList.remove(parameterName);
231 RCP<XpMat> xM = MueLu::TpetraCrs_To_XpetraMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tM);
232 paramList.set<RCP<XpMat> >(parameterName, xM);
234 }
else if (paramList.isType<RCP<tMV> >(parameterName)) {
235 RCP<tMV> tpetra_X = paramList.get<RCP<tMV> >(parameterName);
236 paramList.remove(parameterName);
237 RCP<XpMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpetra_X);
238 paramList.set<RCP<XpMultVec> >(parameterName, X);
241 }
else if (paramList.isType<RCP<tMagMV> >(parameterName)) {
242 RCP<tMagMV> tpetra_X = paramList.get<RCP<tMagMV> >(parameterName);
243 paramList.remove(parameterName);
244 RCP<XpMagMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node>(tpetra_X);
245 paramList.set<RCP<XpMagMultVec> >(parameterName, X);
249 #if defined(MUELU_CAN_USE_MIXED_PRECISION)
250 else if (paramList.isType<RCP<tHalfMagMV> >(parameterName)) {
251 RCP<tHalfMagMV> tpetra_hX = paramList.get<RCP<tHalfMagMV> >(parameterName);
252 paramList.remove(parameterName);
253 RCP<tMagMV> tpetra_X =
rcp(
new tMagMV(tpetra_hX->getMap(), tpetra_hX->getNumVectors()));
255 RCP<XpMagMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node>(tpetra_X);
256 paramList.set<RCP<XpMagMultVec> >(parameterName, X);
261 #ifdef HAVE_MUELU_EPETRA
262 else if (paramList.isType<RCP<Epetra_CrsMatrix> >(parameterName)) {
263 RCP<Epetra_CrsMatrix> eM = paramList.get<RCP<Epetra_CrsMatrix> >(parameterName);
264 paramList.remove(parameterName);
265 RCP<XpEpCrsMat> xeM =
rcp(
new XpEpCrsMat(eM));
266 RCP<XpCrsMat> xCrsM = rcp_dynamic_cast<XpCrsMat>(xeM,
true);
267 RCP<XpCrsMatWrap> xwM =
rcp(
new XpCrsMatWrap(xCrsM));
268 RCP<XpMat> xM = rcp_dynamic_cast<XpMat>(xwM);
269 paramList.set<RCP<XpMat> >(parameterName, xM);
271 }
else if (paramList.isType<RCP<Epetra_MultiVector> >(parameterName)) {
272 RCP<Epetra_MultiVector> epetra_X = Teuchos::null;
273 epetra_X = paramList.get<RCP<Epetra_MultiVector> >(parameterName);
274 paramList.remove(parameterName);
277 RCP<XpMultVec> X = rcp_dynamic_cast<XpMultVec>(xpEpXMult,
true);
278 paramList.set<RCP<XpMultVec> >(parameterName, X);
282 else if (paramList.isType<RCP<const ThyDiagLinOpBase> >(parameterName) ||
283 (paramList.isType<RCP<const ThyLinOpBase> >(parameterName) && !rcp_dynamic_cast<const ThyDiagLinOpBase>(paramList.get<RCP<const ThyLinOpBase> >(parameterName)).
is_null())) {
284 RCP<const ThyDiagLinOpBase> thyM;
285 if (paramList.isType<RCP<const ThyDiagLinOpBase> >(parameterName))
286 thyM = paramList.get<RCP<const ThyDiagLinOpBase> >(parameterName);
288 thyM = rcp_dynamic_cast<
const ThyDiagLinOpBase>(paramList.get<RCP<const ThyLinOpBase> >(parameterName),
true);
289 paramList.remove(parameterName);
290 RCP<const Thyra::VectorBase<Scalar> > diag = thyM->getDiag();
292 RCP<const XpVec> xpDiag;
293 if (!rcp_dynamic_cast<const thyTpV>(diag).
is_null()) {
294 RCP<const tV> tDiag = Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getConstTpetraVector(diag);
295 if (!tDiag.is_null())
298 #ifdef HAVE_MUELU_EPETRA
299 if (xpDiag.is_null()) {
300 RCP<const Epetra_Comm> comm = Thyra::get_Epetra_Comm(*rcp_dynamic_cast<const ThyVSBase>(thyM->range())->getComm());
301 RCP<const Epetra_Map> map = Thyra::get_Epetra_Map(*(thyM->range()), comm);
302 if (!map.is_null()) {
303 RCP<const Epetra_Vector> eDiag = Thyra::get_Epetra_Vector(*map, diag);
304 RCP<Epetra_Vector> nceDiag = rcp_const_cast<
Epetra_Vector>(eDiag);
306 xpDiag = rcp_dynamic_cast<XpVec>(xpEpDiag,
true);
312 paramList.set<RCP<XpMat> >(parameterName, M);
314 }
else if (paramList.isType<RCP<const ThyLinOpBase> >(parameterName)) {
315 RCP<const ThyLinOpBase> thyM = paramList.get<RCP<const ThyLinOpBase> >(parameterName);
316 paramList.remove(parameterName);
318 RCP<XpMat> M = XpThyUtils::toXpetra(Teuchos::rcp_const_cast<ThyLinOpBase>(thyM));
319 paramList.set<RCP<XpMat> >(parameterName, M);
320 }
catch (std::exception& e) {
321 RCP<XpOp> M = XpThyUtils::toXpetraOperator(Teuchos::rcp_const_cast<ThyLinOpBase>(thyM));
325 if (tO->hasDiagonal()) {
326 diag =
rcp(
new tV(tO->getRangeMap()));
327 tO->getLocalDiagCopy(*diag);
331 auto op = rcp_dynamic_cast<XpOp>(tpFOp);
332 paramList.set<RCP<XpOp> >(parameterName, op);
346 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
347 MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MueLuPreconditionerFactory()
348 : paramList_(
rcp(new ParameterList())) {}
350 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
351 MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::~MueLuPreconditionerFactory() =
default;
355 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
356 bool MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isCompatible(
const LinearOpSourceBase<Scalar>& fwdOpSrc)
const {
357 const RCP<const LinearOpBase<Scalar> > fwdOp = fwdOpSrc.getOp();
359 if (Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isTpetra(fwdOp))
return true;
361 #ifdef HAVE_MUELU_EPETRA
362 if (Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isEpetra(fwdOp))
return true;
365 if (Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isBlockedOperator(fwdOp))
return true;
370 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
371 RCP<PreconditionerBase<Scalar> > MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::createPrec()
const {
375 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
376 void MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
377 initializePrec(
const RCP<
const LinearOpSourceBase<Scalar> >& fwdOpSrc, PreconditionerBase<Scalar>* prec,
const ESupportSolveUse supportSolveUse)
const {
380 using Teuchos::rcp_dynamic_cast;
386 typedef Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpThyUtils;
392 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
397 #if defined(MUELU_CAN_USE_MIXED_PRECISION)
398 typedef Xpetra::TpetraHalfPrecisionOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpHalfPrecOp;
399 typedef typename XpHalfPrecOp::HalfScalar HalfScalar;
414 ParameterList paramList = *paramList_;
417 const RCP<const ThyLinOpBase> fwdOp = fwdOpSrc->getOp();
421 bool bIsEpetra = XpThyUtils::isEpetra(fwdOp);
422 bool bIsTpetra = XpThyUtils::isTpetra(fwdOp);
423 bool bIsBlocked = XpThyUtils::isBlockedOperator(fwdOp);
428 RCP<XpMat> A = Teuchos::null;
431 Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<Scalar> >(fwdOp);
441 RCP<XpMat> A00 = XpThyUtils::toXpetra(Teuchos::rcp_const_cast<LinearOpBase<Scalar> >(b00));
444 RCP<const XpMap> rowmap00 = A00->getRowMap();
445 RCP<const Teuchos::Comm<int> > comm = rowmap00->getComm();
448 RCP<XpBlockedCrsMat> bMat =
Teuchos::rcp(
new XpBlockedCrsMat(ThyBlockedOp, comm));
456 A = XpThyUtils::toXpetra(Teuchos::rcp_const_cast<ThyLinOpBase>(fwdOp));
465 RCP<ThyLinOpBase> thyra_precOp = Teuchos::null;
466 thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar> >(defaultPrec->getNonconstUnspecifiedPrecOp(),
true);
471 const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter(
"reuse: type") || paramList.get<std::string>(
"reuse: type") ==
"none");
472 bool useHalfPrecision =
false;
473 if (paramList.isParameter(
"half precision"))
474 useHalfPrecision = paramList.get<
bool>(
"half precision");
475 else if (paramList.isSublist(
"Hierarchy") && paramList.sublist(
"Hierarchy").isParameter(
"half precision"))
476 useHalfPrecision = paramList.sublist(
"Hierarchy").get<
bool>(
"half precision");
477 if (useHalfPrecision)
481 if (startingOver ==
true) {
483 std::list<std::string> convertXpetra = {
"Coordinates",
"Nullspace"};
484 for (
auto it = convertXpetra.begin(); it != convertXpetra.end(); ++it)
485 Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(paramList, *it);
487 for (
int lvlNo = 0; lvlNo < 10; ++lvlNo) {
488 if (paramList.isSublist(
"level " + std::to_string(lvlNo) +
" user data")) {
489 ParameterList& lvlList = paramList.sublist(
"level " + std::to_string(lvlNo) +
" user data");
490 std::list<std::string> convertKeys;
491 for (
auto it = lvlList.begin(); it != lvlList.end(); ++it)
492 convertKeys.push_back(lvlList.name(it));
493 for (
auto it = convertKeys.begin(); it != convertKeys.end(); ++it)
494 Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(lvlList, *it);
498 if (useHalfPrecision) {
499 #if defined(MUELU_CAN_USE_MIXED_PRECISION)
506 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
507 const std::string userName =
"user data";
509 if (userParamList.
isType<RCP<XpmMV> >(
"Coordinates")) {
510 RCP<XpmMV> coords = userParamList.
get<RCP<XpmMV> >(
"Coordinates");
511 userParamList.
remove(
"Coordinates");
512 RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
513 userParamList.
set(
"Coordinates", halfCoords);
515 if (userParamList.
isType<RCP<XpMV> >(
"Nullspace")) {
516 RCP<XpMV> nullspace = userParamList.
get<RCP<XpMV> >(
"Nullspace");
517 userParamList.
remove(
"Nullspace");
518 RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
519 userParamList.
set(
"Nullspace", halfNullspace);
521 if (paramList.isType<RCP<XpmMV> >(
"Coordinates")) {
522 RCP<XpmMV> coords = paramList.
get<RCP<XpmMV> >(
"Coordinates");
523 paramList.remove(
"Coordinates");
524 RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
525 userParamList.
set(
"Coordinates", halfCoords);
527 if (paramList.isType<RCP<XpMV> >(
"Nullspace")) {
528 RCP<XpMV> nullspace = paramList.
get<RCP<XpMV> >(
"Nullspace");
529 paramList.remove(
"Nullspace");
530 RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
531 userParamList.
set(
"Nullspace", halfNullspace);
537 RCP<MueLuHalfXpOp> xpOp =
rcp(
new MueLuHalfXpOp(H));
538 xpPrecOp =
rcp(
new XpHalfPrecOp(xpOp));
543 const std::string userName =
"user data";
545 if (paramList.isType<RCP<XpmMV> >(
"Coordinates")) {
546 RCP<XpmMV> coords = paramList.
get<RCP<XpmMV> >(
"Coordinates");
547 paramList.remove(
"Coordinates");
548 userParamList.
set(
"Coordinates", coords);
550 if (paramList.isType<RCP<XpMV> >(
"Nullspace")) {
551 RCP<XpMV> nullspace = paramList.
get<RCP<XpMV> >(
"Nullspace");
552 paramList.remove(
"Nullspace");
553 userParamList.
set(
"Nullspace", nullspace);
558 xpPrecOp =
rcp(
new MueLuXpOp(H));
562 RCP<ThyXpOp> thyXpOp = rcp_dynamic_cast<ThyXpOp>(thyra_precOp,
true);
563 xpPrecOp = rcp_dynamic_cast<XpOp>(thyXpOp->getXpetraOperator(),
true);
564 #if defined(MUELU_CAN_USE_MIXED_PRECISION)
565 RCP<XpHalfPrecOp> xpHalfPrecOp = rcp_dynamic_cast<XpHalfPrecOp>(xpPrecOp);
566 if (!xpHalfPrecOp.is_null()) {
567 RCP<MueLu::Hierarchy<HalfScalar, LocalOrdinal, GlobalOrdinal, Node> > H = rcp_dynamic_cast<MueLuHalfXpOp>(xpHalfPrecOp->GetHalfPrecisionOperator(),
true)->GetHierarchy();
568 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
571 "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it");
573 "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator");
574 RCP<MueLu::Level> level0 = H->GetLevel(0);
575 RCP<XpHalfOp> O0 = level0->Get<RCP<XpHalfOp> >(
"A");
576 RCP<XphMat> A0 = rcp_dynamic_cast<XphMat>(O0,
true);
582 halfA->SetFixedBlockSize(A0->GetFixedBlockSize());
586 level0->Set(
"A", halfA);
593 RCP<MueLuXpOp> xpOp = rcp_dynamic_cast<MueLuXpOp>(thyXpOp->getXpetraOperator(),
true);
594 RCP<MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node> > H = xpOp->GetHierarchy();
598 "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it");
600 "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator");
601 RCP<MueLu::Level> level0 = H->GetLevel(0);
602 RCP<XpOp> O0 = level0->Get<RCP<XpOp> >(
"A");
603 RCP<XpMat> A0 = rcp_dynamic_cast<XpMat>(O0);
609 A->SetFixedBlockSize(A0->GetFixedBlockSize());
620 RCP<const VectorSpaceBase<Scalar> > thyraRangeSpace = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpPrecOp->getRangeMap());
621 RCP<const VectorSpaceBase<Scalar> > thyraDomainSpace = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpPrecOp->getDomainMap());
623 RCP<ThyLinOpBase> thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace, xpPrecOp);
626 defaultPrec->initializeUnspecified(thyraPrecOp);
629 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
630 void MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
631 uninitializePrec(PreconditionerBase<Scalar>* prec, RCP<
const LinearOpSourceBase<Scalar> >* fwdOp, ESupportSolveUse* supportSolveUse)
const {
640 *fwdOp = Teuchos::null;
643 if (supportSolveUse) {
645 *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
648 defaultPrec->uninitialize();
652 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
653 void MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setParameterList(RCP<ParameterList>
const& paramList) {
655 paramList_ = paramList;
658 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
659 RCP<ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getNonconstParameterList() {
663 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
664 RCP<ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::unsetParameterList() {
665 RCP<ParameterList> savedParamList = paramList_;
666 paramList_ = Teuchos::null;
667 return savedParamList;
670 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
671 RCP<const ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getParameterList()
const {
675 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
676 RCP<const ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getValidParameters()
const {
677 static RCP<const ParameterList> validPL;
687 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
688 std::string MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::description()
const {
689 return "Thyra::MueLuPreconditionerFactory";
693 #endif // HAVE_MUELU_STRATIMIKOS
695 #endif // ifdef THYRA_MUELU_PRECONDITIONER_FACTORY_DEF_HPP
MueLu::DefaultLocalOrdinal LocalOrdinal
bool is_null(const boost::shared_ptr< T > &p)
T & get(const std::string &name, T def_value)
bool nonnull(const std::shared_ptr< T > &p)
bool is_null(const std::shared_ptr< T > &p)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
static RCP< Time > getNewTimer(const std::string &name)
Teuchos::RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > CreateXpetraPreconditioner(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> op, const Teuchos::ParameterList &inParamList)
Helper function to create a MueLu preconditioner that can be used by Xpetra.Given an Xpetra::Matrix...
void deep_copy(MultiVector< DS, DL, DG, DN > &dst, const MultiVector< SS, SL, SG, SN > &src)
bool remove(std::string const &name, bool throwIfNotExists=true)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultScalar Scalar
Concrete Thyra::LinearOpBase subclass for Xpetra::Operator.
RCP< const CrsGraph< int, GlobalOrdinal, Node > > toXpetra(const Epetra_CrsGraph &g)
bool isType(const std::string &name) const
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
static RCP< Matrix > Build(const RCP< const Map > &rowMap, size_t maxNumEntriesPerRow, Xpetra::ProfileType pftype=Xpetra::DynamicProfile)
Exception throws to report errors in the internal logical of the program.
#define TEUCHOS_ASSERT(assertion_test)
static Teuchos::RCP< const Teuchos::ParameterList > List()
Return a "master" list of all valid parameters and their default values.
RCP< Tpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getOperator()
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
Wraps an existing MueLu::Hierarchy as a Xpetra::Operator.