Amesos2 - Direct Sparse Solver Interfaces Version of the Day
Amesos2_Details_LinearSolverFactory_def.hpp
Go to the documentation of this file.
1// @HEADER
2//
3// ***********************************************************************
4//
5// Amesos2: Templated Direct Sparse Solver Package
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 AMESOS2_DETAILS_LINEARSOLVERFACTORY_DEF_HPP
49#define AMESOS2_DETAILS_LINEARSOLVERFACTORY_DEF_HPP
50
51#include "Amesos2_config.h"
52#include "Amesos2_Factory.hpp"
53#include "Amesos2_Solver.hpp"
54#include "Trilinos_Details_LinearSolver.hpp"
55#include "Trilinos_Details_LinearSolverFactory.hpp"
56#include <type_traits>
57
58
59#ifdef HAVE_AMESOS2_EPETRA
60# include "Epetra_CrsMatrix.h"
61#endif // HAVE_AMESOS2_EPETRA
62
63// mfh 23 Jul 2015: Tpetra is currently a required dependency of Amesos2.
64#ifndef HAVE_AMESOS2_TPETRA
65# define HAVE_AMESOS2_TPETRA
66#endif // HAVE_AMESOS2_TPETRA
67
68#ifdef HAVE_AMESOS2_TPETRA
69# include "Tpetra_CrsMatrix.hpp"
70#endif // HAVE_AMESOS2_TPETRA
71
72namespace Amesos2 {
73namespace Details {
74
75// For a given linear algebra implementation's Operator type OP,
76// find the corresponding CrsMatrix type.
77//
78// Amesos2 unfortunately only does ETI for Tpetra::CrsMatrix, even
79// though it could very well take Tpetra::RowMatrix.
80template<class OP>
81struct GetMatrixType {
82 typedef int type; // flag (see below)
83
84
85#ifdef HAVE_AMESOS2_EPETRA
86 static_assert(! std::is_same<OP, Epetra_MultiVector>::value,
87 "Amesos2::Details::GetMatrixType: OP = Epetra_MultiVector. "
88 "This probably means that you mixed up MV and OP.");
89#endif // HAVE_AMESOS2_EPETRA
90
91#ifdef HAVE_AMESOS2_TPETRA
92 static_assert(! std::is_same<OP, Tpetra::MultiVector<typename OP::scalar_type,
93 typename OP::local_ordinal_type, typename OP::global_ordinal_type,
94 typename OP::node_type> >::value,
95 "Amesos2::Details::GetMatrixType: OP = Tpetra::MultiVector. "
96 "This probably means that you mixed up MV and OP.");
97#endif // HAVE_AMESOS2_TPETRA
98};
99
100#ifdef HAVE_AMESOS2_EPETRA
101template<>
102struct GetMatrixType<Epetra_Operator> {
103 typedef Epetra_CrsMatrix type;
104};
105#endif // HAVE_AMESOS2_EPETRA
106
107
108#ifdef HAVE_AMESOS2_TPETRA
109template<class S, class LO, class GO, class NT>
110struct GetMatrixType<Tpetra::Operator<S, LO, GO, NT> > {
111 typedef Tpetra::CrsMatrix<S, LO, GO, NT> type;
112};
113#endif // HAVE_AMESOS2_TPETRA
114
115template<class MV, class OP, class NormType>
116class LinearSolver :
117 public Trilinos::Details::LinearSolver<MV, OP, NormType>,
118 virtual public Teuchos::Describable
119{
120
121#ifdef HAVE_AMESOS2_EPETRA
122 static_assert(! std::is_same<OP, Epetra_MultiVector>::value,
123 "Amesos2::Details::LinearSolver: OP = Epetra_MultiVector. "
124 "This probably means that you mixed up MV and OP.");
125 static_assert(! std::is_same<MV, Epetra_Operator>::value,
126 "Amesos2::Details::LinearSolver: MV = Epetra_Operator. "
127 "This probably means that you mixed up MV and OP.");
128#endif // HAVE_AMESOS2_EPETRA
129
130public:
137 typedef typename GetMatrixType<OP>::type crs_matrix_type;
138 static_assert(! std::is_same<crs_matrix_type, int>::value,
139 "Amesos2::Details::LinearSolver: The given OP type is not "
140 "supported.");
141
143 typedef Amesos2::Solver<crs_matrix_type, MV> solver_type;
144
158 LinearSolver (const std::string& solverName) :
159 solverName_ (solverName)
160 {
161 // FIXME (mfh 25 Aug 2015) Ifpack2::Details::Amesos2Wrapper made
162 // the unfortunate choice to attempt to guess solvers that exist.
163 // Thus, alas, for the sake of backwards compatibility, we must
164 // make an attempt to guess for the user. Furthermore, for strict
165 // backwards compatibility, we should preserve the same (rather
166 // arbitrary) list of choices, in the same order.
167 if (solverName == "") {
168 // KLU2 is the default solver.
169 if (Amesos2::query ("klu2")) {
170 solverName_ = "klu2";
171 }
172 else if (Amesos2::query ("superlu")) {
173 solverName_ = "superlu";
174 }
175 else if (Amesos2::query ("superludist")) {
176 solverName_ = "superludist";
177 }
178 else if (Amesos2::query ("cholmod")) {
179 solverName_ = "cholmod";
180 }
181 else if (Amesos2::query ("cusolver")) {
182 solverName_ = "cusolver";
183 }
184 else if (Amesos2::query ("basker")) {
185 solverName_ = "basker";
186 }
187 else if (Amesos2::query ("shylubasker")) {
188 solverName_ = "shylubasker";
189 }
190 else if (Amesos2::query ("ShyLUBasker")) {
191 solverName_ = "shylubasker";
192 }
193 else if (Amesos2::query ("superlumt")) {
194 solverName_ = "superlumt";
195 }
196 else if (Amesos2::query ("pardiso_mkl")) {
197 solverName_ = "pardiso_mkl";
198 }
199 else if (Amesos2::query ("mumps")) {
200 solverName_ = "mumps";
201 }
202 else if (Amesos2::query ("lapack")) {
203 solverName_ = "lapack";
204 }
205 else if (Amesos2::query ("umfpack")) {
206 solverName_ = "umfpack";
207 }
208 // We don't have to try to rescue the user if their empty solver
209 // name didn't catch any of the above choices.
210 }
211 }
212
214 virtual ~LinearSolver () {}
215
220 void setMatrix (const Teuchos::RCP<const OP>& A) {
221 using Teuchos::null;
222 using Teuchos::RCP;
223 using Teuchos::TypeNameTraits;
224 typedef crs_matrix_type MAT;
225 const char prefix[] = "Amesos2::Details::LinearSolver::setMatrix: ";
226
227 if (A.is_null ()) {
228 solver_ = Teuchos::null;
229 }
230 else {
231 // FIXME (mfh 24 Jul 2015) Even though Amesos2 solvers currently
232 // require a CrsMatrix input, we could add a copy step in order
233 // to let them take a RowMatrix. The issue is that this would
234 // require keeping the original matrix around, and copying again
235 // if it changes.
236 RCP<const MAT> A_mat = Teuchos::rcp_dynamic_cast<const MAT> (A);
237 TEUCHOS_TEST_FOR_EXCEPTION
238 (A_mat.is_null (), std::invalid_argument,
239 "Amesos2::Details::LinearSolver::setMatrix: "
240 "The input matrix A must be a CrsMatrix.");
241 if (solver_.is_null ()) {
242 // Amesos2 solvers must be created with a nonnull matrix.
243 // Thus, we don't actually create the solver until setMatrix
244 // is called for the first time with a nonnull matrix.
245 // Furthermore, Amesos2 solvers do not accept a null matrix
246 // (to setA), so if setMatrix is called with a null matrix, we
247 // free the solver.
248 RCP<solver_type> solver;
249 try {
250 solver = Amesos2::create<MAT, MV> (solverName_, A_mat, null, null);
251 }
252 catch (std::exception& e) {
253 TEUCHOS_TEST_FOR_EXCEPTION
254 (true, std::invalid_argument, prefix << "Failed to create Amesos2 "
255 "solver named \"" << solverName_ << "\". "
256 "Amesos2::create<MAT = " << TypeNameTraits<MAT>::name ()
257 << ", MV = " << TypeNameTraits<MV>::name ()
258 << " threw an exception: " << e.what ());
259 }
260 TEUCHOS_TEST_FOR_EXCEPTION
261 (solver.is_null (), std::invalid_argument, prefix << "Failed to "
262 "create Amesos2 solver named \"" << solverName_ << "\". "
263 "Amesos2::create<MAT = " << TypeNameTraits<MAT>::name ()
264 << ", MV = " << TypeNameTraits<MV>::name ()
265 << " returned null.");
266
267 // Use same parameters as before, if user set parameters.
268 if (! params_.is_null ()) {
269 solver->setParameters (params_);
270 }
271 solver_ = solver;
272 } else if (A_ != A) {
273 solver_->setA (A_mat);
274 }
275 }
276
277 // We also have to keep a pointer to A, so that getMatrix() works.
278 A_ = A;
279 }
280
282 Teuchos::RCP<const OP> getMatrix () const {
283 return A_;
284 }
285
287 void solve (MV& X, const MV& B) {
288 const char prefix[] = "Amesos2::Details::LinearSolver::solve: ";
289 TEUCHOS_TEST_FOR_EXCEPTION
290 (solver_.is_null (), std::runtime_error, prefix << "The solver does not "
291 "exist yet. You must call setMatrix() with a nonnull matrix before you "
292 "may call this method.");
293 TEUCHOS_TEST_FOR_EXCEPTION
294 (A_.is_null (), std::runtime_error, prefix << "The matrix has not been "
295 "set yet. You must call setMatrix() with a nonnull matrix before you "
296 "may call this method.");
297 solver_->solve (&X, &B);
298 }
299
301 void setParameters (const Teuchos::RCP<Teuchos::ParameterList>& params) {
302 if (! solver_.is_null ()) {
303 solver_->setParameters (params);
304 }
305 // Remember them, so that if the solver gets recreated, we'll have
306 // the original parameters.
307 params_ = params;
308 }
309
312 void symbolic () {
313 const char prefix[] = "Amesos2::Details::LinearSolver::symbolic: ";
314 TEUCHOS_TEST_FOR_EXCEPTION
315 (solver_.is_null (), std::runtime_error, prefix << "The solver does not "
316 "exist yet. You must call setMatrix() with a nonnull matrix before you "
317 "may call this method.");
318 TEUCHOS_TEST_FOR_EXCEPTION
319 (A_.is_null (), std::runtime_error, prefix << "The matrix has not been "
320 "set yet. You must call setMatrix() with a nonnull matrix before you "
321 "may call this method.");
322 solver_->symbolicFactorization ();
323 }
324
327 void numeric () {
328 const char prefix[] = "Amesos2::Details::LinearSolver::numeric: ";
329 TEUCHOS_TEST_FOR_EXCEPTION
330 (solver_.is_null (), std::runtime_error, prefix << "The solver does not "
331 "exist yet. You must call setMatrix() with a nonnull matrix before you "
332 "may call this method.");
333 TEUCHOS_TEST_FOR_EXCEPTION
334 (A_.is_null (), std::runtime_error, prefix << "The matrix has not been "
335 "set yet. You must call setMatrix() with a nonnull matrix before you "
336 "may call this method.");
337 solver_->numericFactorization ();
338 }
339
341 std::string description () const {
342 using Teuchos::TypeNameTraits;
343 if (solver_.is_null ()) {
344 std::ostringstream os;
345 os << "\"Amesos2::Details::LinearSolver\": {"
346 << "MV: " << TypeNameTraits<MV>::name ()
347 << ", OP: " << TypeNameTraits<OP>::name ()
348 << ", NormType: " << TypeNameTraits<NormType>::name ()
349 << "}";
350 return os.str ();
351 } else {
352 return solver_->description ();
353 }
354 }
355
357 void
358 describe (Teuchos::FancyOStream& out,
359 const Teuchos::EVerbosityLevel verbLevel =
360 Teuchos::Describable::verbLevel_default) const
361 {
362 using Teuchos::TypeNameTraits;
363 using std::endl;
364 if (solver_.is_null ()) {
365 if (verbLevel > Teuchos::VERB_NONE) {
366 Teuchos::OSTab tab0 (out);
367 out << "\"Amesos2::Details::LinearSolver\":" << endl;
368 Teuchos::OSTab tab1 (out);
369 out << "MV: " << TypeNameTraits<MV>::name () << endl
370 << "OP: " << TypeNameTraits<OP>::name () << endl
371 << "NormType: " << TypeNameTraits<NormType>::name () << endl;
372 }
373 }
374 if (! solver_.is_null ()) {
375 solver_->describe (out, verbLevel);
376 }
377 }
378
379private:
380 std::string solverName_;
381 Teuchos::RCP<solver_type> solver_;
382 Teuchos::RCP<const OP> A_;
383 Teuchos::RCP<Teuchos::ParameterList> params_;
384};
385
386template<class MV, class OP, class NormType>
387Teuchos::RCP<Trilinos::Details::LinearSolver<MV, OP, NormType> >
389getLinearSolver (const std::string& solverName)
390{
391 using Teuchos::rcp;
392 return rcp (new Amesos2::Details::LinearSolver<MV, OP, NormType> (solverName));
393}
394
395template<class MV, class OP, class NormType>
396void
399{
400#ifdef HAVE_TEUCHOSCORE_CXX11
401 typedef std::shared_ptr<Amesos2::Details::LinearSolverFactory<MV, OP, NormType> > ptr_type;
402 //typedef std::shared_ptr<Trilinos::Details::LinearSolverFactory<MV, OP> > base_ptr_type;
403#else
404 typedef Teuchos::RCP<Amesos2::Details::LinearSolverFactory<MV, OP, NormType> > ptr_type;
405 //typedef Teuchos::RCP<Trilinos::Details::LinearSolverFactory<MV, OP> > base_ptr_type;
406#endif // HAVE_TEUCHOSCORE_CXX11
407
409 Trilinos::Details::registerLinearSolverFactory<MV, OP, NormType> ("Amesos2", factory);
410}
411
412} // namespace Details
413} // namespace Amesos2
414
415// Macro for doing explicit instantiation of
416// Amesos2::Details::LinearSolverFactory, for Tpetra objects, with
417// given Tpetra template parameters (SC = Scalar, LO = LocalOrdinal,
418// GO = GlobalOrdinal, NT = Node).
419//
420// We don't have to protect use of Tpetra objects here, or include
421// any header files for them, because this is a macro definition.
422#define AMESOS2_DETAILS_LINEARSOLVERFACTORY_INSTANT(SC, LO, GO, NT) \
423 template class Amesos2::Details::LinearSolverFactory<Tpetra::MultiVector<SC, LO, GO, NT>, \
424 Tpetra::Operator<SC, LO, GO, NT>, \
425 typename Tpetra::MultiVector<SC, LO, GO, NT>::mag_type>;
426
427#endif // AMESOS2_DETAILS_LINEARSOLVERFACTORY_DEF_HPP
Contains declarations for Amesos2::create and Amesos2::query.
Interface for a "factory" that creates Amesos2 solvers.
Definition Amesos2_Details_LinearSolverFactory_decl.hpp:78
static void registerLinearSolverFactory()
Register this LinearSolverFactory with the central registry.
Definition Amesos2_Details_LinearSolverFactory_def.hpp:398
virtual Teuchos::RCP< Trilinos::Details::LinearSolver< MV, OP, NormType > > getLinearSolver(const std::string &solverName)
Get an instance of a Amesos2 solver.
Definition Amesos2_Details_LinearSolverFactory_def.hpp:389