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) 58 using Teuchos::ParameterList;
59 using Teuchos::rcp_dynamic_cast;
60 using Teuchos::rcp_const_cast;
64 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
65 MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::MueLuRefMaxwellPreconditionerFactory() :
66 paramList_(rcp(new ParameterList()))
71 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
72 bool MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isCompatible(
const LinearOpSourceBase<Scalar>& fwdOpSrc)
const {
73 const RCP<const LinearOpBase<Scalar> > fwdOp = fwdOpSrc.getOp();
75 #ifdef HAVE_MUELU_TPETRA 76 if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isTpetra(fwdOp))
return true;
83 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
84 RCP<PreconditionerBase<Scalar> > MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::createPrec()
const {
85 return Teuchos::rcp(
new DefaultPreconditioner<Scalar>);
88 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
89 void MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::
90 initializePrec(
const RCP<
const LinearOpSourceBase<Scalar> >& fwdOpSrc, PreconditionerBase<Scalar>* prec,
const ESupportSolveUse supportSolveUse)
const {
93 typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
94 typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpThyUtils;
95 typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
96 typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
97 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
99 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) 100 typedef Xpetra::TpetraHalfPrecisionOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpHalfPrecOp;
101 typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMV;
102 typedef typename XpHalfPrecOp::HalfScalar HalfScalar;
103 typedef Xpetra::Operator<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XpHalfOp;
104 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
105 typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
106 typedef Xpetra::MultiVector<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XphMV;
107 typedef Xpetra::MultiVector<Magnitude,LocalOrdinal,GlobalOrdinal,Node > XpmMV;
108 typedef Xpetra::MultiVector<HalfMagnitude,LocalOrdinal,GlobalOrdinal,Node > XphmMV;
109 typedef Xpetra::Matrix<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XphMat;
111 Teuchos::TimeMonitor tM(*Teuchos::TimeMonitor::getNewTimer(std::string(
"ThyraMueLuRefMaxwell::initializePrec")));
114 TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
115 TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
116 TEUCHOS_ASSERT(prec);
119 ParameterList paramList = *paramList_;
122 const RCP<const ThyLinOpBase> fwdOp = fwdOpSrc->getOp();
123 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
126 bool bIsEpetra = XpThyUtils::isEpetra(fwdOp);
127 bool bIsTpetra = XpThyUtils::isTpetra(fwdOp);
128 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra ==
true && bIsTpetra ==
true));
130 RCP<const XpCrsMat > xpetraFwdCrsMat = XpThyUtils::toXpetra(fwdOp);
131 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat));
134 RCP<XpCrsMat> xpetraFwdCrsMatNonConst = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat);
135 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst));
138 RCP<XpMat> A = rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst));
139 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A));
142 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(
dynamic_cast<DefaultPreconditioner<Scalar> *
>(prec));
143 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
146 RCP<ThyLinOpBase> thyra_precOp = Teuchos::null;
147 thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar> >(defaultPrec->getNonconstUnspecifiedPrecOp(),
true);
152 const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter(
"refmaxwell: enable reuse") || !paramList.get<
bool>(
"refmaxwell: enable reuse"));
153 const bool useHalfPrecision = paramList.get<
bool>(
"half precision",
false) && bIsTpetra;
156 if (startingOver ==
true) {
159 std::list<std::string> convertXpetra = {
"Coordinates",
"Nullspace",
"M1",
"Ms",
"D0",
"M0inv"};
160 for (
auto it = convertXpetra.begin(); it != convertXpetra.end(); ++it)
161 replaceWithXpetra<Scalar,LocalOrdinal,GlobalOrdinal,Node>(paramList,*it);
163 paramList.set<
bool>(
"refmaxwell: use as preconditioner",
true);
164 if (useHalfPrecision) {
165 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) 168 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
169 if (paramList.isType<RCP<XpmMV> >(
"Coordinates")) {
170 RCP<XpmMV> coords = paramList.get<RCP<XpmMV> >(
"Coordinates");
171 paramList.remove(
"Coordinates");
172 RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
173 paramList.set(
"Coordinates",halfCoords);
175 if (paramList.isType<RCP<XpMV> >(
"Nullspace")) {
176 RCP<XpMV> nullspace = paramList.get<RCP<XpMV> >(
"Nullspace");
177 paramList.remove(
"Nullspace");
178 RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
179 paramList.set(
"Nullspace",halfNullspace);
181 std::list<std::string> convertMat = {
"M1",
"Ms",
"D0",
"M0inv"};
182 for (
auto it = convertMat.begin(); it != convertMat.end(); ++it) {
183 if (paramList.isType<RCP<XpMat> >(*it)) {
184 RCP<XpMat> M = paramList.get<RCP<XpMat> >(*it);
185 paramList.remove(*it);
186 RCP<XphMat> halfM = Xpetra::convertToHalfPrecision(M);
187 paramList.set(*it,halfM);
193 xpPrecOp = rcp(
new XpHalfPrecOp(halfPrec));
195 TEUCHOS_TEST_FOR_EXCEPT(
true);
201 xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
206 RCP<ThyXpOp> thyXpOp = rcp_dynamic_cast<ThyXpOp>(thyra_precOp,
true);
207 RCP<XpOp> xpOp = thyXpOp->getXpetraOperator();
208 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) 209 RCP<XpHalfPrecOp> xpHalfPrecOp = rcp_dynamic_cast<XpHalfPrecOp>(xpOp);
210 if (!xpHalfPrecOp.is_null()) {
212 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
213 preconditioner->resetMatrix(halfA);
214 xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
220 xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
225 RCP<const VectorSpaceBase<Scalar> > thyraRangeSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpPrecOp->getRangeMap());
226 RCP<const VectorSpaceBase<Scalar> > thyraDomainSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpPrecOp->getDomainMap());
228 RCP<ThyLinOpBase > thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace, xpPrecOp);
229 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraPrecOp));
231 defaultPrec->initializeUnspecified(thyraPrecOp);
235 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
236 void MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::
237 uninitializePrec(PreconditionerBase<Scalar>* prec, RCP<
const LinearOpSourceBase<Scalar> >* fwdOp, ESupportSolveUse* supportSolveUse)
const {
238 TEUCHOS_ASSERT(prec);
241 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(
dynamic_cast<DefaultPreconditioner<Scalar> *
>(prec));
242 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
246 *fwdOp = Teuchos::null;
249 if (supportSolveUse) {
251 *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
254 defaultPrec->uninitialize();
259 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
260 void MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::setParameterList(RCP<ParameterList>
const& paramList) {
261 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
262 paramList_ = paramList;
265 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
266 RCP<ParameterList> MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getNonconstParameterList() {
270 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
271 RCP<ParameterList> MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::unsetParameterList() {
272 RCP<ParameterList> savedParamList = paramList_;
273 paramList_ = Teuchos::null;
274 return savedParamList;
277 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
278 RCP<const ParameterList> MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getParameterList()
const {
282 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
283 RCP<const ParameterList> MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getValidParameters()
const {
284 static RCP<const ParameterList> validPL;
286 if (Teuchos::is_null(validPL))
287 validPL = rcp(
new ParameterList());
293 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
294 std::string MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::description()
const {
295 return "Thyra::MueLuRefMaxwellPreconditionerFactory";
299 #endif // HAVE_MUELU_STRATIMIKOS 301 #endif // ifdef THYRA_MUELU_REFMAXWELL_PRECONDITIONER_FACTORY_DEF_HPP Preconditioner (wrapped as a Xpetra::Operator) for Maxwell's equations in curl-curl form...
void resetMatrix(Teuchos::RCP< Matrix > SM_Matrix_new, bool ComputePrec=true)
Reset system matrix.
Concrete Thyra::LinearOpBase subclass for Xpetra::Operator.