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