46 #ifndef THYRA_MUELU_MAXWELL1_PRECONDITIONER_FACTORY_DEF_HPP
47 #define THYRA_MUELU_MAXWELL1_PRECONDITIONER_FACTORY_DEF_HPP
51 #include <Xpetra_CrsMatrixWrap.hpp>
55 #include <MueLu_Maxwell1.hpp>
56 #include <Xpetra_TpetraHalfPrecisionOperator.hpp>
58 #if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA)
61 #if ((defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) && !defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) && !defined(HAVE_TPETRA_INST_COMPLEX_FLOAT)) || \
62 (!defined(HAVE_TPETRA_INST_DOUBLE) && !defined(HAVE_TPETRA_INST_FLOAT) && defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) && defined(HAVE_TPETRA_INST_COMPLEX_FLOAT)) || \
63 (defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) && defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) && defined(HAVE_TPETRA_INST_COMPLEX_FLOAT)))
64 #define MUELU_CAN_USE_MIXED_PRECISION
72 using Teuchos::rcp_const_cast;
73 using Teuchos::rcp_dynamic_cast;
77 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
78 MueLuMaxwell1PreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MueLuMaxwell1PreconditionerFactory()
79 : paramList_(
rcp(new ParameterList())) {}
83 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
84 bool MueLuMaxwell1PreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isCompatible(
const LinearOpSourceBase<Scalar>& fwdOpSrc)
const {
85 const RCP<const LinearOpBase<Scalar>> fwdOp = fwdOpSrc.getOp();
87 if (Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isTpetra(fwdOp))
return true;
89 #ifdef HAVE_MUELU_EPETRA
90 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>> MueLuMaxwell1PreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::createPrec()
const {
101 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
102 void MueLuMaxwell1PreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
103 initializePrec(
const RCP<
const LinearOpSourceBase<Scalar>>& fwdOpSrc, PreconditionerBase<Scalar>* prec,
const ESupportSolveUse supportSolveUse)
const {
106 typedef Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpThyUtils;
108 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
110 #if defined(MUELU_CAN_USE_MIXED_PRECISION)
111 typedef Xpetra::TpetraHalfPrecisionOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpHalfPrecOp;
113 typedef typename XpHalfPrecOp::HalfScalar HalfScalar;
129 ParameterList paramList = *paramList_;
132 const RCP<const ThyLinOpBase> fwdOp = fwdOpSrc->getOp();
136 bool bIsEpetra = XpThyUtils::isEpetra(fwdOp);
137 bool bIsTpetra = XpThyUtils::isTpetra(fwdOp);
142 RCP<XpMat> A = XpThyUtils::toXpetra(Teuchos::rcp_const_cast<ThyLinOpBase>(fwdOp));
150 RCP<ThyLinOpBase> thyra_precOp = Teuchos::null;
151 thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar>>(defaultPrec->getNonconstUnspecifiedPrecOp(),
true);
156 const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter(
"Maxwell1: enable reuse") || !paramList.get<
bool>(
"Maxwell1: enable reuse"));
157 const bool useHalfPrecision = paramList.get<
bool>(
"half precision",
false) && bIsTpetra;
160 if (startingOver ==
true) {
162 std::list<std::string> convertXpetra = {
"Coordinates",
"Nullspace",
"Kn",
"D0"};
163 for (
auto it = convertXpetra.begin(); it != convertXpetra.end(); ++it)
164 Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(paramList, *it);
166 std::list<std::string> sublists = {
"maxwell1: 11list",
"maxwell1: 22list"};
167 for (
auto itSublist = sublists.begin(); itSublist != sublists.end(); ++itSublist)
168 if (paramList.isSublist(*itSublist)) {
169 ParameterList& sublist = paramList.sublist(*itSublist);
170 for (
int lvlNo = 0; lvlNo < 10; ++lvlNo) {
171 if (sublist.isSublist(
"level " + std::to_string(lvlNo) +
" user data")) {
172 ParameterList& lvlList = sublist.sublist(
"level " + std::to_string(lvlNo) +
" user data");
173 std::list<std::string> convertKeys;
174 for (
auto it = lvlList.begin(); it != lvlList.end(); ++it)
175 convertKeys.push_back(lvlList.name(it));
176 for (
auto it = convertKeys.begin(); it != convertKeys.end(); ++it)
177 Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(lvlList, *it);
182 ParameterList& sublist = paramList.sublist(
"maxwell1: 11list");
183 if (sublist.isParameter(
"D0")) {
184 Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(sublist,
"D0");
187 paramList.set<
bool>(
"Maxwell1: use as preconditioner",
true);
188 if (useHalfPrecision) {
189 #if defined(MUELU_CAN_USE_MIXED_PRECISION)
192 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
193 if (paramList.isType<RCP<XpmMV>>(
"Coordinates")) {
194 RCP<XpmMV> coords = paramList.get<RCP<XpmMV>>(
"Coordinates");
195 paramList.remove(
"Coordinates");
196 RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
197 paramList.set(
"Coordinates", halfCoords);
199 if (paramList.isType<RCP<XpMV>>(
"Nullspace")) {
200 RCP<XpMV> nullspace = paramList.get<RCP<XpMV>>(
"Nullspace");
201 paramList.remove(
"Nullspace");
202 RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
203 paramList.set(
"Nullspace", halfNullspace);
205 std::list<std::string> convertMat = {
"Kn",
"D0"};
206 for (
auto it = convertMat.begin(); it != convertMat.end(); ++it) {
207 if (paramList.isType<RCP<XpMat>>(*it)) {
208 RCP<XpMat> M = paramList.get<RCP<XpMat>>(*it);
209 paramList.remove(*it);
210 RCP<XphMat> halfM = Xpetra::convertToHalfPrecision(M);
211 paramList.set(*it, halfM);
217 xpPrecOp =
rcp(
new XpHalfPrecOp(halfPrec));
224 xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
229 RCP<ThyXpOp> thyXpOp = rcp_dynamic_cast<ThyXpOp>(thyra_precOp,
true);
230 RCP<XpOp> xpOp = thyXpOp->getXpetraOperator();
231 #if defined(MUELU_CAN_USE_MIXED_PRECISION)
232 RCP<XpHalfPrecOp> xpHalfPrecOp = rcp_dynamic_cast<XpHalfPrecOp>(xpOp);
233 if (!xpHalfPrecOp.is_null()) {
235 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
236 preconditioner->resetMatrix(halfA);
237 xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
243 xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
248 RCP<const VectorSpaceBase<Scalar>> thyraRangeSpace = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpPrecOp->getRangeMap());
249 RCP<const VectorSpaceBase<Scalar>> thyraDomainSpace = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpPrecOp->getDomainMap());
251 RCP<ThyLinOpBase> thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace, xpPrecOp);
254 defaultPrec->initializeUnspecified(thyraPrecOp);
257 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
258 void MueLuMaxwell1PreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
259 uninitializePrec(PreconditionerBase<Scalar>* prec, RCP<
const LinearOpSourceBase<Scalar>>* fwdOp, ESupportSolveUse* supportSolveUse)
const {
268 *fwdOp = Teuchos::null;
271 if (supportSolveUse) {
273 *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
276 defaultPrec->uninitialize();
280 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
281 void MueLuMaxwell1PreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setParameterList(RCP<ParameterList>
const& paramList) {
283 paramList_ = paramList;
286 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
287 RCP<ParameterList> MueLuMaxwell1PreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getNonconstParameterList() {
291 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
292 RCP<ParameterList> MueLuMaxwell1PreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::unsetParameterList() {
293 RCP<ParameterList> savedParamList = paramList_;
294 paramList_ = Teuchos::null;
295 return savedParamList;
298 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
299 RCP<const ParameterList> MueLuMaxwell1PreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getParameterList()
const {
303 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
304 RCP<const ParameterList> MueLuMaxwell1PreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getValidParameters()
const {
305 static RCP<const ParameterList> validPL;
308 validPL =
rcp(
new ParameterList());
314 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
315 std::string MueLuMaxwell1PreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::description()
const {
316 return "Thyra::MueLuMaxwell1PreconditionerFactory";
320 #endif // HAVE_MUELU_STRATIMIKOS
322 #endif // ifdef THYRA_MUELU_MAXWELL1_PRECONDITIONER_FACTORY_DEF_HPP
void resetMatrix(Teuchos::RCP< Matrix > SM_Matrix_new, bool ComputePrec=true)
Reset system matrix.
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)
Concrete Thyra::LinearOpBase subclass for Xpetra::Operator.
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)