MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_BelosSmoother_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#ifndef MUELU_BELOSSMOOTHER_DEF_HPP
47#define MUELU_BELOSSMOOTHER_DEF_HPP
48
49#include "MueLu_ConfigDefs.hpp"
50
51#if defined(HAVE_MUELU_BELOS)
52
53#include <Teuchos_ParameterList.hpp>
54
55#include <Xpetra_CrsMatrix.hpp>
56#include <Xpetra_Matrix.hpp>
57#include <Xpetra_MultiVectorFactory.hpp>
58#ifdef HAVE_XPETRA_TPETRA
59#include <Xpetra_TpetraMultiVector.hpp>
60#endif
61
63#include "MueLu_Level.hpp"
65#include "MueLu_Utilities.hpp"
66#include "MueLu_Monitor.hpp"
67
68#include <BelosLinearProblem.hpp>
69#include <BelosSolverFactory.hpp>
70
71
72
73namespace MueLu {
74
75 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
76 BelosSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BelosSmoother(const std::string type, const Teuchos::ParameterList& paramList)
77 : type_(type)
78 {
79 bool solverSupported = false;
80#ifdef HAVE_MUELU_TPETRA
81 Belos::SolverFactory<Scalar,tMV,tOP> tFactory;
82 solverSupported = solverSupported || tFactory.isSupported(type);
83#endif
84 // if (!solverSupported) {
85 // Teuchos::Array<std::string> supportedSolverNames = factory.supportedSolverNames();
86
87 // std::ostringstream outString;
88 // outString << "[";
89 // for (auto iter = supportedSolverNames.begin(); iter != supportedSolverNames.end(); ++iter) {
90 // outString << "\"" << *iter << "\"";
91 // if (iter + 1 != supportedSolverNames.end()) {
92 // outString << ", ";
93 // }
94 // }
95 // outString << "]";
96
97 // TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "Belos does not provide the solver '" << type_ << "'.\nSupported Solvers: " << outString.str());
98 // }
99 this->declareConstructionOutcome(!solverSupported, "Belos does not provide the smoother '" + type_ + "'.");
100 if (solverSupported)
101 SetParameterList(paramList);
102 }
103
104 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
106 Factory::SetParameterList(paramList);
107 }
108
109 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
111
112 this->Input(currentLevel, "A");
113
114 }
115
116 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
118 FactoryMonitor m(*this, "Setup Smoother", currentLevel);
119
120 A_ = Factory::Get< RCP<Matrix> >(currentLevel, "A");
121 SetupBelos(currentLevel);
123 this->GetOStream(Statistics1) << description() << std::endl;
124 }
125
126
127 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
129
130 bool useTpetra = A_->getRowMap()->lib() == Xpetra::UseTpetra;
131
132 if (useTpetra) {
133#ifdef HAVE_MUELU_TPETRA
134 tBelosProblem_ = rcp(new Belos::LinearProblem<Scalar, tMV, tOP>());
135 RCP<tOP> tA = Utilities::Op2NonConstTpetraCrs(A_);
136 tBelosProblem_->setOperator(tA);
137
138 Belos::SolverFactory<SC, tMV, tOP> solverFactory;
139 tSolver_ = solverFactory.create(type_, rcpFromRef(const_cast<ParameterList&>(this->GetParameterList())));
140 tSolver_->setProblem(tBelosProblem_);
141#endif
142 } else {
143 TEUCHOS_ASSERT(false);
144 }
145 }
146
147 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
148 void BelosSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Apply(MultiVector& X, const MultiVector& B, bool InitialGuessIsZero) const {
149 TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError, "MueLu::BelosSmoother::Apply(): Setup() has not been called");
150
151 if (A_->getRowMap()->lib() == Xpetra::UseTpetra) {
152#ifdef HAVE_MUELU_TPETRA
153 if (InitialGuessIsZero) {
154 X.putScalar(0.0);
155
156 RCP<Tpetra::MultiVector<SC,LO,GO,NO> > tpX = rcpFromRef(Utilities::MV2NonConstTpetraMV(X));
157 RCP<const Tpetra::MultiVector<SC,LO,GO,NO> > tpB = rcpFromRef(Utilities::MV2TpetraMV(B));
158
159 tBelosProblem_->setInitResVec(tpB);
160 tBelosProblem_->setProblem(tpX, tpB);
161 tSolver_->solve();
162
163 } else {
164 typedef Teuchos::ScalarTraits<Scalar> TST;
165 RCP<MultiVector> Residual = Utilities::Residual(*A_, X, B);
166 RCP<MultiVector> Correction = MultiVectorFactory::Build(A_->getDomainMap(), X.getNumVectors());
167
168 RCP<Tpetra::MultiVector<SC,LO,GO,NO> > tpX = rcpFromRef(Utilities::MV2NonConstTpetraMV(*Correction));
169 RCP<const Tpetra::MultiVector<SC,LO,GO,NO> > tpB = rcpFromRef(Utilities::MV2TpetraMV(*Residual));
170
171 tBelosProblem_->setInitResVec(tpB);
172 tBelosProblem_->setProblem(tpX, tpB);
173 tSolver_->solve();
174
175 X.update(TST::one(), *Correction, TST::one());
176 }
177#endif
178 } else {
179 TEUCHOS_ASSERT(false);
180 }
181
182 }
183
184 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
185 RCP<MueLu::SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node> > BelosSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Copy() const {
186 RCP<BelosSmoother> smoother = rcp(new BelosSmoother(*this) );
187 smoother->SetParameterList(this->GetParameterList());
188 return smoother;
189 }
190
191 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
193 std::ostringstream out;
195 if (A_->getRowMap()->lib() == Xpetra::UseTpetra) {
196#ifdef MUELU_HAVE_TPETRA
197 out << tSolver_->description();
198#endif
199 }
200 } else {
201 out << "BELOS {type = " << type_ << "}";
202 }
203 return out.str();
204 }
205
206 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
207 void BelosSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::print(Teuchos::FancyOStream &out, const VerbLevel verbLevel) const {
209
210 if (verbLevel & Parameters1) {
211 out0 << "Parameter list: " << std::endl;
212 Teuchos::OSTab tab2(out);
213 out << this->GetParameterList();
214 }
215
216 if (verbLevel & External) {
217#ifdef MUELU_HAVE_TPETRA
218 if (tSolver_ != Teuchos::null) {
219 Teuchos::OSTab tab2(out);
220 out << *tSolver_ << std::endl;
221 }
222#endif
223 }
224
225 if (verbLevel & Debug) {
226 if (A_->getRowMap()->lib() == Xpetra::UseTpetra) {
227#ifdef MUELU_HAVE_TPETRA
228 out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl
229 << "-" << std::endl
230 << "RCP<solver_>: " << tSolver_ << std::endl;
231#endif
232 }
233 }
234 }
235
236 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
238 return Teuchos::OrdinalTraits<size_t>::invalid();
239 }
240
241
242} // namespace MueLu
243
244#endif // HAVE_MUELU_BELOS
245#endif // MUELU_BELOSSMOOTHER_DEF_HPP
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
Class that encapsulates Belos smoothers.
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the preconditioner.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
friend class BelosSmoother
Constructor.
void SetParameterList(const Teuchos::ParameterList &paramList)
Set parameters from a parameter list and return with default values.
void DeclareInput(Level &currentLevel) const
Input.
void SetupBelos(Level &currentLevel)
std::string description() const
Return a simple one-line description of this object.
size_t getNodeSmootherComplexity() const
For diagnostic purposes.
void Setup(Level &currentLevel)
Set up the smoother.
RCP< SmootherPrototype > Copy() const
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
virtual void SetParameterList(const Teuchos::ParameterList &paramList)=0
Set parameters from a parameter list and return with default values.
void declareConstructionOutcome(bool fail, std::string msg)
bool IsSetup() const
Get the state of a smoother prototype.
static RCP< Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< const Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2TpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > const vec)
Helper utility to pull out the underlying Tpetra objects from an Xpetra object.
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
static RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &RHS)
Namespace for MueLu classes and methods.
@ Debug
Print additional debugging information.
@ Statistics1
Print more statistics.
@ External
Print external lib objects.
@ Parameters1
Print class parameters (more parameters, more verbose)