MueLu Version of the Day
Loading...
Searching...
No Matches
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"
56#include "MueLu_Ifpack2Smoother.hpp"
57#include "MueLu_BelosSmoother.hpp"
58#include "MueLu_StratimikosSmoother.hpp"
59#include "MueLu_Exceptions.hpp"
60
61namespace 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
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
Class that encapsulates Belos smoothers.
virtual std::string description() const
Return a simple one-line description of this object.
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
Class that encapsulates Ifpack2 smoothers.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
Xpetra::UnderlyingLib lib()
virtual void SetParameterList(const Teuchos::ParameterList &paramList)
Set parameters from a parameter list and return with default values.
bool IsSetup() const
Get the state of a smoother prototype.
Class that encapsulates external library smoothers.
RCP< SmootherPrototype > Copy() const
When this prototype is cloned using Copy(), the clone is an Ifpack or an Ifpack2 smoother.
void Setup(Level &currentLevel)
TrilinosSmoother cannot be turned into a smoother using Setup(). Setup() always returns a RuntimeErro...
LO overlap_
overlap when using the smoother in additive Schwarz mode
static Teuchos::ParameterList Ifpack2ToIfpack1Param(const Teuchos::ParameterList &ifpack2List)
Convert an Ifpack2 parameter list to Ifpack.
static std::string Ifpack2ToIfpack1Type(const std::string &type)
Convert an Ifpack2 preconditioner name to Ifpack.
std::string type_
ifpack1/2-specific key phrase that denote smoother type
friend class TrilinosSmoother
Friend declaration required for clone() functionality.
void DeclareInput(Level &currentLevel) const
Input.
RCP< SmootherPrototype > sEpetra_
Smoother.
void SetFactory(const std::string &varName, const RCP< const FactoryBase > &factory)
Custom SetFactory.
RCP< SmootherPrototype > sStratimikos_
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
TrilinosSmoother cannot be applied. Apply() always returns a RuntimeError exception.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
std::string description() const
Return a simple one-line description of this object.
Namespace for MueLu classes and methods.
@ Warnings0
Important warning messages (one line)
@ Debug
Print additional debugging information.
@ Errors
Errors.
@ Parameters0
Print class parameters.
@ Parameters1
Print class parameters (more parameters, more verbose)