MueLu  Version of the Day
MueLu_DirectSolver_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_DIRECTSOLVER_DEF_HPP
47 #define MUELU_DIRECTSOLVER_DEF_HPP
48 
49 #include <Xpetra_Utils.hpp>
50 #include <Xpetra_Matrix.hpp>
51 
53 
54 #include "MueLu_Exceptions.hpp"
55 #include "MueLu_Level.hpp"
56 
57 #include "MueLu_Amesos2Smoother.hpp"
58 #include "MueLu_AmesosSmoother.hpp"
59 #include "MueLu_BelosSmoother.hpp"
60 #include "MueLu_StratimikosSmoother.hpp"
61 
62 namespace MueLu {
63 
64  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
65  DirectSolver<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DirectSolver(const std::string& type, const Teuchos::ParameterList& paramListIn)
66  : type_(type) {
67  // The original idea behind all smoothers was to use prototype pattern. However, it does not fully work of the dependencies
68  // calculation. Particularly, we need to propagate DeclareInput to proper prototypes. Therefore, both TrilinosSmoother and
69  // DirectSolver do not follow the pattern exactly.
70  // The difference is that in order to propagate the calculation of dependencies, we need to call a DeclareInput on a
71  // constructed object (or objects, as we have two different code branches for Epetra and Tpetra). The only place where we
72  // could construct these objects is the constructor. Thus, we need to store RCPs, and both TrilinosSmoother and DirectSolver
73  // obtain a state: they contain RCP to smoother prototypes.
74  sEpetra_ = Teuchos::null;
75  sTpetra_ = Teuchos::null;
76  sBelos_ = Teuchos::null;
77 
78  ParameterList paramList = paramListIn;
79 
80  // We want DirectSolver to be able to work with both Epetra and Tpetra objects, therefore we try to construct both
81  // Amesos and Amesos2 solver prototypes. The construction really depends on configuration options.
82  triedEpetra_ = triedTpetra_ = false;
83 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_MUELU_AMESOS2)
84  try {
85  sTpetra_ = rcp(new Amesos2Smoother(type_, paramList));
86  if (sTpetra_.is_null())
87  errorTpetra_ = "Unable to construct Amesos2 direct solver";
88  else if (!sTpetra_->constructionSuccessful()) {
89  errorTpetra_ = sTpetra_->constructionErrorMsg();
90  sTpetra_ = Teuchos::null;
91  }
92  } catch (Exceptions::RuntimeError& e) {
93  errorTpetra_ = e.what();
94  } catch (Exceptions::BadCast& e) {
95  errorTpetra_ = e.what();
96  } catch (Teuchos::Exceptions::InvalidParameterName& e) {
97  errorTpetra_ = e.what();
98  }
99  triedTpetra_ = true;
100 #endif
101 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_AMESOS)
102  try {
103  // GetAmesosSmoother masks the template argument matching, and simply throws if template arguments are incompatible with Epetra
104  sEpetra_ = GetAmesosSmoother<SC,LO,GO,NO>(type_, paramList);
105  if (sEpetra_.is_null())
106  errorEpetra_ = "Unable to construct Amesos direct solver";
107  else if (!sEpetra_->constructionSuccessful()) {
108  errorEpetra_ = sEpetra_->constructionErrorMsg();
109  sEpetra_ = Teuchos::null;
110  }
111  } catch (Exceptions::RuntimeError& e) {
112  // AmesosSmoother throws if Scalar != double, LocalOrdinal != int, GlobalOrdinal != int
113  errorEpetra_ = e.what();
114  }
115  triedEpetra_ = true;
116 #endif
117 #if defined(HAVE_MUELU_BELOS)
118  try {
119  sBelos_ = rcp(new BelosSmoother(type_, paramList));
120  if (sBelos_.is_null())
121  errorBelos_ = "Unable to construct Belos solver";
122  else if (!sBelos_->constructionSuccessful()) {
123  errorBelos_ = sBelos_->constructionErrorMsg();
124  sBelos_ = Teuchos::null;
125  }
126  } catch (Exceptions::RuntimeError& e) {
127  errorBelos_ = e.what();
128  } catch (Exceptions::BadCast& e) {
129  errorBelos_ = e.what();
130  }
131  triedBelos_ = true;
132 #endif
133 #if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_TPETRA) && defined(HAVE_MUELU_THYRA)
134  try {
135  sStratimikos_ = rcp(new StratimikosSmoother(type_, paramList));
136  if (sStratimikos_.is_null())
137  errorStratimikos_ = "Unable to construct Stratimikos smoother";
138  else if (!sStratimikos_->constructionSuccessful()) {
139  errorStratimikos_ = sStratimikos_->constructionErrorMsg();
140  sStratimikos_ = Teuchos::null;
141  }
142  } catch (Exceptions::RuntimeError& e){
143  errorStratimikos_ = e.what();
144  }
145  triedStratimikos_ = true;
146 #endif
147 
148  // Check if we were able to construct at least one solver. In many cases that's all we need, for instance if a user
149  // simply wants to use Tpetra only stack, never enables Amesos, and always runs Tpetra objects.
150  TEUCHOS_TEST_FOR_EXCEPTION(!triedEpetra_ && !triedTpetra_ && !triedBelos_ && !triedStratimikos_, Exceptions::RuntimeError, "Unable to construct any direct solver."
151  "Plase enable (TPETRA and AMESOS2) or (EPETRA and AMESOS) or (BELOS) or (STRATIMIKOS)");
152 
153  TEUCHOS_TEST_FOR_EXCEPTION(sEpetra_.is_null() && sTpetra_.is_null() && sBelos_.is_null() && sStratimikos_.is_null(), Exceptions::RuntimeError,
154  "Could not enable any direct solver:\n"
155  << (triedEpetra_ ? "Epetra mode was disabled due to an error:\n" : "")
156  << (triedEpetra_ ? errorEpetra_ : "")
157  << (triedTpetra_ ? "Tpetra mode was disabled due to an error:\n" : "")
158  << (triedTpetra_ ? errorTpetra_ : "")
159  << (triedBelos_ ? "Belos was disabled due to an error:\n" : "")
160  << (triedBelos_ ? errorBelos_ : "")
161  << (triedStratimikos_ ? "Stratimikos was disabled due to an error:\n" : "")
163 
164  this->SetParameterList(paramList);
165  }
166 
167  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
168  void DirectSolver<Scalar, LocalOrdinal, GlobalOrdinal, Node>::SetFactory(const std::string& varName, const RCP<const FactoryBase>& factory) {
169  // We need to propagate SetFactory to proper place
170  if (!sEpetra_.is_null()) sEpetra_->SetFactory(varName, factory);
171  if (!sTpetra_.is_null()) sTpetra_->SetFactory(varName, factory);
172  if (!sBelos_.is_null()) sBelos_->SetFactory(varName, factory);
173  if (!sStratimikos_.is_null()) sStratimikos_->SetFactory(varName, factory);
174  }
175 
176  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
178  if (!sBelos_.is_null())
179  s_ = sBelos_;
180  else if (!sStratimikos_.is_null())
181  s_ = sStratimikos_;
182  else {
183  // Decide whether we are running in Epetra or Tpetra mode
184  //
185  // Theoretically, we could make this decision in the constructor, and create only
186  // one of the smoothers. But we want to be able to reuse, so one can imagine a scenario
187  // where one first runs hierarchy with tpetra matrix, and then with epetra.
188  bool useTpetra = (currentLevel.lib() == Xpetra::UseTpetra);
189  s_ = (useTpetra ? sTpetra_ : sEpetra_);
190  if (s_.is_null()) {
191  if (useTpetra) {
192 #if not defined(HAVE_MUELU_AMESOS2)
193  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError,
194  "Error: running in Tpetra mode, but MueLu with Amesos2 was disabled during the configure stage.\n"
195  "Please make sure that:\n"
196  " - Amesos2 is enabled (Trilinos_ENABLE_Amesos2=ON),\n"
197  " - Amesos2 is available for MueLu to use (MueLu_ENABLE_Amesos2=ON)\n");
198 #else
199  if (triedTpetra_)
200  this->GetOStream(Errors) << "Tpetra mode was disabled due to an error:\n" << errorTpetra_ << std::endl;
201 #endif
202  }
203  if (!useTpetra) {
204 #if not defined(HAVE_MUELU_AMESOS)
205  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError,
206  "Error: running in Epetra mode, but MueLu with Amesos was disabled during the configure stage.\n"
207  "Please make sure that:\n"
208  " - Amesos is enabled (you can do that with Trilinos_ENABLE_Amesos=ON),\n"
209  " - Amesos is available for MueLu to use (MueLu_ENABLE_Amesos=ON)\n");
210 #else
211  if (triedEpetra_)
212  this->GetOStream(Errors) << "Epetra mode was disabled due to an error:\n" << errorEpetra_ << std::endl;
213 #endif
214  }
215  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError,
216  "Direct solver for " << (useTpetra ? "Tpetra" : "Epetra") << " was not constructed");
217  }
218  }
219 
220  s_->DeclareInput(currentLevel);
221  }
222 
223  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
225  if (SmootherPrototype::IsSetup() == true)
226  this->GetOStream(Warnings0) << "MueLu::DirectSolver::Setup(): Setup() has already been called" << std::endl;
227 
228  int oldRank = s_->SetProcRankVerbose(this->GetProcRankVerbose());
229 
230  s_->Setup(currentLevel);
231 
232  s_->SetProcRankVerbose(oldRank);
233 
235 
236  this->SetParameterList(s_->GetParameterList());
237  }
238 
239  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
240  void DirectSolver<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Apply(MultiVector& X, const MultiVector& B, bool InitialGuessIsZero) const {
241  TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError, "MueLu::AmesosSmoother::Apply(): Setup() has not been called");
242 
243  s_->Apply(X, B, InitialGuessIsZero);
244  }
245 
246  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
247  RCP<MueLu::SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node> > DirectSolver<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Copy() const {
248  RCP<DirectSolver> newSmoo = rcp(new DirectSolver(*this));
249 
250  // We need to be quite careful with Copy
251  // We still want DirectSolver to follow Prototype Pattern, so we need to hide the fact that we do have some state
252  if (!sEpetra_.is_null())
253  newSmoo->sEpetra_ = sEpetra_->Copy();
254  if (!sTpetra_.is_null())
255  newSmoo->sTpetra_ = sTpetra_->Copy();
256  if (!sBelos_.is_null())
257  newSmoo->sBelos_ = sBelos_->Copy();
258  if (!sStratimikos_.is_null())
259  newSmoo->sStratimikos_ = sStratimikos_->Copy();
260 
261  // Copy the default mode
262  if (s_.get() == sBelos_.get())
263  newSmoo->s_ = newSmoo->sBelos_;
264  else if (s_.get() == sStratimikos_.get())
265  newSmoo->s_ = newSmoo->sStratimikos_;
266  else if (s_.get() == sTpetra_.get())
267  newSmoo->s_ = newSmoo->sTpetra_;
268  else
269  newSmoo->s_ = newSmoo->sEpetra_;
270  newSmoo->SetParameterList(this->GetParameterList());
271 
272  return newSmoo;
273  }
274 
275  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
277  std::ostringstream out;
278  if (s_ != Teuchos::null) {
279  out << s_->description();
280  } else {
282  out << "{type = " << type_ << "}";
283  }
284  return out.str();
285  }
286 
287  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
288  void DirectSolver<Scalar, LocalOrdinal, GlobalOrdinal, Node>::print(Teuchos::FancyOStream& out, const VerbLevel verbLevel) const {
290 
291  if (verbLevel & Parameters0)
292  out0 << "Prec. type: " << type_ << std::endl;
293 
294  if (verbLevel & Parameters1) {
295  out0 << "Parameter list: " << std::endl;
296  Teuchos::OSTab tab3(out);
297  out << this->GetParameterList();
298  }
299 
300  if (verbLevel & Debug)
301  out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl;
302  }
303 
304 } // namespace MueLu
305 
306 #endif // MUELU_DIRECTSOLVER_DEF_HPP
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Important warning messages (one line)
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
DirectSolver cannot be applied. Apply() always returns a RuntimeError exception.
Exception indicating invalid cast attempted.
void Setup(Level &currentLevel)
DirectSolver cannot be turned into a smoother using Setup(). Setup() always returns a RuntimeError ex...
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
bool IsSetup() const
Get the state of a smoother prototype.
Class that encapsulates Belos smoothers.
Print additional debugging information.
Namespace for MueLu classes and methods.
std::string type_
amesos1/2-specific key phrase that denote smoother type
DirectSolver(const std::string &type="", const Teuchos::ParameterList &paramList=Teuchos::ParameterList())
Constructor Note: only parameters shared by Amesos and Amesos2 should be used for type and paramList ...
Class that encapsulates direct solvers. Autoselection of AmesosSmoother or Amesos2Smoother according ...
virtual void SetParameterList(const Teuchos::ParameterList &paramList)
Set parameters from a parameter list and return with default values.
std::string description() const
Return a simple one-line description of this object.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
Class that encapsulates Amesos2 direct solvers.
RCP< SmootherPrototype > sBelos_
void DeclareInput(Level &currentLevel) const
Input.
Xpetra::UnderlyingLib lib()
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
RCP< SmootherPrototype > sTpetra_
Print class parameters.
void SetFactory(const std::string &varName, const RCP< const FactoryBase > &factory)
Custom SetFactory.
RCP< SmootherPrototype > sEpetra_
Smoother.
Print class parameters (more parameters, more verbose)
Exception throws to report errors in the internal logical of the program.
RCP< SmootherPrototype > sStratimikos_
RCP< SmootherPrototype > Copy() const
When this prototype is cloned using Copy(), the clone is an Amesos or an Amesos2 smoother.
virtual std::string description() const
Return a simple one-line description of this object.