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 
72 namespace Amesos2 {
73 namespace 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.
80 template<class OP>
81 struct 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
101 template<>
102 struct GetMatrixType<Epetra_Operator> {
103  typedef Epetra_CrsMatrix type;
104 };
105 #endif // HAVE_AMESOS2_EPETRA
106 
107 
108 #ifdef HAVE_AMESOS2_TPETRA
109 template<class S, class LO, class GO, class NT>
110 struct GetMatrixType<Tpetra::Operator<S, LO, GO, NT> > {
111  typedef Tpetra::CrsMatrix<S, LO, GO, NT> type;
112 };
113 #endif // HAVE_AMESOS2_TPETRA
114 
115 template<class MV, class OP, class NormType>
116 class 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 
130 public:
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 
379 private:
380  std::string solverName_;
381  Teuchos::RCP<solver_type> solver_;
382  Teuchos::RCP<const OP> A_;
383  Teuchos::RCP<Teuchos::ParameterList> params_;
384 };
385 
386 template<class MV, class OP, class NormType>
387 Teuchos::RCP<Trilinos::Details::LinearSolver<MV, OP, NormType> >
389 getLinearSolver (const std::string& solverName)
390 {
391  using Teuchos::rcp;
392  return rcp (new Amesos2::Details::LinearSolver<MV, OP, NormType> (solverName));
393 }
394 
395 template<class MV, class OP, class NormType>
396 void
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
Interface for a &quot;factory&quot; that creates Amesos2 solvers.
Definition: Amesos2_Details_LinearSolverFactory_decl.hpp:77
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
Interface to Amesos2 solver objects.
Definition: Amesos2_Solver_decl.hpp:78
static void registerLinearSolverFactory()
Register this LinearSolverFactory with the central registry.
Definition: Amesos2_Details_LinearSolverFactory_def.hpp:398
Contains declarations for Amesos2::create and Amesos2::query.