46 #ifndef THYRA_MUELU_PRECONDITIONER_FACTORY_DEF_HPP 47 #define THYRA_MUELU_PRECONDITIONER_FACTORY_DEF_HPP 51 #if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA) 57 using Teuchos::ParameterList;
58 using Teuchos::rcp_dynamic_cast;
59 using Teuchos::rcp_const_cast;
62 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
63 static bool replaceWithXpetra(ParameterList& paramList, std::string parameterName) {
64 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
65 typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
66 typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpThyUtils;
67 typedef Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMatWrap;
68 typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
69 typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
70 typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMultVec;
71 typedef Xpetra::MultiVector<Magnitude,LocalOrdinal,GlobalOrdinal,Node> XpMagMultVec;
72 typedef Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpVec;
74 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
75 typedef Thyra::DiagonalLinearOpBase<Scalar> ThyDiagLinOpBase;
77 typedef Thyra::SpmdVectorSpaceBase<Scalar> ThyVSBase;
79 #ifdef HAVE_MUELU_TPETRA 80 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpCrsMat;
81 typedef Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tV;
82 typedef Thyra::TpetraVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> thyTpV;
83 typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tMV;
84 typedef Tpetra::MultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node> tMagMV;
85 # if defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) 86 typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
87 typedef Tpetra::MultiVector<HalfMagnitude, LocalOrdinal, GlobalOrdinal, Node> tHalfMagMV;
90 #if defined(HAVE_MUELU_EPETRA) 91 typedef Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> XpEpCrsMat;
94 if (paramList.isParameter(parameterName)) {
95 if (paramList.isType<RCP<XpMat> >(parameterName))
97 else if (paramList.isType<RCP<const XpMat> >(parameterName)) {
98 RCP<const XpMat> constM = paramList.get<RCP<const XpMat> >(parameterName);
99 paramList.remove(parameterName);
100 RCP<XpMat> M = rcp_const_cast<XpMat>(constM);
101 paramList.set<RCP<XpMat> >(parameterName, M);
104 else if (paramList.isType<RCP<XpMultVec> >(parameterName))
106 else if (paramList.isType<RCP<const XpMultVec> >(parameterName)) {
107 RCP<const XpMultVec> constX = paramList.get<RCP<const XpMultVec> >(parameterName);
108 paramList.remove(parameterName);
109 RCP<XpMultVec> X = rcp_const_cast<XpMultVec>(constX);
110 paramList.set<RCP<XpMultVec> >(parameterName, X);
113 else if (paramList.isType<RCP<XpMagMultVec> >(parameterName))
115 else if (paramList.isType<RCP<const XpMagMultVec> >(parameterName)) {
116 RCP<const XpMagMultVec> constX = paramList.get<RCP<const XpMagMultVec> >(parameterName);
117 paramList.remove(parameterName);
118 RCP<XpMagMultVec> X = rcp_const_cast<XpMagMultVec>(constX);
119 paramList.set<RCP<XpMagMultVec> >(parameterName, X);
122 #ifdef HAVE_MUELU_TPETRA 123 else if (paramList.isType<RCP<TpCrsMat> >(parameterName)) {
124 RCP<TpCrsMat> tM = paramList.get<RCP<TpCrsMat> >(parameterName);
125 paramList.remove(parameterName);
126 RCP<XpMat> xM = MueLu::TpetraCrs_To_XpetraMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tM);
127 paramList.set<RCP<XpMat> >(parameterName, xM);
129 }
else if (paramList.isType<RCP<tMV> >(parameterName)) {
130 RCP<tMV> tpetra_X = paramList.get<RCP<tMV> >(parameterName);
131 paramList.remove(parameterName);
132 RCP<XpMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tpetra_X);
133 paramList.set<RCP<XpMultVec> >(parameterName, X);
134 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X));
136 }
else if (paramList.isType<RCP<tMagMV> >(parameterName)) {
137 RCP<tMagMV> tpetra_X = paramList.get<RCP<tMagMV> >(parameterName);
138 paramList.remove(parameterName);
139 RCP<XpMagMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Magnitude,LocalOrdinal,GlobalOrdinal,Node>(tpetra_X);
140 paramList.set<RCP<XpMagMultVec> >(parameterName, X);
141 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X));
144 # if defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) 145 else if (paramList.isType<RCP<tHalfMagMV> >(parameterName)) {
146 RCP<tHalfMagMV> tpetra_hX = paramList.get<RCP<tHalfMagMV> >(parameterName);
147 paramList.remove(parameterName);
148 RCP<tMagMV> tpetra_X = rcp(
new tMagMV(tpetra_hX->getMap(),tpetra_hX->getNumVectors()));
149 Tpetra::deep_copy(*tpetra_X,*tpetra_hX);
150 RCP<XpMagMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Magnitude,LocalOrdinal,GlobalOrdinal,Node>(tpetra_X);
151 paramList.set<RCP<XpMagMultVec> >(parameterName, X);
152 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X));
157 #ifdef HAVE_MUELU_EPETRA 158 else if (paramList.isType<RCP<Epetra_CrsMatrix> >(parameterName)) {
159 RCP<Epetra_CrsMatrix> eM = paramList.get<RCP<Epetra_CrsMatrix> >(parameterName);
160 paramList.remove(parameterName);
161 RCP<XpEpCrsMat> xeM = rcp(
new XpEpCrsMat(eM));
162 RCP<XpCrsMat> xCrsM = rcp_dynamic_cast<XpCrsMat>(xeM,
true);
163 RCP<XpCrsMatWrap> xwM = rcp(
new XpCrsMatWrap(xCrsM));
164 RCP<XpMat> xM = rcp_dynamic_cast<XpMat>(xwM);
165 paramList.set<RCP<XpMat> >(parameterName, xM);
167 }
else if (paramList.isType<RCP<Epetra_MultiVector> >(parameterName)) {
168 RCP<Epetra_MultiVector> epetra_X = Teuchos::null;
169 epetra_X = paramList.get<RCP<Epetra_MultiVector> >(parameterName);
170 paramList.remove(parameterName);
171 RCP<Xpetra::EpetraMultiVectorT<int,Node> > xpEpX = rcp(
new Xpetra::EpetraMultiVectorT<int,Node>(epetra_X));
172 RCP<Xpetra::MultiVector<Scalar,int,int,Node> > xpEpXMult = rcp_dynamic_cast<Xpetra::MultiVector<Scalar,int,int,Node> >(xpEpX,
true);
173 RCP<XpMultVec> X = rcp_dynamic_cast<XpMultVec>(xpEpXMult,
true);
174 paramList.set<RCP<XpMultVec> >(parameterName, X);
178 else if (paramList.isType<RCP<const ThyLinOpBase> >(parameterName)) {
179 RCP<const ThyLinOpBase> thyM = paramList.get<RCP<const ThyLinOpBase> >(parameterName);
180 paramList.remove(parameterName);
181 RCP<const XpCrsMat> crsM = XpThyUtils::toXpetra(thyM);
182 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsM));
184 RCP<XpCrsMat> crsMNonConst = rcp_const_cast<XpCrsMat>(crsM);
185 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsMNonConst));
187 RCP<XpMat> M = rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(crsMNonConst));
188 paramList.set<RCP<XpMat> >(parameterName, M);
190 }
else if (paramList.isType<RCP<const ThyDiagLinOpBase> >(parameterName)) {
191 RCP<const ThyDiagLinOpBase> thyM = paramList.get<RCP<const ThyDiagLinOpBase> >(parameterName);
192 paramList.remove(parameterName);
193 RCP<const Thyra::VectorBase<Scalar> > diag = thyM->getDiag();
195 RCP<const XpVec> xpDiag;
196 #ifdef HAVE_MUELU_TPETRA 197 if (!rcp_dynamic_cast<const thyTpV>(diag).is_null()) {
198 RCP<const tV> tDiag = Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getConstTpetraVector(diag);
199 if (!tDiag.is_null())
200 xpDiag = Xpetra::toXpetra(tDiag);
203 #ifdef HAVE_MUELU_EPETRA 204 if (xpDiag.is_null()) {
205 RCP<const Epetra_Comm> comm = Thyra::get_Epetra_Comm(*rcp_dynamic_cast<const ThyVSBase>(thyM->range())->getComm());
206 RCP<const Epetra_Map> map = Thyra::get_Epetra_Map(*(thyM->range()), comm);
207 if (!map.is_null()) {
208 RCP<const Epetra_Vector> eDiag = Thyra::get_Epetra_Vector(*map, diag);
209 RCP<Epetra_Vector> nceDiag = rcp_const_cast<
Epetra_Vector>(eDiag);
210 RCP<Xpetra::EpetraVectorT<int,Node> > xpEpDiag = rcp(
new Xpetra::EpetraVectorT<int,Node>(nceDiag));
211 xpDiag = rcp_dynamic_cast<XpVec>(xpEpDiag,
true);
215 TEUCHOS_ASSERT(!xpDiag.is_null());
216 RCP<XpMat> M = Xpetra::MatrixFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(xpDiag);
217 paramList.set<RCP<XpMat> >(parameterName, M);
230 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
231 MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::MueLuPreconditionerFactory() :
232 paramList_(rcp(new ParameterList()))
237 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
238 bool MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isCompatible(
const LinearOpSourceBase<Scalar>& fwdOpSrc)
const {
239 const RCP<const LinearOpBase<Scalar> > fwdOp = fwdOpSrc.getOp();
241 #ifdef HAVE_MUELU_TPETRA 242 if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isTpetra(fwdOp))
return true;
245 #ifdef HAVE_MUELU_EPETRA 246 if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isEpetra(fwdOp))
return true;
250 if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isBlockedOperator(fwdOp))
return true;
256 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
257 RCP<PreconditionerBase<Scalar> > MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::createPrec()
const {
258 return Teuchos::rcp(
new DefaultPreconditioner<Scalar>);
261 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
262 void MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::
263 initializePrec(
const RCP<
const LinearOpSourceBase<Scalar> >& fwdOpSrc, PreconditionerBase<Scalar>* prec,
const ESupportSolveUse supportSolveUse)
const {
264 Teuchos::TimeMonitor tM(*Teuchos::TimeMonitor::getNewTimer(std::string(
"ThyraMueLu::initializePrec")));
266 using Teuchos::rcp_dynamic_cast;
269 typedef Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> XpMap;
270 typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
272 typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpThyUtils;
273 typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
274 typedef Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpBlockedCrsMat;
275 typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
276 typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMultVec;
278 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
280 typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMV;
281 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
282 typedef Xpetra::MultiVector<Magnitude,LocalOrdinal,GlobalOrdinal,Node > XpmMV;
283 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) 284 typedef Xpetra::TpetraHalfPrecisionOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpHalfPrecOp;
285 typedef typename XpHalfPrecOp::HalfScalar HalfScalar;
286 typedef Xpetra::Operator<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XpHalfOp;
288 typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
289 typedef Xpetra::MultiVector<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XphMV;
290 typedef Xpetra::MultiVector<HalfMagnitude,LocalOrdinal,GlobalOrdinal,Node > XphmMV;
291 typedef Xpetra::Matrix<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XphMat;
296 TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
297 TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
298 TEUCHOS_ASSERT(prec);
301 ParameterList paramList = *paramList_;
304 const RCP<const ThyLinOpBase> fwdOp = fwdOpSrc->getOp();
305 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
308 bool bIsEpetra = XpThyUtils::isEpetra(fwdOp);
309 bool bIsTpetra = XpThyUtils::isTpetra(fwdOp);
310 bool bIsBlocked = XpThyUtils::isBlockedOperator(fwdOp);
311 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra ==
true && bIsTpetra ==
true));
312 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == bIsTpetra) && bIsBlocked ==
false);
313 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra != bIsTpetra) && bIsBlocked ==
true);
315 RCP<XpMat> A = Teuchos::null;
317 Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
318 Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<Scalar> >(fwdOp);
319 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(ThyBlockedOp));
321 TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==
false);
323 Teuchos::RCP<const LinearOpBase<Scalar> > b00 = ThyBlockedOp->getBlock(0,0);
324 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
326 RCP<const XpCrsMat > xpetraFwdCrsMat00 = XpThyUtils::toXpetra(b00);
327 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat00));
330 RCP<XpCrsMat> xpetraFwdCrsMatNonConst00 = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat00);
331 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst00));
334 RCP<XpMat> A00 = rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst00));
335 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A00));
337 RCP<const XpMap> rowmap00 = A00->getRowMap();
338 RCP< const Teuchos::Comm< int > > comm = rowmap00->getComm();
341 RCP<XpBlockedCrsMat> bMat = Teuchos::rcp(
new XpBlockedCrsMat(ThyBlockedOp, comm));
342 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(bMat));
347 RCP<const XpCrsMat > xpetraFwdCrsMat = XpThyUtils::toXpetra(fwdOp);
348 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat));
351 RCP<XpCrsMat> xpetraFwdCrsMatNonConst = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat);
352 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst));
355 A = rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst));
357 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A));
360 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(
dynamic_cast<DefaultPreconditioner<Scalar> *
>(prec));
361 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
364 RCP<ThyLinOpBase> thyra_precOp = Teuchos::null;
365 thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar> >(defaultPrec->getNonconstUnspecifiedPrecOp(),
true);
370 const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter(
"reuse: type") || paramList.get<std::string>(
"reuse: type") ==
"none");
371 bool useHalfPrecision =
false;
372 if (paramList.isParameter(
"half precision"))
373 useHalfPrecision = paramList.get<
bool>(
"half precision");
374 else if (paramList.isSublist(
"Hierarchy") && paramList.sublist(
"Hierarchy").isParameter(
"half precision"))
375 useHalfPrecision = paramList.sublist(
"Hierarchy").get<
bool>(
"half precision");
376 if (useHalfPrecision)
377 TEUCHOS_TEST_FOR_EXCEPTION(!bIsTpetra,
MueLu::Exceptions::RuntimeError,
"The only scalar type Epetra knows is double, so a half precision preconditioner cannot be constructed.");
380 if (startingOver ==
true) {
382 std::list<std::string> convertXpetra = {
"Coordinates",
"Nullspace"};
383 for (
auto it = convertXpetra.begin(); it != convertXpetra.end(); ++it)
384 replaceWithXpetra<Scalar,LocalOrdinal,GlobalOrdinal,Node>(paramList,*it);
386 if (useHalfPrecision) {
387 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) 394 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
395 const std::string userName =
"user data";
396 Teuchos::ParameterList& userParamList = paramList.sublist(userName);
397 if (userParamList.isType<RCP<XpmMV> >(
"Coordinates")) {
398 RCP<XpmMV> coords = userParamList.get<RCP<XpmMV> >(
"Coordinates");
399 userParamList.remove(
"Coordinates");
400 RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
401 userParamList.set(
"Coordinates",halfCoords);
403 if (userParamList.isType<RCP<XpMV> >(
"Nullspace")) {
404 RCP<XpMV> nullspace = userParamList.get<RCP<XpMV> >(
"Nullspace");
405 userParamList.remove(
"Nullspace");
406 RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
407 userParamList.set(
"Nullspace",halfNullspace);
409 if (paramList.isType<RCP<XpmMV> >(
"Coordinates")) {
410 RCP<XpmMV> coords = paramList.get<RCP<XpmMV> >(
"Coordinates");
411 paramList.remove(
"Coordinates");
412 RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
413 userParamList.set(
"Coordinates",halfCoords);
415 if (paramList.isType<RCP<XpMV> >(
"Nullspace")) {
416 RCP<XpMV> nullspace = paramList.get<RCP<XpMV> >(
"Nullspace");
417 paramList.remove(
"Nullspace");
418 RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
419 userParamList.set(
"Nullspace",halfNullspace);
426 RCP<MueLuHalfXpOp> xpOp = rcp(
new MueLuHalfXpOp(H));
427 xpPrecOp = rcp(
new XpHalfPrecOp(xpOp));
429 TEUCHOS_TEST_FOR_EXCEPT(
true);
433 const std::string userName =
"user data";
434 Teuchos::ParameterList& userParamList = paramList.sublist(userName);
435 if (paramList.isType<RCP<XpmMV> >(
"Coordinates")) {
436 RCP<XpmMV> coords = paramList.get<RCP<XpmMV> >(
"Coordinates");
437 paramList.remove(
"Coordinates");
438 userParamList.set(
"Coordinates",coords);
440 if (paramList.isType<RCP<XpMV> >(
"Nullspace")) {
441 RCP<XpMV> nullspace = paramList.get<RCP<XpMV> >(
"Nullspace");
442 paramList.remove(
"Nullspace");
443 userParamList.set(
"Nullspace",nullspace);
448 xpPrecOp = rcp(
new MueLuXpOp(H));
452 RCP<ThyXpOp> thyXpOp = rcp_dynamic_cast<ThyXpOp>(thyra_precOp,
true);
453 xpPrecOp = rcp_dynamic_cast<XpOp>(thyXpOp->getXpetraOperator(),
true);
454 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) 455 RCP<XpHalfPrecOp> xpHalfPrecOp = rcp_dynamic_cast<XpHalfPrecOp>(xpPrecOp);
456 if (!xpHalfPrecOp.is_null()) {
457 RCP<MueLu::Hierarchy<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> > H = rcp_dynamic_cast<MueLuHalfXpOp>(xpHalfPrecOp->GetHalfPrecisionOperator(),
true)->GetHierarchy();
458 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
461 "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it");
463 "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator");
464 RCP<MueLu::Level> level0 = H->GetLevel(0);
465 RCP<XpHalfOp> O0 = level0->Get<RCP<XpHalfOp> >(
"A");
466 RCP<XphMat> A0 = rcp_dynamic_cast<XphMat>(O0,
true);
472 halfA->SetFixedBlockSize(A0->GetFixedBlockSize());
476 level0->Set(
"A", halfA);
483 RCP<MueLuXpOp> xpOp = rcp_dynamic_cast<MueLuXpOp>(thyXpOp->getXpetraOperator(),
true);
484 RCP<MueLu::Hierarchy<Scalar,LocalOrdinal,GlobalOrdinal,Node> > H = xpOp->GetHierarchy();;
487 "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it");
489 "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator");
490 RCP<MueLu::Level> level0 = H->GetLevel(0);
491 RCP<XpOp> O0 = level0->Get<RCP<XpOp> >(
"A");
492 RCP<XpMat> A0 = rcp_dynamic_cast<XpMat>(O0);
498 A->SetFixedBlockSize(A0->GetFixedBlockSize());
509 RCP<const VectorSpaceBase<Scalar> > thyraRangeSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpPrecOp->getRangeMap());
510 RCP<const VectorSpaceBase<Scalar> > thyraDomainSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpPrecOp->getDomainMap());
512 RCP<ThyLinOpBase > thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace, xpPrecOp);
513 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraPrecOp));
515 defaultPrec->initializeUnspecified(thyraPrecOp);
519 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
520 void MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::
521 uninitializePrec(PreconditionerBase<Scalar>* prec, RCP<
const LinearOpSourceBase<Scalar> >* fwdOp, ESupportSolveUse* supportSolveUse)
const {
522 TEUCHOS_ASSERT(prec);
525 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(
dynamic_cast<DefaultPreconditioner<Scalar> *
>(prec));
526 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
530 *fwdOp = Teuchos::null;
533 if (supportSolveUse) {
535 *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
538 defaultPrec->uninitialize();
543 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
544 void MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::setParameterList(RCP<ParameterList>
const& paramList) {
545 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
546 paramList_ = paramList;
549 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
550 RCP<ParameterList> MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getNonconstParameterList() {
554 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
555 RCP<ParameterList> MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::unsetParameterList() {
556 RCP<ParameterList> savedParamList = paramList_;
557 paramList_ = Teuchos::null;
558 return savedParamList;
561 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
562 RCP<const ParameterList> MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getParameterList()
const {
566 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
567 RCP<const ParameterList> MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getValidParameters()
const {
568 static RCP<const ParameterList> validPL;
570 if (Teuchos::is_null(validPL))
571 validPL = rcp(
new ParameterList());
577 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
578 std::string MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::description()
const {
579 return "Thyra::MueLuPreconditionerFactory";
583 #endif // HAVE_MUELU_STRATIMIKOS 585 #endif // ifdef THYRA_MUELU_PRECONDITIONER_FACTORY_DEF_HPP MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
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...
Exception throws to report errors in the internal logical of the program.
Wraps an existing MueLu::Hierarchy as a Xpetra::Operator.