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