Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
cuda_solver_direct.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2018 - 2022 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
17
19
20namespace CUDAWrappers
21{
22 namespace
23 {
24 void
25 cusparsecsr2dense(cusparseHandle_t cusparse_handle,
26 const SparseMatrix<float> &matrix,
27 float * dense_matrix_dev)
28 {
29 auto cusparse_matrix = matrix.get_cusparse_matrix();
30
31 const cusparseStatus_t cusparse_error_code =
32 cusparseScsr2dense(cusparse_handle,
33 matrix.m(),
34 matrix.n(),
35 std::get<3>(cusparse_matrix),
36 std::get<0>(cusparse_matrix),
37 std::get<2>(cusparse_matrix),
38 std::get<1>(cusparse_matrix),
39 dense_matrix_dev,
40 matrix.m());
41 AssertCusparse(cusparse_error_code);
42 }
43
44
45
46 void
47 cusparsecsr2dense(cusparseHandle_t cusparse_handle,
48 const SparseMatrix<double> &matrix,
49 double * dense_matrix_dev)
50 {
51 auto cusparse_matrix = matrix.get_cusparse_matrix();
52
53 const cusparseStatus_t cusparse_error_code =
54 cusparseDcsr2dense(cusparse_handle,
55 matrix.m(),
56 matrix.n(),
57 std::get<3>(cusparse_matrix),
58 std::get<0>(cusparse_matrix),
59 std::get<2>(cusparse_matrix),
60 std::get<1>(cusparse_matrix),
61 dense_matrix_dev,
62 matrix.m());
63 AssertCusparse(cusparse_error_code);
64 }
65
66
67
68 void
69 cusolverDngetrf_buffer_size(cusolverDnHandle_t cusolver_dn_handle,
70 int m,
71 int n,
72 float * dense_matrix_dev,
73 int & workspace_size)
74 {
75 const cusolverStatus_t cusolver_error_code = cusolverDnSgetrf_bufferSize(
76 cusolver_dn_handle, m, n, dense_matrix_dev, m, &workspace_size);
77 AssertCusolver(cusolver_error_code);
78 }
79
80
81
82 void
83 cusolverDngetrf_buffer_size(cusolverDnHandle_t cusolver_dn_handle,
84 int m,
85 int n,
86 double * dense_matrix_dev,
87 int & workspace_size)
88 {
89 const cusolverStatus_t cusolver_error_code = cusolverDnDgetrf_bufferSize(
90 cusolver_dn_handle, m, n, dense_matrix_dev, m, &workspace_size);
91 AssertCusolver(cusolver_error_code);
92 }
93
94
95
96 void
97 cusolverDngetrf(cusolverDnHandle_t cusolver_dn_handle,
98 int m,
99 int n,
100 float * dense_matrix_dev,
101 float * workspace_dev,
102 int * pivot_dev,
103 int * info_dev)
104 {
105 const cusolverStatus_t cusolver_error_code =
106 cusolverDnSgetrf(cusolver_dn_handle,
107 m,
108 n,
109 dense_matrix_dev,
110 m,
111 workspace_dev,
112 pivot_dev,
113 info_dev);
114 AssertCusolver(cusolver_error_code);
115 }
116
117
118
119 void
120 cusolverDngetrf(cusolverDnHandle_t cusolver_dn_handle,
121 int m,
122 int n,
123 double * dense_matrix_dev,
124 double * workspace_dev,
125 int * pivot_dev,
126 int * info_dev)
127 {
128 const cusolverStatus_t cusolver_error_code =
129 cusolverDnDgetrf(cusolver_dn_handle,
130 m,
131 n,
132 dense_matrix_dev,
133 m,
134 workspace_dev,
135 pivot_dev,
136 info_dev);
137 AssertCusolver(cusolver_error_code);
138 }
139
140
141
142 void
143 cusolverDngetrs(cusolverDnHandle_t cusolver_dn_handle,
144 int m,
145 float * dense_matrix_dev,
146 int * pivot_dev,
147 float * b,
148 int * info_dev)
149 {
150 const int n_rhs = 1;
151 const cusolverStatus_t cusolver_error_code =
152 cusolverDnSgetrs(cusolver_dn_handle,
153 CUBLAS_OP_N,
154 m,
155 n_rhs,
156 dense_matrix_dev,
157 m,
158 pivot_dev,
159 b,
160 m,
161 info_dev);
162 AssertCusolver(cusolver_error_code);
163 }
164
165
166
167 void
168 cusolverDngetrs(cusolverDnHandle_t cusolver_dn_handle,
169 int m,
170 double * dense_matrix_dev,
171 int * pivot_dev,
172 double * b,
173 int * info_dev)
174 {
175 const int n_rhs = 1;
176 const cusolverStatus_t cusolver_error_code =
177 cusolverDnDgetrs(cusolver_dn_handle,
178 CUBLAS_OP_N,
179 m,
180 n_rhs,
181 dense_matrix_dev,
182 m,
183 pivot_dev,
184 b,
185 m,
186 info_dev);
187 AssertCusolver(cusolver_error_code);
188 }
189
190
191
192 void
193 cusolverSpcsrlsvluHost(cusolverSpHandle_t cusolver_sp_handle,
194 const unsigned int n_rows,
195 const unsigned int nnz,
196 cusparseMatDescr_t descr,
197 const float * val_host,
198 const int * row_ptr_host,
199 const int * column_index_host,
200 const float * b_host,
201 float * x_host)
202 {
203 int singularity = 0;
204 const cusolverStatus_t cusolver_error_code =
205 cusolverSpScsrlsvluHost(cusolver_sp_handle,
206 n_rows,
207 nnz,
208 descr,
209 val_host,
210 row_ptr_host,
211 column_index_host,
212 b_host,
213 0.,
214 1,
215 x_host,
216 &singularity);
217 AssertCusolver(cusolver_error_code);
218 Assert(singularity == -1, ExcMessage("Coarse matrix is singular"));
219 }
220
221
222
223 void
224 cusolverSpcsrlsvluHost(cusolverSpHandle_t cusolver_sp_handle,
225 const unsigned int n_rows,
226 unsigned int nnz,
227 cusparseMatDescr_t descr,
228 const double * val_host,
229 const int * row_ptr_host,
230 const int * column_index_host,
231 const double * b_host,
232 double * x_host)
233 {
234 int singularity = 0;
235 const cusolverStatus_t cusolver_error_code =
236 cusolverSpDcsrlsvluHost(cusolver_sp_handle,
237 n_rows,
238 nnz,
239 descr,
240 val_host,
241 row_ptr_host,
242 column_index_host,
243 b_host,
244 0.,
245 1,
246 x_host,
247 &singularity);
248 AssertCusolver(cusolver_error_code);
249 Assert(singularity == -1, ExcMessage("Coarse matrix is singular"));
250 }
251
252
253
254 void
255 cholesky_factorization(cusolverSpHandle_t cusolver_sp_handle,
256 const SparseMatrix<float> &matrix,
257 const float * b,
258 float * x)
259 {
260 auto cusparse_matrix = matrix.get_cusparse_matrix();
261 int singularity = 0;
262
263 const cusolverStatus_t cusolver_error_code =
264 cusolverSpScsrlsvchol(cusolver_sp_handle,
265 matrix.m(),
266 matrix.n_nonzero_elements(),
267 std::get<3>(cusparse_matrix),
268 std::get<0>(cusparse_matrix),
269 std::get<2>(cusparse_matrix),
270 std::get<1>(cusparse_matrix),
271 b,
272 0.,
273 0,
274 x,
275 &singularity);
276 AssertCusolver(cusolver_error_code);
277 Assert(singularity == -1, ExcMessage("Coarse matrix is not SPD"));
278 }
279
280
281
282 void
283 cholesky_factorization(cusolverSpHandle_t cusolver_sp_handle,
284 const SparseMatrix<double> &matrix,
285 const double * b,
286 double * x)
287 {
288 auto cusparse_matrix = matrix.get_cusparse_matrix();
289 int singularity = 0;
290
291 const cusolverStatus_t cusolver_error_code =
292 cusolverSpDcsrlsvchol(cusolver_sp_handle,
293 matrix.m(),
294 matrix.n_nonzero_elements(),
295 std::get<3>(cusparse_matrix),
296 std::get<0>(cusparse_matrix),
297 std::get<2>(cusparse_matrix),
298 std::get<1>(cusparse_matrix),
299 b,
300 0.,
301 0,
302 x,
303 &singularity);
304 AssertCusolver(cusolver_error_code);
305 Assert(singularity == -1, ExcMessage("Coarse matrix is not SPD"));
306 }
307
308
309
310 template <typename Number>
311 void
312 lu_factorization(cusparseHandle_t cusparse_handle,
313 cusolverDnHandle_t cusolver_dn_handle,
314 const SparseMatrix<Number> &matrix,
315 const Number * b_dev,
316 Number * x_dev)
317 {
318 // Change the format of the matrix from sparse to dense
319 unsigned int const m = matrix.m();
320 unsigned int const n = matrix.n();
321 Assert(m == n, ExcMessage("The matrix is not square"));
322 Number *dense_matrix_dev;
323 Utilities::CUDA::malloc(dense_matrix_dev, m * n);
324
325 // Change the format of matrix to dense
326 cusparsecsr2dense(cusparse_handle, matrix, dense_matrix_dev);
327
328 // Create the working space
329 int workspace_size = 0;
330 cusolverDngetrf_buffer_size(
331 cusolver_dn_handle, m, n, dense_matrix_dev, workspace_size);
332 Assert(workspace_size > 0, ExcMessage("No workspace was allocated"));
333 Number *workspace_dev;
334 Utilities::CUDA::malloc(workspace_dev, workspace_size);
335
336 // LU factorization
337 int *pivot_dev;
338 Utilities::CUDA::malloc(pivot_dev, m);
339 int *info_dev;
340 Utilities::CUDA::malloc(info_dev, 1);
341
342 cusolverDngetrf(cusolver_dn_handle,
343 m,
344 n,
345 dense_matrix_dev,
346 workspace_dev,
347 pivot_dev,
348 info_dev);
349
350#ifdef DEBUG
351 int info = 0;
352 cudaError_t cuda_error_code_debug =
353 cudaMemcpy(&info, info_dev, sizeof(int), cudaMemcpyDeviceToHost);
354 AssertCuda(cuda_error_code_debug);
355 Assert(info == 0,
356 ExcMessage("There was a problem during the LU factorization"));
357#endif
358
359 // Solve Ax = b
360 cudaError_t cuda_error_code =
361 cudaMemcpy(x_dev, b_dev, m * sizeof(Number), cudaMemcpyDeviceToDevice);
362 AssertCuda(cuda_error_code);
363 cusolverDngetrs(
364 cusolver_dn_handle, m, dense_matrix_dev, pivot_dev, x_dev, info_dev);
365#ifdef DEBUG
366 cuda_error_code =
367 cudaMemcpy(&info, info_dev, sizeof(int), cudaMemcpyDeviceToHost);
368 AssertCuda(cuda_error_code);
369 Assert(info == 0, ExcMessage("There was a problem during the LU solve"));
370#endif
371
372 // Free the memory allocated
373 Utilities::CUDA::free(dense_matrix_dev);
374 Utilities::CUDA::free(workspace_dev);
375 Utilities::CUDA::free(pivot_dev);
376 Utilities::CUDA::free(info_dev);
377 }
378
379
380
381 template <typename Number>
382 void
383 lu_factorization(cusolverSpHandle_t cusolver_sp_handle,
384 const SparseMatrix<Number> &matrix,
385 const Number * b_dev,
386 Number * x_dev)
387 {
388 // cuSOLVER does not support LU factorization of sparse matrix on the
389 // device, so we need to move everything to the host first and then back
390 // to the host.
391 const unsigned int nnz = matrix.n_nonzero_elements();
392 const unsigned int n_rows = matrix.m();
393 std::vector<Number> val_host(nnz);
394 std::vector<int> column_index_host(nnz);
395 std::vector<int> row_ptr_host(n_rows + 1);
396 auto cusparse_matrix = matrix.get_cusparse_matrix();
397 Utilities::CUDA::copy_to_host(std::get<0>(cusparse_matrix), val_host);
398 Utilities::CUDA::copy_to_host(std::get<1>(cusparse_matrix),
399 column_index_host);
400 Utilities::CUDA::copy_to_host(std::get<2>(cusparse_matrix), row_ptr_host);
401 std::vector<Number> b_host(n_rows);
402 Utilities::CUDA::copy_to_host(b_dev, b_host);
403 std::vector<Number> x_host(n_rows);
404 Utilities::CUDA::copy_to_host(x_dev, x_host);
405
406 cusolverSpcsrlsvluHost(cusolver_sp_handle,
407 n_rows,
408 nnz,
409 std::get<3>(cusparse_matrix),
410 val_host.data(),
411 row_ptr_host.data(),
412 column_index_host.data(),
413 b_host.data(),
414 x_host.data());
415
416 // Move the solution back to the device
417 Utilities::CUDA::copy_to_dev(x_host, x_dev);
418 }
419 } // namespace
420
421
422
423 template <typename Number>
425 const std::string &solver_type)
426 : solver_type(solver_type)
427 {}
428
429
430
431 template <typename Number>
433 SolverControl & cn,
434 const AdditionalData & data)
435 : cuda_handle(handle)
436 , solver_control(cn)
437 , additional_data(data.solver_type)
438 {}
439
440
441
442 template <typename Number>
445 {
446 return solver_control;
447 }
448
449
450
451 template <typename Number>
452 void
454 const SparseMatrix<Number> & A,
457 {
458 if (additional_data.solver_type == "Cholesky")
459 cholesky_factorization(cuda_handle.cusolver_sp_handle,
460 A,
461 b.get_values(),
462 x.get_values());
463 else if (additional_data.solver_type == "LU_dense")
464 lu_factorization(cuda_handle.cusparse_handle,
465 cuda_handle.cusolver_dn_handle,
466 A,
467 b.get_values(),
468 x.get_values());
469 else if (additional_data.solver_type == "LU_host")
470 lu_factorization(cuda_handle.cusolver_sp_handle,
471 A,
472 b.get_values(),
473 x.get_values());
474 else
475 AssertThrow(false,
476 ExcMessage("The provided solver name " +
477 additional_data.solver_type + " is invalid."));
478
479 // Force the SolverControl object to report convergence
480 solver_control.check(0, 0);
481 }
482
483
484 // Explicit Instanationation
485 template class SolverDirect<float>;
486 template class SolverDirect<double>;
487} // namespace CUDAWrappers
488
const Utilities::CUDA::Handle & cuda_handle
SolverDirect(const Utilities::CUDA::Handle &handle, SolverControl &cn, const AdditionalData &data=AdditionalData())
const AdditionalData additional_data
void solve(const SparseMatrix< Number > &A, LinearAlgebra::CUDAWrappers::Vector< Number > &x, const LinearAlgebra::CUDAWrappers::Vector< Number > &b)
SolverControl & control() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define AssertCusparse(error_code)
#define Assert(cond, exc)
#define AssertCusolver(error_code)
#define AssertCuda(error_code)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
@ matrix
Contents is actually a matrix.
void malloc(T *&pointer, const unsigned int n_elements)
Definition cuda.h:85
void copy_to_host(const ArrayView< const T, MemorySpace::CUDA > &in, ArrayView< T, MemorySpace::Host > &out)
Definition cuda.h:132
void copy_to_dev(const ArrayView< const T, MemorySpace::Host > &in, ArrayView< T, MemorySpace::CUDA > &out)
Definition cuda.h:148
void free(T *&pointer)
Definition cuda.h:97
AdditionalData(const std::string &solver_type="LU_dense")