MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_Details_LinearSolverFactory_def.hpp
Go to the documentation of this file.
1// @HEADER
2//
3// ***********************************************************************
4//
5// MueLu: A package for multigrid based preconditioning
6// Copyright 2011 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 Michael A. Heroux (maherou@sandia.gov)
39//
40// ***********************************************************************
41//
42// @HEADER
43
47
48#ifndef MUELU_DETAILS_LINEARSOLVERFACTORY_DEF_HPP
49#define MUELU_DETAILS_LINEARSOLVERFACTORY_DEF_HPP
50
51#include "MueLu_config.hpp"
52#include "Trilinos_Details_LinearSolver.hpp"
53#include "Trilinos_Details_LinearSolverFactory.hpp"
54#include <type_traits>
55
56#ifdef HAVE_MUELU_EPETRA
57# include "Epetra_CrsMatrix.h"
59#endif // HAVE_MUELU_EPETRA
60
61// Tpetra is not a required dependency of MueLu.
62#ifdef HAVE_MUELU_TPETRA
63# include "Tpetra_Operator.hpp"
65#endif // HAVE_MUELU_TPETRA
66
67namespace MueLu {
68namespace Details {
69
70template<class MV, class OP, class NormType>
72 public Trilinos::Details::LinearSolver<MV, OP, NormType>,
73 virtual public Teuchos::Describable
74{
75
76public:
77
80
82 virtual ~LinearSolver () {}
83
88 void setMatrix (const Teuchos::RCP<const OP>& A);
89
91 Teuchos::RCP<const OP> getMatrix () const {
92 return A_;
93 }
94
96 void solve (MV& X, const MV& B);
97
99 void setParameters (const Teuchos::RCP<Teuchos::ParameterList>& params);
100
103 void symbolic () {}
104
107 void numeric ();
108
110 std::string description () const;
111
113 void
114 describe (Teuchos::FancyOStream& out,
115 const Teuchos::EVerbosityLevel verbLevel =
116 Teuchos::Describable::verbLevel_default) const;
117
118private:
119 Teuchos::RCP<const OP> A_;
120 Teuchos::RCP<Teuchos::ParameterList> params_;
121};
122
123// Why does MueLu_EpetraOperator insist on HAVE_MUELU_SERIAL?
124#if defined(HAVE_MUELU_SERIAL) and defined(HAVE_MUELU_EPETRA)
125template<>
127 public Trilinos::Details::LinearSolver<Epetra_MultiVector, Epetra_Operator, double>,
128 virtual public Teuchos::Describable
129{
130
131public:
132
134 LinearSolver () :
135 changedA_(false),
136 changedParams_(false)
137 {}
138
140 virtual ~LinearSolver () {}
141
146 void setMatrix (const Teuchos::RCP<const Epetra_Operator>& A)
147 {
148 const char prefix[] = "MueLu::Details::LinearSolver::setMatrix: ";
149
150 if(A != A_)
151 {
152 if(solver_ != Teuchos::null)
153 changedA_ = true;
154
155 A_ = rcp_dynamic_cast<const Epetra_CrsMatrix>(A);
156 TEUCHOS_TEST_FOR_EXCEPTION
157 (A_.is_null(), std::runtime_error, prefix << "MueLu requires "
158 "an Epetra_CrsMatrix, but the matrix you provided is of a "
159 "different type. Please provide an Epetra_CrsMatrix instead.");
160 }
161 }
162
164 Teuchos::RCP<const Epetra_Operator> getMatrix () const {
165 return A_;
166 }
167
170 {
171 // TODO amk: Do we assume the user has called numeric before solve, or should we call it for them?
172 const char prefix[] = "MueLu::Details::LinearSolver::solve: ";
173 TEUCHOS_TEST_FOR_EXCEPTION
174 (solver_.is_null (), std::runtime_error, prefix << "The solver does not "
175 "exist yet. You must call numeric() before you may call this method.");
176 TEUCHOS_TEST_FOR_EXCEPTION
177 (changedA_, std::runtime_error, prefix << "The matrix A has been reset "
178 "since the last call to numeric(). Please call numeric() again.");
179 TEUCHOS_TEST_FOR_EXCEPTION
180 (changedParams_, std::runtime_error, prefix << "The parameters have been reset "
181 "since the last call to numeric(). Please call numeric() again.");
182
183 int err = solver_->ApplyInverse(B, X);
184
185 TEUCHOS_TEST_FOR_EXCEPTION
186 (err != 0, std::runtime_error, prefix << "EpetraOperator::ApplyInverse returned "
187 "nonzero error code " << err);
188 }
189
191 void setParameters (const Teuchos::RCP<Teuchos::ParameterList>& params)
192 {
193 if(solver_ != Teuchos::null && params != params_)
194 changedParams_ = true;
195
196 params_ = params;
197 }
198
201 void symbolic () {}
202
205 void numeric ()
206 {
207 const char prefix[] = "MueLu::Details::LinearSolver::numeric: ";
208
209 // If the solver is up-to-date, leave it alone
210 if(solver_ == Teuchos::null || changedA_ || changedParams_)
211 {
212 changedA_ = false;
213 changedParams_ = false;
214
215 TEUCHOS_TEST_FOR_EXCEPTION
216 (A_ == Teuchos::null, std::runtime_error, prefix << "The matrix has not been "
217 "set yet. You must call setMatrix() with a nonnull matrix before you may "
218 "call this method.");
219
220 // TODO: We should not have to cast away the constness here
221 // TODO: See bug 6462
222 if(params_ != Teuchos::null)
223 solver_ = CreateEpetraPreconditioner(rcp_const_cast<Epetra_CrsMatrix>(A_), *params_);
224 else
225 solver_ = CreateEpetraPreconditioner(rcp_const_cast<Epetra_CrsMatrix>(A_));
226 }
227 }
228
230 std::string description () const
231 {
232 if (solver_.is_null()) {
233 return "\"MueLu::Details::LinearSolver\": {MV: Epetra_MultiVector, OP: Epetra_Operator, NormType: double}";
234 }
235 else {
236 return solver_->GetHierarchy()->description ();
237 }
238 }
239
241 void
242 describe (Teuchos::FancyOStream& out,
243 const Teuchos::EVerbosityLevel verbLevel =
244 Teuchos::Describable::verbLevel_default) const
245 {
246 using std::endl;
247 if (solver_.is_null()) {
248 if(verbLevel > Teuchos::VERB_NONE) {
249 Teuchos::OSTab tab0 (out);
250 out << "\"MueLu::Details::LinearSolver\":" << endl;
251 Teuchos::OSTab tab1 (out);
252 out << "MV: Epetra_MultiVector" << endl
253 << "OP: Epetra_Operator" << endl
254 << "NormType: double" << endl;
255 }
256 }
257 else {
258 solver_->GetHierarchy()->describe (out, verbLevel);
259 }
260 }
261
262private:
263 Teuchos::RCP<const Epetra_CrsMatrix> A_;
264 Teuchos::RCP<Teuchos::ParameterList> params_;
265 Teuchos::RCP<EpetraOperator> solver_;
266 bool changedA_;
267 bool changedParams_;
268};
269#endif // HAVE_MUELU_EPETRA
270
271#ifdef HAVE_MUELU_TPETRA
272template<class Scalar, class LO, class GO, class Node>
273class LinearSolver<Tpetra::MultiVector<Scalar,LO,GO,Node>,
274 Tpetra::Operator<Scalar,LO,GO,Node>,
275 typename Teuchos::ScalarTraits<Scalar>::magnitudeType> :
276 public Trilinos::Details::LinearSolver<Tpetra::MultiVector<Scalar,LO,GO,Node>,
277 Tpetra::Operator<Scalar,LO,GO,Node>,
278 typename Teuchos::ScalarTraits<Scalar>::magnitudeType>,
279 virtual public Teuchos::Describable
280{
281
282public:
283
286 changedA_(false),
287 changedParams_(false)
288 {}
289
291 virtual ~LinearSolver () {}
292
297 void setMatrix (const Teuchos::RCP<const Tpetra::Operator<Scalar,LO,GO,Node> >& A)
298 {
299 if(A != A_)
300 {
301 if(solver_ != Teuchos::null)
302 changedA_ = true;
303
304 A_ = A;
305 }
306 }
307
309 Teuchos::RCP<const Tpetra::Operator<Scalar,LO,GO,Node> > getMatrix () const {
310 return A_;
311 }
312
314 void solve (Tpetra::MultiVector<Scalar,LO,GO,Node>& X, const Tpetra::MultiVector<Scalar,LO,GO,Node>& B)
315 {
316 // TODO amk: Do we assume the user has called numeric before solve, or should we call it for them?
317 const char prefix[] = "MueLu::Details::LinearSolver::solve: ";
318 TEUCHOS_TEST_FOR_EXCEPTION
319 (solver_.is_null (), std::runtime_error, prefix << "The solver does not "
320 "exist yet. You must call numeric() before you may call this method.");
321 TEUCHOS_TEST_FOR_EXCEPTION
322 (changedA_, std::runtime_error, prefix << "The matrix A has been reset "
323 "since the last call to numeric(). Please call numeric() again.");
324 TEUCHOS_TEST_FOR_EXCEPTION
325 (changedParams_, std::runtime_error, prefix << "The parameters have been reset "
326 "since the last call to numeric(). Please call numeric() again.");
327
328 solver_->apply(B, X);
329 }
330
332 void setParameters (const Teuchos::RCP<Teuchos::ParameterList>& params)
333 {
334 if(solver_ != Teuchos::null && params != params_)
335 changedParams_ = true;
336
337 params_ = params;
338 }
339
342 void symbolic () {}
343
346 void numeric ()
347 {
348 const char prefix[] = "MueLu::Details::LinearSolver::numeric: ";
349
350 // If the solver is up-to-date, leave it alone
351 if(solver_ == Teuchos::null || changedParams_)
352 {
353 TEUCHOS_TEST_FOR_EXCEPTION
354 (A_ == Teuchos::null, std::runtime_error, prefix << "The matrix has not been "
355 "set yet. You must call setMatrix() with a nonnull matrix before you may "
356 "call this method.");
357
358 // TODO: We should not have to cast away the constness here
359 // TODO: See bug 6462
360 if(params_ != Teuchos::null)
361 solver_ = CreateTpetraPreconditioner(rcp_const_cast<Tpetra::Operator<Scalar,LO,GO,Node> >(A_), *params_);
362 else
363 solver_ = CreateTpetraPreconditioner(rcp_const_cast<Tpetra::Operator<Scalar,LO,GO,Node> >(A_));
364 }
365 else if(changedA_)
366 {
367 TEUCHOS_TEST_FOR_EXCEPTION
368 (A_ == Teuchos::null, std::runtime_error, prefix << "The matrix has not been "
369 "set yet. You must call setMatrix() with a nonnull matrix before you may "
370 "call this method.");
371
372 // TODO: We should not have to cast away the constness here
373 // TODO: See bug 6462
374 RCP<const Tpetra::CrsMatrix<Scalar,LO,GO,Node> > helperMat;
375 helperMat = rcp_dynamic_cast<const Tpetra::CrsMatrix<Scalar,LO,GO,Node> >(A_);
376 TEUCHOS_TEST_FOR_EXCEPTION
377 (helperMat.is_null(), std::runtime_error, prefix << "MueLu requires "
378 "a Tpetra::CrsMatrix, but the matrix you provided is of a "
379 "different type. Please provide a Tpetra::CrsMatrix instead.");
380 ReuseTpetraPreconditioner(rcp_const_cast<Tpetra::CrsMatrix<Scalar,LO,GO,Node> >(helperMat), *solver_);
381 }
382
383 changedA_ = false;
384 changedParams_ = false;
385 }
386
388 std::string description () const
389 {
390 using Teuchos::TypeNameTraits;
391 if (solver_.is_null()) {
392 std::ostringstream os;
393 os << "\"MueLu::Details::LinearSolver\": {"
394 << "MV: " << TypeNameTraits<Tpetra::MultiVector<Scalar,LO,GO,Node> >::name()
395 << "OP: " << TypeNameTraits<Tpetra::Operator<Scalar,LO,GO,Node> >::name()
396 << "NormType: " << TypeNameTraits<typename Teuchos::ScalarTraits<Scalar>::magnitudeType>::name()
397 << "}";
398 return os.str ();
399 }
400 else {
401 return solver_->GetHierarchy()->description ();
402 }
403 }
404
406 void
407 describe (Teuchos::FancyOStream& out,
408 const Teuchos::EVerbosityLevel verbLevel =
409 Teuchos::Describable::verbLevel_default) const
410 {
411 using Teuchos::TypeNameTraits;
412 using std::endl;
413 if (solver_.is_null()) {
414 if(verbLevel > Teuchos::VERB_NONE) {
415 Teuchos::OSTab tab0 (out);
416 out << "\"MueLu::Details::LinearSolver\":" << endl;
417 Teuchos::OSTab tab1 (out);
418 out << "MV: " << TypeNameTraits<Tpetra::MultiVector<Scalar,LO,GO,Node> >::name() << endl
419 << "OP: " << TypeNameTraits<Tpetra::Operator<Scalar,LO,GO,Node> >::name() << endl
420 << "NormType: " << TypeNameTraits<typename Teuchos::ScalarTraits<Scalar>::magnitudeType>::name() << endl;
421 }
422 }
423 else {
424 solver_->GetHierarchy()->describe (out, verbLevel);
425 }
426 }
427
428private:
429 Teuchos::RCP<const Tpetra::Operator<Scalar,LO,GO,Node> > A_;
430 Teuchos::RCP<Teuchos::ParameterList> params_;
431 Teuchos::RCP<TpetraOperator<Scalar,LO,GO,Node> > solver_;
434};
435#endif // HAVE_MUELU_TPETRA
436
437template<class MV, class OP, class NormType>
438Teuchos::RCP<Trilinos::Details::LinearSolver<MV, OP, NormType> >
440getLinearSolver (const std::string& solverName)
441{
442 using Teuchos::rcp;
444}
445
446template<class MV, class OP, class NormType>
447void
450{
451#ifdef HAVE_TEUCHOSCORE_CXX11
452 typedef std::shared_ptr<MueLu::Details::LinearSolverFactory<MV, OP, NormType> > ptr_type;
453 //typedef std::shared_ptr<Trilinos::Details::LinearSolverFactory<MV, OP> > base_ptr_type;
454#else
455 typedef Teuchos::RCP<MueLu::Details::LinearSolverFactory<MV, OP, NormType> > ptr_type;
456 //typedef Teuchos::RCP<Trilinos::Details::LinearSolverFactory<MV, OP> > base_ptr_type;
457#endif // HAVE_TEUCHOSCORE_CXX11
458
460 Trilinos::Details::registerLinearSolverFactory<MV, OP, NormType> ("MueLu", factory);
461}
462
463} // namespace Details
464} // namespace MueLu
465
466// Macro for doing explicit instantiation of
467// MueLu::Details::LinearSolverFactory, for Tpetra objects, with
468// given Tpetra template parameters (SC = Scalar, LO = LocalOrdinal,
469// GO = GlobalOrdinal, NT = Node).
470//
471// We don't have to protect use of Tpetra objects here, or include
472// any header files for them, because this is a macro definition.
473#define MUELU_DETAILS_LINEARSOLVERFACTORY_INSTANT(SC, LO, GO, NT) \
474 template class MueLu::Details::LinearSolverFactory<Tpetra::MultiVector<SC, LO, GO, NT>, \
475 Tpetra::Operator<SC, LO, GO, NT>, \
476 typename Tpetra::MultiVector<SC, LO, GO, NT>::mag_type>;
477
478#endif // MUELU_DETAILS_LINEARSOLVERFACTORY_DEF_HPP
Various adapters that will create a MueLu preconditioner that is an Epetra_Operator.
Various adapters that will create a MueLu preconditioner that is a Tpetra::Operator.
Interface for a "factory" that creates MueLu solvers.
static void registerLinearSolverFactory()
Register this LinearSolverFactory with the central registry.
virtual Teuchos::RCP< Trilinos::Details::LinearSolver< MV, OP, NormType > > getLinearSolver(const std::string &solverName)
Get an instance of a MueLu solver.
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Implementation of Teuchos::Describable::describe.
void solve(Tpetra::MultiVector< Scalar, LO, GO, Node > &X, const Tpetra::MultiVector< Scalar, LO, GO, Node > &B)
Solve the linear system(s) AX=B.
Teuchos::RCP< Teuchos::ParameterList > params_
std::string description() const
Implementation of Teuchos::Describable::description.
void symbolic()
Set up any part of the solve that depends on the structure of the input matrix, but not its numerical...
virtual ~LinearSolver()
Destructor (virtual for memory safety).
Teuchos::RCP< const OP > getMatrix() const
Get a pointer to this Solver's matrix.
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > &params)
Set this solver's parameters.
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Implementation of Teuchos::Describable::describe.
void solve(MV &X, const MV &B)
Solve the linear system(s) AX=B.
void numeric()
Set up any part of the solve that depends on both the structure and the numerical values of the input...
void setMatrix(const Teuchos::RCP< const OP > &A)
Set the Solver's matrix.
Namespace for MueLu classes and methods.
Teuchos::RCP< MueLu::TpetraOperator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > CreateTpetraPreconditioner(const Teuchos::RCP< Tpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &inA, Teuchos::ParameterList &inParamList)
Helper function to create a MueLu or AMGX preconditioner that can be used by Tpetra....
void ReuseTpetraPreconditioner(const Teuchos::RCP< Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &inA, MueLu::TpetraOperator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op)
Helper function to reuse an existing MueLu preconditioner.
Teuchos::RCP< MueLu::EpetraOperator > CreateEpetraPreconditioner(const Teuchos::RCP< Epetra_CrsMatrix > &inA, Teuchos::ParameterList &paramListIn)
Helper function to create a MueLu preconditioner that can be used by Epetra.Given a EpetraCrs_Matrix,...