Amesos2 - Direct Sparse Solver Interfaces  Version of the Day
Amesos2_Superlu_decl.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 
53 #ifndef AMESOS2_SUPERLU_DECL_HPP
54 #define AMESOS2_SUPERLU_DECL_HPP
55 
56 #include "Amesos2_SolverTraits.hpp"
57 #include "Amesos2_SolverCore.hpp"
59 
60 #if defined(KOKKOSKERNELS_ENABLE_SUPERNODAL_SPTRSV) && defined(KOKKOSKERNELS_ENABLE_TPL_SUPERLU)
61 #include "KokkosKernels_Handle.hpp"
62 #endif
63 
64 namespace Amesos2 {
65 
66 
74 template <class Matrix,
75  class Vector>
76 class Superlu : public SolverCore<Amesos2::Superlu, Matrix, Vector>
77 {
78  friend class SolverCore<Amesos2::Superlu,Matrix,Vector>; // Give our base access
79  // to our private
80  // implementation funcs
81 public:
82 
84  static const char* name; // declaration. Initialization outside.
85 
86  typedef Superlu<Matrix,Vector> type;
87  typedef SolverCore<Amesos2::Superlu,Matrix,Vector> super_type;
88 
89  // Since typedef's are not inheritted, go grab them
90  typedef typename super_type::scalar_type scalar_type;
91  typedef typename super_type::local_ordinal_type local_ordinal_type;
92  typedef typename super_type::global_ordinal_type global_ordinal_type;
93  typedef typename super_type::global_size_type global_size_type;
94 
95  typedef TypeMap<Amesos2::Superlu,scalar_type> type_map;
96 
97  /*
98  * The SuperLU interface will need two other typedef's, which are:
99  * - the superlu type that corresponds to scalar_type and
100  * - the corresponding type to use for magnitude
101  */
102  typedef typename type_map::type slu_type;
103  typedef typename type_map::convert_type slu_convert_type;
104  typedef typename type_map::magnitude_type magnitude_type;
105 
106  typedef FunctionMap<Amesos2::Superlu,slu_type> function_map;
107 
108 #ifdef HAVE_AMESOS2_SUPERLU5_API
109  typedef typename function_map::GlobalLU_type GlobalLU_t;
110 #endif
111 
113 
114 
121  Superlu(Teuchos::RCP<const Matrix> A,
122  Teuchos::RCP<Vector> X,
123  Teuchos::RCP<const Vector> B);
124 
125 
127  ~Superlu( );
128 
130 
132  std::string description() const;
133 
134 private:
135 
141  int preOrdering_impl();
142 
143 
152 
153 
160 
161 
173  int solve_impl(const Teuchos::Ptr<MultiVecAdapter<Vector> > X,
174  const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const;
175 
176 
180  bool matrixShapeOK_impl() const;
181 
182 
216  void setParameters_impl(
217  const Teuchos::RCP<Teuchos::ParameterList> & parameterList );
218 
219 
226  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters_impl() const;
227 
228 
237  bool loadA_impl(EPhase current_phase);
238 
239  typedef Kokkos::DefaultHostExecutionSpace HostExecSpaceType;
240 
241  // struct holds all data necessary to make a superlu factorization or solve call
242  mutable struct SLUData {
243  SLU::SuperMatrix A, B, X, L, U; // matrix A in NCformat
244  SLU::SuperMatrix AC; // permuted matrix A in NCPformat
245 
246  SLU::superlu_options_t options;
247  SLU::mem_usage_t mem_usage;
248 #ifdef HAVE_AMESOS2_SUPERLU5_API
249  GlobalLU_t lu; // Use for gssvx and gsisx in SuperLU 5.0
250 #endif
251  SLU::SuperLUStat_t stat;
252 
253 
254 
255  typedef Kokkos::View<magnitude_type*, HostExecSpaceType> host_mag_array;
256  typedef Kokkos::View<int*, HostExecSpaceType> host_int_array;
257  host_mag_array berr;
258  host_mag_array ferr;
259  host_int_array perm_r;
260  host_int_array perm_c;
261  host_int_array etree;
262  host_mag_array R;
263  host_mag_array C;
264 
265 #if defined(KOKKOSKERNELS_ENABLE_SUPERNODAL_SPTRSV) && defined(KOKKOSKERNELS_ENABLE_TPL_SUPERLU)
266  host_int_array parents;
267 #endif
268 
269  char equed;
270  bool rowequ, colequ; // flags what type of equilibration
271  // has been performed
272 
273  int relax;
274  int panel_size;
275  } data_;
276 
277  typedef int size_type;
278  typedef int ordinal_type;
279  typedef Kokkos::View<size_type*, HostExecSpaceType> host_size_type_array;
280  typedef Kokkos::View<ordinal_type*, HostExecSpaceType> host_ordinal_type_array;
281  typedef Kokkos::View<slu_type*, HostExecSpaceType> host_value_type_array;
282 
283  // The following Arrays are persisting storage arrays for A, X, and B
285  host_value_type_array host_nzvals_view_;
286  Teuchos::Array<slu_convert_type> convert_nzvals_; // copy to SuperLU native array before calling SuperLU
287 
289  host_size_type_array host_rows_view_;
291  host_ordinal_type_array host_col_ptr_view_;
292 
293  typedef typename Kokkos::View<slu_type**, Kokkos::LayoutLeft, HostExecSpaceType>
294  host_solve_array_t;
295 
297  mutable host_solve_array_t host_xValues_;
298  mutable Teuchos::Array<slu_convert_type> convert_xValues_; // copy to SuperLU native array before calling SuperLU
299 
301  mutable host_solve_array_t host_bValues_;
302  mutable Teuchos::Array<slu_convert_type> convert_bValues_; // copy to SuperLU native array before calling SuperLU
303 
304 #if defined(KOKKOSKERNELS_ENABLE_SUPERNODAL_SPTRSV) && defined(KOKKOSKERNELS_ENABLE_TPL_SUPERLU)
305  typedef Kokkos::DefaultExecutionSpace DeviceExecSpaceType;
306 
307  #ifdef KOKKOS_ENABLE_CUDA
308  // solver will be UVM off even though Tpetra is CudaUVMSpace
309  typedef typename Kokkos::CudaSpace DeviceMemSpaceType;
310  #else
311  typedef typename DeviceExecSpaceType::memory_space DeviceMemSpaceType;
312  #endif
313 
314  typedef Kokkos::View<slu_type**, Kokkos::LayoutLeft, DeviceMemSpaceType>
315  device_solve_array_t;
316  // For triangular solves we have both host and device versions of xValues and
317  // bValues because a parameter can turn it on or off.
318  mutable device_solve_array_t device_xValues_;
319  mutable device_solve_array_t device_bValues_;
320  typedef Kokkos::View<int*, DeviceMemSpaceType> device_int_array;
321  typedef Kokkos::View<magnitude_type*, DeviceMemSpaceType> device_mag_array;
322  device_int_array device_trsv_perm_r_;
323  device_int_array device_trsv_perm_c_;
324  device_mag_array device_trsv_R_;
325  device_mag_array device_trsv_C_;
326  mutable device_solve_array_t device_trsv_rhs_;
327  mutable device_solve_array_t device_trsv_sol_;
328  typedef KokkosKernels::Experimental::KokkosKernelsHandle <size_type, ordinal_type, slu_type,
329  DeviceExecSpaceType, DeviceMemSpaceType, DeviceMemSpaceType> kernel_handle_type;
330  mutable kernel_handle_type device_khL_;
331  mutable kernel_handle_type device_khU_;
332  /* parameters for SpTRSV */
333  bool sptrsv_invert_diag_;
334  bool sptrsv_invert_offdiag_;
335  bool sptrsv_u_in_csr_;
336  bool sptrsv_merge_supernodes_;
337  bool sptrsv_use_spmv_;
338 #endif
339 
340  /* Note: In the above, must use "Amesos2::Superlu" rather than
341  * "Superlu" because otherwise the compiler references the
342  * specialized type of the class, and not the templated type that is
343  * required for Amesos2::TypeMap
344  */
345 
346  /* SuperLU can accept input in either compressed-row or
347  * compressed-column storage. We will store and pass matrices in
348  * *compressed-column* format.
349  */
350 
351  /*
352  * Internal flag that is used for the numericFactorization_impl
353  * routine. If true, then the superlu gstrf routine should have
354  * SamePattern_SameRowPerm in its options. Otherwise, it should
355  * factor from scratch.
356  *
357  * This is somewhat of a kludge to get around the fact that the
358  * superlu routines all expect something different from the options
359  * struct. The big issue is that we don't want gstrf doing the
360  * symbolic factorization if it doesn't need to. On the other hand,
361  * we can't leave options.Fact set to SamePattern_SameRowPerm
362  * because the solver driver needs it to be set at FACTORED. But
363  * having it set at FACTORED upon re-entrance into
364  * numericFactorization prompts gstrf to redo the symbolic
365  * factorization.
366  */
367  bool same_symbolic_;
368  bool ILU_Flag_;
369 
370  bool is_contiguous_;
371  bool use_triangular_solves_;
372 
373  void triangular_solve_factor();
374 
375  /* call metis before SuperLU */
376  bool use_metis_;
377  bool symmetrize_metis_;
378 
379  public: // for GPU
380  void triangular_solve() const; // Only for internal use - public to support kernels
381 }; // End class Superlu
382 
383 
384 // Specialize solver_traits struct for SuperLU
385 template <>
386 struct solver_traits<Superlu> {
387 #ifdef HAVE_TEUCHOS_COMPLEX
388  typedef Meta::make_list6<float, double,
389  std::complex<float>, std::complex<double>,
390  Kokkos::complex<float>, Kokkos::complex<double>>
391  supported_scalars;
392 #else
393  typedef Meta::make_list2<float, double> supported_scalars;
394 #endif
395 };
396 
397 template <typename Scalar, typename LocalOrdinal, typename ExecutionSpace>
398 struct solver_supports_matrix<Superlu,
399  KokkosSparse::CrsMatrix<Scalar, LocalOrdinal, ExecutionSpace>> {
400  static const bool value = true;
401 };
402 
403 } // end namespace Amesos2
404 
405 #endif // AMESOS2_SUPERLU_DECL_HPP
Amesos2::SolverCore: A templated interface for interaction with third-party direct sparse solvers...
Definition: Amesos2_SolverCore_decl.hpp:105
int solve_impl(const Teuchos::Ptr< MultiVecAdapter< Vector > > X, const Teuchos::Ptr< const MultiVecAdapter< Vector > > B) const
Superlu specific solve.
Definition: Amesos2_Superlu_def.hpp:462
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal structures.
Definition: Amesos2_Superlu_def.hpp:924
int numericFactorization_impl()
Superlu specific numeric factorization.
Definition: Amesos2_Superlu_def.hpp:313
Map types to solver-specific data-types and enums.
Definition: Amesos2_TypeMap.hpp:82
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_Superlu_def.hpp:780
host_value_type_array host_nzvals_view_
Stores the values of the nonzero entries for SuperLU.
Definition: Amesos2_Superlu_decl.hpp:285
Provides a mechanism to map function calls to the correct Solver function based on the scalar type of...
Provides traits about solvers.
Definition: Amesos2_SolverTraits.hpp:71
host_size_type_array host_rows_view_
Stores the location in Ai_ and Aval_ that starts row j.
Definition: Amesos2_Superlu_decl.hpp:289
Definition: Amesos2_AbstractConcreteMatrixAdapter.hpp:48
host_solve_array_t host_xValues_
Persisting 1D store for X.
Definition: Amesos2_Superlu_decl.hpp:297
host_ordinal_type_array host_col_ptr_view_
Stores the row indices of the nonzero entries.
Definition: Amesos2_Superlu_decl.hpp:291
host_solve_array_t host_bValues_
Persisting 1D store for B.
Definition: Amesos2_Superlu_decl.hpp:301
Interface to Amesos2 solver objects.
Definition: Amesos2_Solver_decl.hpp:78
std::string description() const
Returns a short description of this Solver.
Definition: Amesos2_Superlu_def.hpp:151
std::string name() const
Return the name of this solver.
Definition: Amesos2_SolverCore_def.hpp:509
Passes functions to TPL functions based on type.
Definition: Amesos2_FunctionMap.hpp:76
Provides access to interesting solver traits.
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_Superlu_def.hpp:664
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > &parameterList)
Definition: Amesos2_Superlu_def.hpp:675
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_Superlu_def.hpp:205
int symbolicFactorization_impl()
Perform symbolic factorization of the matrix using Superlu.
Definition: Amesos2_Superlu_def.hpp:230
Amesos2 interface to the SuperLU package.
Definition: Amesos2_Superlu_decl.hpp:76