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
28using Teuchos::ParameterList;
31using Teuchos::rcp_const_cast;
32using Teuchos::rcp_dynamic_cast;
34template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
35bool Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(ParameterList& paramList, std::string parameterName) {
36 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
37 typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
38 typedef Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpThyUtils;
41 typedef Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpMat;
42 typedef Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpMultVec;
43 typedef Xpetra::MultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node> XpMagMultVec;
44 typedef Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpVec;
46 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
47 typedef Thyra::DiagonalLinearOpBase<Scalar> ThyDiagLinOpBase;
51 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpCrsMat;
52 typedef Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> tOp;
53 typedef Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tV;
54 typedef Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> thyTpV;
55 typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tMV;
56 typedef Tpetra::MultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node> tMagMV;
57#if defined(MUELU_CAN_USE_MIXED_PRECISION)
58 typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
59 typedef Tpetra::MultiVector<HalfMagnitude, LocalOrdinal, GlobalOrdinal, Node> tHalfMagMV;
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 = Xpetra::toXpetra(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 = Xpetra::toXpetra(tpetra_X);
97 paramList.set<RCP<XpMultVec> >(parameterName, X);
98 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(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 = Xpetra::toXpetra(tpetra_X);
104 paramList.set<RCP<XpMagMultVec> >(parameterName, X);
105 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(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()));
113 Tpetra::deep_copy(*tpetra_X, *tpetra_hX);
114 RCP<XpMagMultVec> X = Xpetra::toXpetra(tpetra_X);
115 paramList.set<RCP<XpMagMultVec> >(parameterName, X);
116 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(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())
134 xpDiag = Xpetra::toXpetra(tDiag);
136 TEUCHOS_ASSERT(!xpDiag.is_null());
137 RCP<XpMat> M = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(xpDiag);
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));
148 RCP<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpOp = rcp_dynamic_cast<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(M,
true);
149 RCP<tOp> tO = tpOp->getOperator();
151 if (tO->hasDiagonal()) {
152 diag = rcp(
new tV(tO->getRangeMap()));
153 tO->getLocalDiagCopy(*diag);
156 RCP<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpFOp = rcp(
new Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>(fTpRow));
157 auto op = rcp_dynamic_cast<XpOp>(tpFOp);
158 paramList.set<RCP<XpOp> >(parameterName, op);
169#ifdef HAVE_MUELU_EPETRA
170template <
class GlobalOrdinal>
171bool Converters<double, int, GlobalOrdinal, Tpetra::KokkosCompat::KokkosSerialWrapperNode>::replaceWithXpetra(ParameterList& paramList, std::string parameterName) {
174 typedef Tpetra::KokkosCompat::KokkosSerialWrapperNode
Node;
175 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
176 typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
177 typedef Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpThyUtils;
178 typedef Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpCrsMatWrap;
179 typedef Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpCrsMat;
180 typedef Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpMat;
181 typedef Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpMultVec;
182 typedef Xpetra::MultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node> XpMagMultVec;
183 typedef Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpVec;
185 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
186 typedef Thyra::DiagonalLinearOpBase<Scalar> ThyDiagLinOpBase;
187 typedef Thyra::SpmdVectorSpaceBase<Scalar> ThyVSBase;
189 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpCrsMat;
190 typedef Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> tOp;
191 typedef Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tV;
192 typedef Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> thyTpV;
193 typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tMV;
194 typedef Tpetra::MultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node> tMagMV;
195#if defined(MUELU_CAN_USE_MIXED_PRECISION)
196 typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
197 typedef Tpetra::MultiVector<HalfMagnitude, LocalOrdinal, GlobalOrdinal, Node> tHalfMagMV;
199#if defined(HAVE_MUELU_EPETRA)
200 typedef Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node> XpEpCrsMat;
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 = Xpetra::toXpetra(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 = Xpetra::toXpetra(tpetra_X);
238 paramList.set<RCP<XpMultVec> >(parameterName, X);
239 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(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 = Xpetra::toXpetra(tpetra_X);
245 paramList.set<RCP<XpMagMultVec> >(parameterName, X);
246 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(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()));
254 Tpetra::deep_copy(*tpetra_X, *tpetra_hX);
255 RCP<XpMagMultVec> X = Xpetra::toXpetra(tpetra_X);
256 paramList.set<RCP<XpMagMultVec> >(parameterName, X);
257 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(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);
275 RCP<Xpetra::EpetraMultiVectorT<int, Node> > xpEpX = rcp(
new Xpetra::EpetraMultiVectorT<int, Node>(epetra_X));
276 RCP<Xpetra::MultiVector<Scalar, int, int, Node> > xpEpXMult = rcp_dynamic_cast<Xpetra::MultiVector<Scalar, int, int, Node> >(xpEpX,
true);
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())
296 xpDiag = Xpetra::toXpetra(tDiag);
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);
305 RCP<Xpetra::EpetraVectorT<int, Node> > xpEpDiag = rcp(
new Xpetra::EpetraVectorT<int, Node>(nceDiag));
306 xpDiag = rcp_dynamic_cast<XpVec>(xpEpDiag,
true);
310 TEUCHOS_ASSERT(!xpDiag.is_null());
311 RCP<XpMat> M = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(xpDiag);
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));
322 RCP<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpOp = rcp_dynamic_cast<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(M,
true);
323 RCP<tOp> tO = tpOp->getOperator();
325 if (tO->hasDiagonal()) {
326 diag = rcp(
new tV(tO->getRangeMap()));
327 tO->getLocalDiagCopy(*diag);
330 RCP<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpFOp = rcp(
new Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>(fTpRow));
331 auto op = rcp_dynamic_cast<XpOp>(tpFOp);
332 paramList.set<RCP<XpOp> >(parameterName, op);
346template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
347MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MueLuPreconditionerFactory()
348 : paramList_(rcp(new ParameterList())) {}
350template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
351MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::~MueLuPreconditionerFactory() =
default;
355template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
356bool 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;
370template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
371RCP<PreconditionerBase<Scalar> > MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::createPrec()
const {
372 return Teuchos::rcp(
new DefaultPreconditioner<Scalar>);
375template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
376void MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
377 initializePrec(
const RCP<
const LinearOpSourceBase<Scalar> >& fwdOpSrc, PreconditionerBase<Scalar>* prec,
const ESupportSolveUse supportSolveUse)
const {
378 Teuchos::TimeMonitor tM(*Teuchos::TimeMonitor::getNewTimer(std::string(
"ThyraMueLu::initializePrec")));
380 using Teuchos::rcp_dynamic_cast;
383 typedef Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> XpMap;
384 typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
386 typedef Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpThyUtils;
388 typedef Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpBlockedCrsMat;
389 typedef Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpMat;
392 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
394 typedef Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpMV;
395 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
396 typedef Xpetra::MultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node> XpmMV;
397#if defined(MUELU_CAN_USE_MIXED_PRECISION)
398 typedef Xpetra::TpetraHalfPrecisionOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpHalfPrecOp;
399 typedef typename XpHalfPrecOp::HalfScalar HalfScalar;
400 typedef Xpetra::Operator<HalfScalar, LocalOrdinal, GlobalOrdinal, Node> XpHalfOp;
402 typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
403 typedef Xpetra::MultiVector<HalfScalar, LocalOrdinal, GlobalOrdinal, Node> XphMV;
404 typedef Xpetra::MultiVector<HalfMagnitude, LocalOrdinal, GlobalOrdinal, Node> XphmMV;
405 typedef Xpetra::Matrix<HalfScalar, LocalOrdinal, GlobalOrdinal, Node> XphMat;
409 TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
410 TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
411 TEUCHOS_ASSERT(prec);
414 ParameterList paramList = *paramList_;
417 const RCP<const ThyLinOpBase> fwdOp = fwdOpSrc->getOp();
418 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
421 bool bIsEpetra = XpThyUtils::isEpetra(fwdOp);
422 bool bIsTpetra = XpThyUtils::isTpetra(fwdOp);
423 bool bIsBlocked = XpThyUtils::isBlockedOperator(fwdOp);
424 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra ==
true && bIsTpetra ==
true));
425 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == bIsTpetra) && bIsBlocked ==
false);
426 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra != bIsTpetra) && bIsBlocked ==
true);
428 RCP<XpMat> A = Teuchos::null;
430 Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
431 Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(fwdOp);
432 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(ThyBlockedOp));
434 TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0, 0) ==
false);
436 Teuchos::RCP<const LinearOpBase<Scalar> > b00 = ThyBlockedOp->getBlock(0, 0);
437 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
441 RCP<XpMat> A00 = XpThyUtils::toXpetra(Teuchos::rcp_const_cast<LinearOpBase<Scalar> >(b00));
442 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A00));
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));
449 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(bMat));
456 A = XpThyUtils::toXpetra(Teuchos::rcp_const_cast<ThyLinOpBase>(fwdOp));
458 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A));
461 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(
dynamic_cast<DefaultPreconditioner<Scalar>*
>(prec));
462 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
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)
478 TEUCHOS_TEST_FOR_EXCEPTION(!bIsTpetra,
MueLu::Exceptions::RuntimeError,
"The only scalar type Epetra knows is double, so a half precision preconditioner cannot be constructed.");
481 if (startingOver ==
true) {
483 std::list<std::string> convertXpetra = {
"Coordinates",
"Nullspace",
"Material",
"BlockNumber"};
484 for (
auto it = convertXpetra.begin(); it != convertXpetra.end(); ++it)
485 Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(paramList, *it);
487 if (paramList.isSublist(
"user data")) {
488 auto& userData = paramList.sublist(
"user data");
489 for (
auto it = convertXpetra.begin(); it != convertXpetra.end(); ++it)
490 Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(userData, *it);
493 for (
int lvlNo = 0; lvlNo < 10; ++lvlNo) {
494 if (paramList.isSublist(
"level " + std::to_string(lvlNo) +
" user data")) {
495 ParameterList& lvlList = paramList.sublist(
"level " + std::to_string(lvlNo) +
" user data");
496 std::list<std::string> convertKeys;
497 for (
auto it = lvlList.begin(); it != lvlList.end(); ++it)
498 convertKeys.push_back(lvlList.name(it));
499 for (
auto it = convertKeys.begin(); it != convertKeys.end(); ++it)
500 Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(lvlList, *it);
504 if (useHalfPrecision) {
505#if defined(MUELU_CAN_USE_MIXED_PRECISION)
512 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
513 const std::string userName =
"user data";
514 Teuchos::ParameterList& userParamList = paramList.sublist(userName);
515 if (userParamList.isType<RCP<XpmMV> >(
"Coordinates")) {
516 RCP<XpmMV> coords = userParamList.get<RCP<XpmMV> >(
"Coordinates");
517 userParamList.remove(
"Coordinates");
518 RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
519 userParamList.set(
"Coordinates", halfCoords);
521 if (userParamList.isType<RCP<XpMV> >(
"Nullspace")) {
522 RCP<XpMV> nullspace = userParamList.get<RCP<XpMV> >(
"Nullspace");
523 userParamList.remove(
"Nullspace");
524 RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
525 userParamList.set(
"Nullspace", halfNullspace);
527 if (paramList.isType<RCP<XpmMV> >(
"Coordinates")) {
528 RCP<XpmMV> coords = paramList.get<RCP<XpmMV> >(
"Coordinates");
529 paramList.remove(
"Coordinates");
530 RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
531 userParamList.set(
"Coordinates", halfCoords);
533 if (paramList.isType<RCP<XpMV> >(
"Nullspace")) {
534 RCP<XpMV> nullspace = paramList.get<RCP<XpMV> >(
"Nullspace");
535 paramList.remove(
"Nullspace");
536 RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
537 userParamList.set(
"Nullspace", halfNullspace);
543 RCP<MueLuHalfXpOp> xpOp = rcp(
new MueLuHalfXpOp(H));
544 xpPrecOp = rcp(
new XpHalfPrecOp(xpOp));
546 TEUCHOS_TEST_FOR_EXCEPT(
true);
549 const std::string userName =
"user data";
550 Teuchos::ParameterList& userParamList = paramList.sublist(userName);
551 if (paramList.isType<RCP<XpmMV> >(
"Coordinates")) {
552 RCP<XpmMV> coords = paramList.get<RCP<XpmMV> >(
"Coordinates");
553 paramList.remove(
"Coordinates");
554 userParamList.set(
"Coordinates", coords);
556 if (paramList.isType<RCP<XpMV> >(
"Nullspace")) {
557 RCP<XpMV> nullspace = paramList.get<RCP<XpMV> >(
"Nullspace");
558 paramList.remove(
"Nullspace");
559 userParamList.set(
"Nullspace", nullspace);
564 xpPrecOp = rcp(
new MueLuXpOp(H));
568 RCP<ThyXpOp> thyXpOp = rcp_dynamic_cast<ThyXpOp>(thyra_precOp,
true);
569 xpPrecOp = rcp_dynamic_cast<XpOp>(thyXpOp->getXpetraOperator(),
true);
570#if defined(MUELU_CAN_USE_MIXED_PRECISION)
571 RCP<XpHalfPrecOp> xpHalfPrecOp = rcp_dynamic_cast<XpHalfPrecOp>(xpPrecOp);
572 if (!xpHalfPrecOp.is_null()) {
573 RCP<MueLu::Hierarchy<HalfScalar, LocalOrdinal, GlobalOrdinal, Node> > H = rcp_dynamic_cast<MueLuHalfXpOp>(xpHalfPrecOp->GetHalfPrecisionOperator(),
true)->GetHierarchy();
574 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
577 "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it");
579 "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator");
580 RCP<MueLu::Level> level0 = H->GetLevel(0);
581 RCP<XpHalfOp> O0 = level0->Get<RCP<XpHalfOp> >(
"A");
582 RCP<XphMat> A0 = rcp_dynamic_cast<XphMat>(O0,
true);
588 halfA->SetFixedBlockSize(A0->GetFixedBlockSize());
592 level0->Set(
"A", halfA);
599 RCP<MueLuXpOp> xpOp = rcp_dynamic_cast<MueLuXpOp>(thyXpOp->getXpetraOperator(),
true);
600 RCP<MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node> > H = xpOp->GetHierarchy();
604 "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it");
606 "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator");
607 RCP<MueLu::Level> level0 = H->GetLevel(0);
608 RCP<XpOp> O0 = level0->Get<RCP<XpOp> >(
"A");
609 RCP<XpMat> A0 = rcp_dynamic_cast<XpMat>(O0);
615 A->SetFixedBlockSize(A0->GetFixedBlockSize());
626 RCP<const VectorSpaceBase<Scalar> > thyraRangeSpace = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpPrecOp->getRangeMap());
627 RCP<const VectorSpaceBase<Scalar> > thyraDomainSpace = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpPrecOp->getDomainMap());
629 RCP<ThyLinOpBase> thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace, xpPrecOp);
630 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraPrecOp));
632 defaultPrec->initializeUnspecified(thyraPrecOp);
635template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
636void MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
637 uninitializePrec(PreconditionerBase<Scalar>* prec, RCP<
const LinearOpSourceBase<Scalar> >* fwdOp, ESupportSolveUse* supportSolveUse)
const {
638 TEUCHOS_ASSERT(prec);
641 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(
dynamic_cast<DefaultPreconditioner<Scalar>*
>(prec));
642 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
646 *fwdOp = Teuchos::null;
649 if (supportSolveUse) {
651 *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
654 defaultPrec->uninitialize();
658template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
659void MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setParameterList(RCP<ParameterList>
const& paramList) {
660 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
661 paramList_ = paramList;
664template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
665RCP<ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getNonconstParameterList() {
669template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
670RCP<ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::unsetParameterList() {
671 RCP<ParameterList> savedParamList = paramList_;
672 paramList_ = Teuchos::null;
673 return savedParamList;
676template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
677RCP<const ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getParameterList()
const {
681template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
682RCP<const ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getValidParameters()
const {
683 static RCP<const ParameterList> validPL;
685 if (Teuchos::is_null(validPL)) {
693template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
694std::string MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::description()
const {
695 return "Thyra::MueLuPreconditionerFactory";
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultScalar Scalar
Exception throws to report errors in the internal logical of the program.
static Teuchos::RCP< const Teuchos::ParameterList > List()
Return a "master" list of all valid parameters and their default values.
Wraps an existing MueLu::Hierarchy as a Xpetra::Operator.
Concrete Thyra::LinearOpBase subclass for Xpetra::Operator.
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,...