MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
6 #include "Xpetra_CrsMatrixWrap.hpp"
8 
9 #include "MueLu_config.hpp" // for HAVE_MUELU_DEBUG
10 #include "MueLu_RefMaxwell.hpp"
11 #include "MueLu_Exceptions.hpp"
12 
14 
15 #if defined(HAVE_MUELU_SERIAL) and defined(HAVE_MUELU_EPETRA)
16 
17 namespace MueLu {
18 
19 int AztecEpetraOperator::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
20  try {
21  // There is no rcpFromRef(const T&), so we need to do const_cast
22  const Xpetra::EpetraMultiVectorT<GO, NO> eX(Teuchos::rcpFromRef(const_cast<Epetra_MultiVector&>(X)));
23  Xpetra::EpetraMultiVectorT<GO, NO> eY(Teuchos::rcpFromRef(Y));
24  // Generally, we assume two different vectors, but AztecOO uses a single vector
25  if (X.Values() == Y.Values()) {
26  // X and Y point to the same memory, use an additional vector
28  tmpY->putScalar(0.0);
29  xOp_->apply(eX, *tmpY);
30  // deep copy solution from MueLu
31  eY.update(1.0, *tmpY, 0.0);
32  } else {
33  // X and Y point to different memory, pass the vectors through
34  eY.putScalar(0.0);
35  xOp_->apply(eX, eY);
36  }
37 
38  } catch (std::exception& e) {
39  // TODO: error msg directly on std::cerr?
40  std::cerr << "Caught an exception in MueLu::AztecEpetraOperator::ApplyInverse():" << std::endl
41  << e.what() << std::endl;
42  return -1;
43  }
44  return 0;
45 }
46 
47 const Epetra_Comm& AztecEpetraOperator::Comm() const {
48  const Epetra_Map emap = Xpetra::toEpetra(xOp_->getDomainMap());
49  return emap.Comm();
50 }
51 
52 const Epetra_Map& AztecEpetraOperator::OperatorDomainMap() const {
53  if (Teuchos::rcp_dynamic_cast<MueLu::RefMaxwell<SC, LO, GO, NO>>(xOp_) != Teuchos::null) {
54  RCP<Xpetra::Matrix<SC, LO, GO, NO>> A = Teuchos::rcp_dynamic_cast<MueLu::RefMaxwell<SC, LO, GO, NO>>(xOp_)->getJacobian();
55  RCP<const Xpetra::CrsMatrixWrap<SC, LO, GO, NO>> crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<SC, LO, GO, NO>>(A);
56  if (crsOp == Teuchos::null)
57  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
58  const RCP<const Xpetra::EpetraCrsMatrixT<GO, NO>>& tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GO, NO>>(crsOp->getCrsMatrix());
59  if (tmp_ECrsMtx == Teuchos::null)
60  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
61  return tmp_ECrsMtx->getEpetra_CrsMatrixNonConst()->DomainMap();
62  }
63 
64  // TAW there is a problem with the Map version of toEeptra leading to code crashes
65  // we probably need to create a new copy of the Epetra_Map
66  Teuchos::RCP<const Map> map = xOp_->getDomainMap();
67  return Xpetra::toEpetra(map);
68 }
69 
70 const Epetra_Map& AztecEpetraOperator::OperatorRangeMap() const {
71  if (Teuchos::rcp_dynamic_cast<MueLu::RefMaxwell<SC, LO, GO, NO>>(xOp_) != Teuchos::null) {
72  RCP<Xpetra::Matrix<SC, LO, GO, NO>> A = Teuchos::rcp_dynamic_cast<MueLu::RefMaxwell<SC, LO, GO, NO>>(xOp_)->getJacobian();
73  RCP<const Xpetra::CrsMatrixWrap<SC, LO, GO, NO>> crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<SC, LO, GO, NO>>(A);
74  if (crsOp == Teuchos::null)
75  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
76  const RCP<const Xpetra::EpetraCrsMatrixT<GO, NO>>& tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GO, NO>>(crsOp->getCrsMatrix());
77  if (tmp_ECrsMtx == Teuchos::null)
78  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
79  return tmp_ECrsMtx->getEpetra_CrsMatrixNonConst()->RangeMap();
80  }
81 
82  // TAW there is a problem with the Map version of toEeptra leading to code crashes
83  // we probably need to create a new copy of the Epetra_Map
84  Teuchos::RCP<const Map> map = xOp_->getRangeMap();
85  return Xpetra::toEpetra(map);
86 }
87 
88 } // namespace MueLu
89 
90 #endif /*#if defined(HAVE_MUELU_SERIAL) and defined(HAVE_MUELU_EPETRA)*/
91 
92 #endif /* PACKAGES_MUELU_ADAPTERS_AZTECOO_MUELU_AZTECEPETRAOPERATOR_CPP_ */
void update(const Scalar &alpha, const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Scalar &beta)
void putScalar(const Scalar &value)
const Epetra_CrsGraph & toEpetra(const RCP< const CrsGraph< int, GlobalOrdinal, Node > > &graph)
Preconditioner (wrapped as a Xpetra::Operator) for Maxwell&#39;s equations in curl-curl form...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
size_t getNumVectors() const
const Epetra_Comm & Comm() const
Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const