MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_ShiftedLaplacianOperator_def.hpp
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
47#ifndef MUELU_SHIFTEDLAPLACIANOPERATOR_DEF_HPP
48#define MUELU_SHIFTEDLAPLACIANOPERATOR_DEF_HPP
49
50#include "MueLu_ConfigDefs.hpp"
51
52#ifdef HAVE_MUELU_TPETRA
53
54#include <Xpetra_Matrix.hpp>
55#include <Xpetra_CrsMatrixWrap.hpp>
56#include <Xpetra_BlockedCrsMatrix.hpp>
57#include <Xpetra_TpetraMultiVector.hpp>
58#include <Xpetra_MultiVectorFactory.hpp>
59
61#include "MueLu_Hierarchy.hpp"
62#include "MueLu_Utilities.hpp"
63
64
65namespace MueLu {
66
67// ------------- getDomainMap -----------------------
68
69template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
70Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >
72getDomainMap () const
73{
74 typedef Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> XMatrix;
75
76 RCP<MueLu::Level> L0 = Hierarchy_->GetLevel (0);
77 RCP<XMatrix> A = L0->Get<RCP<XMatrix> > ("A");
78
79 RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpbA =
80 Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(A);
81 if (tpbA != Teuchos::null) {
82 return Xpetra::toTpetraNonZero (tpbA->getDomainMap ());
83 }
84
85 RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpA =
87 return tpA->getDomainMap ();
88}
89
90// ------------- getRangeMap -----------------------
91
92template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
93Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >
95getRangeMap () const
96{
97 typedef Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> XMatrix;
98
99 RCP<MueLu::Level> L0 = Hierarchy_->GetLevel(0);
100 RCP<XMatrix> A = L0->Get< RCP<XMatrix> >("A");
101
102 RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpbA =
103 Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(A);
104 if(tpbA != Teuchos::null)
105 return Xpetra::toTpetraNonZero(tpbA->getRangeMap());
106
107 RCP< Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpA =
109 return tpA->getRangeMap();
110}
111
112// ------------- apply -----------------------
113
114template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
115void ShiftedLaplacianOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node>::apply(const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X,
116 Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y,
117 Teuchos::ETransp /* mode */, Scalar /* alpha */, Scalar /* beta */) const {
118
119 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TMV;
120 typedef Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XTMV;
121 // typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XMV; // unused
122
123 TMV & temp_x = const_cast<TMV &>(X);
124 const XTMV tX(rcpFromRef(temp_x));
125 XTMV tY(rcpFromRef(Y));
126
127 try {
128 tY.putScalar(0.0);
129 Hierarchy_->Iterate(tX, tY, cycles_, true);
130 }
131
132 catch(std::exception& e) {
133 //FIXME add message and rethrow
134 std::cerr << "Caught an exception in MueLu::ShiftedLaplacianOperator::ApplyInverse():" << std::endl
135 << e.what() << std::endl;
136 }
137
138 // update solution with 2-grid error correction
139 /*if(option_==1) {
140 for(int j=0; j<cycles_; j++) {
141 RCP<XMV> residual = MueLu::Utilities<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Residual(*A_, tY, tX);
142 RCP<XMV> coarseResidual = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(R_->getRangeMap(), tX.getNumVectors());
143 RCP<XMV> coarseError = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(R_->getRangeMap(), tX.getNumVectors());
144 R_ -> apply(*residual, *coarseResidual, Teuchos::NO_TRANS, (Scalar) 1.0, (Scalar) 0.0);
145 RCP<TMV> tcoarseR = MueLu::Utilities<Scalar,LocalOrdinal,GlobalOrdinal,Node>::MV2NonConstTpetraMV(coarseResidual);
146 RCP<TMV> tcoarseE = MueLu::Utilities<Scalar,LocalOrdinal,GlobalOrdinal,Node>::MV2NonConstTpetraMV(coarseError);
147 BelosLP_ -> setProblem(tcoarseE,tcoarseR);
148 BelosSM_ -> solve();
149 RCP<XMV> fineError = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(P_->getRangeMap(), tX.getNumVectors());
150 XTMV tmpcoarseE(rcpFromRef(*tcoarseE));
151 P_ -> apply(tmpcoarseE, *fineError, Teuchos::NO_TRANS, (Scalar) 1.0, (Scalar) 0.0);
152 tY.update((Scalar) 1.0, *fineError, (Scalar) 1.0);
153 }
154 }
155
156 try {
157 Hierarchy_->Iterate(tX, tY, 1, false);
158 }
159
160 catch(std::exception& e) {
161 //FIXME add message and rethrow
162 std::cerr << "Caught an exception in MueLu::ShiftedLaplacianOperator::ApplyInverse():" << std::endl
163 << e.what() << std::endl;
164 }*/
165
166}
167
168// ------------- apply -----------------------
169template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
173
174} // namespace
175#endif //ifdef HAVE_MUELU_TPETRA
176
177#endif //ifdef MUELU_SHIFTEDLAPLACIANOPERATOR_DEF_HPP
MueLu::DefaultScalar Scalar
Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getRangeMap() const
Returns the Tpetra::Map object associated with the range of this operator.
bool hasTransposeApply() const
Indicates whether this operator supports applying the adjoint operator.
Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getDomainMap() const
Returns the Tpetra::Map object associated with the domain of this operator.
void apply(const Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=Teuchos::ScalarTraits< Scalar >::one(), Scalar beta=Teuchos::ScalarTraits< Scalar >::one()) const
Returns in Y the result of a Tpetra::Operator applied to a Tpetra::MultiVector X.
static RCP< Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
Namespace for MueLu classes and methods.