47 #ifndef THYRA_MUELU_REFMAXWELL_PRECONDITIONER_FACTORY_DECL_HPP 48 #define THYRA_MUELU_REFMAXWELL_PRECONDITIONER_FACTORY_DECL_HPP 52 #if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA) 55 #include "Thyra_DefaultPreconditioner.hpp" 56 #include "Thyra_BlockedLinearOpBase.hpp" 57 #include "Thyra_DiagonalLinearOpBase.hpp" 60 #ifdef HAVE_MUELU_TPETRA 61 #include "Thyra_TpetraLinearOp.hpp" 62 #include "Thyra_TpetraThyraWrappers.hpp" 64 #ifdef HAVE_MUELU_EPETRA 65 #include "Thyra_EpetraLinearOp.hpp" 66 #include "Thyra_EpetraThyraWrappers.hpp" 69 #include "Teuchos_Ptr.hpp" 70 #include "Teuchos_TestForException.hpp" 71 #include "Teuchos_Assert.hpp" 72 #include "Teuchos_Time.hpp" 74 #include <Xpetra_CrsMatrixWrap.hpp> 75 #include <Xpetra_CrsMatrix.hpp> 76 #include <Xpetra_Matrix.hpp> 77 #include <Xpetra_ThyraUtils.hpp> 80 #include <MueLu_RefMaxwell.hpp> 81 #ifdef HAVE_MUELU_TPETRA 82 #include <MueLu_TpetraOperator.hpp> 83 #include <Xpetra_TpetraOperator.hpp> 84 #include <Xpetra_TpetraHalfPrecisionOperator.hpp> 86 #ifdef HAVE_MUELU_EPETRA 88 #include <Xpetra_EpetraOperator.hpp> 91 #include "Thyra_PreconditionerFactoryBase.hpp" 93 #include "Kokkos_DefaultNode.hpp" 107 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node = KokkosClassic::DefaultNode::DefaultNodeType>
108 class MueLuRefMaxwellPreconditionerFactory :
public PreconditionerFactoryBase<Scalar> {
115 MueLuRefMaxwellPreconditionerFactory();
122 bool isCompatible(
const LinearOpSourceBase<Scalar>& fwdOp)
const;
124 Teuchos::RCP<PreconditionerBase<Scalar> > createPrec()
const;
126 void initializePrec(
const Teuchos::RCP<
const LinearOpSourceBase<Scalar> >& fwdOp,
127 PreconditionerBase<Scalar>* prec,
128 const ESupportSolveUse supportSolveUse
131 void uninitializePrec(PreconditionerBase<Scalar>* prec,
132 Teuchos::RCP<
const LinearOpSourceBase<Scalar> >* fwdOp,
133 ESupportSolveUse* supportSolveUse
142 void setParameterList(
const Teuchos::RCP<Teuchos::ParameterList>& paramList);
144 Teuchos::RCP<Teuchos::ParameterList> unsetParameterList();
146 Teuchos::RCP<Teuchos::ParameterList> getNonconstParameterList();
148 Teuchos::RCP<const Teuchos::ParameterList> getParameterList()
const;
157 std::string description()
const;
165 Teuchos::RCP<Teuchos::ParameterList> paramList_;
169 #ifdef HAVE_MUELU_EPETRA 178 class MueLuRefMaxwellPreconditionerFactory<double,int,int,Xpetra::
EpetraNode> :
public PreconditionerFactoryBase<double> {
189 MueLuRefMaxwellPreconditionerFactory() : paramList_(rcp(new ParameterList())) { }
197 bool isCompatible(
const LinearOpSourceBase<Scalar>& fwdOpSrc)
const {
198 const RCP<const LinearOpBase<Scalar> > fwdOp = fwdOpSrc.getOp();
200 #ifdef HAVE_MUELU_TPETRA 201 if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isTpetra(fwdOp))
return true;
204 #ifdef HAVE_MUELU_EPETRA 205 if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isEpetra(fwdOp))
return true;
212 Teuchos::RCP<PreconditionerBase<Scalar> > createPrec()
const {
213 return Teuchos::rcp(
new DefaultPreconditioner<Scalar>);
217 void initializePrec(
const Teuchos::RCP<
const LinearOpSourceBase<Scalar> >& fwdOpSrc,
218 PreconditionerBase<Scalar>* prec,
219 const ESupportSolveUse
221 using Teuchos::rcp_dynamic_cast;
224 typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
225 typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpThyUtils;
226 typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
227 typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
228 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
230 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) 231 typedef Xpetra::TpetraHalfPrecisionOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpHalfPrecOp;
232 typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMV;
233 typedef typename XpHalfPrecOp::HalfScalar HalfScalar;
234 typedef Xpetra::Operator<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XpHalfOp;
235 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
236 typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
237 typedef Xpetra::MultiVector<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XphMV;
238 typedef Xpetra::MultiVector<Magnitude,LocalOrdinal,GlobalOrdinal,Node > XpmMV;
239 typedef Xpetra::MultiVector<HalfMagnitude,LocalOrdinal,GlobalOrdinal,Node > XphmMV;
240 typedef Xpetra::Matrix<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XphMat;
242 Teuchos::TimeMonitor tM(*Teuchos::TimeMonitor::getNewTimer(std::string(
"ThyraMueLuRefMaxwell::initializePrec")));
245 TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
246 TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
247 TEUCHOS_ASSERT(prec);
250 ParameterList paramList = *paramList_;
253 const RCP<const ThyLinOpBase> fwdOp = fwdOpSrc->getOp();
254 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
257 bool bIsEpetra = XpThyUtils::isEpetra(fwdOp);
258 bool bIsTpetra = XpThyUtils::isTpetra(fwdOp);
259 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra ==
true && bIsTpetra ==
true));
261 RCP<const XpCrsMat > xpetraFwdCrsMat = XpThyUtils::toXpetra(fwdOp);
262 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat));
265 RCP<XpCrsMat> xpetraFwdCrsMatNonConst = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat);
266 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst));
269 RCP<XpMat> A = rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst));
270 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A));
273 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(
dynamic_cast<DefaultPreconditioner<Scalar> *
>(prec));
274 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
277 RCP<ThyLinOpBase> thyra_precOp = Teuchos::null;
278 thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar> >(defaultPrec->getNonconstUnspecifiedPrecOp(),
true);
283 const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter(
"refmaxwell: enable reuse") || !paramList.get<
bool>(
"refmaxwell: enable reuse"));
284 const bool useHalfPrecision = paramList.get<
bool>(
"half precision",
false) && bIsTpetra;
287 if (startingOver ==
true) {
290 std::list<std::string> convertXpetra = {
"Coordinates",
"Nullspace",
"M1",
"Ms",
"D0",
"M0inv"};
291 for (
auto it = convertXpetra.begin(); it != convertXpetra.end(); ++it)
292 replaceWithXpetra<Scalar,LocalOrdinal,GlobalOrdinal,Node>(paramList,*it);
294 paramList.set<
bool>(
"refmaxwell: use as preconditioner",
true);
295 if (useHalfPrecision) {
296 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) 299 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
300 if (paramList.isType<RCP<XpmMV> >(
"Coordinates")) {
301 RCP<XpmMV> coords = paramList.get<RCP<XpmMV> >(
"Coordinates");
302 paramList.remove(
"Coordinates");
303 RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
304 paramList.set(
"Coordinates",halfCoords);
306 if (paramList.isType<RCP<XpMV> >(
"Nullspace")) {
307 RCP<XpMV> nullspace = paramList.get<RCP<XpMV> >(
"Nullspace");
308 paramList.remove(
"Nullspace");
309 RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
310 paramList.set(
"Nullspace",halfNullspace);
312 std::list<std::string> convertMat = {
"M1",
"Ms",
"D0",
"M0inv"};
313 for (
auto it = convertMat.begin(); it != convertMat.end(); ++it) {
314 if (paramList.isType<RCP<XpMat> >(*it)) {
315 RCP<XpMat> M = paramList.get<RCP<XpMat> >(*it);
316 paramList.remove(*it);
317 RCP<XphMat> halfM = Xpetra::convertToHalfPrecision(M);
318 paramList.set(*it,halfM);
324 xpPrecOp = rcp(
new XpHalfPrecOp(halfPrec));
326 TEUCHOS_TEST_FOR_EXCEPT(
true);
332 xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
337 RCP<ThyXpOp> thyXpOp = rcp_dynamic_cast<ThyXpOp>(thyra_precOp,
true);
338 RCP<XpOp> xpOp = thyXpOp->getXpetraOperator();
339 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) 340 RCP<XpHalfPrecOp> xpHalfPrecOp = rcp_dynamic_cast<XpHalfPrecOp>(xpOp);
341 if (!xpHalfPrecOp.is_null()) {
343 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
344 preconditioner->resetMatrix(halfA);
345 xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
351 xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
356 RCP<const VectorSpaceBase<Scalar> > thyraRangeSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpPrecOp->getRangeMap());
357 RCP<const VectorSpaceBase<Scalar> > thyraDomainSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpPrecOp->getDomainMap());
359 RCP<ThyLinOpBase > thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace, xpPrecOp);
360 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraPrecOp));
362 defaultPrec->initializeUnspecified(thyraPrecOp);
366 void uninitializePrec(PreconditionerBase<Scalar>* prec,
367 Teuchos::RCP<
const LinearOpSourceBase<Scalar> >* fwdOp,
368 ESupportSolveUse* supportSolveUse
370 TEUCHOS_ASSERT(prec);
373 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(
dynamic_cast<DefaultPreconditioner<Scalar> *
>(prec));
374 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
378 *fwdOp = Teuchos::null;
381 if (supportSolveUse) {
383 *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
386 defaultPrec->uninitialize();
395 void setParameterList(
const Teuchos::RCP<Teuchos::ParameterList>& paramList) {
396 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
397 paramList_ = paramList;
400 Teuchos::RCP<Teuchos::ParameterList> unsetParameterList() {
401 RCP<ParameterList> savedParamList = paramList_;
402 paramList_ = Teuchos::null;
403 return savedParamList;
406 Teuchos::RCP<Teuchos::ParameterList> getNonconstParameterList() {
return paramList_; }
408 Teuchos::RCP<const Teuchos::ParameterList> getParameterList()
const {
return paramList_; }
411 static RCP<const ParameterList> validPL;
413 if (Teuchos::is_null(validPL))
414 validPL = rcp(
new ParameterList());
424 std::string description()
const {
return "Thyra::MueLuRefMaxwellPreconditionerFactory"; }
431 Teuchos::RCP<Teuchos::ParameterList> paramList_;
434 #endif // HAVE_MUELU_EPETRA 438 #endif // #ifdef HAVE_MUELU_STRATIMIKOS 440 #endif // THYRA_MUELU_REFMAXWELL_PRECONDITIONER_FACTORY_DECL_HPP MueLu::DefaultLocalOrdinal LocalOrdinal
Preconditioner (wrapped as a Xpetra::Operator) for Maxwell's equations in curl-curl form...
MueLu::DefaultScalar Scalar
void resetMatrix(Teuchos::RCP< Matrix > SM_Matrix_new, bool ComputePrec=true)
Reset system matrix.
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Concrete Thyra::LinearOpBase subclass for Xpetra::Operator.
Kokkos::Compat::KokkosSerialWrapperNode EpetraNode
void getValidParameters(Teuchos::ParameterList ¶ms)