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
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
64namespace Amesos2 {
65
66
74template <class Matrix,
75 class Vector>
76class 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
81public:
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
109
110
117 Superlu(Teuchos::RCP<const Matrix> A,
118 Teuchos::RCP<Vector> X,
119 Teuchos::RCP<const Vector> B);
120
121
123 ~Superlu( );
124
126
128 std::string description() const;
129
130private:
131
137 int preOrdering_impl();
138
139
148
149
156
157
169 int solve_impl(const Teuchos::Ptr<MultiVecAdapter<Vector> > X,
170 const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const;
171
172
176 bool matrixShapeOK_impl() const;
177
178
213 const Teuchos::RCP<Teuchos::ParameterList> & parameterList );
214
215
222 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters_impl() const;
223
224
233 bool loadA_impl(EPhase current_phase);
234
235 typedef Kokkos::DefaultHostExecutionSpace HostExecSpaceType;
236
237 // struct holds all data necessary to make a superlu factorization or solve call
238 mutable struct SLUData {
239 SLU::SuperMatrix A, B, X, L, U; // matrix A in NCformat
240 SLU::SuperMatrix AC; // permuted matrix A in NCPformat
241
242 SLU::superlu_options_t options;
243 SLU::mem_usage_t mem_usage;
244#ifdef HAVE_AMESOS2_SUPERLU5_API
245 SLU::GlobalLU_t lu; // Use for gssvx and gsisx in SuperLU 5.0
246#endif
247 SLU::SuperLUStat_t stat;
248
249
250
251 typedef Kokkos::View<magnitude_type*, HostExecSpaceType> host_mag_array;
252 typedef Kokkos::View<int*, HostExecSpaceType> host_int_array;
253 host_mag_array berr;
254 host_mag_array ferr;
255 host_int_array perm_r;
256 host_int_array perm_c;
257 host_int_array etree;
258 host_mag_array R;
259 host_mag_array C;
260
261#if defined(KOKKOSKERNELS_ENABLE_SUPERNODAL_SPTRSV) && defined(KOKKOSKERNELS_ENABLE_TPL_SUPERLU)
262 host_int_array parents;
263#endif
264
265 char equed;
266 bool rowequ, colequ; // flags what type of equilibration
267 // has been performed
268 magnitude_type anorm, rcond; // condition number estimate
269
270 int relax;
271 int panel_size;
272 } data_;
273
274 typedef int size_type;
275 typedef int ordinal_type;
276 typedef Kokkos::View<size_type*, HostExecSpaceType> host_size_type_array;
277 typedef Kokkos::View<ordinal_type*, HostExecSpaceType> host_ordinal_type_array;
278 typedef Kokkos::View<slu_type*, HostExecSpaceType> host_value_type_array;
279
280 // The following Arrays are persisting storage arrays for A, X, and B
282 host_value_type_array host_nzvals_view_;
283 Teuchos::Array<slu_convert_type> convert_nzvals_; // copy to SuperLU native array before calling SuperLU
284
286 host_size_type_array host_rows_view_;
288 host_ordinal_type_array host_col_ptr_view_;
289
290 typedef typename Kokkos::View<slu_type**, Kokkos::LayoutLeft, HostExecSpaceType>
291 host_solve_array_t;
292
294 mutable host_solve_array_t host_xValues_;
295 mutable Teuchos::Array<slu_convert_type> convert_xValues_; // copy to SuperLU native array before calling SuperLU
296
298 mutable host_solve_array_t host_bValues_;
299 mutable Teuchos::Array<slu_convert_type> convert_bValues_; // copy to SuperLU native array before calling SuperLU
300
301#if defined(KOKKOSKERNELS_ENABLE_SUPERNODAL_SPTRSV) && defined(KOKKOSKERNELS_ENABLE_TPL_SUPERLU)
302 typedef Kokkos::DefaultExecutionSpace DeviceExecSpaceType;
303
304 #ifdef KOKKOS_ENABLE_CUDA
305 // solver will be UVM off even though Tpetra is CudaUVMSpace
306 typedef typename Kokkos::CudaSpace DeviceMemSpaceType;
307 #else
308 typedef typename DeviceExecSpaceType::memory_space DeviceMemSpaceType;
309 #endif
310
311 typedef Kokkos::View<slu_type**, Kokkos::LayoutLeft, DeviceMemSpaceType>
312 device_solve_array_t;
313 // For triangular solves we have both host and device versions of xValues and
314 // bValues because a parameter can turn it on or off.
315 mutable device_solve_array_t device_xValues_;
316 mutable device_solve_array_t device_bValues_;
317 typedef Kokkos::View<int*, DeviceMemSpaceType> device_int_array;
318 typedef Kokkos::View<magnitude_type*, DeviceMemSpaceType> device_mag_array;
319 device_int_array device_trsv_perm_r_;
320 device_int_array device_trsv_perm_c_;
321 device_mag_array device_trsv_R_;
322 device_mag_array device_trsv_C_;
323 mutable device_solve_array_t device_trsv_rhs_;
324 mutable device_solve_array_t device_trsv_sol_;
325 typedef KokkosKernels::Experimental::KokkosKernelsHandle <size_type, ordinal_type, slu_type,
326 DeviceExecSpaceType, DeviceMemSpaceType, DeviceMemSpaceType> kernel_handle_type;
327 mutable kernel_handle_type device_khL_;
328 mutable kernel_handle_type device_khU_;
329 /* parameters for SpTRSV */
330 bool sptrsv_invert_diag_;
331 bool sptrsv_invert_offdiag_;
332 bool sptrsv_u_in_csr_;
333 bool sptrsv_merge_supernodes_;
334 bool sptrsv_use_spmv_;
335#endif
336
337 /* Note: In the above, must use "Amesos2::Superlu" rather than
338 * "Superlu" because otherwise the compiler references the
339 * specialized type of the class, and not the templated type that is
340 * required for Amesos2::TypeMap
341 */
342
343 /* SuperLU can accept input in either compressed-row or
344 * compressed-column storage. We will store and pass matrices in
345 * *compressed-column* format.
346 */
347
348 /*
349 * Internal flag that is used for the numericFactorization_impl
350 * routine. If true, then the superlu gstrf routine should have
351 * SamePattern_SameRowPerm in its options. Otherwise, it should
352 * factor from scratch.
353 *
354 * This is somewhat of a kludge to get around the fact that the
355 * superlu routines all expect something different from the options
356 * struct. The big issue is that we don't want gstrf doing the
357 * symbolic factorization if it doesn't need to. On the other hand,
358 * we can't leave options.Fact set to SamePattern_SameRowPerm
359 * because the solver driver needs it to be set at FACTORED. But
360 * having it set at FACTORED upon re-entrance into
361 * numericFactorization prompts gstrf to redo the symbolic
362 * factorization.
363 */
364 bool same_symbolic_;
365 bool ILU_Flag_;
366
367 bool is_contiguous_;
368 bool use_triangular_solves_;
369
370 void triangular_solve_factor();
371
372 /* call metis before SuperLU */
373 bool use_metis_;
374 bool symmetrize_metis_;
375
376 public: // for GPU
377 void triangular_solve() const; // Only for internal use - public to support kernels
378}; // End class Superlu
379
380
381// Specialize solver_traits struct for SuperLU
382template <>
383struct solver_traits<Superlu> {
384#ifdef HAVE_TEUCHOS_COMPLEX
385 typedef Meta::make_list6<float, double,
386 std::complex<float>, std::complex<double>,
387 Kokkos::complex<float>, Kokkos::complex<double>>
388 supported_scalars;
389#else
390 typedef Meta::make_list2<float, double> supported_scalars;
391#endif
392};
393
394template <typename Scalar, typename LocalOrdinal, typename ExecutionSpace>
395struct solver_supports_matrix<Superlu,
396 KokkosSparse::CrsMatrix<Scalar, LocalOrdinal, ExecutionSpace>> {
397 static const bool value = true;
398};
399
400} // end namespace Amesos2
401
402#endif // AMESOS2_SUPERLU_DECL_HPP
Provides access to interesting solver traits.
Provides a mechanism to map function calls to the correct Solver function based on the scalar type of...
Amesos2::SolverCore: A templated interface for interaction with third-party direct sparse solvers.
Definition Amesos2_SolverCore_decl.hpp:106
Amesos2 interface to the SuperLU package.
Definition Amesos2_Superlu_decl.hpp:77
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition Amesos2_Superlu_def.hpp:699
int numericFactorization_impl()
Superlu specific numeric factorization.
Definition Amesos2_Superlu_def.hpp:317
host_solve_array_t host_xValues_
Persisting 1D store for X.
Definition Amesos2_Superlu_decl.hpp:294
host_ordinal_type_array host_col_ptr_view_
Stores the row indices of the nonzero entries.
Definition Amesos2_Superlu_decl.hpp:288
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition Amesos2_Superlu_def.hpp:205
static const char * name
Name of this solver interface.
Definition Amesos2_Superlu_decl.hpp:84
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition Amesos2_Superlu_def.hpp:818
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:497
host_size_type_array host_rows_view_
Stores the location in Ai_ and Aval_ that starts row j.
Definition Amesos2_Superlu_decl.hpp:286
int symbolicFactorization_impl()
Perform symbolic factorization of the matrix using Superlu.
Definition Amesos2_Superlu_def.hpp:230
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > &parameterList)
Definition Amesos2_Superlu_def.hpp:710
host_value_type_array host_nzvals_view_
Stores the values of the nonzero entries for SuperLU.
Definition Amesos2_Superlu_decl.hpp:282
host_solve_array_t host_bValues_
Persisting 1D store for B.
Definition Amesos2_Superlu_decl.hpp:298
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal structures.
Definition Amesos2_Superlu_def.hpp:962
std::string description() const
Returns a short description of this Solver.
Definition Amesos2_Superlu_def.hpp:151
Passes functions to TPL functions based on type.
Definition Amesos2_FunctionMap.hpp:77
Map types to solver-specific data-types and enums.
Definition Amesos2_TypeMap.hpp:82
Provides traits about solvers.
Definition Amesos2_SolverTraits.hpp:71