46#ifndef THYRA_MUELU_MAXWELL1_PRECONDITIONER_FACTORY_DEF_HPP
47#define THYRA_MUELU_MAXWELL1_PRECONDITIONER_FACTORY_DEF_HPP
51#include <Xpetra_CrsMatrixWrap.hpp>
52#include <Xpetra_CrsMatrix.hpp>
53#include <Xpetra_Matrix.hpp>
54#include <Xpetra_ThyraUtils.hpp>
55#include <MueLu_Maxwell1.hpp>
56#ifdef HAVE_MUELU_TPETRA
57# include <Xpetra_TpetraHalfPrecisionOperator.hpp>
61#if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA)
64#if (defined(HAVE_MUELU_TPETRA) && \
65 ((defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) && !defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) && !defined(HAVE_TPETRA_INST_COMPLEX_FLOAT)) || \
66 (!defined(HAVE_TPETRA_INST_DOUBLE) && !defined(HAVE_TPETRA_INST_FLOAT) && defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) && defined(HAVE_TPETRA_INST_COMPLEX_FLOAT)) || \
67 (defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) && defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) && defined(HAVE_TPETRA_INST_COMPLEX_FLOAT))))
68# define MUELU_CAN_USE_MIXED_PRECISION
75 using Teuchos::ParameterList;
76 using Teuchos::rcp_dynamic_cast;
77 using Teuchos::rcp_const_cast;
81 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
82 MueLuMaxwell1PreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::MueLuMaxwell1PreconditionerFactory() :
83 paramList_(rcp(new ParameterList()))
88 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
89 bool MueLuMaxwell1PreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isCompatible(
const LinearOpSourceBase<Scalar>& fwdOpSrc)
const {
90 const RCP<const LinearOpBase<Scalar> > fwdOp = fwdOpSrc.getOp();
92#ifdef HAVE_MUELU_TPETRA
93 if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isTpetra(fwdOp))
return true;
96#ifdef HAVE_MUELU_EPETRA
97 if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isEpetra(fwdOp))
return true;
104 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
105 RCP<PreconditionerBase<Scalar> > MueLuMaxwell1PreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::createPrec()
const {
106 return Teuchos::rcp(
new DefaultPreconditioner<Scalar>);
109 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
110 void MueLuMaxwell1PreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::
111 initializePrec(
const RCP<
const LinearOpSourceBase<Scalar> >& fwdOpSrc, PreconditionerBase<Scalar>* prec,
const ESupportSolveUse supportSolveUse)
const {
114 typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
115 typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpThyUtils;
116 typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
117 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
119#if defined(MUELU_CAN_USE_MIXED_PRECISION)
120 typedef Xpetra::TpetraHalfPrecisionOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpHalfPrecOp;
121 typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMV;
122 typedef typename XpHalfPrecOp::HalfScalar HalfScalar;
123 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
124 typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
125 typedef Xpetra::MultiVector<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XphMV;
126 typedef Xpetra::MultiVector<Magnitude,LocalOrdinal,GlobalOrdinal,Node > XpmMV;
127 typedef Xpetra::MultiVector<HalfMagnitude,LocalOrdinal,GlobalOrdinal,Node > XphmMV;
128 typedef Xpetra::Matrix<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XphMat;
130 Teuchos::TimeMonitor tM(*Teuchos::TimeMonitor::getNewTimer(std::string(
"ThyraMueLuMaxwell1::initializePrec")));
133 TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
134 TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
135 TEUCHOS_ASSERT(prec);
138 ParameterList paramList = *paramList_;
141 const RCP<const ThyLinOpBase> fwdOp = fwdOpSrc->getOp();
142 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
145 bool bIsEpetra = XpThyUtils::isEpetra(fwdOp);
146 bool bIsTpetra = XpThyUtils::isTpetra(fwdOp);
147 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra ==
true && bIsTpetra ==
true));
151 RCP<XpMat> A = XpThyUtils::toXpetra(Teuchos::rcp_const_cast<ThyLinOpBase>(fwdOp));
152 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A));
155 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(
dynamic_cast<DefaultPreconditioner<Scalar> *
>(prec));
156 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
159 RCP<ThyLinOpBase> thyra_precOp = Teuchos::null;
160 thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar> >(defaultPrec->getNonconstUnspecifiedPrecOp(),
true);
165 const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter(
"Maxwell1: enable reuse") || !paramList.get<
bool>(
"Maxwell1: enable reuse"));
166 const bool useHalfPrecision = paramList.get<
bool>(
"half precision",
false) && bIsTpetra;
169 if (startingOver ==
true) {
172 std::list<std::string> convertXpetra = {
"Coordinates",
"Nullspace",
"Kn",
"D0"};
173 for (
auto it = convertXpetra.begin(); it != convertXpetra.end(); ++it)
174 Converters<Scalar,LocalOrdinal,GlobalOrdinal,Node>::replaceWithXpetra(paramList,*it);
176 std::list<std::string> sublists = {
"maxwell1: 11list",
"maxwell1: 22list"};
177 for (
auto itSublist = sublists.begin(); itSublist != sublists.end(); ++itSublist)
178 if (paramList.isSublist(*itSublist)) {
179 ParameterList& sublist = paramList.sublist(*itSublist);
180 for (
int lvlNo=0; lvlNo < 10; ++lvlNo) {
181 if (sublist.isSublist(
"level " + std::to_string(lvlNo) +
" user data")) {
182 ParameterList& lvlList = sublist.sublist(
"level " + std::to_string(lvlNo) +
" user data");
183 std::list<std::string> convertKeys;
184 for (
auto it = lvlList.begin(); it != lvlList.end(); ++it)
185 convertKeys.push_back(lvlList.name(it));
186 for (
auto it = convertKeys.begin(); it != convertKeys.end(); ++it)
187 Converters<Scalar,LocalOrdinal,GlobalOrdinal,Node>::replaceWithXpetra(lvlList,*it);
192 ParameterList& sublist = paramList.sublist(
"maxwell1: 11list");
193 if (sublist.isParameter(
"D0")) {
194 Converters<Scalar,LocalOrdinal,GlobalOrdinal,Node>::replaceWithXpetra(sublist,
"D0");
197 paramList.set<
bool>(
"Maxwell1: use as preconditioner",
true);
198 if (useHalfPrecision) {
199#if defined(MUELU_CAN_USE_MIXED_PRECISION)
202 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
203 if (paramList.isType<RCP<XpmMV> >(
"Coordinates")) {
204 RCP<XpmMV> coords = paramList.get<RCP<XpmMV> >(
"Coordinates");
205 paramList.remove(
"Coordinates");
206 RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
207 paramList.set(
"Coordinates",halfCoords);
209 if (paramList.isType<RCP<XpMV> >(
"Nullspace")) {
210 RCP<XpMV> nullspace = paramList.get<RCP<XpMV> >(
"Nullspace");
211 paramList.remove(
"Nullspace");
212 RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
213 paramList.set(
"Nullspace",halfNullspace);
215 std::list<std::string> convertMat = {
"Kn",
"D0"};
216 for (
auto it = convertMat.begin(); it != convertMat.end(); ++it) {
217 if (paramList.isType<RCP<XpMat> >(*it)) {
218 RCP<XpMat> M = paramList.get<RCP<XpMat> >(*it);
219 paramList.remove(*it);
220 RCP<XphMat> halfM = Xpetra::convertToHalfPrecision(M);
221 paramList.set(*it,halfM);
227 xpPrecOp = rcp(
new XpHalfPrecOp(halfPrec));
229 TEUCHOS_TEST_FOR_EXCEPT(
true);
235 xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
240 RCP<ThyXpOp> thyXpOp = rcp_dynamic_cast<ThyXpOp>(thyra_precOp,
true);
241 RCP<XpOp> xpOp = thyXpOp->getXpetraOperator();
242#if defined(MUELU_CAN_USE_MIXED_PRECISION)
243 RCP<XpHalfPrecOp> xpHalfPrecOp = rcp_dynamic_cast<XpHalfPrecOp>(xpOp);
244 if (!xpHalfPrecOp.is_null()) {
245 RCP<MueLu::Maxwell1<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> > preconditioner = rcp_dynamic_cast<MueLu::Maxwell1<HalfScalar,LocalOrdinal,GlobalOrdinal,Node>>(xpHalfPrecOp->GetHalfPrecisionOperator(),
true);
246 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
247 preconditioner->resetMatrix(halfA);
248 xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
252 RCP<MueLu::Maxwell1<Scalar,LocalOrdinal,GlobalOrdinal,Node> > preconditioner = rcp_dynamic_cast<MueLu::Maxwell1<Scalar,LocalOrdinal,GlobalOrdinal,Node>>(xpOp,
true);
253 preconditioner->resetMatrix(A);
254 xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
259 RCP<const VectorSpaceBase<Scalar> > thyraRangeSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpPrecOp->getRangeMap());
260 RCP<const VectorSpaceBase<Scalar> > thyraDomainSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpPrecOp->getDomainMap());
262 RCP<ThyLinOpBase > thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace, xpPrecOp);
263 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraPrecOp));
265 defaultPrec->initializeUnspecified(thyraPrecOp);
269 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
270 void MueLuMaxwell1PreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::
271 uninitializePrec(PreconditionerBase<Scalar>* prec, RCP<
const LinearOpSourceBase<Scalar> >* fwdOp, ESupportSolveUse* supportSolveUse)
const {
272 TEUCHOS_ASSERT(prec);
275 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(
dynamic_cast<DefaultPreconditioner<Scalar> *
>(prec));
276 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
280 *fwdOp = Teuchos::null;
283 if (supportSolveUse) {
285 *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
288 defaultPrec->uninitialize();
293 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
294 void MueLuMaxwell1PreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::setParameterList(RCP<ParameterList>
const& paramList) {
295 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
296 paramList_ = paramList;
299 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
300 RCP<ParameterList> MueLuMaxwell1PreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getNonconstParameterList() {
304 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
305 RCP<ParameterList> MueLuMaxwell1PreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::unsetParameterList() {
306 RCP<ParameterList> savedParamList = paramList_;
307 paramList_ = Teuchos::null;
308 return savedParamList;
311 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
312 RCP<const ParameterList> MueLuMaxwell1PreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getParameterList()
const {
316 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
317 RCP<const ParameterList> MueLuMaxwell1PreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getValidParameters()
const {
318 static RCP<const ParameterList> validPL;
320 if (Teuchos::is_null(validPL))
321 validPL = rcp(
new ParameterList());
327 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
328 std::string MueLuMaxwell1PreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::description()
const {
329 return "Thyra::MueLuMaxwell1PreconditionerFactory";
Preconditioner (wrapped as a Xpetra::Operator) for Maxwell's equations in curl-curl form.
Concrete Thyra::LinearOpBase subclass for Xpetra::Operator.