MueLu  Version of the Day
MueLu_TrilinosSmoother_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_TRILINOSSMOOTHER_DEF_HPP
47 #define MUELU_TRILINOSSMOOTHER_DEF_HPP
48 
49 #include <Xpetra_Map.hpp>
50 #include <Xpetra_Matrix.hpp>
51 
53 
54 #include "MueLu_Level.hpp"
55 #include "MueLu_IfpackSmoother.hpp"
56 #include "MueLu_Ifpack2Smoother.hpp"
57 #include "MueLu_BelosSmoother.hpp"
58 #include "MueLu_StratimikosSmoother.hpp"
59 #include "MueLu_Exceptions.hpp"
60 
61 namespace MueLu {
62 
63  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
64  TrilinosSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TrilinosSmoother(const std::string& type, const Teuchos::ParameterList& paramListIn, const LO& overlap)
65  : type_(type), overlap_(overlap)
66  {
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  sStratimikos_ = Teuchos::null;
78 
79  TEUCHOS_TEST_FOR_EXCEPTION(overlap_ < 0, Exceptions::RuntimeError, "Overlap parameter is negative (" << overlap << ")");
80 
81  // paramList contains only parameters which are understood by Ifpack/Ifpack2
82  ParameterList paramList = paramListIn;
83 
84 
85  // We want TrilinosSmoother to be able to work with both Epetra and Tpetra objects, therefore we try to construct both
86  // Ifpack and Ifpack2 smoother prototypes. The construction really depends on configuration options.
88 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_MUELU_IFPACK2)
89  try {
90  sTpetra_ = rcp(new Ifpack2Smoother(type_, paramList, overlap_));
91  if (sTpetra_.is_null())
92  errorTpetra_ = "Unable to construct Ifpack2 smoother";
93  else if (!sTpetra_->constructionSuccessful()) {
94  errorTpetra_ = sTpetra_->constructionErrorMsg();
95  sTpetra_ = Teuchos::null;
96  }
97  } catch (Exceptions::RuntimeError& e) {
98  errorTpetra_ = e.what();
99  } catch (Exceptions::BadCast& e) {
100  errorTpetra_ = e.what();
101  } catch (Teuchos::Exceptions::InvalidParameterName& e) {
102  errorTpetra_ = e.what();
103  }
104  triedTpetra_ = true;
105 #endif
106 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_IFPACK)
107  try {
108  // GetIfpackSmoother masks the template argument matching, and simply throws if template arguments are incompatible with Epetra
110  if (sEpetra_.is_null())
111  errorEpetra_ = "Unable to construct Ifpack smoother";
112  else if (!sEpetra_->constructionSuccessful()) {
113  errorEpetra_ = sEpetra_->constructionErrorMsg();
114  sEpetra_ = Teuchos::null;
115  }
116  } catch (Exceptions::RuntimeError& e) {
117  // IfpackSmoother throws if Scalar != double, LocalOrdinal != int, GlobalOrdinal != int
118  errorEpetra_ = e.what();
119  }
120  triedEpetra_ = true;
121 #endif
122 #if defined(HAVE_MUELU_BELOS)
123  try {
124  sBelos_ = rcp(new BelosSmoother(type_, paramList));
125  if (sBelos_.is_null())
126  errorBelos_ = "Unable to construct Belos smoother";
127  else if (!sBelos_->constructionSuccessful()) {
128  errorBelos_ = sBelos_->constructionErrorMsg();
129  sBelos_ = Teuchos::null;
130  }
131  } catch (Exceptions::RuntimeError& e){
132  errorBelos_ = e.what();
133  } catch (Exceptions::BadCast& e) {
134  errorBelos_ = e.what();
135  }
136  triedBelos_ = true;
137 #endif
138 #if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_TPETRA) && defined(HAVE_MUELU_THYRA)
139  try {
140  sStratimikos_ = rcp(new StratimikosSmoother(type_, paramList));
141  if (sStratimikos_.is_null())
142  errorStratimikos_ = "Unable to construct Stratimikos smoother";
143  else if (!sStratimikos_->constructionSuccessful()) {
144  errorStratimikos_ = sStratimikos_->constructionErrorMsg();
145  sStratimikos_ = Teuchos::null;
146  }
147  } catch (Exceptions::RuntimeError& e){
148  errorStratimikos_ = e.what();
149  }
150  triedStratimikos_ = true;
151 #endif
152 
153  // Check if we were able to construct at least one smoother. In many cases that's all we need, for instance if a user
154  // simply wants to use Tpetra only stack, never enables Ifpack, and always runs Tpetra objects.
155  TEUCHOS_TEST_FOR_EXCEPTION(!triedEpetra_ && !triedTpetra_ && !triedBelos_ && !triedStratimikos_, Exceptions::RuntimeError, "Unable to construct any smoother."
156  "Please enable (TPETRA and IFPACK2) or (EPETRA and IFPACK) or (BELOS) or (STRATIMIKOS)");
157 
158  TEUCHOS_TEST_FOR_EXCEPTION(sEpetra_.is_null() && sTpetra_.is_null() && sBelos_.is_null() && sStratimikos_.is_null(), Exceptions::RuntimeError,
159  "Could not construct any smoother:\n"
160  << (triedEpetra_ ? "=> Failed to build an Epetra smoother due to the following exception:\n" : "=> Epetra and/or Ifpack are not enabled.\n")
161  << (triedEpetra_ ? errorEpetra_ + "\n" : "")
162  << (triedTpetra_ ? "=> Failed to build a Tpetra smoother due to the following exception:\n" : "=> Tpetra and/or Ifpack2 are not enabled.\n")
163  << (triedTpetra_ ? errorTpetra_ + "\n" : "")
164  << (triedBelos_ ? "=> Failed to build a Belos smoother due to the following exception:\n" : "=> Belos not enabled.\n")
165  << (triedBelos_ ? errorBelos_ + "\n" : "")
166  << (triedStratimikos_ ? "=> Failed to build a Stratimikos smoother due to the following exception:\n" : "=> Stratimikos not enabled.\n")
167  << (triedStratimikos_ ? errorStratimikos_ + "\n" : ""));
168 
169  this->SetParameterList(paramList);
170  }
171 
172  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
173  void TrilinosSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::SetFactory(const std::string& varName, const RCP<const FactoryBase>& factory) {
174  // We need to propagate SetFactory to proper place
175  if (!sEpetra_.is_null()) sEpetra_->SetFactory(varName, factory);
176  if (!sTpetra_.is_null()) sTpetra_->SetFactory(varName, factory);
177  if (!sBelos_.is_null()) sBelos_->SetFactory(varName, factory);
178  if (!sStratimikos_.is_null()) sStratimikos_->SetFactory(varName, factory);
179  }
180 
181  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
183  if (!sBelos_.is_null())
184  s_ = sBelos_;
185  else if (!sStratimikos_.is_null())
186  s_ = sStratimikos_;
187  else {
188  // Decide whether we are running in Epetra or Tpetra mode
189  //
190  // Theoretically, we could make this decision in the constructor, and create only
191  // one of the smoothers. But we want to be able to reuse, so one can imagine a scenario
192  // where one first runs hierarchy with tpetra matrix, and then with epetra.
193 
194  bool useTpetra = (currentLevel.lib() == Xpetra::UseTpetra);
195  s_ = (useTpetra ? sTpetra_ : sEpetra_);
196  if (s_.is_null()) {
197  if (useTpetra) {
198 #if not defined(HAVE_MUELU_IFPACK2)
199  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError,
200  "Error: running in Tpetra mode, but MueLu with Ifpack2 was disabled during the configure stage.\n"
201  "Please make sure that:\n"
202  " - Ifpack2 is enabled (Trilinos_ENABLE_Ifpack2=ON),\n"
203  " - Ifpack2 is available for MueLu to use (MueLu_ENABLE_Ifpack2=ON)\n");
204 #else
205  if (triedTpetra_)
206  this->GetOStream(Errors) << "Tpetra mode was disabled due to an error:\n" << errorTpetra_ << std::endl;
207 #endif
208  }
209  if (!useTpetra) {
210 #if not defined(HAVE_MUELU_IFPACK)
211  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError,
212  "Error: running in Epetra mode, but MueLu with Ifpack was disabled during the configure stage.\n"
213  "Please make sure that:\n"
214  " - Ifpack is enabled (you can do that with Trilinos_ENABLE_Ifpack=ON),\n"
215  " - Ifpack is available for MueLu to use (MueLu_ENABLE_Ifpack=ON)\n");
216 #else
217  if (triedEpetra_)
218  this->GetOStream(Errors) << "Epetra mode was disabled due to an error:\n" << errorEpetra_ << std::endl;
219 #endif
220  }
221  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError,
222  "Smoother for " << (useTpetra ? "Tpetra" : "Epetra") << " was not constructed");
223  }
224  }
225 
226  s_->DeclareInput(currentLevel);
227  }
228 
229  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
231  if (SmootherPrototype::IsSetup() == true)
232  this->GetOStream(Warnings0) << "MueLu::TrilinosSmoother::Setup(): Setup() has already been called" << std::endl;
233 
234  int oldRank = s_->SetProcRankVerbose(this->GetProcRankVerbose());
235 
236  s_->Setup(currentLevel);
237 
238  s_->SetProcRankVerbose(oldRank);
239 
241 
242  this->SetParameterList(s_->GetParameterList());
243  }
244 
245  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
246  void TrilinosSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Apply(MultiVector& X, const MultiVector& B, bool InitialGuessIsZero) const {
247  TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError, "MueLu::TrilinosSmoother::Apply(): Setup() has not been called");
248 
249  s_->Apply(X, B, InitialGuessIsZero);
250  }
251 
252  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
253  RCP<MueLu::SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
255  RCP<TrilinosSmoother> newSmoo = rcp(new TrilinosSmoother(type_, this->GetParameterList(), overlap_));
256 
257  // We need to be quite careful with Copy
258  // We still want TrilinosSmoother to follow Prototype Pattern, so we need to hide the fact that we do have some state
259  if (!sEpetra_.is_null())
260  newSmoo->sEpetra_ = sEpetra_->Copy();
261  if (!sTpetra_.is_null())
262  newSmoo->sTpetra_ = sTpetra_->Copy();
263  if (!sBelos_.is_null())
264  newSmoo->sBelos_ = sBelos_->Copy();
265  if (!sStratimikos_.is_null())
266  newSmoo->sStratimikos_ = sStratimikos_->Copy();
267 
268  // Copy the default mode
269  if (s_.get() == sBelos_.get())
270  newSmoo->s_ = newSmoo->sBelos_;
271  else if (s_.get() == sStratimikos_.get())
272  newSmoo->s_ = newSmoo->sStratimikos_;
273  else if (s_.get() == sTpetra_.get())
274  newSmoo->s_ = newSmoo->sTpetra_;
275  else
276  newSmoo->s_ = newSmoo->sEpetra_;
277  newSmoo->SetParameterList(this->GetParameterList());
278 
279  return newSmoo;
280  }
281 
282  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
284  if (type == "RELAXATION") { return "point relaxation stand-alone"; }
285  if (type == "CHEBYSHEV") { return "Chebyshev"; }
286  if (type == "ILUT") { return "ILUT"; }
287  if (type == "RILUK") { return "ILU"; }
288  if (type == "ILU") { return "ILU"; }
289  if (type == "Amesos") { return "Amesos"; }
290 
291  // special types
292 
293  // Note: for Ifpack there is no distinction between block and banded relaxation as there is no
294  // BandedContainer or TridiagonalContainer.
295  if (type == "LINESMOOTHING_BLOCKRELAXATION") { return "LINESMOOTHING_BLOCKRELAXATION"; }
296  if (type == "LINESMOOTHING_BLOCK RELAXATION") { return "LINESMOOTHING_BLOCKRELAXATION"; }
297  if (type == "LINESMOOTHING_BLOCK_RELAXATION") { return "LINESMOOTHING_BLOCKRELAXATION"; }
298  if (type == "LINESMOOTHING_BANDEDRELAXATION") { return "LINESMOOTHING_BLOCKRELAXATION"; }
299  if (type == "LINESMOOTHING_BANDED RELAXATION") { return "LINESMOOTHING_BLOCKRELAXATION"; }
300  if (type == "LINESMOOTHING_BANDED_RELAXATION") { return "LINESMOOTHING_BLOCKRELAXATION"; }
301  if (type == "LINESMOOTHING_TRIDIRELAXATION") { return "LINESMOOTHING_BLOCKRELAXATION"; }
302  if (type == "LINESMOOTHING_TRIDI RELAXATION") { return "LINESMOOTHING_BLOCKRELAXATION"; }
303  if (type == "LINESMOOTHING_TRIDI_RELAXATION") { return "LINESMOOTHING_BLOCKRELAXATION"; }
304  if (type == "LINESMOOTHING_TRIDIAGONALRELAXATION") { return "LINESMOOTHING_BLOCKRELAXATION"; }
305  if (type == "LINESMOOTHING_TRIDIAGONAL RELAXATION") { return "LINESMOOTHING_BLOCKRELAXATION"; }
306  if (type == "LINESMOOTHING_TRIDIAGONAL_RELAXATION") { return "LINESMOOTHING_BLOCKRELAXATION"; }
307  if (type == "AGGREGATE") { return "AGGREGATE"; }
308  if(type == "BLOCK_RELAXATION" ||
309  type == "BLOCK RELAXATION" ||
310  type == "BLOCKRELAXATION" ||
311  // Banded
312  type == "BANDED_RELAXATION" ||
313  type == "BANDED RELAXATION" ||
314  type == "BANDEDRELAXATION" ||
315  // Tridiagonal
316  type == "TRIDI_RELAXATION" ||
317  type == "TRIDI RELAXATION" ||
318  type == "TRIDIRELAXATION" ||
319  type == "TRIDIAGONAL_RELAXATION" ||
320  type == "TRIDIAGONAL RELAXATION" ||
321  type == "TRIDIAGONALRELAXATION") {
322  return "block relaxation stand-alone";
323  }
324 
325  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "Cannot convert Ifpack2 preconditioner name to Ifpack: unknown type: \"" + type + "\"");
326  }
327 
328  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
329  Teuchos::ParameterList TrilinosSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Ifpack2ToIfpack1Param(const Teuchos::ParameterList& ifpack2List) {
330  Teuchos::ParameterList ifpack1List = ifpack2List;
331 
332  if (ifpack2List.isParameter("relaxation: type") && ifpack2List.get<std::string>("relaxation: type") == "Symmetric Gauss-Seidel")
333  ifpack1List.set("relaxation: type", "symmetric Gauss-Seidel");
334 
335  if (ifpack2List.isParameter("fact: iluk level-of-fill")) {
336  ifpack1List.remove("fact: iluk level-of-fill");
337  ifpack1List.set("fact: level-of-fill", ifpack2List.get<int>("fact: iluk level-of-fill"));
338  }
339 
340  return ifpack1List;
341  }
342 
343  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
345  std::ostringstream out;
346  if (s_ != Teuchos::null) {
347  out << s_->description();
348  } else {
350  out << "{type = " << type_ << "}";
351  }
352  return out.str();
353  }
354 
355  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
356  void TrilinosSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::print(Teuchos::FancyOStream& out, const VerbLevel verbLevel) const {
358 
359  if (verbLevel & Parameters0)
360  out0 << "Prec. type: " << type_ << std::endl;
361 
362  if (verbLevel & Parameters1) {
363  out0 << "PrecType: " << type_ << std::endl;
364  out0 << "Parameter list: " << std::endl;
365  Teuchos::OSTab tab2(out);
366  out << this->GetParameterList();
367  out0 << "Overlap: " << overlap_ << std::endl;
368  }
369 
370  if (verbLevel & Debug) {
371  out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl
372  << "-" << std::endl
373  << "Epetra PrecType: " << Ifpack2ToIfpack1Type(type_) << std::endl
374  << "Epetra Parameter list: " << std::endl;
375  Teuchos::OSTab tab2(out);
376  out << Ifpack2ToIfpack1Param(this->GetParameterList());;
377  }
378  }
379 
380 } // namespace MueLu
381 
382 #endif // MUELU_TRILINOSSMOOTHER_DEF_HPP
Important warning messages (one line)
Exception indicating invalid cast attempted.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
RCP< SmootherPrototype > Copy() const
When this prototype is cloned using Copy(), the clone is an Ifpack or an Ifpack2 smoother.
RCP< SmootherPrototype > sStratimikos_
bool IsSetup() const
Get the state of a smoother prototype.
Class that encapsulates external library smoothers.
static std::string Ifpack2ToIfpack1Type(const std::string &type)
Convert an Ifpack2 preconditioner name to Ifpack.
static Teuchos::ParameterList Ifpack2ToIfpack1Param(const Teuchos::ParameterList &ifpack2List)
Convert an Ifpack2 parameter list to Ifpack.
Class that encapsulates Belos smoothers.
Print additional debugging information.
void DeclareInput(Level &currentLevel) const
Input.
Namespace for MueLu classes and methods.
void Setup(Level &currentLevel)
TrilinosSmoother cannot be turned into a smoother using Setup(). Setup() always returns a RuntimeErro...
friend class TrilinosSmoother
Friend declaration required for clone() functionality.
std::string description() const
Return a simple one-line description of this object.
virtual void SetParameterList(const Teuchos::ParameterList &paramList)
Set parameters from a parameter list and return with default values.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
std::string type_
ifpack1/2-specific key phrase that denote smoother type
Xpetra::UnderlyingLib lib()
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
Print class parameters.
RCP< SmootherPrototype > sEpetra_
Smoother.
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.
LO overlap_
overlap when using the smoother in additive Schwarz mode
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
TrilinosSmoother cannot be applied. Apply() always returns a RuntimeError exception.
Class that encapsulates Ifpack2 smoothers.
virtual std::string description() const
Return a simple one-line description of this object.
void SetFactory(const std::string &varName, const RCP< const FactoryBase > &factory)
Custom SetFactory.