22#ifdef DEAL_II_WITH_CUDA
32 template <
typename Number>
39 threadIdx.x + blockIdx.x * blockDim.x;
50 const float * A_val_dev,
51 const int * A_row_ptr_dev,
52 const int * A_column_index_dev,
53 cusparseSpMatDescr_t &sp_descr)
55 cusparseStatus_t error_code = cusparseCreateCsr(
60 reinterpret_cast<void *
>(
const_cast<int *
>(A_row_ptr_dev)),
61 reinterpret_cast<void *
>(
const_cast<int *
>(A_column_index_dev)),
62 reinterpret_cast<void *
>(
const_cast<float *
>(A_val_dev)),
65 CUSPARSE_INDEX_BASE_ZERO,
76 const double * A_val_dev,
77 const int * A_row_ptr_dev,
78 const int * A_column_index_dev,
79 cusparseSpMatDescr_t &sp_descr)
81 cusparseStatus_t error_code = cusparseCreateCsr(
86 reinterpret_cast<void *
>(
const_cast<int *
>(A_row_ptr_dev)),
87 reinterpret_cast<void *
>(
const_cast<int *
>(A_column_index_dev)),
88 reinterpret_cast<void *
>(
const_cast<double *
>(A_val_dev)),
91 CUSPARSE_INDEX_BASE_ZERO,
103 const cusparseSpMatDescr_t sp_descr,
109 float beta = add ? 1. : 0.;
110 cusparseOperation_t cusparse_operation =
111 transpose ? CUSPARSE_OPERATION_TRANSPOSE :
112 CUSPARSE_OPERATION_NON_TRANSPOSE;
115 cusparseDnVecDescr_t x_cuvec;
116 cusparseStatus_t error_code =
117 cusparseCreateDnVec(&x_cuvec,
119 reinterpret_cast<void *
>(
const_cast<float *
>(x)),
123 cusparseDnVecDescr_t y_cuvec;
125 cusparseCreateDnVec(&y_cuvec,
127 reinterpret_cast<void *
>(
const_cast<float *
>(y)),
132 size_t buffer_size = 0;
133 error_code = cusparseSpMV_bufferSize(handle,
141 CUSPARSE_MV_ALG_DEFAULT,
144 void * buffer =
nullptr;
145 cudaError_t cuda_error_code = cudaMalloc(&buffer, buffer_size);
149 error_code = cusparseSpMV(handle,
157 CUSPARSE_MV_ALG_DEFAULT,
161 cuda_error_code = cudaFree(buffer);
163 error_code = cusparseDestroyDnVec(x_cuvec);
165 error_code = cusparseDestroyDnVec(y_cuvec);
176 const cusparseSpMatDescr_t sp_descr,
182 double beta = add ? 1. : 0.;
183 cusparseOperation_t cusparse_operation =
184 transpose ? CUSPARSE_OPERATION_TRANSPOSE :
185 CUSPARSE_OPERATION_NON_TRANSPOSE;
188 cusparseDnVecDescr_t x_cuvec;
189 cusparseStatus_t error_code =
190 cusparseCreateDnVec(&x_cuvec,
192 reinterpret_cast<void *
>(
const_cast<double *
>(x)),
196 cusparseDnVecDescr_t y_cuvec;
198 cusparseCreateDnVec(&y_cuvec,
200 reinterpret_cast<void *
>(
const_cast<double *
>(y)),
205 size_t buffer_size = 0;
206 error_code = cusparseSpMV_bufferSize(handle,
214 CUSPARSE_MV_ALG_DEFAULT,
217 void * buffer =
nullptr;
218 cudaError_t cuda_error_code = cudaMalloc(&buffer, buffer_size);
222 error_code = cusparseSpMV(handle,
230 CUSPARSE_MV_ALG_DEFAULT,
234 cuda_error_code = cudaFree(buffer);
236 error_code = cusparseDestroyDnVec(x_cuvec);
238 error_code = cusparseDestroyDnVec(y_cuvec);
244 template <
typename Number>
247 const Number * val_dev,
248 const int * column_index_dev,
249 const int * row_ptr_dev,
253 threadIdx.x + blockIdx.x * blockDim.x;
257 for (
int j = row_ptr_dev[row]; j < row_ptr_dev[row + 1]; ++j)
258 atomicAdd(&sums[column_index_dev[j]],
std::abs(val_dev[j]));
264 template <
typename Number>
267 const Number * val_dev,
269 const int *row_ptr_dev,
273 threadIdx.x + blockIdx.x * blockDim.x;
277 sums[row] = (Number)0.;
278 for (
int j = row_ptr_dev[row]; j < row_ptr_dev[row + 1]; ++j)
286 template <
typename Number>
290 , val_dev(nullptr,
Utilities::CUDA::delete_device_data<Number>)
291 , column_index_dev(nullptr,
Utilities::CUDA::delete_device_data<
int>)
292 , row_ptr_dev(nullptr,
Utilities::CUDA::delete_device_data<
int>)
299 template <
typename Number>
302 const ::SparseMatrix<Number> &sparse_matrix_host)
303 : val_dev(nullptr,
Utilities::CUDA::delete_device_data<Number>)
304 , column_index_dev(nullptr,
Utilities::CUDA::delete_device_data<
int>)
305 , row_ptr_dev(nullptr,
Utilities::CUDA::delete_device_data<
int>)
309 reinit(handle, sparse_matrix_host);
314 template <
typename Number>
316 : cusparse_handle(other.cusparse_handle)
318 , n_rows(other.n_rows)
319 , n_cols(other.n_cols)
320 , val_dev(
std::move(other.val_dev))
321 , column_index_dev(
std::move(other.column_index_dev))
322 , row_ptr_dev(
std::move(other.row_ptr_dev))
324 , sp_descr(other.sp_descr)
329 other.descr =
nullptr;
330 other.sp_descr =
nullptr;
335 template <
typename Number>
338 if (descr !=
nullptr)
340 const cusparseStatus_t cusparse_error_code =
341 cusparseDestroyMatDescr(descr);
346 if (sp_descr !=
nullptr)
348 const cusparseStatus_t cusparse_error_code =
349 cusparseDestroySpMat(sp_descr);
360 template <
typename Number>
366 n_rows = other.n_rows;
367 n_cols = other.n_cols;
368 val_dev = std::move(other.val_dev);
369 column_index_dev = std::move(other.column_index_dev);
370 row_ptr_dev = std::move(other.row_ptr_dev);
372 sp_descr = other.sp_descr;
377 other.descr =
nullptr;
378 other.sp_descr =
nullptr;
385 template <
typename Number>
389 const ::SparseMatrix<Number> &sparse_matrix_host)
392 nnz = sparse_matrix_host.n_nonzero_elements();
393 n_rows = sparse_matrix_host.m();
394 n_cols = sparse_matrix_host.n();
395 unsigned int const row_ptr_size = n_rows + 1;
396 std::vector<Number> val;
398 std::vector<int> column_index;
399 column_index.reserve(nnz);
400 std::vector<int> row_ptr(row_ptr_size, 0);
404 for (
int row = 0; row < n_rows; ++row)
406 auto p_end = sparse_matrix_host.end(row);
407 unsigned int counter = 0;
408 for (
auto p = sparse_matrix_host.begin(row); p != p_end; ++p)
410 val.emplace_back(p->value());
411 column_index.emplace_back(p->column());
414 row_ptr[row + 1] = row_ptr[row] + counter;
417 unsigned int const offset = row_ptr[row];
418 int const diag_index = column_index[offset];
419 Number diag_elem = sparse_matrix_host.diag_element(row);
420 unsigned int pos = 1;
421 while ((column_index[offset + pos] < row) && (pos < counter))
423 val[offset + pos - 1] = val[offset + pos];
424 column_index[offset + pos - 1] = column_index[offset + pos];
427 val[offset + pos - 1] = diag_elem;
428 column_index[offset + pos - 1] = diag_index;
432 val_dev.reset(Utilities::CUDA::allocate_device_data<Number>(nnz));
433 cudaError_t error_code = cudaMemcpy(val_dev.get(),
435 nnz *
sizeof(Number),
436 cudaMemcpyHostToDevice);
440 column_index_dev.reset(Utilities::CUDA::allocate_device_data<int>(nnz));
442 error_code = cudaMemcpy(column_index_dev.get(),
445 cudaMemcpyHostToDevice);
449 row_ptr_dev.reset(Utilities::CUDA::allocate_device_data<int>(row_ptr_size));
451 error_code = cudaMemcpy(row_ptr_dev.get(),
453 row_ptr_size *
sizeof(
int),
454 cudaMemcpyHostToDevice);
458 cusparseStatus_t cusparse_error_code = cusparseCreateMatDescr(&descr);
460 cusparse_error_code =
461 cusparseSetMatType(descr, CUSPARSE_MATRIX_TYPE_GENERAL);
463 cusparse_error_code =
464 cusparseSetMatIndexBase(descr, CUSPARSE_INDEX_BASE_ZERO);
473 column_index_dev.get(),
479 template <
typename Number>
484 const int n_blocks = 1 + (nnz - 1) /
block_size;
485 internal::scale<Number>
486 <<<n_blocks,
block_size>>>(val_dev.get(), factor, nnz);
494 template <
typename Number>
500 const int n_blocks = 1 + (nnz - 1) /
block_size;
501 internal::scale<Number>
502 <<<n_blocks,
block_size>>>(val_dev.get(), 1. / factor, nnz);
510 template <
typename Number>
528 template <
typename Number>
546 template <
typename Number>
564 template <
typename Number>
582 template <
typename Number>
595 template <
typename Number>
609 template <
typename Number>
617 dst.
sadd(-1., 1., b);
624 template <
typename Number>
629 const int n_blocks = 1 + (nnz - 1) /
block_size;
630 internal::l1_norm<Number>
633 column_index_dev.get(),
643 template <
typename Number>
648 const int n_blocks = 1 + (nnz - 1) /
block_size;
649 internal::linfty_norm<Number>
652 column_index_dev.get(),
662 template <
typename Number>
667 cudaError_t cuda_error = cudaMemcpy(matrix_values.
get_values(),
669 nnz *
sizeof(Number),
670 cudaMemcpyDeviceToDevice);
673 return matrix_values.
l2_norm();
678 template <
typename Number>
679 std::tuple<Number *, int *, int *, cusparseMatDescr_t, cusparseSpMatDescr_t>
682 return std::make_tuple(val_dev.get(),
683 column_index_dev.get(),
void vmult_add(LinearAlgebra::CUDAWrappers::Vector< Number > &dst, const LinearAlgebra::CUDAWrappers::Vector< Number > &src) const
std::tuple< Number *, int *, int *, cusparseMatDescr_t, cusparseSpMatDescr_t > get_cusparse_matrix() const
void Tvmult_add(LinearAlgebra::CUDAWrappers::Vector< Number > &dst, const LinearAlgebra::CUDAWrappers::Vector< Number > &src) const
cusparseHandle_t cusparse_handle
void Tvmult(LinearAlgebra::CUDAWrappers::Vector< Number > &dst, const LinearAlgebra::CUDAWrappers::Vector< Number > &src) const
void reinit(Utilities::CUDA::Handle &handle, const ::SparseMatrix< Number > &sparse_matrix_host)
Number frobenius_norm() const
Number matrix_norm_square(const LinearAlgebra::CUDAWrappers::Vector< Number > &v) const
SparseMatrix & operator*=(const Number factor)
Number residual(LinearAlgebra::CUDAWrappers::Vector< Number > &dst, const LinearAlgebra::CUDAWrappers::Vector< Number > &x, const LinearAlgebra::CUDAWrappers::Vector< Number > &b) const
SparseMatrix & operator/=(const Number factor)
SparseMatrix & operator=(CUDAWrappers::SparseMatrix< Number > &&)
void vmult(LinearAlgebra::CUDAWrappers::Vector< Number > &dst, const LinearAlgebra::CUDAWrappers::Vector< Number > &src) const
Number linfty_norm() const
Number matrix_scalar_product(const LinearAlgebra::CUDAWrappers::Vector< Number > &u, const LinearAlgebra::CUDAWrappers::Vector< Number > &v) const
virtual void sadd(const Number s, const Number a, const VectorSpaceVector< Number > &V) override
virtual real_type l2_norm() const override
Number * get_values() const
virtual real_type linfty_norm() const override
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#define AssertCusparse(error_code)
#define AssertCudaKernel()
static ::ExceptionBase & ExcZero()
#define Assert(cond, exc)
#define AssertIsFinite(number)
#define AssertNothrowCusparse(error_code)
#define AssertCuda(error_code)
void create_sp_mat_descr(int m, int n, int nnz, const float *A_val_dev, const int *A_row_ptr_dev, const int *A_column_index_dev, cusparseSpMatDescr_t &sp_descr)
__global__ void linfty_norm(const typename SparseMatrix< Number >::size_type n_rows, const Number *val_dev, const int *, const int *row_ptr_dev, Number *sums)
__global__ void l1_norm(const typename SparseMatrix< Number >::size_type n_rows, const Number *val_dev, const int *column_index_dev, const int *row_ptr_dev, Number *sums)
__global__ void scale(Number *val, const Number a, const typename SparseMatrix< Number >::size_type N)
void csrmv(cusparseHandle_t handle, bool transpose, int m, int n, const cusparseSpMatDescr_t sp_descr, const float *x, bool add, float *y)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
cusparseHandle_t cusparse_handle