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