46#ifndef THYRA_MUELU_PRECONDITIONER_FACTORY_DEF_HPP
47#define THYRA_MUELU_PRECONDITIONER_FACTORY_DEF_HPP
52#if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA)
55#if (defined(HAVE_MUELU_TPETRA) && \
56 ((defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) && !defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) && !defined(HAVE_TPETRA_INST_COMPLEX_FLOAT)) || \
57 (!defined(HAVE_TPETRA_INST_DOUBLE) && !defined(HAVE_TPETRA_INST_FLOAT) && defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) && defined(HAVE_TPETRA_INST_COMPLEX_FLOAT)) || \
58 (defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) && defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) && defined(HAVE_TPETRA_INST_COMPLEX_FLOAT))))
59# define MUELU_CAN_USE_MIXED_PRECISION
66 using Teuchos::ParameterList;
67 using Teuchos::rcp_dynamic_cast;
68 using Teuchos::rcp_const_cast;
71 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
72 bool Converters<Scalar,LocalOrdinal,GlobalOrdinal,Node>::replaceWithXpetra(ParameterList& paramList, std::string parameterName) {
73 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
74 typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
75 typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpThyUtils;
78 typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
79 typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMultVec;
80 typedef Xpetra::MultiVector<Magnitude,LocalOrdinal,GlobalOrdinal,Node> XpMagMultVec;
81 typedef Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpVec;
83 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
84 typedef Thyra::DiagonalLinearOpBase<Scalar> ThyDiagLinOpBase;
88#ifdef HAVE_MUELU_TPETRA
89 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpCrsMat;
90 typedef Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> tOp;
91 typedef Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tV;
92 typedef Thyra::TpetraVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> thyTpV;
93 typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tMV;
94 typedef Tpetra::MultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node> tMagMV;
95# if defined(MUELU_CAN_USE_MIXED_PRECISION)
96 typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
97 typedef Tpetra::MultiVector<HalfMagnitude, LocalOrdinal, GlobalOrdinal, Node> tHalfMagMV;
101 if (paramList.isParameter(parameterName)) {
102 if (paramList.isType<RCP<XpMat> >(parameterName))
104 else if (paramList.isType<RCP<const XpMat> >(parameterName)) {
105 RCP<const XpMat> constM = paramList.get<RCP<const XpMat> >(parameterName);
106 paramList.remove(parameterName);
107 RCP<XpMat> M = rcp_const_cast<XpMat>(constM);
108 paramList.set<RCP<XpMat> >(parameterName, M);
111 else if (paramList.isType<RCP<XpMultVec> >(parameterName))
113 else if (paramList.isType<RCP<const XpMultVec> >(parameterName)) {
114 RCP<const XpMultVec> constX = paramList.get<RCP<const XpMultVec> >(parameterName);
115 paramList.remove(parameterName);
116 RCP<XpMultVec> X = rcp_const_cast<XpMultVec>(constX);
117 paramList.set<RCP<XpMultVec> >(parameterName, X);
120 else if (paramList.isType<RCP<XpMagMultVec> >(parameterName))
122 else if (paramList.isType<RCP<const XpMagMultVec> >(parameterName)) {
123 RCP<const XpMagMultVec> constX = paramList.get<RCP<const XpMagMultVec> >(parameterName);
124 paramList.remove(parameterName);
125 RCP<XpMagMultVec> X = rcp_const_cast<XpMagMultVec>(constX);
126 paramList.set<RCP<XpMagMultVec> >(parameterName, X);
129#ifdef HAVE_MUELU_TPETRA
130 else if (paramList.isType<RCP<TpCrsMat> >(parameterName)) {
131 RCP<TpCrsMat> tM = paramList.get<RCP<TpCrsMat> >(parameterName);
132 paramList.remove(parameterName);
133 RCP<XpMat> xM = MueLu::TpetraCrs_To_XpetraMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tM);
134 paramList.set<RCP<XpMat> >(parameterName, xM);
136 }
else if (paramList.isType<RCP<tMV> >(parameterName)) {
137 RCP<tMV> tpetra_X = paramList.get<RCP<tMV> >(parameterName);
138 paramList.remove(parameterName);
139 RCP<XpMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tpetra_X);
140 paramList.set<RCP<XpMultVec> >(parameterName, X);
141 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X));
143 }
else if (paramList.isType<RCP<tMagMV> >(parameterName)) {
144 RCP<tMagMV> tpetra_X = paramList.get<RCP<tMagMV> >(parameterName);
145 paramList.remove(parameterName);
146 RCP<XpMagMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Magnitude,LocalOrdinal,GlobalOrdinal,Node>(tpetra_X);
147 paramList.set<RCP<XpMagMultVec> >(parameterName, X);
148 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X));
151#if defined(MUELU_CAN_USE_MIXED_PRECISION)
152 else if (paramList.isType<RCP<tHalfMagMV> >(parameterName)) {
153 RCP<tHalfMagMV> tpetra_hX = paramList.get<RCP<tHalfMagMV> >(parameterName);
154 paramList.remove(parameterName);
155 RCP<tMagMV> tpetra_X = rcp(
new tMagMV(tpetra_hX->getMap(),tpetra_hX->getNumVectors()));
156 Tpetra::deep_copy(*tpetra_X,*tpetra_hX);
157 RCP<XpMagMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Magnitude,LocalOrdinal,GlobalOrdinal,Node>(tpetra_X);
158 paramList.set<RCP<XpMagMultVec> >(parameterName, X);
159 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X));
164 else if (paramList.isType<RCP<const ThyDiagLinOpBase> >(parameterName) ||
165 (paramList.isType<RCP<const ThyLinOpBase> >(parameterName) && !
166 rcp_dynamic_cast<const ThyDiagLinOpBase>(paramList.get<RCP<const ThyLinOpBase> >(parameterName)).is_null())) {
167 RCP<const ThyDiagLinOpBase> thyM;
168 if (paramList.isType<RCP<const ThyDiagLinOpBase> >(parameterName))
169 thyM = paramList.get<RCP<const ThyDiagLinOpBase> >(parameterName);
171 thyM = rcp_dynamic_cast<const ThyDiagLinOpBase>(paramList.get<RCP<const ThyLinOpBase> >(parameterName),
true);
172 paramList.remove(parameterName);
173 RCP<const Thyra::VectorBase<Scalar> > diag = thyM->getDiag();
175 RCP<const XpVec> xpDiag;
176#ifdef HAVE_MUELU_TPETRA
177 if (!rcp_dynamic_cast<const thyTpV>(diag).is_null()) {
178 RCP<const tV> tDiag = Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getConstTpetraVector(diag);
179 if (!tDiag.is_null())
180 xpDiag = Xpetra::toXpetra(tDiag);
183 TEUCHOS_ASSERT(!xpDiag.is_null());
184 RCP<XpMat> M = Xpetra::MatrixFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(xpDiag);
185 paramList.set<RCP<XpMat> >(parameterName, M);
188 else if (paramList.isType<RCP<const ThyLinOpBase> >(parameterName)) {
189 RCP<const ThyLinOpBase> thyM = paramList.get<RCP<const ThyLinOpBase> >(parameterName);
190 paramList.remove(parameterName);
192 RCP<XpMat> M = XpThyUtils::toXpetra(Teuchos::rcp_const_cast<ThyLinOpBase>(thyM));
193 paramList.set<RCP<XpMat> >(parameterName, M);
194 }
catch (std::exception& e) {
195 RCP<XpOp> M = XpThyUtils::toXpetraOperator(Teuchos::rcp_const_cast<ThyLinOpBase>(thyM));
196 RCP<Xpetra::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tpOp = rcp_dynamic_cast<Xpetra::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(M,
true);
197 RCP<tOp> tO = tpOp->getOperator();
199 if (tO->hasDiagonal()) {
200 diag = rcp(
new tV(tO->getRangeMap()));
201 tO->getLocalDiagCopy(*diag);
204 RCP<Xpetra::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tpFOp = rcp(
new Xpetra::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> (fTpRow));
205 auto op = rcp_dynamic_cast<XpOp>(tpFOp);
206 paramList.set<RCP<XpOp> >(parameterName, op);
219#ifdef HAVE_MUELU_EPETRA
220 template <
class GlobalOrdinal>
221 bool Converters<double, int, GlobalOrdinal, Kokkos::Compat::KokkosSerialWrapperNode>::replaceWithXpetra(ParameterList& paramList, std::string parameterName) {
224 typedef Kokkos::Compat::KokkosSerialWrapperNode
Node;
225 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
226 typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
227 typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpThyUtils;
228 typedef Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMatWrap;
229 typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
230 typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
231 typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMultVec;
232 typedef Xpetra::MultiVector<Magnitude,LocalOrdinal,GlobalOrdinal,Node> XpMagMultVec;
233 typedef Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpVec;
235 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
236 typedef Thyra::DiagonalLinearOpBase<Scalar> ThyDiagLinOpBase;
237 typedef Thyra::SpmdVectorSpaceBase<Scalar> ThyVSBase;
239#ifdef HAVE_MUELU_TPETRA
240 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpCrsMat;
241 typedef Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> tOp;
242 typedef Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tV;
243 typedef Thyra::TpetraVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> thyTpV;
244 typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tMV;
245 typedef Tpetra::MultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node> tMagMV;
246#if defined(MUELU_CAN_USE_MIXED_PRECISION)
247 typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
248 typedef Tpetra::MultiVector<HalfMagnitude, LocalOrdinal, GlobalOrdinal, Node> tHalfMagMV;
251#if defined(HAVE_MUELU_EPETRA)
252 typedef Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> XpEpCrsMat;
255 if (paramList.isParameter(parameterName)) {
256 if (paramList.isType<RCP<XpMat> >(parameterName))
258 else if (paramList.isType<RCP<const XpMat> >(parameterName)) {
259 RCP<const XpMat> constM = paramList.get<RCP<const XpMat> >(parameterName);
260 paramList.remove(parameterName);
261 RCP<XpMat> M = rcp_const_cast<XpMat>(constM);
262 paramList.set<RCP<XpMat> >(parameterName, M);
265 else if (paramList.isType<RCP<XpMultVec> >(parameterName))
267 else if (paramList.isType<RCP<const XpMultVec> >(parameterName)) {
268 RCP<const XpMultVec> constX = paramList.get<RCP<const XpMultVec> >(parameterName);
269 paramList.remove(parameterName);
270 RCP<XpMultVec> X = rcp_const_cast<XpMultVec>(constX);
271 paramList.set<RCP<XpMultVec> >(parameterName, X);
274 else if (paramList.isType<RCP<XpMagMultVec> >(parameterName))
276 else if (paramList.isType<RCP<const XpMagMultVec> >(parameterName)) {
277 RCP<const XpMagMultVec> constX = paramList.get<RCP<const XpMagMultVec> >(parameterName);
278 paramList.remove(parameterName);
279 RCP<XpMagMultVec> X = rcp_const_cast<XpMagMultVec>(constX);
280 paramList.set<RCP<XpMagMultVec> >(parameterName, X);
283#ifdef HAVE_MUELU_TPETRA
284 else if (paramList.isType<RCP<TpCrsMat> >(parameterName)) {
285 RCP<TpCrsMat> tM = paramList.get<RCP<TpCrsMat> >(parameterName);
286 paramList.remove(parameterName);
287 RCP<XpMat> xM = MueLu::TpetraCrs_To_XpetraMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tM);
288 paramList.set<RCP<XpMat> >(parameterName, xM);
290 }
else if (paramList.isType<RCP<tMV> >(parameterName)) {
291 RCP<tMV> tpetra_X = paramList.get<RCP<tMV> >(parameterName);
292 paramList.remove(parameterName);
293 RCP<XpMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tpetra_X);
294 paramList.set<RCP<XpMultVec> >(parameterName, X);
295 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X));
297 }
else if (paramList.isType<RCP<tMagMV> >(parameterName)) {
298 RCP<tMagMV> tpetra_X = paramList.get<RCP<tMagMV> >(parameterName);
299 paramList.remove(parameterName);
300 RCP<XpMagMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Magnitude,LocalOrdinal,GlobalOrdinal,Node>(tpetra_X);
301 paramList.set<RCP<XpMagMultVec> >(parameterName, X);
302 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X));
305#if defined(MUELU_CAN_USE_MIXED_PRECISION)
306 else if (paramList.isType<RCP<tHalfMagMV> >(parameterName)) {
307 RCP<tHalfMagMV> tpetra_hX = paramList.get<RCP<tHalfMagMV> >(parameterName);
308 paramList.remove(parameterName);
309 RCP<tMagMV> tpetra_X = rcp(
new tMagMV(tpetra_hX->getMap(),tpetra_hX->getNumVectors()));
310 Tpetra::deep_copy(*tpetra_X,*tpetra_hX);
311 RCP<XpMagMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Magnitude,LocalOrdinal,GlobalOrdinal,Node>(tpetra_X);
312 paramList.set<RCP<XpMagMultVec> >(parameterName, X);
313 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X));
318#ifdef HAVE_MUELU_EPETRA
319 else if (paramList.isType<RCP<Epetra_CrsMatrix> >(parameterName)) {
320 RCP<Epetra_CrsMatrix> eM = paramList.get<RCP<Epetra_CrsMatrix> >(parameterName);
321 paramList.remove(parameterName);
322 RCP<XpEpCrsMat> xeM = rcp(
new XpEpCrsMat(eM));
323 RCP<XpCrsMat> xCrsM = rcp_dynamic_cast<XpCrsMat>(xeM,
true);
324 RCP<XpCrsMatWrap> xwM = rcp(
new XpCrsMatWrap(xCrsM));
325 RCP<XpMat> xM = rcp_dynamic_cast<XpMat>(xwM);
326 paramList.set<RCP<XpMat> >(parameterName, xM);
328 }
else if (paramList.isType<RCP<Epetra_MultiVector> >(parameterName)) {
329 RCP<Epetra_MultiVector> epetra_X = Teuchos::null;
330 epetra_X = paramList.get<RCP<Epetra_MultiVector> >(parameterName);
331 paramList.remove(parameterName);
332 RCP<Xpetra::EpetraMultiVectorT<int,Node> > xpEpX = rcp(
new Xpetra::EpetraMultiVectorT<int,Node>(epetra_X));
333 RCP<Xpetra::MultiVector<Scalar,int,int,Node> > xpEpXMult = rcp_dynamic_cast<Xpetra::MultiVector<Scalar,int,int,Node> >(xpEpX,
true);
334 RCP<XpMultVec> X = rcp_dynamic_cast<XpMultVec>(xpEpXMult,
true);
335 paramList.set<RCP<XpMultVec> >(parameterName, X);
339 else if (paramList.isType<RCP<const ThyDiagLinOpBase> >(parameterName) ||
340 (paramList.isType<RCP<const ThyLinOpBase> >(parameterName) && !
341 rcp_dynamic_cast<const ThyDiagLinOpBase>(paramList.get<RCP<const ThyLinOpBase> >(parameterName)).is_null())) {
342 RCP<const ThyDiagLinOpBase> thyM;
343 if (paramList.isType<RCP<const ThyDiagLinOpBase> >(parameterName))
344 thyM = paramList.get<RCP<const ThyDiagLinOpBase> >(parameterName);
346 thyM = rcp_dynamic_cast<const ThyDiagLinOpBase>(paramList.get<RCP<const ThyLinOpBase> >(parameterName),
true);
347 paramList.remove(parameterName);
348 RCP<const Thyra::VectorBase<Scalar> > diag = thyM->getDiag();
350 RCP<const XpVec> xpDiag;
351#ifdef HAVE_MUELU_TPETRA
352 if (!rcp_dynamic_cast<const thyTpV>(diag).is_null()) {
353 RCP<const tV> tDiag = Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getConstTpetraVector(diag);
354 if (!tDiag.is_null())
355 xpDiag = Xpetra::toXpetra(tDiag);
358#ifdef HAVE_MUELU_EPETRA
359 if (xpDiag.is_null()) {
360 RCP<const Epetra_Comm> comm = Thyra::get_Epetra_Comm(*rcp_dynamic_cast<const ThyVSBase>(thyM->range())->getComm());
361 RCP<const Epetra_Map> map = Thyra::get_Epetra_Map(*(thyM->range()), comm);
362 if (!map.is_null()) {
363 RCP<const Epetra_Vector> eDiag = Thyra::get_Epetra_Vector(*map, diag);
364 RCP<Epetra_Vector> nceDiag = rcp_const_cast<Epetra_Vector>(eDiag);
365 RCP<Xpetra::EpetraVectorT<int,Node> > xpEpDiag = rcp(
new Xpetra::EpetraVectorT<int,Node>(nceDiag));
366 xpDiag = rcp_dynamic_cast<XpVec>(xpEpDiag,
true);
370 TEUCHOS_ASSERT(!xpDiag.is_null());
371 RCP<XpMat> M = Xpetra::MatrixFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(xpDiag);
372 paramList.set<RCP<XpMat> >(parameterName, M);
375 else if (paramList.isType<RCP<const ThyLinOpBase> >(parameterName)) {
376 RCP<const ThyLinOpBase> thyM = paramList.get<RCP<const ThyLinOpBase> >(parameterName);
377 paramList.remove(parameterName);
379 RCP<XpMat> M = XpThyUtils::toXpetra(Teuchos::rcp_const_cast<ThyLinOpBase>(thyM));
380 paramList.set<RCP<XpMat> >(parameterName, M);
381 }
catch (std::exception& e) {
382 RCP<XpOp> M = XpThyUtils::toXpetraOperator(Teuchos::rcp_const_cast<ThyLinOpBase>(thyM));
383 RCP<Xpetra::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tpOp = rcp_dynamic_cast<Xpetra::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(M,
true);
384 RCP<tOp> tO = tpOp->getOperator();
386 if (tO->hasDiagonal()) {
387 diag = rcp(
new tV(tO->getRangeMap()));
388 tO->getLocalDiagCopy(*diag);
391 RCP<Xpetra::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tpFOp = rcp(
new Xpetra::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> (fTpRow));
392 auto op = rcp_dynamic_cast<XpOp>(tpFOp);
393 paramList.set<RCP<XpOp> >(parameterName, op);
408 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
409 MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::MueLuPreconditionerFactory() :
410 paramList_(rcp(new ParameterList()))
415 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
416 bool MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isCompatible(
const LinearOpSourceBase<Scalar>& fwdOpSrc)
const {
417 const RCP<const LinearOpBase<Scalar> > fwdOp = fwdOpSrc.getOp();
419#ifdef HAVE_MUELU_TPETRA
420 if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isTpetra(fwdOp))
return true;
423#ifdef HAVE_MUELU_EPETRA
424 if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isEpetra(fwdOp))
return true;
428 if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isBlockedOperator(fwdOp))
return true;
434 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
435 RCP<PreconditionerBase<Scalar> > MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::createPrec()
const {
436 return Teuchos::rcp(
new DefaultPreconditioner<Scalar>);
439 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
440 void MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::
441 initializePrec(
const RCP<
const LinearOpSourceBase<Scalar> >& fwdOpSrc, PreconditionerBase<Scalar>* prec,
const ESupportSolveUse supportSolveUse)
const {
442 Teuchos::TimeMonitor tM(*Teuchos::TimeMonitor::getNewTimer(std::string(
"ThyraMueLu::initializePrec")));
444 using Teuchos::rcp_dynamic_cast;
447 typedef Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> XpMap;
448 typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
450 typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpThyUtils;
452 typedef Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpBlockedCrsMat;
453 typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
456 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
458 typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMV;
459 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
460 typedef Xpetra::MultiVector<Magnitude,LocalOrdinal,GlobalOrdinal,Node > XpmMV;
461#if defined(MUELU_CAN_USE_MIXED_PRECISION)
462 typedef Xpetra::TpetraHalfPrecisionOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpHalfPrecOp;
463 typedef typename XpHalfPrecOp::HalfScalar HalfScalar;
464 typedef Xpetra::Operator<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XpHalfOp;
466 typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
467 typedef Xpetra::MultiVector<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XphMV;
468 typedef Xpetra::MultiVector<HalfMagnitude,LocalOrdinal,GlobalOrdinal,Node > XphmMV;
469 typedef Xpetra::Matrix<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XphMat;
474 TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
475 TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
476 TEUCHOS_ASSERT(prec);
479 ParameterList paramList = *paramList_;
482 const RCP<const ThyLinOpBase> fwdOp = fwdOpSrc->getOp();
483 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
486 bool bIsEpetra = XpThyUtils::isEpetra(fwdOp);
487 bool bIsTpetra = XpThyUtils::isTpetra(fwdOp);
488 bool bIsBlocked = XpThyUtils::isBlockedOperator(fwdOp);
489 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra ==
true && bIsTpetra ==
true));
490 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == bIsTpetra) && bIsBlocked ==
false);
491 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra != bIsTpetra) && bIsBlocked ==
true);
493 RCP<XpMat> A = Teuchos::null;
495 Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
496 Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(fwdOp);
497 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(ThyBlockedOp));
499 TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==
false);
501 Teuchos::RCP<const LinearOpBase<Scalar> > b00 = ThyBlockedOp->getBlock(0,0);
502 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
506 RCP<XpMat> A00 = XpThyUtils::toXpetra(Teuchos::rcp_const_cast<LinearOpBase<Scalar> >(b00));
507 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A00));
509 RCP<const XpMap> rowmap00 = A00->getRowMap();
510 RCP< const Teuchos::Comm< int > > comm = rowmap00->getComm();
513 RCP<XpBlockedCrsMat> bMat = Teuchos::rcp(
new XpBlockedCrsMat(ThyBlockedOp, comm));
514 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(bMat));
521 A = XpThyUtils::toXpetra(Teuchos::rcp_const_cast<ThyLinOpBase>(fwdOp));
523 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A));
526 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(
dynamic_cast<DefaultPreconditioner<Scalar> *
>(prec));
527 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
530 RCP<ThyLinOpBase> thyra_precOp = Teuchos::null;
531 thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar> >(defaultPrec->getNonconstUnspecifiedPrecOp(),
true);
536 const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter(
"reuse: type") || paramList.get<std::string>(
"reuse: type") ==
"none");
537 bool useHalfPrecision =
false;
538 if (paramList.isParameter(
"half precision"))
539 useHalfPrecision = paramList.get<
bool>(
"half precision");
540 else if (paramList.isSublist(
"Hierarchy") && paramList.sublist(
"Hierarchy").isParameter(
"half precision"))
541 useHalfPrecision = paramList.sublist(
"Hierarchy").get<
bool>(
"half precision");
542 if (useHalfPrecision)
543 TEUCHOS_TEST_FOR_EXCEPTION(!bIsTpetra,
MueLu::Exceptions::RuntimeError,
"The only scalar type Epetra knows is double, so a half precision preconditioner cannot be constructed.");
546 if (startingOver ==
true) {
548 std::list<std::string> convertXpetra = {
"Coordinates",
"Nullspace"};
549 for (
auto it = convertXpetra.begin(); it != convertXpetra.end(); ++it)
550 Converters<Scalar,LocalOrdinal,GlobalOrdinal,Node>::replaceWithXpetra(paramList,*it);
552 for (
int lvlNo=0; lvlNo < 10; ++lvlNo) {
553 if (paramList.isSublist(
"level " + std::to_string(lvlNo) +
" user data")) {
554 ParameterList& lvlList = paramList.sublist(
"level " + std::to_string(lvlNo) +
" user data");
555 std::list<std::string> convertKeys;
556 for (
auto it = lvlList.begin(); it != lvlList.end(); ++it)
557 convertKeys.push_back(lvlList.name(it));
558 for (
auto it = convertKeys.begin(); it != convertKeys.end(); ++it)
559 Converters<Scalar,LocalOrdinal,GlobalOrdinal,Node>::replaceWithXpetra(lvlList,*it);
563 if (useHalfPrecision) {
564#if defined(MUELU_CAN_USE_MIXED_PRECISION)
571 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
572 const std::string userName =
"user data";
573 Teuchos::ParameterList& userParamList = paramList.sublist(userName);
574 if (userParamList.isType<RCP<XpmMV> >(
"Coordinates")) {
575 RCP<XpmMV> coords = userParamList.get<RCP<XpmMV> >(
"Coordinates");
576 userParamList.remove(
"Coordinates");
577 RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
578 userParamList.set(
"Coordinates",halfCoords);
580 if (userParamList.isType<RCP<XpMV> >(
"Nullspace")) {
581 RCP<XpMV> nullspace = userParamList.get<RCP<XpMV> >(
"Nullspace");
582 userParamList.remove(
"Nullspace");
583 RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
584 userParamList.set(
"Nullspace",halfNullspace);
586 if (paramList.isType<RCP<XpmMV> >(
"Coordinates")) {
587 RCP<XpmMV> coords = paramList.get<RCP<XpmMV> >(
"Coordinates");
588 paramList.remove(
"Coordinates");
589 RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
590 userParamList.set(
"Coordinates",halfCoords);
592 if (paramList.isType<RCP<XpMV> >(
"Nullspace")) {
593 RCP<XpMV> nullspace = paramList.get<RCP<XpMV> >(
"Nullspace");
594 paramList.remove(
"Nullspace");
595 RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
596 userParamList.set(
"Nullspace",halfNullspace);
603 RCP<MueLuHalfXpOp> xpOp = rcp(
new MueLuHalfXpOp(H));
604 xpPrecOp = rcp(
new XpHalfPrecOp(xpOp));
606 TEUCHOS_TEST_FOR_EXCEPT(
true);
610 const std::string userName =
"user data";
611 Teuchos::ParameterList& userParamList = paramList.sublist(userName);
612 if (paramList.isType<RCP<XpmMV> >(
"Coordinates")) {
613 RCP<XpmMV> coords = paramList.get<RCP<XpmMV> >(
"Coordinates");
614 paramList.remove(
"Coordinates");
615 userParamList.set(
"Coordinates",coords);
617 if (paramList.isType<RCP<XpMV> >(
"Nullspace")) {
618 RCP<XpMV> nullspace = paramList.get<RCP<XpMV> >(
"Nullspace");
619 paramList.remove(
"Nullspace");
620 userParamList.set(
"Nullspace",nullspace);
625 xpPrecOp = rcp(
new MueLuXpOp(H));
629 RCP<ThyXpOp> thyXpOp = rcp_dynamic_cast<ThyXpOp>(thyra_precOp,
true);
630 xpPrecOp = rcp_dynamic_cast<XpOp>(thyXpOp->getXpetraOperator(),
true);
631#if defined(MUELU_CAN_USE_MIXED_PRECISION)
632 RCP<XpHalfPrecOp> xpHalfPrecOp = rcp_dynamic_cast<XpHalfPrecOp>(xpPrecOp);
633 if (!xpHalfPrecOp.is_null()) {
634 RCP<MueLu::Hierarchy<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> > H = rcp_dynamic_cast<MueLuHalfXpOp>(xpHalfPrecOp->GetHalfPrecisionOperator(),
true)->GetHierarchy();
635 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
638 "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it");
640 "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator");
641 RCP<MueLu::Level> level0 = H->GetLevel(0);
642 RCP<XpHalfOp> O0 = level0->Get<RCP<XpHalfOp> >(
"A");
643 RCP<XphMat> A0 = rcp_dynamic_cast<XphMat>(O0,
true);
649 halfA->SetFixedBlockSize(A0->GetFixedBlockSize());
653 level0->Set(
"A", halfA);
660 RCP<MueLuXpOp> xpOp = rcp_dynamic_cast<MueLuXpOp>(thyXpOp->getXpetraOperator(),
true);
661 RCP<MueLu::Hierarchy<Scalar,LocalOrdinal,GlobalOrdinal,Node> > H = xpOp->GetHierarchy();;
664 "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it");
666 "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator");
667 RCP<MueLu::Level> level0 = H->GetLevel(0);
668 RCP<XpOp> O0 = level0->Get<RCP<XpOp> >(
"A");
669 RCP<XpMat> A0 = rcp_dynamic_cast<XpMat>(O0);
675 A->SetFixedBlockSize(A0->GetFixedBlockSize());
686 RCP<const VectorSpaceBase<Scalar> > thyraRangeSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpPrecOp->getRangeMap());
687 RCP<const VectorSpaceBase<Scalar> > thyraDomainSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpPrecOp->getDomainMap());
689 RCP<ThyLinOpBase > thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace, xpPrecOp);
690 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraPrecOp));
692 defaultPrec->initializeUnspecified(thyraPrecOp);
696 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
697 void MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::
698 uninitializePrec(PreconditionerBase<Scalar>* prec, RCP<
const LinearOpSourceBase<Scalar> >* fwdOp, ESupportSolveUse* supportSolveUse)
const {
699 TEUCHOS_ASSERT(prec);
702 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(
dynamic_cast<DefaultPreconditioner<Scalar> *
>(prec));
703 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
707 *fwdOp = Teuchos::null;
710 if (supportSolveUse) {
712 *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
715 defaultPrec->uninitialize();
720 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
721 void MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::setParameterList(RCP<ParameterList>
const& paramList) {
722 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
723 paramList_ = paramList;
726 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
727 RCP<ParameterList> MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getNonconstParameterList() {
731 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
732 RCP<ParameterList> MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::unsetParameterList() {
733 RCP<ParameterList> savedParamList = paramList_;
734 paramList_ = Teuchos::null;
735 return savedParamList;
738 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
739 RCP<const ParameterList> MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getParameterList()
const {
743 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
744 RCP<const ParameterList> MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getValidParameters()
const {
745 static RCP<const ParameterList> validPL;
747 if (Teuchos::is_null(validPL))
748 validPL = rcp(
new ParameterList());
754 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
755 std::string MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::description()
const {
756 return "Thyra::MueLuPreconditionerFactory";
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultScalar Scalar
Exception throws to report errors in the internal logical of the program.
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,...