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 ("superlumt")) {
193  solverName_ = "superlumt";
194  }
195  else if (Amesos2::query ("pardiso_mkl")) {
196  solverName_ = "pardiso_mkl";
197  }
198  else if (Amesos2::query ("mumps")) {
199  solverName_ = "mumps";
200  }
201  else if (Amesos2::query ("lapack")) {
202  solverName_ = "lapack";
203  }
204  // We don't have to try to rescue the user if their empty solver
205  // name didn't catch any of the above choices.
206  }
207  }
208 
210  virtual ~LinearSolver () {}
211 
216  void setMatrix (const Teuchos::RCP<const OP>& A) {
217  using Teuchos::null;
218  using Teuchos::RCP;
219  using Teuchos::TypeNameTraits;
220  typedef crs_matrix_type MAT;
221  const char prefix[] = "Amesos2::Details::LinearSolver::setMatrix: ";
222 
223  if (A.is_null ()) {
224  solver_ = Teuchos::null;
225  }
226  else {
227  // FIXME (mfh 24 Jul 2015) Even though Amesos2 solvers currently
228  // require a CrsMatrix input, we could add a copy step in order
229  // to let them take a RowMatrix. The issue is that this would
230  // require keeping the original matrix around, and copying again
231  // if it changes.
232  RCP<const MAT> A_mat = Teuchos::rcp_dynamic_cast<const MAT> (A);
233  TEUCHOS_TEST_FOR_EXCEPTION
234  (A_mat.is_null (), std::invalid_argument,
235  "Amesos2::Details::LinearSolver::setMatrix: "
236  "The input matrix A must be a CrsMatrix.");
237  if (solver_.is_null ()) {
238  // Amesos2 solvers must be created with a nonnull matrix.
239  // Thus, we don't actually create the solver until setMatrix
240  // is called for the first time with a nonnull matrix.
241  // Furthermore, Amesos2 solvers do not accept a null matrix
242  // (to setA), so if setMatrix is called with a null matrix, we
243  // free the solver.
244  RCP<solver_type> solver;
245  try {
246  solver = Amesos2::create<MAT, MV> (solverName_, A_mat, null, null);
247  }
248  catch (std::exception& e) {
249  TEUCHOS_TEST_FOR_EXCEPTION
250  (true, std::invalid_argument, prefix << "Failed to create Amesos2 "
251  "solver named \"" << solverName_ << "\". "
252  "Amesos2::create<MAT = " << TypeNameTraits<MAT>::name ()
253  << ", MV = " << TypeNameTraits<MV>::name ()
254  << " threw an exception: " << e.what ());
255  }
256  TEUCHOS_TEST_FOR_EXCEPTION
257  (solver.is_null (), std::invalid_argument, prefix << "Failed to "
258  "create Amesos2 solver named \"" << solverName_ << "\". "
259  "Amesos2::create<MAT = " << TypeNameTraits<MAT>::name ()
260  << ", MV = " << TypeNameTraits<MV>::name ()
261  << " returned null.");
262 
263  // Use same parameters as before, if user set parameters.
264  if (! params_.is_null ()) {
265  solver->setParameters (params_);
266  }
267  solver_ = solver;
268  } else if (A_ != A) {
269  solver_->setA (A_mat);
270  }
271  }
272 
273  // We also have to keep a pointer to A, so that getMatrix() works.
274  A_ = A;
275  }
276 
278  Teuchos::RCP<const OP> getMatrix () const {
279  return A_;
280  }
281 
283  void solve (MV& X, const MV& B) {
284  const char prefix[] = "Amesos2::Details::LinearSolver::solve: ";
285  TEUCHOS_TEST_FOR_EXCEPTION
286  (solver_.is_null (), std::runtime_error, prefix << "The solver does not "
287  "exist yet. You must call setMatrix() with a nonnull matrix before you "
288  "may call this method.");
289  TEUCHOS_TEST_FOR_EXCEPTION
290  (A_.is_null (), std::runtime_error, prefix << "The matrix has not been "
291  "set yet. You must call setMatrix() with a nonnull matrix before you "
292  "may call this method.");
293  solver_->solve (&X, &B);
294  }
295 
297  void setParameters (const Teuchos::RCP<Teuchos::ParameterList>& params) {
298  if (! solver_.is_null ()) {
299  solver_->setParameters (params);
300  }
301  // Remember them, so that if the solver gets recreated, we'll have
302  // the original parameters.
303  params_ = params;
304  }
305 
308  void symbolic () {
309  const char prefix[] = "Amesos2::Details::LinearSolver::symbolic: ";
310  TEUCHOS_TEST_FOR_EXCEPTION
311  (solver_.is_null (), std::runtime_error, prefix << "The solver does not "
312  "exist yet. You must call setMatrix() with a nonnull matrix before you "
313  "may call this method.");
314  TEUCHOS_TEST_FOR_EXCEPTION
315  (A_.is_null (), std::runtime_error, prefix << "The matrix has not been "
316  "set yet. You must call setMatrix() with a nonnull matrix before you "
317  "may call this method.");
318  solver_->symbolicFactorization ();
319  }
320 
323  void numeric () {
324  const char prefix[] = "Amesos2::Details::LinearSolver::numeric: ";
325  TEUCHOS_TEST_FOR_EXCEPTION
326  (solver_.is_null (), std::runtime_error, prefix << "The solver does not "
327  "exist yet. You must call setMatrix() with a nonnull matrix before you "
328  "may call this method.");
329  TEUCHOS_TEST_FOR_EXCEPTION
330  (A_.is_null (), std::runtime_error, prefix << "The matrix has not been "
331  "set yet. You must call setMatrix() with a nonnull matrix before you "
332  "may call this method.");
333  solver_->numericFactorization ();
334  }
335 
337  std::string description () const {
338  using Teuchos::TypeNameTraits;
339  if (solver_.is_null ()) {
340  std::ostringstream os;
341  os << "\"Amesos2::Details::LinearSolver\": {"
342  << "MV: " << TypeNameTraits<MV>::name ()
343  << ", OP: " << TypeNameTraits<OP>::name ()
344  << ", NormType: " << TypeNameTraits<NormType>::name ()
345  << "}";
346  return os.str ();
347  } else {
348  return solver_->description ();
349  }
350  }
351 
353  void
354  describe (Teuchos::FancyOStream& out,
355  const Teuchos::EVerbosityLevel verbLevel =
356  Teuchos::Describable::verbLevel_default) const
357  {
358  using Teuchos::TypeNameTraits;
359  using std::endl;
360  if (solver_.is_null ()) {
361  if (verbLevel > Teuchos::VERB_NONE) {
362  Teuchos::OSTab tab0 (out);
363  out << "\"Amesos2::Details::LinearSolver\":" << endl;
364  Teuchos::OSTab tab1 (out);
365  out << "MV: " << TypeNameTraits<MV>::name () << endl
366  << "OP: " << TypeNameTraits<OP>::name () << endl
367  << "NormType: " << TypeNameTraits<NormType>::name () << endl;
368  }
369  }
370  if (! solver_.is_null ()) {
371  solver_->describe (out, verbLevel);
372  }
373  }
374 
375 private:
376  std::string solverName_;
377  Teuchos::RCP<solver_type> solver_;
378  Teuchos::RCP<const OP> A_;
379  Teuchos::RCP<Teuchos::ParameterList> params_;
380 };
381 
382 template<class MV, class OP, class NormType>
383 Teuchos::RCP<Trilinos::Details::LinearSolver<MV, OP, NormType> >
385 getLinearSolver (const std::string& solverName)
386 {
387  using Teuchos::rcp;
388  return rcp (new Amesos2::Details::LinearSolver<MV, OP, NormType> (solverName));
389 }
390 
391 template<class MV, class OP, class NormType>
392 void
395 {
396 #ifdef HAVE_TEUCHOSCORE_CXX11
397  typedef std::shared_ptr<Amesos2::Details::LinearSolverFactory<MV, OP, NormType> > ptr_type;
398  //typedef std::shared_ptr<Trilinos::Details::LinearSolverFactory<MV, OP> > base_ptr_type;
399 #else
400  typedef Teuchos::RCP<Amesos2::Details::LinearSolverFactory<MV, OP, NormType> > ptr_type;
401  //typedef Teuchos::RCP<Trilinos::Details::LinearSolverFactory<MV, OP> > base_ptr_type;
402 #endif // HAVE_TEUCHOSCORE_CXX11
403 
405  Trilinos::Details::registerLinearSolverFactory<MV, OP, NormType> ("Amesos2", factory);
406 }
407 
408 } // namespace Details
409 } // namespace Amesos2
410 
411 // Macro for doing explicit instantiation of
412 // Amesos2::Details::LinearSolverFactory, for Tpetra objects, with
413 // given Tpetra template parameters (SC = Scalar, LO = LocalOrdinal,
414 // GO = GlobalOrdinal, NT = Node).
415 //
416 // We don't have to protect use of Tpetra objects here, or include
417 // any header files for them, because this is a macro definition.
418 #define AMESOS2_DETAILS_LINEARSOLVERFACTORY_INSTANT(SC, LO, GO, NT) \
419  template class Amesos2::Details::LinearSolverFactory<Tpetra::MultiVector<SC, LO, GO, NT>, \
420  Tpetra::Operator<SC, LO, GO, NT>, \
421  typename Tpetra::MultiVector<SC, LO, GO, NT>::mag_type>;
422 
423 #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:385
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:394
Contains declarations for Amesos2::create and Amesos2::query.