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