46#ifndef THYRA_MUELU_REFMAXWELL_PRECONDITIONER_FACTORY_DEF_HPP
47#define THYRA_MUELU_REFMAXWELL_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
67 using Teuchos::ParameterList;
68 using Teuchos::rcp_dynamic_cast;
69 using Teuchos::rcp_const_cast;
73 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
74 MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::MueLuRefMaxwellPreconditionerFactory() :
75 paramList_(rcp(new ParameterList()))
80 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
81 bool MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isCompatible(
const LinearOpSourceBase<Scalar>& fwdOpSrc)
const {
82 const RCP<const LinearOpBase<Scalar> > fwdOp = fwdOpSrc.getOp();
84#ifdef HAVE_MUELU_TPETRA
85 if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isTpetra(fwdOp))
return true;
88#ifdef HAVE_MUELU_EPETRA
89 if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isEpetra(fwdOp))
return true;
96 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
97 RCP<PreconditionerBase<Scalar> > MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::createPrec()
const {
98 return Teuchos::rcp(
new DefaultPreconditioner<Scalar>);
101 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
102 void MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::
103 initializePrec(
const RCP<
const LinearOpSourceBase<Scalar> >& fwdOpSrc, PreconditionerBase<Scalar>* prec,
const ESupportSolveUse supportSolveUse)
const {
106 typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
107 typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpThyUtils;
108 typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
109 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
111#if defined(MUELU_CAN_USE_MIXED_PRECISION)
112 typedef Xpetra::TpetraHalfPrecisionOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpHalfPrecOp;
113 typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMV;
114 typedef typename XpHalfPrecOp::HalfScalar HalfScalar;
115 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
116 typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
117 typedef Xpetra::MultiVector<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XphMV;
118 typedef Xpetra::MultiVector<Magnitude,LocalOrdinal,GlobalOrdinal,Node > XpmMV;
119 typedef Xpetra::MultiVector<HalfMagnitude,LocalOrdinal,GlobalOrdinal,Node > XphmMV;
120 typedef Xpetra::Matrix<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XphMat;
122 Teuchos::TimeMonitor tM(*Teuchos::TimeMonitor::getNewTimer(std::string(
"ThyraMueLuRefMaxwell::initializePrec")));
125 TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
126 TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
127 TEUCHOS_ASSERT(prec);
130 ParameterList paramList = *paramList_;
133 const RCP<const ThyLinOpBase> fwdOp = fwdOpSrc->getOp();
134 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
137 bool bIsEpetra = XpThyUtils::isEpetra(fwdOp);
138 bool bIsTpetra = XpThyUtils::isTpetra(fwdOp);
139 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra ==
true && bIsTpetra ==
true));
143 RCP<XpMat> A = XpThyUtils::toXpetra(Teuchos::rcp_const_cast<ThyLinOpBase>(fwdOp));
144 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A));
147 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(
dynamic_cast<DefaultPreconditioner<Scalar> *
>(prec));
148 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
151 RCP<ThyLinOpBase> thyra_precOp = Teuchos::null;
152 thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar> >(defaultPrec->getNonconstUnspecifiedPrecOp(),
true);
157 const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter(
"refmaxwell: enable reuse") || !paramList.get<
bool>(
"refmaxwell: enable reuse"));
158 const bool useHalfPrecision = paramList.get<
bool>(
"half precision",
false) && bIsTpetra;
161 if (startingOver ==
true) {
164 std::list<std::string> convertXpetra = {
"Coordinates",
"Nullspace",
"M1",
"Ms",
"D0",
"M0inv"};
165 for (
auto it = convertXpetra.begin(); it != convertXpetra.end(); ++it)
166 Converters<Scalar,LocalOrdinal,GlobalOrdinal,Node>::replaceWithXpetra(paramList,*it);
168 paramList.set<
bool>(
"refmaxwell: use as preconditioner",
true);
169 if (useHalfPrecision) {
170#if defined(MUELU_CAN_USE_MIXED_PRECISION)
173 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
174 if (paramList.isType<RCP<XpmMV> >(
"Coordinates")) {
175 RCP<XpmMV> coords = paramList.get<RCP<XpmMV> >(
"Coordinates");
176 paramList.remove(
"Coordinates");
177 RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
178 paramList.set(
"Coordinates",halfCoords);
180 if (paramList.isType<RCP<XpMV> >(
"Nullspace")) {
181 RCP<XpMV> nullspace = paramList.get<RCP<XpMV> >(
"Nullspace");
182 paramList.remove(
"Nullspace");
183 RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
184 paramList.set(
"Nullspace",halfNullspace);
186 std::list<std::string> convertMat = {
"M1",
"Ms",
"D0",
"M0inv"};
187 for (
auto it = convertMat.begin(); it != convertMat.end(); ++it) {
188 if (paramList.isType<RCP<XpMat> >(*it)) {
189 RCP<XpMat> M = paramList.get<RCP<XpMat> >(*it);
190 paramList.remove(*it);
191 RCP<XphMat> halfM = Xpetra::convertToHalfPrecision(M);
192 paramList.set(*it,halfM);
198 xpPrecOp = rcp(
new XpHalfPrecOp(halfPrec));
200 TEUCHOS_TEST_FOR_EXCEPT(
true);
206 xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
211 RCP<ThyXpOp> thyXpOp = rcp_dynamic_cast<ThyXpOp>(thyra_precOp,
true);
212 RCP<XpOp> xpOp = thyXpOp->getXpetraOperator();
213#if defined(MUELU_CAN_USE_MIXED_PRECISION)
214 RCP<XpHalfPrecOp> xpHalfPrecOp = rcp_dynamic_cast<XpHalfPrecOp>(xpOp);
215 if (!xpHalfPrecOp.is_null()) {
216 RCP<MueLu::RefMaxwell<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> > preconditioner = rcp_dynamic_cast<MueLu::RefMaxwell<HalfScalar,LocalOrdinal,GlobalOrdinal,Node>>(xpHalfPrecOp->GetHalfPrecisionOperator(),
true);
217 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
218 preconditioner->resetMatrix(halfA);
219 xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
223 RCP<MueLu::RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node> > preconditioner = rcp_dynamic_cast<MueLu::RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node>>(xpOp,
true);
224 preconditioner->resetMatrix(A);
225 xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
230 RCP<const VectorSpaceBase<Scalar> > thyraRangeSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpPrecOp->getRangeMap());
231 RCP<const VectorSpaceBase<Scalar> > thyraDomainSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpPrecOp->getDomainMap());
233 RCP<ThyLinOpBase > thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace, xpPrecOp);
234 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraPrecOp));
236 defaultPrec->initializeUnspecified(thyraPrecOp);
240 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
241 void MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::
242 uninitializePrec(PreconditionerBase<Scalar>* prec, RCP<
const LinearOpSourceBase<Scalar> >* fwdOp, ESupportSolveUse* supportSolveUse)
const {
243 TEUCHOS_ASSERT(prec);
246 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(
dynamic_cast<DefaultPreconditioner<Scalar> *
>(prec));
247 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
251 *fwdOp = Teuchos::null;
254 if (supportSolveUse) {
256 *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
259 defaultPrec->uninitialize();
264 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
265 void MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::setParameterList(RCP<ParameterList>
const& paramList) {
266 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
267 paramList_ = paramList;
270 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
271 RCP<ParameterList> MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getNonconstParameterList() {
275 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
276 RCP<ParameterList> MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::unsetParameterList() {
277 RCP<ParameterList> savedParamList = paramList_;
278 paramList_ = Teuchos::null;
279 return savedParamList;
282 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
283 RCP<const ParameterList> MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getParameterList()
const {
287 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
288 RCP<const ParameterList> MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getValidParameters()
const {
289 static RCP<const ParameterList> validPL;
291 if (Teuchos::is_null(validPL))
292 validPL = rcp(
new ParameterList());
298 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
299 std::string MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::description()
const {
300 return "Thyra::MueLuRefMaxwellPreconditionerFactory";
Preconditioner (wrapped as a Xpetra::Operator) for Maxwell's equations in curl-curl form.
Concrete Thyra::LinearOpBase subclass for Xpetra::Operator.