Amesos2 - Direct Sparse Solver Interfaces Version of the Day
Amesos2_cuSOLVER_def.hpp
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
44#ifndef AMESOS2_CUSOLVER_DEF_HPP
45#define AMESOS2_CUSOLVER_DEF_HPP
46
47#include <Teuchos_Tuple.hpp>
48#include <Teuchos_ParameterList.hpp>
49#include <Teuchos_StandardParameterEntryValidators.hpp>
50
52#include "Amesos2_cuSOLVER_decl.hpp"
53
54namespace Amesos2 {
55
56template <class Matrix, class Vector>
58 Teuchos::RCP<const Matrix> A,
59 Teuchos::RCP<Vector> X,
60 Teuchos::RCP<const Vector> B )
61 : SolverCore<Amesos2::cuSOLVER,Matrix,Vector>(A, X, B)
62{
63 auto status = cusolverSpCreate(&data_.handle);
64 TEUCHOS_TEST_FOR_EXCEPTION( status != CUSOLVER_STATUS_SUCCESS,
65 std::runtime_error, "cusolverSpCreate failed");
66
67 status = cusolverSpCreateCsrcholInfo(&data_.chol_info);
68 TEUCHOS_TEST_FOR_EXCEPTION( status != CUSOLVER_STATUS_SUCCESS,
69 std::runtime_error, "cusolverSpCreateCsrcholInfo failed");
70
71 auto sparse_status = cusparseCreateMatDescr(&data_.desc);
72 TEUCHOS_TEST_FOR_EXCEPTION( sparse_status != CUSPARSE_STATUS_SUCCESS,
73 std::runtime_error, "cusparseCreateMatDescr failed");
74}
75
76template <class Matrix, class Vector>
78{
79 cusparseDestroyMatDescr(data_.desc);
80 cusolverSpDestroyCsrcholInfo(data_.chol_info);
81 cusolverSpDestroy(data_.handle);
82}
83
84template<class Matrix, class Vector>
85int
87{
88#ifdef HAVE_AMESOS2_TIMERS
89 Teuchos::TimeMonitor preOrderTimer(this->timers_.preOrderTime_);
90#endif
91 if(do_optimization()) {
92 this->matrixA_->returnRowPtr_kokkos_view(device_row_ptr_view_);
93 this->matrixA_->returnColInd_kokkos_view(device_cols_view_);
94
95 // reorder to optimize cuSolver
96 if(data_.bReorder) {
97 Amesos2::Util::reorder(
98 device_row_ptr_view_, device_cols_view_,
99 device_perm_, device_peri_, sorted_nnz,
100 true);
101 }
102 }
103
104 return(0);
105}
106
107template <class Matrix, class Vector>
108int
110{
111#ifdef HAVE_AMESOS2_TIMERS
112 Teuchos::TimeMonitor symFactTimer(this->timers_.symFactTime_);
113#endif
114
115 int err = 0;
116 if ( this->root_ ) {
117 const int size = this->globalNumRows_;
118 const int nnz = device_cols_view_.size(); // reorder may have changed this
119 const int * colIdx = device_cols_view_.data();
120 const int * rowPtr = device_row_ptr_view_.data();
121 auto status = cusolverSpXcsrcholAnalysis(
122 data_.handle, size, nnz, data_.desc, rowPtr, colIdx, data_.chol_info);
123 err = (status != CUSOLVER_STATUS_SUCCESS) ? 1 : 0;
124 }
125
126 Teuchos::broadcast(*(this->getComm()), 0, &err);
127 TEUCHOS_TEST_FOR_EXCEPTION(err != 0,
128 std::runtime_error, "Amesos2 cuSolver symbolic failed.");
129
130 return err;
131}
132
133template <class Matrix, class Vector>
134int
136{
137 int err = 0;
138 if(do_optimization()) { // just supporting one rank right now
139 this->matrixA_->returnValues_kokkos_view(device_nzvals_view_);
140
141 // reorder to optimize cuSolver
142 if(data_.bReorder) {
143 // must have original row and cols - maybe cache this from 1st symbiolic setup
144 // this setup exists to support the refactor option
145 device_size_type_array orig_device_row_ptr_view;
146 device_ordinal_type_array orig_device_cols_view;
147 this->matrixA_->returnRowPtr_kokkos_view(orig_device_row_ptr_view);
148 this->matrixA_->returnColInd_kokkos_view(orig_device_cols_view);
149 Amesos2::Util::reorder_values(
150 device_nzvals_view_, orig_device_row_ptr_view, device_row_ptr_view_, orig_device_cols_view,
151 device_perm_, device_peri_, sorted_nnz);
152 }
153
154 const int size = this->globalNumRows_;
155 const int nnz = device_cols_view_.size(); // reorder may have changed this
156 const cusolver_type * values = device_nzvals_view_.data();
157 const int * colIdx = device_cols_view_.data();
158 const int * rowPtr = device_row_ptr_view_.data();
159
160 size_t internalDataInBytes, workspaceInBytes;
161 auto status = function_map::bufferInfo(data_.handle, size, nnz, data_.desc,
162 values, rowPtr, colIdx, data_.chol_info,
163 &internalDataInBytes, &workspaceInBytes);
164
165 if(status == CUSOLVER_STATUS_SUCCESS) {
166 const size_t buffer_size = workspaceInBytes / sizeof(cusolver_type);
167 if(buffer_size > buffer_.extent(0)) {
168 buffer_ = device_value_type_array(
169 Kokkos::ViewAllocateWithoutInitializing("cusolver buf"), buffer_size);
170 }
171 status = function_map::numeric(data_.handle, size, nnz, data_.desc,
172 values, rowPtr, colIdx, data_.chol_info, buffer_.data());
173 }
174 err = (status != CUSOLVER_STATUS_SUCCESS) ? 1 : 0;
175 }
176
177 Teuchos::broadcast(*(this->getComm()), 0, &err);
178 TEUCHOS_TEST_FOR_EXCEPTION(err != 0,
179 std::runtime_error, "Amesos2 cuSolver numeric failed.");
180
181 return err;
182}
183
184template <class Matrix, class Vector>
185int
187 const Teuchos::Ptr<MultiVecAdapter<Vector> > X,
188 const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const
189{
190 const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
191 const ordinal_type nrhs = X->getGlobalNumVectors();
192
193 bool bAssignedX;
194 { // Get values from RHS B
195#ifdef HAVE_AMESOS2_TIMERS
196 Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
197 Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
198#endif
199
200 const bool initialize_data = true;
201 const bool do_not_initialize_data = false;
202 Util::get_1d_copy_helper_kokkos_view<MultiVecAdapter<Vector>,
203 device_solve_array_t>::do_get(initialize_data, B, this->bValues_, Teuchos::as<size_t>(ld_rhs),
204 ROOTED, this->rowIndexBase_);
205
206 // In general we may want to write directly to the x space without a copy.
207 // So we 'get' x which may be a direct view assignment to the MV.
208 bAssignedX = Util::get_1d_copy_helper_kokkos_view<MultiVecAdapter<Vector>,
209 device_solve_array_t>::do_get(do_not_initialize_data, X, this->xValues_, Teuchos::as<size_t>(ld_rhs),
210 ROOTED, this->rowIndexBase_);
211 }
212
213 int err = 0;
214
215 if ( this->root_ ) { // Do solve!
216#ifdef HAVE_AMESOS2_TIMER
217 Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
218#endif
219
220 const int size = this->globalNumRows_;
221
222 if(data_.bReorder) {
223 Amesos2::Util::apply_reorder_permutation(
224 this->bValues_, this->permute_result_, this->device_perm_);
225 }
226 else {
227 this->permute_result_ = this->bValues_; // no permutation
228 }
229
230 for(ordinal_type n = 0; n < nrhs; ++n) {
231 const cusolver_type * b = this->permute_result_.data() + n * size;
232 cusolver_type * x = this->xValues_.data() + n * size;
233 auto status = function_map::solve(
234 data_.handle, size, b, x, data_.chol_info, buffer_.data());
235 err = (status != CUSOLVER_STATUS_SUCCESS) ? 1 : 0;
236 if(err != 0) {
237 break;
238 }
239 }
240
241 if(data_.bReorder && err == 0) {
242 Amesos2::Util::apply_reorder_permutation(
243 this->xValues_, this->permute_result_, this->device_peri_);
244 Kokkos::deep_copy(this->xValues_, this->permute_result_); // full copy since permute_result_ is reused
245 }
246 }
247
248 /* Update X's global values */
249
250 // if bDidAssignX, then we solved straight to the adapter's X memory space without
251 // requiring additional memory allocation, so the x data is already in place.
252 if(!bAssignedX) {
253#ifdef HAVE_AMESOS2_TIMERS
254 Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
255#endif
256
257 Util::template put_1d_data_helper_kokkos_view<
258 MultiVecAdapter<Vector>,device_solve_array_t>::do_put(X, xValues_,
259 Teuchos::as<size_t>(ld_rhs), ROOTED, this->rowIndexBase_);
260 }
261
262 Teuchos::broadcast(*(this->getComm()), 0, &err);
263 TEUCHOS_TEST_FOR_EXCEPTION(err != 0,
264 std::runtime_error, "Amesos2 cuSolver solve failed.");
265
266 return err;
267}
268
269template <class Matrix, class Vector>
270bool
272{
273 return( this->matrixA_->getGlobalNumRows() == this->matrixA_->getGlobalNumCols() );
274}
275
276template <class Matrix, class Vector>
277void
278cuSOLVER<Matrix,Vector>::setParameters_impl(const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
279{
280 using Teuchos::RCP;
281 using Teuchos::ParameterEntryValidator;
282
283 RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
284
285 if( parameterList->isParameter("Reorder") ){
286 RCP<const ParameterEntryValidator> reorder_validator = valid_params->getEntry("Reorder").validator();
287 parameterList->getEntry("Reorder").setValidator(reorder_validator);
288 }
289
290 data_.bReorder = parameterList->get<bool>("Reorder", true);
291}
292
293template <class Matrix, class Vector>
294Teuchos::RCP<const Teuchos::ParameterList>
296{
297 static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
298
299 if( is_null(valid_params) ){
300 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
301
302 pl->set("Reorder", true, "Whether GIDs contiguous");
303
304 valid_params = pl;
305 }
306
307 return valid_params;
308}
309
310template <class Matrix, class Vector>
311bool
313 return (this->root_ && (this->matrixA_->getComm()->getSize() == 1));
314}
315
316template <class Matrix, class Vector>
317bool
319{
320 if(current_phase == SOLVE) {
321 return(false);
322 }
323
324 if(!do_optimization()) { // we're only doing serial right now for cuSolver
325 TEUCHOS_TEST_FOR_EXCEPTION( true, std::runtime_error,
326 "cuSolver is only implemented for serial.");
327 }
328
329 return true;
330}
331
332template<class Matrix, class Vector>
333const char* cuSOLVER<Matrix,Vector>::name = "cuSOLVER";
334
335} // end namespace Amesos2
336
337#endif // AMESOS2_CUSOLVER_DEF_HPP
@ ROOTED
Definition: Amesos2_TypeDecl.hpp:127
Amesos2::SolverCore: A templated interface for interaction with third-party direct sparse solvers.
Definition: Amesos2_SolverCore_decl.hpp:106
Amesos2 interface to cuSOLVER.
Definition: Amesos2_cuSOLVER_decl.hpp:60
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > &parameterList)
Definition: Amesos2_cuSOLVER_def.hpp:278
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_cuSOLVER_def.hpp:86
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal structures.
Definition: Amesos2_cuSOLVER_def.hpp:318
cuSOLVER(Teuchos::RCP< const Matrix > A, Teuchos::RCP< Vector > X, Teuchos::RCP< const Vector > B)
Initialize from Teuchos::RCP.
Definition: Amesos2_cuSOLVER_def.hpp:57
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_cuSOLVER_def.hpp:271
int solve_impl(const Teuchos::Ptr< MultiVecAdapter< Vector > > X, const Teuchos::Ptr< const MultiVecAdapter< Vector > > B) const
cuSOLVER specific solve.
Definition: Amesos2_cuSOLVER_def.hpp:186
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_cuSOLVER_def.hpp:295
int symbolicFactorization_impl()
Perform symbolic factorization of the matrix using cuSOLVER.
Definition: Amesos2_cuSOLVER_def.hpp:109
bool do_optimization() const
can we optimize size_type and ordinal_type for straight pass through
Definition: Amesos2_cuSOLVER_def.hpp:312
int numericFactorization_impl()
cuSOLVER specific numeric factorization
Definition: Amesos2_cuSOLVER_def.hpp:135
~cuSOLVER()
Destructor.
Definition: Amesos2_cuSOLVER_def.hpp:77
EPhase
Used to indicate a phase in the direct solution.
Definition: Amesos2_TypeDecl.hpp:65
A templated MultiVector class adapter for Amesos2.
Definition: Amesos2_MultiVecAdapter_decl.hpp:176