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 
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 
14 namespace MueLu {
15 
16 int 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
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 
44 const Epetra_Comm& AztecEpetraOperator::Comm() const {
45  const Epetra_Map emap = Xpetra::toEpetra(xOp_->getDomainMap());
46  return emap.Comm();
47 }
48 
49 const 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::EpetraCrsMatrix> &tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::EpetraCrsMatrix>(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 
67 const 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::EpetraCrsMatrix> &tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::EpetraCrsMatrix>(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_ */
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