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_Exceptions.hpp"
58 
59 namespace MueLu {
60 
61  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
62  TrilinosSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TrilinosSmoother(const std::string& type, const Teuchos::ParameterList& paramListIn, const LO& overlap)
63  : type_(type), overlap_(overlap)
64  {
65  // The original idea behind all smoothers was to use prototype pattern. However, it does not fully work of the dependencies
66  // calculation. Particularly, we need to propagate DeclareInput to proper prototypes. Therefore, both TrilinosSmoother and
67  // DirectSolver do not follow the pattern exactly.
68  // The difference is that in order to propagate the calculation of dependencies, we need to call a DeclareInput on a
69  // constructed object (or objects, as we have two different code branches for Epetra and Tpetra). The only place where we
70  // could construct these objects is the constructor. Thus, we need to store RCPs, and both TrilinosSmoother and DirectSolver
71  // obtain a state: they contain RCP to smoother prototypes.
72  sEpetra_ = Teuchos::null;
73  sTpetra_ = Teuchos::null;
74 
75  TEUCHOS_TEST_FOR_EXCEPTION(overlap_ < 0, Exceptions::RuntimeError, "Overlap parameter is negative (" << overlap << ")");
76 
77  // paramList contains only parameters which are understood by Ifpack/Ifpack2
78  ParameterList paramList = paramListIn;
79 
80 
81  // We want TrilinosSmoother to be able to work with both Epetra and Tpetra objects, therefore we try to construct both
82  // Ifpack and Ifpack2 smoother prototypes. The construction really depends on configuration options.
83  triedEpetra_ = triedTpetra_ = false;
84 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_MUELU_IFPACK2)
85  try {
86  sTpetra_ = rcp(new Ifpack2Smoother(type_, paramList, overlap_));
87  TEUCHOS_TEST_FOR_EXCEPTION(sTpetra_.is_null(), Exceptions::RuntimeError, "Unable to construct Ifpack2 smoother");
88  } catch (Exceptions::RuntimeError& e) {
89  errorTpetra_ = e.what();
90  } catch (Exceptions::BadCast& e) {
91  errorTpetra_ = e.what();
92  }
93  triedTpetra_ = true;
94 #endif
95 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_IFPACK)
96  try {
97  // GetIfpackSmoother masks the template argument matching, and simply throws if template arguments are incompatible with Epetra
99  TEUCHOS_TEST_FOR_EXCEPTION(sEpetra_.is_null(), Exceptions::RuntimeError, "Unable to construct Ifpack smoother");
100  } catch (Exceptions::RuntimeError& e) {
101  // IfpackSmoother throws if Scalar != double, LocalOrdinal != int, GlobalOrdinal != int
102  errorEpetra_ = e.what();
103  }
104  triedEpetra_ = true;
105 #endif
106 
107  // Check if we were able to construct at least one smoother. In many cases that's all we need, for instance if a user
108  // simply wants to use Tpetra only stack, never enables Ifpack, and always runs Tpetra objects.
109  TEUCHOS_TEST_FOR_EXCEPTION(!triedEpetra_ && !triedTpetra_, Exceptions::RuntimeError, "Unable to construct any smoother."
110  "Please enable (TPETRA and IFPACK2) or (EPETRA and IFPACK)");
111 
112  TEUCHOS_TEST_FOR_EXCEPTION(sEpetra_.is_null() && sTpetra_.is_null(), Exceptions::RuntimeError,
113  "Could not construct any smoother:\n"
114  << (triedEpetra_ ? "=> Failed to build an Epetra smoother due to the following exception:\n" : "=> Epetra and/or Ifpack are not enabled.\n")
115  << (triedEpetra_ ? errorEpetra_ + "\n" : "")
116  << (triedTpetra_ ? "=> Failed to build a Tpetra smoother due to the following exception:\n" : "=> Tpetra and/or Ifpack2 are not enabled.\n")
117  << (triedTpetra_ ? errorTpetra_ + "\n" : ""));
118 
119  this->SetParameterList(paramList);
120  }
121 
122  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
123  void TrilinosSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::SetFactory(const std::string& varName, const RCP<const FactoryBase>& factory) {
124  // We need to propagate SetFactory to proper place
125  if (!sEpetra_.is_null()) sEpetra_->SetFactory(varName, factory);
126  if (!sTpetra_.is_null()) sTpetra_->SetFactory(varName, factory);
127  }
128 
129  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
131  // Decide whether we are running in Epetra or Tpetra mode
132  //
133  // Theoretically, we could make this decision in the constructor, and create only
134  // one of the smoothers. But we want to be able to reuse, so one can imagine a scenario
135  // where one first runs hierarchy with tpetra matrix, and then with epetra.
136  bool useTpetra = (currentLevel.lib() == Xpetra::UseTpetra);
137  s_ = (useTpetra ? sTpetra_ : sEpetra_);
138  if (s_.is_null()) {
139  if (useTpetra) {
140 #if not defined(HAVE_MUELU_IFPACK22)
141  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError,
142  "Error: running in Tpetra mode, but MueLu with Ifpack2 was disabled during the configure stage.\n"
143  "Please make sure that:\n"
144  " - Ifpack2 is enabled (Trilinos_ENABLE_Ifpack2=ON),\n"
145  " - Ifpack2 is available for MueLu to use (MueLu_ENABLE_Ifpack2=ON)\n");
146 #else
147  if (triedTpetra_)
148  this->GetOStream(Errors) << "Tpetra mode was disabled due to an error:\n" << errorTpetra_ << std::endl;
149 #endif
150  }
151  if (!useTpetra) {
152 #if not defined(HAVE_MUELU_IFPACK)
153  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError,
154  "Error: running in Epetra mode, but MueLu with Ifpack was disabled during the configure stage.\n"
155  "Please make sure that:\n"
156  " - Ifpack is enabled (you can do that with Trilinos_ENABLE_Ifpack=ON),\n"
157  " - Ifpack is available for MueLu to use (MueLu_ENABLE_Ifpack=ON)\n");
158 #else
159  if (triedEpetra_)
160  this->GetOStream(Errors) << "Epetra mode was disabled due to an error:\n" << errorEpetra_ << std::endl;
161 #endif
162  }
163  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError,
164  "Smoother for " << (useTpetra ? "Tpetra" : "Epetra") << " was not constructed");
165  }
166 
167  s_->DeclareInput(currentLevel);
168  }
169 
170  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
172  if (SmootherPrototype::IsSetup() == true)
173  this->GetOStream(Warnings0) << "MueLu::TrilinosSmoother::Setup(): Setup() has already been called" << std::endl;
174 
175  int oldRank = s_->SetProcRankVerbose(this->GetProcRankVerbose());
176 
177  s_->Setup(currentLevel);
178 
179  s_->SetProcRankVerbose(oldRank);
180 
182 
183  this->SetParameterList(s_->GetParameterList());
184  }
185 
186  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
187  void TrilinosSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Apply(MultiVector& X, const MultiVector& B, bool InitialGuessIsZero) const {
188  TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError, "MueLu::TrilinosSmoother::Apply(): Setup() has not been called");
189 
190  s_->Apply(X, B, InitialGuessIsZero);
191  }
192 
193  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
194  RCP<MueLu::SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
196  RCP<TrilinosSmoother> newSmoo = rcp(new TrilinosSmoother(type_, this->GetParameterList(), overlap_));
197 
198  // We need to be quite careful with Copy
199  // We still want TrilinosSmoother to follow Prototype Pattern, so we need to hide the fact that we do have some state
200  if (!sEpetra_.is_null())
201  newSmoo->sEpetra_ = sEpetra_->Copy();
202  if (!sTpetra_.is_null())
203  newSmoo->sTpetra_ = sTpetra_->Copy();
204 
205  // Copy the default mode
206  newSmoo->s_ = (s_.get() == sTpetra_.get() ? newSmoo->sTpetra_ : newSmoo->sEpetra_);
207  newSmoo->SetParameterList(this->GetParameterList());
208 
209  return newSmoo;
210  }
211 
212  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
214  if (type == "RELAXATION") { return "point relaxation stand-alone"; }
215  if (type == "CHEBYSHEV") { return "Chebyshev"; }
216  if (type == "ILUT") { return "ILUT"; }
217  if (type == "RILUK") { return "ILU"; }
218  if (type == "ILU") { return "ILU"; }
219  if (type == "Amesos") { return "Amesos"; }
220 
221  // special types
222 
223  // Note: for Ifpack there is no distinction between block and banded relaxation as there is no BandedContainer.
224  if (type == "LINESMOOTHING_BLOCKRELAXATION") { return "LINESMOOTHING_BLOCKRELAXATION"; }
225  if (type == "LINESMOOTHING_BLOCK RELAXATION") { return "LINESMOOTHING_BLOCKRELAXATION"; }
226  if (type == "LINESMOOTHING_BLOCK_RELAXATION") { return "LINESMOOTHING_BLOCKRELAXATION"; }
227  if (type == "LINESMOOTHING_BANDEDRELAXATION") { return "LINESMOOTHING_BLOCKRELAXATION"; }
228  if (type == "LINESMOOTHING_BANDED RELAXATION") { return "LINESMOOTHING_BLOCKRELAXATION"; }
229  if (type == "LINESMOOTHING_BANDED_RELAXATION") { return "LINESMOOTHING_BLOCKRELAXATION"; }
230 
231  TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "Cannot convert Ifpack2 preconditioner name to Ifpack: unknown type: \"" + type + "\"");
232  }
233 
234  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
235  Teuchos::ParameterList TrilinosSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Ifpack2ToIfpack1Param(const Teuchos::ParameterList& ifpack2List) {
236  Teuchos::ParameterList ifpack1List = ifpack2List;
237 
238  if (ifpack2List.isParameter("relaxation: type") && ifpack2List.get<std::string>("relaxation: type") == "Symmetric Gauss-Seidel")
239  ifpack1List.set("relaxation: type", "symmetric Gauss-Seidel");
240 
241  if (ifpack2List.isParameter("fact: iluk level-of-fill")) {
242  ifpack1List.remove("fact: iluk level-of-fill");
243  ifpack1List.set("fact: level-of-fill", ifpack2List.get<int>("fact: iluk level-of-fill"));
244  }
245 
246  return ifpack1List;
247  }
248 
249  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
251  std::ostringstream out;
252  if (s_ != Teuchos::null) {
253  out << s_->description();
254  } else {
256  out << "{type = " << type_ << "}";
257  }
258  return out.str();
259  }
260 
261  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
262  void TrilinosSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::print(Teuchos::FancyOStream& out, const VerbLevel verbLevel) const {
264 
265  if (verbLevel & Parameters0)
266  out0 << "Prec. type: " << type_ << std::endl;
267 
268  if (verbLevel & Parameters1) {
269  out0 << "PrecType: " << type_ << std::endl;
270  out0 << "Parameter list: " << std::endl;
271  Teuchos::OSTab tab2(out);
272  out << this->GetParameterList();
273  out0 << "Overlap: " << overlap_ << std::endl;
274  }
275 
276  if (verbLevel & Debug) {
277  out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl
278  << "-" << std::endl
279  << "Epetra PrecType: " << Ifpack2ToIfpack1Type(type_) << std::endl
280  << "Epetra Parameter list: " << std::endl;
281  Teuchos::OSTab tab2(out);
282  out << Ifpack2ToIfpack1Param(this->GetParameterList());;
283  }
284  }
285 
286 } // namespace MueLu
287 
288 #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.
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.
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 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.