MueLu  Version of the Day
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 
73 namespace 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
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.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Print external lib objects.
bool IsSetup() const
Get the state of a smoother prototype.
size_t getNodeSmootherComplexity() const
For diagnostic purposes.
std::string description() const
Return a simple one-line description of this object.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that encapsulates Belos smoothers.
Print more statistics.
Print additional debugging information.
void Setup(Level &currentLevel)
Set up the smoother.
Namespace for MueLu classes and methods.
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)
virtual void SetParameterList(const Teuchos::ParameterList &paramList)
Set parameters from a parameter list and return with default values.
void declareConstructionOutcome(bool fail, std::string msg)
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
static RCP< Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
void SetupBelos(Level &currentLevel)
friend class BelosSmoother
Constructor.
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
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.
Print class parameters (more parameters, more verbose)
Exception throws to report errors in the internal logical of the program.
RCP< SmootherPrototype > Copy() const
void SetParameterList(const Teuchos::ParameterList &paramList)
Set parameters from a parameter list and return with default values.
void DeclareInput(Level &currentLevel) const
Input.