From: Bruno Turcksin Date: Mon, 19 Feb 2018 01:56:13 +0000 (-0500) Subject: Add wrappers for cuSPARSE csr sparse matrix X-Git-Tag: v9.0.0-rc1~377^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=988a56a8f3e24383e4fe4aab31eda72ae5bd73a4;p=dealii.git Add wrappers for cuSPARSE csr sparse matrix --- diff --git a/include/deal.II/lac/cuda_sparse_matrix.h b/include/deal.II/lac/cuda_sparse_matrix.h new file mode 100644 index 0000000000..2d28bd6dfd --- /dev/null +++ b/include/deal.II/lac/cuda_sparse_matrix.h @@ -0,0 +1,326 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#ifndef dealii_cuda_sparse_matrix_h +#define dealii_cuda_sparse_matrix_h + +#include +#include + +#ifdef DEAL_II_WITH_CUDA +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace CUDAWrappers +{ + /** + * This class is a wrapper around cuSPARSE csr sparse matrix. Unlike deal.II's + * own SparseMatrix all elements within each row are stored in increasing + * column index order. + * + * @note Instantiations for this template are provided for @ and + * @. + * + * @ingroup Matrix1 + * @author Bruno Turcksin + * @date 2018 + */ + template + class SparseMatrix: public virtual Subscriptor + { + public: + /** + * Declare type for container size. + */ + typedef unsigned int size_type; + + /** + * Type of the matrix entries. + */ + typedef Number value_type; + + /** + * Declare a type that holds real-valued numbers with the same precision + * as the template argument to this class. + */ + typedef Number real_type; + + /** + * @name Constructors and initialization + */ + //@{ + /** + * Constructor. Initialize the matrix to be empty, without any structure, + * i.e., the matrix is not usable at all. This constructor is therefore + * only useful for matrices which are members of a class. + * + * You have to initialize the matrix before usage with reinit. + */ + SparseMatrix(); + + /** + * Constructor. Takes a cuSPARSE handle and a sparse matrix on the host. + * The sparse matrix on the host is copied on the device and the elements + * are reordered according to the format supported by cuSPARSE. + */ + SparseMatrix(cusparseHandle_t handle, + const ::dealii::SparseMatrix &sparse_matrix_host); + + /** + * Move constructor. Create a new SparseMatrix by stealing the internal + * data. + */ + SparseMatrix(CUDAWrappers::SparseMatrix &&); + + /** + * Copy constructor is deleted. + */ + SparseMatrix(const CUDAWrappers::SparseMatrix &) = delete; + + /** + * Destructor. Free all memory. + */ + ~SparseMatrix(); + + /** + * Reinitialize the sparse matrix. The sparse matrix on the host is copied + * to the device and the elementes are reordered according to the format + * supported by cuSPARSE. + */ + void reinit(cusparseHandle_t handle, + const ::dealii::SparseMatrix &sparse_matrix_host); + //@} + + /** + * @name Information on the matrix + */ + //@{ + /** + * Return the dimension of the codomain (or range) space. Note that the + * matrix is of dimension $m \times n$. + */ + size_type m() const; + + /** + * Return the dimension of the domain space. Note that the matrix is of + * dimension $m \times n$. + */ + size_type n() const; + + /** + * Return the number of nonzero elements of this matrix. Actually, it + * returns the number of entries in the sparsity pattern; if any of the + * entries should happen to be zero, it is counted anyway. + */ + std::size_t n_nonzero_elements() const; + //@} + + /** + * @name Modifying entries + */ + //@{ + /** + * Multiply the entire matrix by a fixed factor. + */ + SparseMatrix &operator*= (const Number factor); + + /** + * Divide the entrie matrix by a fixed factor. + */ + SparseMatrix &operator/= (const Number factor); + //@} + + /** + * @name Multiplications + */ + //@{ + /** + * Matrix-vector multiplication: let $dst = M \cdot src$ with $M$ + * being this matrix. + */ + void vmult(LinearAlgebra::CUDAWrappers::Vector &dst, + const LinearAlgebra::CUDAWrappers::Vector &src) const; + + /** + * Matrix-vector multiplication: let $dst = M^T \cdot src$ with + * $M$ being this matrix. This function does the same as vmult() but + * takes thes transposed matrix. + */ + void Tvmult(LinearAlgebra::CUDAWrappers::Vector &dst, + const LinearAlgebra::CUDAWrappers::Vector &src) const; + + /** + * Adding matrix-vector multiplication. Add $M \cdot src$ on $dst$ + * with $M$ being this matrix. + */ + void vmult_add(LinearAlgebra::CUDAWrappers::Vector &dst, + const LinearAlgebra::CUDAWrappers::Vector &src) const; + + /** + * Adding matrix-vector multiplication. Add $M^T \cdot src$ to + * $dst$ with $M$ being this matrix. This function foes the same + * as vmult_add() but takes the transposed matrix. + */ + void Tvmult_add(LinearAlgebra::CUDAWrappers::Vector &dst, + const LinearAlgebra::CUDAWrappers::Vector &src) const; + + /** + * Return the square of the norm of the vector $v$ with respect to the the + * norm induced by this matrix, i.e., $\left(v,Mv\right)$. This is useful, + * e.g., in the finite context, where the $L_2$ norm of a function equals + * the matrix norm with respect to the mass matrix of the vector + * representing the nodal values of the finite element function. + * + * Obviously, the matrix needs to be quadratic for this operation. + */ + Number matrix_norm_square(const LinearAlgebra::CUDAWrappers::Vector &v) const; + + /** + * Compute the matrix scalar product $\left(u,Mv\right)$. + */ + Number matrix_scalar_product(const LinearAlgebra::CUDAWrappers::Vector &u, + const LinearAlgebra::CUDAWrappers::Vector &v) const; + + /** + * Compute the residual of an equation $M \cdot x=b$, where the residual is + * defined to be $r=b-M \cdot x$. Write the residual into $dst$. The + * $l_2$ norm of the residual vector is returned. + * + * Source $x$ and destination $dst$ must not be the same vector. + */ + Number residual(LinearAlgebra::CUDAWrappers::Vector &dst, + const LinearAlgebra::CUDAWrappers::Vector &x, + const LinearAlgebra::CUDAWrappers::Vector &b) const; + //@} + + /** + * @name Matrix norms + */ + //@{ + /** + * Return the $l_1$-norm of the matrix, that is $|M|_1=\max_{\mathrm{all\ + * columns\ }j}\sum_{\mathrm{all\ rows\ }i} |M_{ij}|$, (max. sum of + * columns). This is the natural matrix norm that is compatible to the + * $l_1$-norm for vectors, i.e., $|Mv|_1\leq |M|_1 |v|_1$. + */ + Number l1_norm() const; + + /** + * Return the $l_\infty$-norm of the matrix, that is + * $|M|_\infty=\max_{\mathrm{all\ rows\ }i}\sum_{\mathrm{all\ columns\ }j} + * |M_{ij}|$, (max. sum of rows). This is the natural norm that is + * compatible to the $l_\infty$-norm of vectors, i.e., $|Mv|_\infty \leq + * |M|_\infty |v|_\infty$. + */ + Number linfty_norm() const; + + /** + * Return the frobenius norm of the matrix, i.e., the square root of the + * sum of squares of all entries in the matrix. + */ + Number frobenius_norm() const; + //@} + + /** + *@name Access to underlying CUDA data + */ + //@{ + /** + * Return a tuple containing the pointer to the values of matrix, the + * pointer to the columns indices, the pointer to the rows pointer, and + * the cuSPARSE matrix description. + */ + std::tuple + get_cusparse_matrix(); + //*} + + private: + /** + * cuSPARSE used to call cuSPARSE function. The cuSPARSE handle needs to + * be mutable to be called in a const function. + */ + mutable cusparseHandle_t cusparse_handle; + + /** + * Number of non-zero elements in the sparse matrix. + */ + int nnz; + + /** + * Number of rows of the sparse matrix. + */ + int n_rows; + + /** + * Number of columns of the sparse matrix. + */ + int n_cols; + + /** + * Pointer to the values (on the device) of the sparse matrix. + */ + Number *val_dev; + + /** + * Pointer to the column indices (on the device) of the sparse matrix. + */ + int *column_index_dev; + + /** + * Pointer to the row pointer (on the device) of the sparse matrix. + */ + int *row_ptr_dev; + + /** + * cuSPARSE description of the sparse matrix. + */ + cusparseMatDescr_t descr; + }; + + + + template + inline + unsigned int SparseMatrix::m() const + { + return n_rows; + } + + + + template + inline + unsigned int SparseMatrix::n() const + { + return n_cols; + } + + + + template + inline + std::size_t SparseMatrix::n_nonzero_elements() const + { + return nnz; + } +} + +DEAL_II_NAMESPACE_CLOSE + +#endif +#endif diff --git a/source/lac/CMakeLists.txt b/source/lac/CMakeLists.txt index 8464b2426b..dea7cf02f8 100644 --- a/source/lac/CMakeLists.txt +++ b/source/lac/CMakeLists.txt @@ -147,6 +147,7 @@ IF(DEAL_II_WITH_CUDA) SET(_separate_src ${_separate_src} cuda_vector.cu + cuda_sparse_matrix.cu ) ENDIF() diff --git a/source/lac/cuda_sparse_matrix.cu b/source/lac/cuda_sparse_matrix.cu new file mode 100644 index 0000000000..5008d6f4ec --- /dev/null +++ b/source/lac/cuda_sparse_matrix.cu @@ -0,0 +1,469 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#include +#include +#include +#include + +#ifdef DEAL_II_WITH_CUDA + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace CUDAWrappers +{ + using ::dealii::CUDAWrappers::block_size; + + namespace internal + { + template + __global__ void scale(Number *val, + const Number a, + const typename SparseMatrix::size_type N) + { + const typename SparseMatrix::size_type idx = threadIdx.x + + blockIdx.x * blockDim.x; + if (idx + __global__ void l1_norm(const typename SparseMatrix::size_type n_rows, + const Number *val_dev, const int *column_index_dev, + const int *row_ptr_dev, Number *sums) + { + const typename SparseMatrix::size_type row = threadIdx.x + + blockIdx.x * blockDim.x; + + if (row + __global__ void linfty_norm(const typename SparseMatrix::size_type n_rows, + const Number *val_dev, const int *column_index_dev, + const int *row_ptr_dev, Number *sums) + { + const typename SparseMatrix::size_type row = threadIdx.x + + blockIdx.x * blockDim.x; + + if (row + SparseMatrix::SparseMatrix() + : + nnz(0), + n_rows(0), + val_dev(nullptr), + column_index_dev(nullptr), + row_ptr_dev(nullptr), + descr(nullptr) + {} + + + + template + SparseMatrix::SparseMatrix(cusparseHandle_t handle, + const ::dealii::SparseMatrix &sparse_matrix_host) + : + val_dev(nullptr), + column_index_dev(nullptr), + row_ptr_dev(nullptr), + descr(nullptr) + { + reinit(handle, sparse_matrix_host); + } + + + + template + SparseMatrix::SparseMatrix(CUDAWrappers::SparseMatrix &&other) + { + cusparse_handle = other.cusparse_handle; + nnz = other.nnz; + n_rows = other.n_rows; + n_cols = other.n_cols; + val_dev = other.val_dev; + column_index_dev = other.column_index_dev; + row_ptr_dev = other.row_ptr_dev; + descr = other.descr; + + other.nnz = 0; + other.n_rows = 0; + other.n_cols = 0; + other.val_dev = nullptr; + other.column_index_dev = nullptr; + other.row_ptr_dev = nullptr; + other.descr = nullptr; + } + + + + template + SparseMatrix::~SparseMatrix() + { + if (val_dev != nullptr) + { + cudaError_t error_code = cudaFree(val_dev); + AssertCuda(error_code); + val_dev = nullptr; + } + + if (column_index_dev != nullptr) + { + cudaError_t error_code = cudaFree(column_index_dev); + AssertCuda(error_code); + column_index_dev = nullptr; + } + + if (row_ptr_dev != nullptr) + { + cudaError_t error_code = cudaFree(row_ptr_dev); + AssertCuda(error_code); + row_ptr_dev = nullptr; + } + + if (descr != nullptr) + { + cusparseStatus_t cusparse_error_code = cusparseDestroyMatDescr(descr); + AssertCusparse(cusparse_error_code); + descr = nullptr; + } + + nnz = 0; + n_rows = 0; + } + + + + template + void SparseMatrix::reinit(cusparseHandle_t handle, + const ::dealii::SparseMatrix &sparse_matrix_host) + { + cusparse_handle = handle; + nnz = sparse_matrix_host.n_nonzero_elements(); + n_rows = sparse_matrix_host.m(); + n_cols = sparse_matrix_host.n(); + unsigned int const row_ptr_size = n_rows + 1; + std::vector val; + val.reserve(nnz); + std::vector column_index; + column_index.reserve(nnz); + std::vector row_ptr(row_ptr_size, 0); + + // dealii::SparseMatrix stores the diagonal first in each row so we need to do some + // reordering + for (int row = 0; row < n_rows; ++row) + { + auto p_end = sparse_matrix_host.end(row); + unsigned int counter = 0; + for (auto p = sparse_matrix_host.begin(row); p != p_end; ++p) + { + val.emplace_back(p->value()); + column_index.emplace_back(p->column()); + ++counter; + } + row_ptr[row + 1] = row_ptr[row] + counter; + + // Sort the elements in the row + unsigned int const offset = row_ptr[row]; + int const diag_index = column_index[offset]; + Number diag_elem = sparse_matrix_host.diag_element(row); + unsigned int pos = 1; + while ((column_index[offset + pos] < row) && (pos < counter)) + { + val[offset + pos - 1] = val[offset + pos]; + column_index[offset + pos - 1] = column_index[offset + pos]; + ++pos; + } + val[offset + pos - 1] = diag_elem; + column_index[offset + pos - 1] = diag_index; + } + + // Copy the elements to the gpu + cudaError_t error_code = cudaMalloc(&val_dev, nnz * sizeof(Number)); + AssertCuda(error_code); + error_code = cudaMemcpy(val_dev, &val[0], nnz * sizeof(Number), + cudaMemcpyHostToDevice); + AssertCuda(error_code); + + // Copy the column indices to the gpu + error_code = cudaMalloc(&column_index_dev, nnz * sizeof(int)); + AssertCuda(error_code); + error_code = cudaMemcpy(column_index_dev, &column_index[0], nnz * sizeof(int), + cudaMemcpyHostToDevice); + AssertCuda(error_code); + + // Copy the row pointer to the gpu + error_code = cudaMalloc(&row_ptr_dev, row_ptr_size * sizeof(int)); + AssertCuda(error_code); + error_code = cudaMemcpy(row_ptr_dev, &row_ptr[0], row_ptr_size * sizeof(int), + cudaMemcpyHostToDevice); + AssertCuda(error_code); + + // Create the matrix descriptor + cusparseStatus_t cusparse_error_code = cusparseCreateMatDescr(&descr); + AssertCusparse(cusparse_error_code); + cusparse_error_code = cusparseSetMatType(descr, CUSPARSE_MATRIX_TYPE_GENERAL); + AssertCusparse(cusparse_error_code); + cusparse_error_code = cusparseSetMatIndexBase(descr, CUSPARSE_INDEX_BASE_ZERO); + AssertCusparse(cusparse_error_code); + } + + + + template + SparseMatrix &SparseMatrix::operator*= (const Number factor) + { + AssertIsFinite(factor); + const int n_blocks = 1 + (nnz-1)/block_size; + internal::scale <<>>(val_dev, factor, nnz); + + // Check that the kernel was launched correctly + AssertCuda(cudaGetLastError()); + // Check that there was no problem during the execution of the kernel + AssertCuda(cudaDeviceSynchronize()); + + return *this; + } + + + + template + SparseMatrix &SparseMatrix::operator/= (const Number factor) + { + AssertIsFinite(factor); + Assert(factor!=Number(0.), ExcZero()); + const int n_blocks = 1 + (nnz-1)/block_size; + internal::scale <<>>(val_dev, 1./factor, nnz); + + // Check that the kernel was launched correctly + AssertCuda(cudaGetLastError()); + // Check that there was no problem during the execution of the kernel + AssertCuda(cudaDeviceSynchronize()); + + return *this; + } + + + + template + void SparseMatrix::vmult(LinearAlgebra::CUDAWrappers::Vector &dst, + const + LinearAlgebra::CUDAWrappers::Vector &src) const + { + internal::csrmv(cusparse_handle, false, n_rows, n_cols, nnz, descr, val_dev, + row_ptr_dev, column_index_dev, src.get_values(), false, + dst.get_values()); + } + + + + template + void SparseMatrix::Tvmult(LinearAlgebra::CUDAWrappers::Vector &dst, + const LinearAlgebra::CUDAWrappers::Vector &src) const + { + internal::csrmv(cusparse_handle, true, n_rows, n_cols, nnz, descr, val_dev, + row_ptr_dev, column_index_dev, src.get_values(), false, + dst.get_values()); + } + + + + template + void SparseMatrix::vmult_add(LinearAlgebra::CUDAWrappers::Vector &dst, + const LinearAlgebra::CUDAWrappers::Vector &src) const + { + internal::csrmv(cusparse_handle, false, n_rows, n_cols, nnz, descr, val_dev, + row_ptr_dev, column_index_dev, src.get_values(), true, + dst.get_values()); + } + + + + template + void SparseMatrix::Tvmult_add(LinearAlgebra::CUDAWrappers::Vector &dst, + const LinearAlgebra::CUDAWrappers::Vector &src) const + { + internal::csrmv(cusparse_handle, true, n_rows, n_cols, nnz, descr, val_dev, + row_ptr_dev, column_index_dev, src.get_values(), true, + dst.get_values()); + } + + + + template + Number SparseMatrix::matrix_norm_square(const LinearAlgebra::CUDAWrappers::Vector &v) const + { + LinearAlgebra::CUDAWrappers::Vector tmp = v; + vmult(tmp, v); + + return v*tmp; + } + + + + template + Number SparseMatrix::matrix_scalar_product(const LinearAlgebra::CUDAWrappers::Vector &u, + const LinearAlgebra::CUDAWrappers::Vector &v) const + { + LinearAlgebra::CUDAWrappers::Vector tmp = v; + vmult(tmp, v); + + return u*tmp; + } + + + + template + Number SparseMatrix::residual(LinearAlgebra::CUDAWrappers::Vector &dst, + const LinearAlgebra::CUDAWrappers::Vector &x, + const LinearAlgebra::CUDAWrappers::Vector &b) const + { + vmult(dst, x); + dst.sadd(-1., 1., b); + + return dst.l2_norm(); + } + + + + template + Number SparseMatrix::l1_norm() const + { + LinearAlgebra::CUDAWrappers::Vector column_sums(n_cols); + const int n_blocks = 1 + (nnz-1)/block_size; + internal::l1_norm <<>>(n_rows, val_dev, + column_index_dev, + row_ptr_dev, + column_sums.get_values()); + // Check that the kernel was launched correctly + AssertCuda(cudaGetLastError()); + // Check that there was no problem during the execution of the kernel + AssertCuda(cudaDeviceSynchronize()); + + return column_sums.linfty_norm(); + } + + + + template + Number SparseMatrix::linfty_norm() const + { + LinearAlgebra::CUDAWrappers::Vector row_sums(n_rows); + const int n_blocks = 1 + (nnz-1)/block_size; + internal::linfty_norm <<>>(n_rows, val_dev, + column_index_dev, + row_ptr_dev, + row_sums.get_values()); + // Check that the kernel was launched correctly + AssertCuda(cudaGetLastError()); + // Check that there was no problem during the execution of the kernel + AssertCuda(cudaDeviceSynchronize()); + + return row_sums.linfty_norm(); + } + + + + template + Number SparseMatrix::frobenius_norm() const + { + LinearAlgebra::CUDAWrappers::Vector matrix_values(nnz); + cudaError_t cuda_error = cudaMemcpy(matrix_values.get_values(), val_dev, + nnz*sizeof(Number), + cudaMemcpyDeviceToDevice); + + return matrix_values.l2_norm(); + } + + + + template + std::tuple + SparseMatrix::get_cusparse_matrix() + { + return std::make_tuple(val_dev, column_index_dev, row_ptr_dev, descr); + } + + + + template class SparseMatrix; + template class SparseMatrix; +} +DEAL_II_NAMESPACE_CLOSE + +#endif