Reference documentation for deal.II version GIT edc7d5c3ce 2023-09-25 07:10: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\}}\)
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 
20 namespace CUDAWrappers
21 {
22  namespace
23  {
24  void
25  cusparsecsr2dense(cusparseHandle_t cusparse_handle,
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,
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,
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,
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,
315  const Number *b_dev,
316  Number *x_dev)
317  {
318  // Change the format of the matrix from sparse to dense
319  const unsigned int m = matrix.m();
320  const unsigned int 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,
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>
443  SolverControl &
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:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
#define AssertCusparse(error_code)
Definition: exceptions.h:2027
#define Assert(cond, exc)
Definition: exceptions.h:1616
#define AssertCusolver(error_code)
Definition: exceptions.h:2088
#define AssertCuda(error_code)
Definition: exceptions.h:1942
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
@ matrix
Contents is actually a matrix.
static const char A
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
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")