MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
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>
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>());
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 
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 
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>
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>
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>
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.
Print external lib objects.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Class that encapsulates Belos smoothers.
Print more statistics.
Print additional debugging information.
void Setup(Level &currentLevel)
Set up the smoother.
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
size_t getNodeSmootherComplexity() const
For diagnostic purposes.
RCP< SmootherPrototype > Copy() const
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.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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.
bool IsSetup() const
Get the state of a smoother prototype.
#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.
Print class parameters (more parameters, more verbose)
Exception throws to report errors in the internal logical of the program.
#define TEUCHOS_ASSERT(assertion_test)
void SetParameterList(const Teuchos::ParameterList &paramList)
Set parameters from a parameter list and return with default values.
void DeclareInput(Level &currentLevel) const
Input.
std::string toString(const T &t)
std::string description() const
Return a simple one-line description of this object.