MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_AztecEpetraOperator.cpp
Go to the documentation of this file.
1
2#ifndef PACKAGES_MUELU_ADAPTERS_AZTECOO_MUELU_AZTECEPETRAOPERATOR_CPP_
3#define PACKAGES_MUELU_ADAPTERS_AZTECOO_MUELU_AZTECEPETRAOPERATOR_CPP_
4
5#include "Xpetra_EpetraMultiVector.hpp"
6
7#include "MueLu_config.hpp" // for HAVE_MUELU_DEBUG
8#include "MueLu_RefMaxwell.hpp"
9
11
12#if defined(HAVE_MUELU_SERIAL) and defined(HAVE_MUELU_EPETRA)
13
14namespace MueLu {
15
16int AztecEpetraOperator::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
17 try {
18 // There is no rcpFromRef(const T&), so we need to do const_cast
19 const Xpetra::EpetraMultiVectorT<GO,NO> eX(Teuchos::rcpFromRef(const_cast<Epetra_MultiVector&>(X)));
20 Xpetra::EpetraMultiVectorT<GO,NO> eY(Teuchos::rcpFromRef(Y));
21 // Generally, we assume two different vectors, but AztecOO uses a single vector
22 if (X.Values() == Y.Values()) {
23 // X and Y point to the same memory, use an additional vector
24 Teuchos::RCP<Xpetra::EpetraMultiVectorT<GO,NO> > tmpY = Teuchos::rcp(new Xpetra::EpetraMultiVectorT<GO,NO>(eY.getMap(), eY.getNumVectors()));
25 tmpY->putScalar(0.0);
26 xOp_->apply(eX, *tmpY);
27 // deep copy solution from MueLu
28 eY.update(1.0, *tmpY, 0.0);
29 } else {
30 // X and Y point to different memory, pass the vectors through
31 eY.putScalar(0.0);
32 xOp_->apply(eX, eY);
33 }
34
35 } catch (std::exception& e) {
36 //TODO: error msg directly on std::cerr?
37 std::cerr << "Caught an exception in MueLu::AztecEpetraOperator::ApplyInverse():" << std::endl
38 << e.what() << std::endl;
39 return -1;
40 }
41 return 0;
42}
43
44const Epetra_Comm& AztecEpetraOperator::Comm() const {
45 const Epetra_Map emap = Xpetra::toEpetra(xOp_->getDomainMap());
46 return emap.Comm();
47}
48
49const Epetra_Map& AztecEpetraOperator::OperatorDomainMap() const {
50 if(Teuchos::rcp_dynamic_cast<MueLu::RefMaxwell<SC,LO,GO,NO> >(xOp_) != Teuchos::null) {
51 RCP<Xpetra::Matrix<SC,LO,GO,NO> > A = Teuchos::rcp_dynamic_cast<MueLu::RefMaxwell<SC,LO,GO,NO> >(xOp_)->getJacobian();
52 RCP<const Xpetra::CrsMatrixWrap<SC,LO,GO,NO> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<SC,LO,GO,NO> >(A);
53 if (crsOp == Teuchos::null)
54 throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
55 const RCP<const Xpetra::EpetraCrsMatrixT<GO,NO>> &tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GO,NO>>(crsOp->getCrsMatrix());
56 if (tmp_ECrsMtx == Teuchos::null)
57 throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
58 return tmp_ECrsMtx->getEpetra_CrsMatrixNonConst()->DomainMap();
59 }
60
61 // TAW there is a problem with the Map version of toEeptra leading to code crashes
62 // we probably need to create a new copy of the Epetra_Map
63 Teuchos::RCP<const Map> map = xOp_->getDomainMap();
64 return Xpetra::toEpetra(map);
65}
66
67const Epetra_Map & AztecEpetraOperator::OperatorRangeMap() const {
68
69 if(Teuchos::rcp_dynamic_cast<MueLu::RefMaxwell<SC,LO,GO,NO> >(xOp_) != Teuchos::null) {
70 RCP<Xpetra::Matrix<SC,LO,GO,NO> > A = Teuchos::rcp_dynamic_cast<MueLu::RefMaxwell<SC,LO,GO,NO> >(xOp_)->getJacobian();
71 RCP<const Xpetra::CrsMatrixWrap<SC,LO,GO,NO> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<SC,LO,GO,NO> >(A);
72 if (crsOp == Teuchos::null)
73 throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
74 const RCP<const Xpetra::EpetraCrsMatrixT<GO,NO>> &tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GO,NO>>(crsOp->getCrsMatrix());
75 if (tmp_ECrsMtx == Teuchos::null)
76 throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
77 return tmp_ECrsMtx->getEpetra_CrsMatrixNonConst()->RangeMap();
78 }
79
80 // TAW there is a problem with the Map version of toEeptra leading to code crashes
81 // we probably need to create a new copy of the Epetra_Map
82 Teuchos::RCP<const Map> map = xOp_->getRangeMap();
83 return Xpetra::toEpetra(map);
84}
85
86}
87
88#endif /*#if defined(HAVE_MUELU_SERIAL) and defined(HAVE_MUELU_EPETRA)*/
89
90#endif /* PACKAGES_MUELU_ADAPTERS_AZTECOO_MUELU_AZTECEPETRAOPERATOR_CPP_ */
const Epetra_Comm & Comm() const
double * Values() const
Preconditioner (wrapped as a Xpetra::Operator) for Maxwell's equations in curl-curl form.
Namespace for MueLu classes and methods.