MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_EpetraOperator.cpp
Go to the documentation of this file.
1// @HEADER
2//
3// ***********************************************************************
4//
5// MueLu: A package for multigrid based preconditioning
6// Copyright 2012 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact
39// Jonathan Hu (jhu@sandia.gov)
40// Andrey Prokopenko (aprokop@sandia.gov)
41// Ray Tuminaro (rstumin@sandia.gov)
42//
43// ***********************************************************************
44//
45// @HEADER
46#include <Xpetra_Matrix.hpp>
47#include <Xpetra_CrsMatrixWrap.hpp>
48#include <Xpetra_BlockedCrsMatrix.hpp>
49#include <Xpetra_EpetraMultiVector.hpp>
50
52#include "MueLu_Level.hpp"
53
54#if defined(HAVE_MUELU_SERIAL) and defined(HAVE_MUELU_EPETRA)
55
56namespace MueLu {
57
58int EpetraOperator::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
59 try {
60 // There is no rcpFromRef(const T&), so we need to do const_cast
61 const Xpetra::EpetraMultiVectorT<GO,NO> eX(rcpFromRef(const_cast<Epetra_MultiVector&>(X)));
62 Xpetra::EpetraMultiVectorT<GO,NO> eY(rcpFromRef(Y));
63 // Generally, we assume two different vectors, but AztecOO uses a single vector
64 if (X.Values() == Y.Values()) {
65 // X and Y point to the same memory, use an additional vector
66 RCP<Xpetra::EpetraMultiVectorT<GO,NO> > tmpY = Teuchos::rcp(new Xpetra::EpetraMultiVectorT<GO,NO>(eY.getMap(), eY.getNumVectors()));
67 // InitialGuessIsZero in MueLu::Hierarchy.Iterate() does not zero out components, it
68 // only assumes that user provided an already zeroed out vector
69 bool initialGuessZero = true;
70 tmpY->putScalar(0.0);
71 // apply one V-cycle as preconditioner
72 Hierarchy_->Iterate(eX, *tmpY, 1, initialGuessZero);
73 // deep copy solution from MueLu
74 eY.update(1.0, *tmpY, 0.0);
75 } else {
76 // X and Y point to different memory, pass the vectors through
77
78 // InitialGuessIsZero in MueLu::Hierarchy.Iterate() does not zero out components, it
79 // only assumes that user provided an already zeroed out vector
80 bool initialGuessZero = true;
81 eY.putScalar(0.0);
82 Hierarchy_->Iterate(eX, eY, 1, initialGuessZero);
83 }
84
85 } catch (std::exception& e) {
86 //TODO: error msg directly on std::cerr?
87 std::cerr << "Caught an exception in MueLu::EpetraOperator::ApplyInverse():" << std::endl
88 << e.what() << std::endl;
89 return -1;
90 }
91 return 0;
92}
93
94const Epetra_Comm& EpetraOperator::Comm() const {
95 RCP<Matrix> A = Hierarchy_->GetLevel(0)->Get<RCP<Matrix> >("A");
96
97 //TODO: This code is not pretty
98 RCP<Xpetra::BlockedCrsMatrix<SC, LO, GO, NO> > epbA = Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<SC, LO, GO, NO> >(A);
99 if (epbA != Teuchos::null) {
100 RCP<const Xpetra::Matrix<SC, LO, GO, NO> > blockMat = epbA->getMatrix(0,0);
101 RCP<const Xpetra::CrsMatrixWrap<SC, LO, GO, NO> > blockCrsWrap = Teuchos::rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<SC, LO, GO, NO> >(blockMat);
102 if (blockCrsWrap == Teuchos::null)
103 throw Exceptions::BadCast("MueLu::EpetraOperator::Comm(): Cast from block (0,0) to CrsMatrixWrap failed. Could be a block matrix. TODO implement recursive support for block matrices.");
104 RCP<const Xpetra::EpetraCrsMatrixT<GO,NO>> tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GO,NO> >(blockCrsWrap->getCrsMatrix());
105 if (tmp_ECrsMtx == Teuchos::null)
106 throw Exceptions::BadCast("MueLu::EpetraOperator::Comm(): Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
107 RCP<Epetra_CrsMatrix> epA = tmp_ECrsMtx->getEpetra_CrsMatrixNonConst();
108 return epA->Comm();
109 }
110
111 RCP<const Xpetra::CrsMatrixWrap<SC,LO,GO,NO> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<SC,LO,GO,NO> >(A);
112 if (crsOp == Teuchos::null)
113 throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
114 const RCP<const Xpetra::EpetraCrsMatrixT<GO,NO>> &tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GO,NO>>(crsOp->getCrsMatrix());
115 if (tmp_ECrsMtx == Teuchos::null)
116 throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
117 return tmp_ECrsMtx->getEpetra_CrsMatrixNonConst()->Comm();
118}
119
120const Epetra_Map& EpetraOperator::OperatorDomainMap() const {
121 RCP<Xpetra::Matrix<SC,LO,GO,NO> > A = Hierarchy_->GetLevel(0)->Get<RCP<Matrix> >("A");
122
123 RCP<Xpetra::BlockedCrsMatrix<SC, LO, GO, NO> > epbA = Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<SC, LO, GO, NO> >(A);
124 if (epbA != Teuchos::null)
125 return Xpetra::toEpetra(epbA->getFullDomainMap()); // TODO check me
126
127 RCP<const Xpetra::CrsMatrixWrap<SC,LO,GO,NO> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<SC,LO,GO,NO> >(A);
128 if (crsOp == Teuchos::null)
129 throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
130 const RCP<const Xpetra::EpetraCrsMatrixT<GO,NO>> &tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GO,NO>>(crsOp->getCrsMatrix());
131 if (tmp_ECrsMtx == Teuchos::null)
132 throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
133 return tmp_ECrsMtx->getEpetra_CrsMatrixNonConst()->DomainMap();
134}
135
136const Epetra_Map & EpetraOperator::OperatorRangeMap() const {
137 RCP<Xpetra::Matrix<SC,LO,GO,NO> > A = Hierarchy_->GetLevel(0)->Get<RCP<Matrix> >("A");
138
139 RCP<Xpetra::BlockedCrsMatrix<SC, LO, GO, NO> > epbA = Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<SC, LO, GO, NO> >(A);
140 if (epbA != Teuchos::null)
141 return Xpetra::toEpetra(epbA->getFullRangeMap());
142
143 RCP<const Xpetra::CrsMatrixWrap<SC,LO,GO,NO> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<SC,LO,GO,NO> >(A);
144 if (crsOp == Teuchos::null)
145 throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
146 const RCP<const Xpetra::EpetraCrsMatrixT<GO,NO>> &tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GO,NO>>(crsOp->getCrsMatrix());
147 if (tmp_ECrsMtx == Teuchos::null)
148 throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
149 return tmp_ECrsMtx->getEpetra_CrsMatrixNonConst()->RangeMap();
150}
151
152} // namespace
153
154#endif // #if defined(HAVE_MUELU_SERIAL) and defined(HAVE_MUELU_EPETRA)
double * Values() const
Namespace for MueLu classes and methods.