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_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) && !defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) && !defined(HAVE_TPETRA_INST_COMPLEX_FLOAT)) || \
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 #define MUELU_CAN_USE_MIXED_PRECISION
66 using Teuchos::rcp_const_cast;
67 using Teuchos::rcp_dynamic_cast;
71 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
72 MueLuRefMaxwellPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MueLuRefMaxwellPreconditionerFactory()
73 : paramList_(
rcp(new ParameterList())) {}
77 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
78 bool MueLuRefMaxwellPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isCompatible(
const LinearOpSourceBase<Scalar>& fwdOpSrc)
const {
79 const RCP<const LinearOpBase<Scalar>> fwdOp = fwdOpSrc.getOp();
81 if (Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isTpetra(fwdOp))
return true;
83 #ifdef HAVE_MUELU_EPETRA
84 if (Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isEpetra(fwdOp))
return true;
90 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
91 RCP<PreconditionerBase<Scalar>> MueLuRefMaxwellPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::createPrec()
const {
95 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
96 void MueLuRefMaxwellPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
97 initializePrec(
const RCP<
const LinearOpSourceBase<Scalar>>& fwdOpSrc, PreconditionerBase<Scalar>* prec,
const ESupportSolveUse supportSolveUse)
const {
100 typedef Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpThyUtils;
102 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
104 #if defined(MUELU_CAN_USE_MIXED_PRECISION)
105 typedef Xpetra::TpetraHalfPrecisionOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpHalfPrecOp;
107 typedef typename XpHalfPrecOp::HalfScalar HalfScalar;
123 ParameterList paramList = *paramList_;
126 const RCP<const ThyLinOpBase> fwdOp = fwdOpSrc->getOp();
130 bool bIsEpetra = XpThyUtils::isEpetra(fwdOp);
131 bool bIsTpetra = XpThyUtils::isTpetra(fwdOp);
136 RCP<XpMat> A = XpThyUtils::toXpetra(Teuchos::rcp_const_cast<ThyLinOpBase>(fwdOp));
144 RCP<ThyLinOpBase> thyra_precOp = Teuchos::null;
145 thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar>>(defaultPrec->getNonconstUnspecifiedPrecOp(),
true);
150 const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter(
"refmaxwell: enable reuse") || !paramList.get<
bool>(
"refmaxwell: enable reuse"));
151 const bool useHalfPrecision = paramList.get<
bool>(
"half precision",
false) && bIsTpetra;
154 if (startingOver ==
true) {
156 std::list<std::string> convertMat = {
157 "Dk_1",
"Dk_2",
"D0",
158 "Mk_one",
"Mk_1_one",
"M1_beta",
"M1_alpha",
159 "invMk_1_invBeta",
"invMk_2_invAlpha",
161 "M1",
"Ms",
"M0inv"};
162 std::list<std::string> convertMV = {
"Coordinates",
"Nullspace"};
163 std::list<std::string> convertXpetra;
164 convertXpetra.insert(convertXpetra.end(), convertMV.begin(), convertMV.end());
165 convertXpetra.insert(convertXpetra.end(), convertMat.begin(), convertMat.end());
166 for (
auto it = convertXpetra.begin(); it != convertXpetra.end(); ++it)
167 Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(paramList, *it);
169 paramList.set<
bool>(
"refmaxwell: use as preconditioner",
true);
170 if (useHalfPrecision) {
171 #if defined(MUELU_CAN_USE_MIXED_PRECISION)
174 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
175 if (paramList.isType<RCP<XpmMV>>(
"Coordinates")) {
176 RCP<XpmMV> coords = paramList.get<RCP<XpmMV>>(
"Coordinates");
177 paramList.remove(
"Coordinates");
178 RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
179 paramList.set(
"Coordinates", halfCoords);
181 if (paramList.isType<RCP<XpMV>>(
"Nullspace")) {
182 RCP<XpMV> nullspace = paramList.get<RCP<XpMV>>(
"Nullspace");
183 paramList.remove(
"Nullspace");
184 RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
185 paramList.set(
"Nullspace", halfNullspace);
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));
205 xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
210 RCP<ThyXpOp> thyXpOp = rcp_dynamic_cast<ThyXpOp>(thyra_precOp,
true);
211 RCP<XpOp> xpOp = thyXpOp->getXpetraOperator();
212 #if defined(MUELU_CAN_USE_MIXED_PRECISION)
213 RCP<XpHalfPrecOp> xpHalfPrecOp = rcp_dynamic_cast<XpHalfPrecOp>(xpOp);
214 if (!xpHalfPrecOp.is_null()) {
216 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
217 preconditioner->resetMatrix(halfA);
218 xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
224 xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
229 RCP<const VectorSpaceBase<Scalar>> thyraRangeSpace = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpPrecOp->getRangeMap());
230 RCP<const VectorSpaceBase<Scalar>> thyraDomainSpace = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpPrecOp->getDomainMap());
232 RCP<ThyLinOpBase> thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace, xpPrecOp);
235 defaultPrec->initializeUnspecified(thyraPrecOp);
238 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
239 void MueLuRefMaxwellPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
240 uninitializePrec(PreconditionerBase<Scalar>* prec, RCP<
const LinearOpSourceBase<Scalar>>* fwdOp, ESupportSolveUse* supportSolveUse)
const {
249 *fwdOp = Teuchos::null;
252 if (supportSolveUse) {
254 *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
257 defaultPrec->uninitialize();
261 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
262 void MueLuRefMaxwellPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setParameterList(RCP<ParameterList>
const& paramList) {
264 paramList_ = paramList;
267 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
268 RCP<ParameterList> MueLuRefMaxwellPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getNonconstParameterList() {
272 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
273 RCP<ParameterList> MueLuRefMaxwellPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::unsetParameterList() {
274 RCP<ParameterList> savedParamList = paramList_;
275 paramList_ = Teuchos::null;
276 return savedParamList;
279 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
280 RCP<const ParameterList> MueLuRefMaxwellPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getParameterList()
const {
284 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
285 RCP<const ParameterList> MueLuRefMaxwellPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getValidParameters()
const {
286 static RCP<const ParameterList> validPL;
289 validPL =
rcp(
new ParameterList());
295 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
296 std::string MueLuRefMaxwellPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::description()
const {
297 return "Thyra::MueLuRefMaxwellPreconditionerFactory";
301 #endif // HAVE_MUELU_STRATIMIKOS
303 #endif // ifdef THYRA_MUELU_REFMAXWELL_PRECONDITIONER_FACTORY_DEF_HPP
bool nonnull(const std::shared_ptr< T > &p)
bool is_null(const std::shared_ptr< T > &p)
Preconditioner (wrapped as a Xpetra::Operator) for Maxwell's equations in curl-curl form...
static RCP< Time > getNewTimer(const std::string &name)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void resetMatrix(Teuchos::RCP< Matrix > SM_Matrix_new, bool ComputePrec=true)
Reset system matrix.
Concrete Thyra::LinearOpBase subclass for Xpetra::Operator.
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)