47 #ifndef THYRA_MUELU_REFMAXWELL_PRECONDITIONER_FACTORY_DECL_HPP 48 #define THYRA_MUELU_REFMAXWELL_PRECONDITIONER_FACTORY_DECL_HPP 52 #ifdef HAVE_MUELU_STRATIMIKOS 55 #include "Thyra_DefaultPreconditioner.hpp" 56 #include "Thyra_BlockedLinearOpBase.hpp" 57 #include "Thyra_DiagonalLinearOpBase.hpp" 59 #ifdef HAVE_MUELU_TPETRA 60 #include "Thyra_TpetraLinearOp.hpp" 61 #include "Thyra_TpetraThyraWrappers.hpp" 63 #ifdef HAVE_MUELU_EPETRA 64 #include "Thyra_EpetraLinearOp.hpp" 65 #include "Thyra_EpetraThyraWrappers.hpp" 68 #include "Teuchos_Ptr.hpp" 69 #include "Teuchos_TestForException.hpp" 70 #include "Teuchos_Assert.hpp" 71 #include "Teuchos_Time.hpp" 73 #include <Xpetra_CrsMatrixWrap.hpp> 74 #include <Xpetra_CrsMatrix.hpp> 75 #include <Xpetra_Matrix.hpp> 76 #include <Xpetra_ThyraUtils.hpp> 78 #include <MueLu_Hierarchy.hpp> 80 #include <MueLu_HierarchyUtils.hpp> 81 #include <MueLu_Utilities.hpp> 82 #include <MueLu_ParameterListInterpreter.hpp> 83 #include <MueLu_MLParameterListInterpreter.hpp> 86 #include <MueLu_RefMaxwell.hpp> 87 #ifdef HAVE_MUELU_TPETRA 88 #include <MueLu_TpetraOperator.hpp> 90 #ifdef HAVE_MUELU_EPETRA 92 #include <Xpetra_EpetraOperator.hpp> 95 #include "Thyra_PreconditionerFactoryBase.hpp" 97 #include "Kokkos_DefaultNode.hpp" 110 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node = KokkosClassic::DefaultNode::DefaultNodeType>
125 bool isCompatible(
const LinearOpSourceBase<Scalar>& fwdOp)
const;
127 Teuchos::RCP<PreconditionerBase<Scalar> >
createPrec()
const;
129 void initializePrec(
const Teuchos::RCP<
const LinearOpSourceBase<Scalar> >& fwdOp,
130 PreconditionerBase<Scalar>* prec,
131 const ESupportSolveUse supportSolveUse
135 Teuchos::RCP<
const LinearOpSourceBase<Scalar> >* fwdOp,
136 ESupportSolveUse* supportSolveUse
145 void setParameterList(
const Teuchos::RCP<Teuchos::ParameterList>& paramList);
172 #ifdef HAVE_MUELU_EPETRA 201 const RCP<const LinearOpBase<Scalar> > fwdOp = fwdOpSrc.getOp();
203 #ifdef HAVE_MUELU_TPETRA 204 if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isTpetra(fwdOp))
return true;
207 #ifdef HAVE_MUELU_EPETRA 208 if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isEpetra(fwdOp))
return true;
211 if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isBlockedOperator(fwdOp))
return true;
217 Teuchos::RCP<PreconditionerBase<Scalar> >
createPrec()
const {
218 return Teuchos::rcp(
new DefaultPreconditioner<Scalar>);
222 void initializePrec(
const Teuchos::RCP<
const LinearOpSourceBase<Scalar> >& fwdOpSrc,
223 PreconditionerBase<Scalar>* prec,
224 const ESupportSolveUse supportSolveUse
226 using Teuchos::rcp_dynamic_cast;
229 typedef Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> XpMap;
230 typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
231 typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpThyUtils;
232 typedef Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMatWrap;
233 typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
234 typedef Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpBlockedCrsMat;
235 typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
236 typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMultVec;
237 typedef Xpetra::MultiVector<double,LocalOrdinal,GlobalOrdinal,Node> XpMultVecDouble;
238 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
239 typedef Thyra::DiagonalLinearOpBase<Scalar> ThyDiagLinOpBase;
240 #ifdef HAVE_MUELU_TPETRA 242 #if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT))) || \ 243 (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT)))) 244 typedef Thyra::TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyTpLinOp;
247 #if defined(HAVE_MUELU_EPETRA) 248 typedef Thyra::EpetraLinearOp ThyEpLinOp;
249 typedef Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> XpEpCrsMat;
257 TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
259 TEUCHOS_ASSERT(prec);
265 const RCP<const ThyLinOpBase> fwdOp = fwdOpSrc->getOp();
266 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
269 bool bIsEpetra = XpThyUtils::isEpetra(fwdOp);
270 bool bIsTpetra = XpThyUtils::isTpetra(fwdOp);
271 bool bIsBlocked = XpThyUtils::isBlockedOperator(fwdOp);
272 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra ==
true && bIsTpetra ==
true));
273 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == bIsTpetra) && bIsBlocked ==
false);
274 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra != bIsTpetra) && bIsBlocked ==
true);
276 RCP<XpMat> A = Teuchos::null;
278 Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
279 Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<Scalar> >(fwdOp);
280 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(ThyBlockedOp));
282 TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==
false);
284 Teuchos::RCP<const LinearOpBase<Scalar> > b00 = ThyBlockedOp->getBlock(0,0);
285 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
287 RCP<const XpCrsMat > xpetraFwdCrsMat00 = XpThyUtils::toXpetra(b00);
288 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat00));
291 RCP<XpCrsMat> xpetraFwdCrsMatNonConst00 = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat00);
292 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst00));
295 RCP<XpMat> A00 = rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst00));
296 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A00));
298 RCP<const XpMap> rowmap00 = A00->getRowMap();
299 RCP< const Teuchos::Comm< int > > comm = rowmap00->getComm();
302 RCP<XpBlockedCrsMat> bMat = Teuchos::rcp(
new XpBlockedCrsMat(ThyBlockedOp, comm));
303 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(bMat));
308 RCP<const XpCrsMat > xpetraFwdCrsMat = XpThyUtils::toXpetra(fwdOp);
309 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat));
312 RCP<XpCrsMat> xpetraFwdCrsMatNonConst = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat);
313 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst));
316 A = rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst));
318 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A));
321 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(
dynamic_cast<DefaultPreconditioner<Scalar> *
>(prec));
322 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
325 RCP<ThyLinOpBase> thyra_precOp = Teuchos::null;
326 thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar> >(defaultPrec->getNonconstUnspecifiedPrecOp(),
true);
329 RCP<MueLu::RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node> > preconditioner = Teuchos::null;
334 const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter(
"reuse: type") || paramList.get<std::string>(
"reuse: type") ==
"none");
336 if (startingOver ==
true) {
338 Teuchos::RCP<XpMultVecDouble> coordinates = Teuchos::null;
340 paramList.set<RCP<XpMultVecDouble> >(
"Coordinates", coordinates);
343 #ifdef HAVE_MUELU_TPETRA 345 #if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT))) || \ 346 (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT)))) 347 typedef Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tV;
348 typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tMV;
349 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpCrsMat;
350 if (paramList.isType<Teuchos::RCP<tMV> >(
"Nullspace")) {
351 RCP<tMV> tpetra_nullspace = paramList.get<RCP<tMV> >(
"Nullspace");
352 paramList.remove(
"Nullspace");
353 RCP<XpMultVec> nullspace = MueLu::TpetraMultiVector_To_XpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tpetra_nullspace);
354 paramList.set<RCP<XpMultVec> >(
"Nullspace", nullspace);
355 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(nullspace));
358 if (paramList.isParameter(
"M1")) {
359 if (paramList.isType<Teuchos::RCP<TpCrsMat> >(
"M1")) {
360 RCP<TpCrsMat> tM1 = paramList.get<RCP<TpCrsMat> >(
"M1");
361 paramList.remove(
"M1");
362 RCP<XpCrsMat> xM1 = rcp_dynamic_cast<XpCrsMat>(tM1,
true);
363 paramList.set<RCP<XpCrsMat> >(
"M1", xM1);
364 }
else if (paramList.isType<Teuchos::RCP<const ThyLinOpBase> >(
"M1")) {
365 RCP<const ThyLinOpBase> thyM1 = paramList.get<RCP<const ThyLinOpBase> >(
"M1");
366 paramList.remove(
"M1");
367 RCP<const XpCrsMat> crsM1 = XpThyUtils::toXpetra(thyM1);
368 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsM1));
370 RCP<XpCrsMat> crsM1NonConst = Teuchos::rcp_const_cast<XpCrsMat>(crsM1);
371 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsM1NonConst));
373 RCP<XpMat> M1 = rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(crsM1NonConst));
374 paramList.set<RCP<XpMat> >(
"M1", M1);
375 }
else if (paramList.isType<Teuchos::RCP<XpMat> >(
"M1")) {
382 if (paramList.isParameter(
"D0")) {
383 if(paramList.isType<Teuchos::RCP<TpCrsMat> >(
"D0")) {
384 RCP<TpCrsMat> tD0 = paramList.get<RCP<TpCrsMat> >(
"D0");
385 paramList.remove(
"D0");
386 RCP<XpCrsMat> xD0 = rcp_dynamic_cast<XpCrsMat>(tD0,
true);
387 paramList.set<RCP<XpCrsMat> >(
"D0", xD0);
388 }
else if (paramList.isType<Teuchos::RCP<const ThyLinOpBase> >(
"D0")) {
389 RCP<const ThyLinOpBase> thyD0 = paramList.get<RCP<const ThyLinOpBase> >(
"D0");
390 paramList.remove(
"D0");
391 RCP<const XpCrsMat> crsD0 = XpThyUtils::toXpetra(thyD0);
392 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsD0));
394 RCP<XpCrsMat> crsD0NonConst = Teuchos::rcp_const_cast<XpCrsMat>(crsD0);
395 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsD0NonConst));
397 RCP<XpMat> D0 = rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(crsD0NonConst));
398 paramList.set<RCP<XpMat> >(
"D0", D0);
399 }
else if (paramList.isType<Teuchos::RCP<XpMat> >(
"D0")) {
406 if (paramList.isParameter(
"M0inv")) {
407 if (paramList.isType<Teuchos::RCP<TpCrsMat> >(
"M0inv")) {
408 RCP<TpCrsMat> tM0inv = paramList.get<RCP<TpCrsMat> >(
"M0inv");
409 paramList.remove(
"M0inv");
410 RCP<XpCrsMat> xM0inv = rcp_dynamic_cast<XpCrsMat>(tM0inv,
true);
411 paramList.set<RCP<XpCrsMat> >(
"M0inv", xM0inv);
412 }
else if (paramList.isType<Teuchos::RCP<const ThyDiagLinOpBase> >(
"M0inv")) {
413 RCP<const ThyDiagLinOpBase> thyM0inv = paramList.get<RCP<const ThyDiagLinOpBase> >(
"M0inv");
414 paramList.remove(
"M0inv");
415 RCP<const Thyra::VectorBase<Scalar> > diag = thyM0inv->getDiag();
416 RCP<const Thyra::TpetraVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > ttDiag = rcp_dynamic_cast<
const Thyra::TpetraVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(diag);
417 RCP<const tV> tDiag = Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getConstTpetraVector(diag);
418 RCP<XpMat> M0inv = Xpetra::MatrixFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(Xpetra::toXpetra(tDiag));
419 paramList.set<RCP<XpMat> >(
"M0inv", M0inv);
420 }
else if (paramList.isType<Teuchos::RCP<const ThyLinOpBase> >(
"M0inv")) {
421 RCP<const ThyLinOpBase> thyM0inv = paramList.get<RCP<const ThyLinOpBase> >(
"M0inv");
422 paramList.remove(
"M0inv");
423 RCP<const XpCrsMat> crsM0inv = XpThyUtils::toXpetra(thyM0inv);
424 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsM0inv));
426 RCP<XpCrsMat> crsM0invNonConst = Teuchos::rcp_const_cast<XpCrsMat>(crsM0inv);
427 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsM0invNonConst));
429 RCP<XpMat> M0inv = rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(crsM0invNonConst));
430 paramList.set<RCP<XpMat> >(
"M0inv", M0inv);
431 }
else if (paramList.isType<Teuchos::RCP<XpMat> >(
"M0inv")) {
439 "Thyra::MueLuRefMaxwellPreconditionerFactory: Tpetra does not support GO=int and or EpetraNode.");
443 #ifdef HAVE_MUELU_EPETRA 445 if (paramList.isType<RCP<Epetra_MultiVector> >(
"Nullspace")) {
446 RCP<Epetra_MultiVector> epetra_nullspace = Teuchos::null;
447 epetra_nullspace = paramList.get<RCP<Epetra_MultiVector> >(
"Nullspace");
448 paramList.remove(
"Nullspace");
449 RCP<Xpetra::EpetraMultiVectorT<int,Node> > xpEpNullspace = Teuchos::rcp(
new Xpetra::EpetraMultiVectorT<int,Node>(epetra_nullspace));
450 RCP<Xpetra::MultiVector<double,int,int,Node> > xpEpNullspaceMult = rcp_dynamic_cast<Xpetra::MultiVector<double,int,int,Node> >(xpEpNullspace,
true);
451 RCP<XpMultVec> nullspace = rcp_dynamic_cast<XpMultVec>(xpEpNullspaceMult,
true);
452 paramList.set<RCP<XpMultVec> >(
"Nullspace", nullspace);
455 if (paramList.isParameter(
"M1")) {
456 if (paramList.isType<Teuchos::RCP<Epetra_CrsMatrix> >(
"M1")) {
457 RCP<Epetra_CrsMatrix> eM1 = paramList.get<RCP<Epetra_CrsMatrix> >(
"M1");
458 paramList.remove(
"M1");
459 RCP<XpEpCrsMat> xeM1 = Teuchos::rcp(
new XpEpCrsMat(eM1));
460 RCP<XpCrsMat> xCrsM1 = rcp_dynamic_cast<XpCrsMat>(xeM1,
true);
461 RCP<XpCrsMatWrap> xwM1 = Teuchos::rcp(
new XpCrsMatWrap(xCrsM1));
462 RCP<XpMat> xM1 = rcp_dynamic_cast<XpMat>(xwM1);
463 paramList.set<RCP<XpMat> >(
"M1", xM1);
465 else if (paramList.isType<Teuchos::RCP<const ThyLinOpBase> >(
"M1")) {
466 RCP<const ThyLinOpBase> thyM1 = paramList.get<RCP<const ThyLinOpBase> >(
"M1");
467 paramList.remove(
"M1");
468 RCP<const XpCrsMat> crsM1 = XpThyUtils::toXpetra(thyM1);
469 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsM1));
471 RCP<XpCrsMat> crsM1NonConst = Teuchos::rcp_const_cast<XpCrsMat>(crsM1);
472 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsM1NonConst));
474 RCP<XpMat> M1 = rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(crsM1NonConst));
475 paramList.set<RCP<XpMat> >(
"M1", M1);
476 }
else if (paramList.isType<Teuchos::RCP<XpMat> >(
"M1")) {
483 if (paramList.isParameter(
"D0")) {
484 if (paramList.isType<Teuchos::RCP<Epetra_CrsMatrix> >(
"D0")) {
485 RCP<Epetra_CrsMatrix> eD0 = paramList.get<RCP<Epetra_CrsMatrix> >(
"D0");
486 paramList.remove(
"D0");
487 RCP<XpEpCrsMat> xeD0 = Teuchos::rcp(
new XpEpCrsMat(eD0));
488 RCP<XpCrsMat> xCrsD0 = rcp_dynamic_cast<XpCrsMat>(xeD0,
true);
489 RCP<XpCrsMatWrap> xwD0 = Teuchos::rcp(
new XpCrsMatWrap(xCrsD0));
490 RCP<XpMat> xD0 = rcp_dynamic_cast<XpMat>(xwD0);
491 paramList.set<RCP<XpMat> >(
"D0", xD0);
493 else if (paramList.isType<Teuchos::RCP<const ThyLinOpBase> >(
"D0")) {
494 RCP<const ThyLinOpBase> thyD0 = paramList.get<RCP<const ThyLinOpBase> >(
"D0");
495 paramList.remove(
"D0");
496 RCP<const XpCrsMat> crsD0 = XpThyUtils::toXpetra(thyD0);
497 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsD0));
499 RCP<XpCrsMat> crsD0NonConst = Teuchos::rcp_const_cast<XpCrsMat>(crsD0);
500 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsD0NonConst));
502 RCP<XpMat> D0 = rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(crsD0NonConst));
503 paramList.set<RCP<XpMat> >(
"D0", D0);
504 }
else if (paramList.isType<Teuchos::RCP<XpMat> >(
"D0")) {
511 if (paramList.isParameter(
"M0inv")) {
512 if (paramList.isType<Teuchos::RCP<Epetra_CrsMatrix> >(
"M0inv")) {
513 RCP<Epetra_CrsMatrix> eM0inv = paramList.get<RCP<Epetra_CrsMatrix> >(
"M0inv");
514 paramList.remove(
"M0inv");
515 RCP<XpEpCrsMat> xeM0inv = Teuchos::rcp(
new XpEpCrsMat(eM0inv));
516 RCP<XpCrsMat> xCrsM0inv = rcp_dynamic_cast<XpCrsMat>(xeM0inv,
true);
517 RCP<XpCrsMatWrap> xwM0inv = Teuchos::rcp(
new XpCrsMatWrap(xCrsM0inv));
518 RCP<XpMat> xM0inv = rcp_dynamic_cast<XpMat>(xwM0inv);
519 paramList.set<RCP<XpMat> >(
"M0inv", xM0inv);
521 else if (paramList.isType<Teuchos::RCP<const ThyDiagLinOpBase> >(
"M0inv")) {
522 RCP<const ThyDiagLinOpBase> thyM0inv = paramList.get<RCP<const ThyDiagLinOpBase> >(
"M0inv");
523 paramList.remove(
"M0inv");
525 RCP<const Teuchos::Comm<int> > comm = A->getDomainMap()->getComm();
526 RCP<const Epetra_Map> map = Thyra::get_Epetra_Map(*(thyM0inv->range()), Xpetra::toEpetra(comm));
528 RCP<const Thyra::VectorBase<double> > diag = thyM0inv->getDiag();
529 RCP<const Epetra_Vector> eDiag = Thyra::get_Epetra_Vector(*map, diag);
530 RCP<Epetra_Vector> nceDiag = Teuchos::rcp_const_cast<Epetra_Vector>(eDiag);
531 RCP<Xpetra::EpetraVectorT<int,Node> > xpEpDiag = Teuchos::rcp(
new Xpetra::EpetraVectorT<int,Node>(nceDiag));
532 RCP<const Xpetra::Vector<double,int,int,Node> > xpDiag = rcp_dynamic_cast<
const Xpetra::Vector<double,int,int,Node> >(xpEpDiag,
true);
533 RCP<XpMat> M0inv = Xpetra::MatrixFactory<double,int,int,Node>::Build(xpDiag);
534 paramList.set<RCP<XpMat> >(
"M0inv", M0inv);
535 }
else if (paramList.isType<Teuchos::RCP<const ThyLinOpBase> >(
"M0inv")) {
536 RCP<const ThyLinOpBase> thyM0inv = paramList.get<RCP<const ThyLinOpBase> >(
"M0inv");
537 paramList.remove(
"M0inv");
538 RCP<const XpCrsMat> crsM0inv = XpThyUtils::toXpetra(thyM0inv);
539 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsM0inv));
541 RCP<XpCrsMat> crsM0invNonConst = Teuchos::rcp_const_cast<XpCrsMat>(crsM0inv);
542 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsM0invNonConst));
544 RCP<XpMat> M0inv = rcp(
new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(crsM0invNonConst));
545 paramList.set<RCP<XpMat> >(
"M0inv", M0inv);
546 }
else if (paramList.isType<Teuchos::RCP<XpMat> >(
"M0inv")) {
555 paramList.set<
bool>(
"refmaxwell: use as preconditioner",
true);
562 #if defined(HAVE_MUELU_TPETRA) 564 #if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT))) || \ 565 (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT)))) 566 RCP<ThyTpLinOp> tpetr_precOp = rcp_dynamic_cast<ThyTpLinOp>(thyra_precOp);
570 "Thyra::MueLuRefMaxwellPreconditionerFactory: Tpetra does not support GO=int and or EpetraNode.");
574 #if defined(HAVE_MUELU_EPETRA)// && defined(HAVE_MUELU_SERIAL) 576 RCP<ThyEpLinOp> epetr_precOp = rcp_dynamic_cast<ThyEpLinOp>(thyra_precOp);
585 RCP<ThyLinOpBase > thyraPrecOp = Teuchos::null;
586 RCP<const VectorSpaceBase<Scalar> > thyraRangeSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(preconditioner->getRangeMap());
587 RCP<const VectorSpaceBase<Scalar> > thyraDomainSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(preconditioner->getDomainMap());
589 RCP<XpOp> xpOp = Teuchos::rcp_dynamic_cast<XpOp>(preconditioner);
590 thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace,xpOp);
592 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraPrecOp));
594 defaultPrec->initializeUnspecified(thyraPrecOp);
599 Teuchos::RCP<
const LinearOpSourceBase<Scalar> >* fwdOp,
600 ESupportSolveUse* supportSolveUse
602 TEUCHOS_ASSERT(prec);
605 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(
dynamic_cast<DefaultPreconditioner<Scalar> *
>(prec));
606 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
610 *fwdOp = Teuchos::null;
613 if (supportSolveUse) {
615 *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
618 defaultPrec->uninitialize();
628 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
633 RCP<ParameterList> savedParamList =
paramList_;
635 return savedParamList;
643 static RCP<const ParameterList> validPL;
645 if (Teuchos::is_null(validPL))
646 validPL = rcp(
new ParameterList());
656 std::string
description()
const {
return "Thyra::MueLuRefMaxwellPreconditionerFactory"; }
666 #endif // HAVE_MUELU_EPETRA 670 #endif // #ifdef HAVE_MUELU_STRATIMIKOS 672 #endif // THYRA_MUELU_REFMAXWELL_PRECONDITIONER_FACTORY_DECL_HPP
void uninitializePrec(PreconditionerBase< Scalar > *prec, Teuchos::RCP< const LinearOpSourceBase< Scalar > > *fwdOp, ESupportSolveUse *supportSolveUse) const
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
Teuchos::RCP< const Teuchos::ParameterList > getParameterList() const
void initializePrec(const Teuchos::RCP< const LinearOpSourceBase< Scalar > > &fwdOpSrc, PreconditionerBase< Scalar > *prec, const ESupportSolveUse supportSolveUse) const
MueLuRefMaxwellPreconditionerFactory()
Concrete preconditioner factory subclass for Thyra based on MueLu.Add support for MueLu preconditione...
std::string description() const
Teuchos::RCP< Teuchos::ParameterList > paramList_
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
static RCP< Xpetra::MultiVector< double, LocalOrdinal, GlobalOrdinal, Node > > ExtractCoordinatesFromParameterList(ParameterList ¶mList)
bool isCompatible(const LinearOpSourceBase< Scalar > &fwdOpSrc) const
Preconditioner (wrapped as a Xpetra::Operator) for Maxwell's equations in curl-curl form...
std::string description() const
MueLuRefMaxwellPreconditionerFactory()
Teuchos::RCP< Teuchos::ParameterList > paramList_
void initializePrec(const Teuchos::RCP< const LinearOpSourceBase< Scalar > > &fwdOp, PreconditionerBase< Scalar > *prec, const ESupportSolveUse supportSolveUse) const
Teuchos::RCP< PreconditionerBase< Scalar > > createPrec() 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
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > ¶mList)
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
Exception throws to report errors in the internal logical of the program.
Kokkos::Compat::KokkosSerialWrapperNode EpetraNode
Teuchos::RCP< const Teuchos::ParameterList > getParameterList() const
void uninitializePrec(PreconditionerBase< Scalar > *prec, Teuchos::RCP< const LinearOpSourceBase< Scalar > > *fwdOp, ESupportSolveUse *supportSolveUse) const