46 #ifndef THYRA_MUELU_REFMAXWELL_PRECONDITIONER_FACTORY_DEF_HPP 47 #define THYRA_MUELU_REFMAXWELL_PRECONDITIONER_FACTORY_DEF_HPP 51 #ifdef HAVE_MUELU_STRATIMIKOS 57 using Teuchos::ParameterList;
62 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
64 paramList_(rcp(new ParameterList()))
69 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
71 const RCP<const LinearOpBase<Scalar> > fwdOp = fwdOpSrc.getOp();
73 #ifdef HAVE_MUELU_TPETRA 74 if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isTpetra(fwdOp))
return true;
77 if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isBlockedOperator(fwdOp))
return true;
83 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
85 return Teuchos::rcp(
new DefaultPreconditioner<Scalar>);
88 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
90 initializePrec(
const RCP<
const LinearOpSourceBase<Scalar> >& fwdOpSrc, PreconditionerBase<Scalar>* prec,
const ESupportSolveUse supportSolveUse)
const {
91 using Teuchos::rcp_dynamic_cast;
94 typedef Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> XpMap;
95 typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
96 typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpThyUtils;
97 typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
98 typedef Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpBlockedCrsMat;
99 typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
100 typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMultVec;
101 typedef Xpetra::MultiVector<double,LocalOrdinal,GlobalOrdinal,Node> XpMultVecDouble;
102 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
103 typedef Thyra::DiagonalLinearOpBase<Scalar> ThyDiagLinOpBase;
104 #ifdef HAVE_MUELU_TPETRA 105 typedef Thyra::TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyTpLinOp;
107 Teuchos::TimeMonitor tM(*Teuchos::TimeMonitor::getNewTimer(std::string(
"ThyraMueLuRefMaxwell::initializePrec")));
110 TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
111 TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
112 TEUCHOS_ASSERT(prec);
115 ParameterList paramList = *paramList_;
118 const RCP<const ThyLinOpBase> fwdOp = fwdOpSrc->getOp();
119 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
122 bool bIsEpetra = XpThyUtils::isEpetra(fwdOp);
123 bool bIsTpetra = XpThyUtils::isTpetra(fwdOp);
124 bool bIsBlocked = XpThyUtils::isBlockedOperator(fwdOp);
125 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra ==
true && bIsTpetra ==
true));
126 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == bIsTpetra) && bIsBlocked ==
false);
127 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra != bIsTpetra) && bIsBlocked ==
true);
129 RCP<XpMat> A = Teuchos::null;
131 Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
132 Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<Scalar> >(fwdOp);
133 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(ThyBlockedOp));
135 TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==
false);
137 Teuchos::RCP<const LinearOpBase<Scalar> > b00 = ThyBlockedOp->getBlock(0,0);
138 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
140 RCP<const XpCrsMat > xpetraFwdCrsMat00 = XpThyUtils::toXpetra(b00);
141 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat00));
144 RCP<XpCrsMat> xpetraFwdCrsMatNonConst00 = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat00);
145 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst00));
148 RCP<XpMat> A00 = rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst00));
149 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A00));
151 RCP<const XpMap> rowmap00 = A00->getRowMap();
152 RCP< const Teuchos::Comm< int > > comm = rowmap00->getComm();
155 RCP<XpBlockedCrsMat> bMat = Teuchos::rcp(
new XpBlockedCrsMat(ThyBlockedOp, comm));
156 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(bMat));
161 RCP<const XpCrsMat > xpetraFwdCrsMat = XpThyUtils::toXpetra(fwdOp);
162 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat));
165 RCP<XpCrsMat> xpetraFwdCrsMatNonConst = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat);
166 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst));
169 A = rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst));
171 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A));
174 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(
dynamic_cast<DefaultPreconditioner<Scalar> *
>(prec));
175 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
178 RCP<ThyLinOpBase> thyra_precOp = Teuchos::null;
179 thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar> >(defaultPrec->getNonconstUnspecifiedPrecOp(),
true);
182 RCP<MueLu::RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node> > preconditioner = Teuchos::null;
187 const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter(
"reuse: type") || paramList.get<std::string>(
"reuse: type") ==
"none");
189 if (startingOver ==
true) {
191 Teuchos::RCP<XpMultVecDouble> coordinates = Teuchos::null;
193 Teuchos::TimeMonitor tM_coords(*Teuchos::TimeMonitor::getNewTimer(std::string(
"ThyraMueLuRefMaxwell::initializePrec get coords")));
195 paramList.set<RCP<XpMultVecDouble> >(
"Coordinates", coordinates);
199 #ifdef HAVE_MUELU_TPETRA 201 typedef Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tV;
202 typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tMV;
203 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpCrsMat;
204 Teuchos::TimeMonitor tMwrap(*Teuchos::TimeMonitor::getNewTimer(std::string(
"ThyraMueLuRefMaxwell::initializePrec wrap objects")));
205 if (paramList.isType<Teuchos::RCP<tMV> >(
"Nullspace")) {
206 Teuchos::TimeMonitor tM_nullspace(*Teuchos::TimeMonitor::getNewTimer(std::string(
"ThyraMueLuRefMaxwell::initializePrec wrap nullspace")));
207 RCP<tMV> tpetra_nullspace = paramList.get<RCP<tMV> >(
"Nullspace");
208 paramList.remove(
"Nullspace");
209 RCP<XpMultVec> nullspace = MueLu::TpetraMultiVector_To_XpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tpetra_nullspace);
210 paramList.set<RCP<XpMultVec> >(
"Nullspace", nullspace);
211 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(nullspace));
214 if (paramList.isParameter(
"M1")) {
215 if (paramList.isType<Teuchos::RCP<TpCrsMat> >(
"M1")) {
216 Teuchos::TimeMonitor tM_M1(*Teuchos::TimeMonitor::getNewTimer(std::string(
"ThyraMueLuRefMaxwell::initializePrec wrap M1")));
217 RCP<TpCrsMat> tM1 = paramList.get<RCP<TpCrsMat> >(
"M1");
218 paramList.remove(
"M1");
219 RCP<XpCrsMat> xM1 = rcp_dynamic_cast<XpCrsMat>(tM1,
true);
220 paramList.set<RCP<XpCrsMat> >(
"M1", xM1);
221 }
else if (paramList.isType<Teuchos::RCP<const ThyLinOpBase> >(
"M1")) {
222 Teuchos::TimeMonitor tM_M1(*Teuchos::TimeMonitor::getNewTimer(std::string(
"ThyraMueLuRefMaxwell::initializePrec wrap M1")));
223 RCP<const ThyLinOpBase> thyM1 = paramList.get<RCP<const ThyLinOpBase> >(
"M1");
224 paramList.remove(
"M1");
225 RCP<const XpCrsMat> crsM1 = XpThyUtils::toXpetra(thyM1);
226 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsM1));
228 RCP<XpCrsMat> crsM1NonConst = Teuchos::rcp_const_cast<XpCrsMat>(crsM1);
229 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsM1NonConst));
231 RCP<XpMat> M1 = rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(crsM1NonConst));
232 paramList.set<RCP<XpMat> >(
"M1", M1);
233 }
else if (paramList.isType<Teuchos::RCP<XpMat> >(
"M1")) {
240 if (paramList.isParameter(
"D0")) {
241 if (paramList.isType<Teuchos::RCP<TpCrsMat> >(
"D0")) {
242 Teuchos::TimeMonitor tM_D0(*Teuchos::TimeMonitor::getNewTimer(std::string(
"ThyraMueLuRefMaxwell::initializePrec wrap D0")));
243 RCP<TpCrsMat> tD0 = paramList.get<RCP<TpCrsMat> >(
"D0");
244 paramList.remove(
"D0");
245 RCP<XpCrsMat> xD0 = rcp_dynamic_cast<XpCrsMat>(tD0,
true);
246 paramList.set<RCP<XpCrsMat> >(
"D0", xD0);
247 }
else if (paramList.isType<Teuchos::RCP<const ThyLinOpBase> >(
"D0")) {
248 Teuchos::TimeMonitor tM_D0(*Teuchos::TimeMonitor::getNewTimer(std::string(
"ThyraMueLuRefMaxwell::initializePrec wrap D0")));
249 RCP<const ThyLinOpBase> thyD0 = paramList.get<RCP<const ThyLinOpBase> >(
"D0");
250 paramList.remove(
"D0");
251 RCP<const XpCrsMat> crsD0 = XpThyUtils::toXpetra(thyD0);
252 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsD0));
254 RCP<XpCrsMat> crsD0NonConst = Teuchos::rcp_const_cast<XpCrsMat>(crsD0);
255 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsD0NonConst));
257 RCP<XpMat> D0 = rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(crsD0NonConst));
258 paramList.set<RCP<XpMat> >(
"D0", D0);
259 }
else if (paramList.isType<Teuchos::RCP<XpMat> >(
"D0")) {
266 if (paramList.isParameter(
"M0inv")) {
267 if (paramList.isType<Teuchos::RCP<TpCrsMat> >(
"M0inv")) {
268 Teuchos::TimeMonitor tM_M0inv(*Teuchos::TimeMonitor::getNewTimer(std::string(
"ThyraMueLuRefMaxwell::initializePrec wrap M0inv")));
269 RCP<TpCrsMat> tM0inv = paramList.get<RCP<TpCrsMat> >(
"M0inv");
270 paramList.remove(
"M0inv");
271 RCP<XpCrsMat> xM0inv = rcp_dynamic_cast<XpCrsMat>(tM0inv,
true);
272 paramList.set<RCP<XpCrsMat> >(
"M0inv", xM0inv);
273 }
else if (paramList.isType<Teuchos::RCP<const ThyDiagLinOpBase> >(
"M0inv")) {
274 Teuchos::TimeMonitor tM_M0inv(*Teuchos::TimeMonitor::getNewTimer(std::string(
"ThyraMueLuRefMaxwell::initializePrec wrap M0inv")));
275 RCP<const ThyDiagLinOpBase> thyM0inv = paramList.get<RCP<const ThyDiagLinOpBase> >(
"M0inv");
276 paramList.remove(
"M0inv");
277 RCP<const Thyra::VectorBase<Scalar> > diag = thyM0inv->getDiag();
278 RCP<const tV> tDiag = Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getConstTpetraVector(diag);
279 RCP<XpMat> M0inv = Xpetra::MatrixFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(Xpetra::toXpetra(tDiag));
280 paramList.set<RCP<XpMat> >(
"M0inv", M0inv);
281 }
else if (paramList.isType<Teuchos::RCP<const ThyLinOpBase> >(
"M0inv")) {
282 Teuchos::TimeMonitor tM_M0inv(*Teuchos::TimeMonitor::getNewTimer(std::string(
"ThyraMueLuRefMaxwell::initializePrec wrap M0inv")));
283 RCP<const ThyLinOpBase> thyM0inv = paramList.get<RCP<const ThyLinOpBase> >(
"M0inv");
284 paramList.remove(
"M0inv");
285 RCP<const XpCrsMat> crsM0inv = XpThyUtils::toXpetra(thyM0inv);
286 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsM0inv));
288 RCP<XpCrsMat> crsM0invNonConst = Teuchos::rcp_const_cast<XpCrsMat>(crsM0inv);
289 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsM0invNonConst));
291 RCP<XpMat> M0inv = rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(crsM0invNonConst));
292 paramList.set<RCP<XpMat> >(
"M0inv", M0inv);
293 }
else if (paramList.isType<Teuchos::RCP<XpMat> >(
"M0inv")) {
305 Teuchos::TimeMonitor tMbuild(*Teuchos::TimeMonitor::getNewTimer(std::string(
"ThyraMueLuRefMaxwell::initializePrec build prec")));
306 paramList.set<
bool>(
"refmaxwell: use as preconditioner",
true);
314 #if defined(HAVE_MUELU_TPETRA) 317 RCP<ThyTpLinOp> tpetr_precOp = rcp_dynamic_cast<ThyTpLinOp>(thyra_precOp);
327 RCP<ThyLinOpBase > thyraPrecOp = Teuchos::null;
328 RCP<const VectorSpaceBase<Scalar> > thyraRangeSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(preconditioner->getRangeMap());
329 RCP<const VectorSpaceBase<Scalar> > thyraDomainSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(preconditioner->getDomainMap());
331 RCP<XpOp> xpOp = Teuchos::rcp_dynamic_cast<XpOp>(preconditioner);
332 thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace,xpOp);
334 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraPrecOp));
336 defaultPrec->initializeUnspecified(thyraPrecOp);
340 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
342 uninitializePrec(PreconditionerBase<Scalar>* prec, RCP<
const LinearOpSourceBase<Scalar> >* fwdOp, ESupportSolveUse* supportSolveUse)
const {
343 TEUCHOS_ASSERT(prec);
346 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(
dynamic_cast<DefaultPreconditioner<Scalar> *
>(prec));
347 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
351 *fwdOp = Teuchos::null;
354 if (supportSolveUse) {
356 *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
359 defaultPrec->uninitialize();
364 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
366 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
367 paramList_ = paramList;
370 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
375 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
377 RCP<ParameterList> savedParamList = paramList_;
378 paramList_ = Teuchos::null;
379 return savedParamList;
382 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
387 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
389 static RCP<const ParameterList> validPL;
391 if (Teuchos::is_null(validPL))
392 validPL = rcp(
new ParameterList());
398 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
400 return "Thyra::MueLuRefMaxwellPreconditionerFactory";
404 #endif // HAVE_MUELU_STRATIMIKOS 406 #endif // ifdef THYRA_MUELU_REFMAXWELL_PRECONDITIONER_FACTORY_DEF_HPP Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
Teuchos::RCP< const Teuchos::ParameterList > getParameterList() const
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
static RCP< Xpetra::MultiVector< double, LocalOrdinal, GlobalOrdinal, Node > > ExtractCoordinatesFromParameterList(ParameterList ¶mList)
Preconditioner (wrapped as a Xpetra::Operator) for Maxwell's equations in curl-curl form...
std::string description() const
MueLuRefMaxwellPreconditionerFactory()
void initializePrec(const Teuchos::RCP< const LinearOpSourceBase< Scalar > > &fwdOp, PreconditionerBase< Scalar > *prec, const ESupportSolveUse supportSolveUse) const
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > ¶mList)
Teuchos::RCP< PreconditionerBase< Scalar > > createPrec() const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
bool isCompatible(const LinearOpSourceBase< Scalar > &fwdOp) const
Exception throws to report errors in the internal logical of the program.
void uninitializePrec(PreconditionerBase< Scalar > *prec, Teuchos::RCP< const LinearOpSourceBase< Scalar > > *fwdOp, ESupportSolveUse *supportSolveUse) const