]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add wrappers for cuSPARSE csr sparse matrix
authorBruno Turcksin <bruno.turcksin@gmail.com>
Mon, 19 Feb 2018 01:56:13 +0000 (20:56 -0500)
committerBruno Turcksin <bruno.turcksin@gmail.com>
Wed, 28 Feb 2018 02:02:13 +0000 (21:02 -0500)
include/deal.II/lac/cuda_sparse_matrix.h [new file with mode: 0644]
source/lac/CMakeLists.txt
source/lac/cuda_sparse_matrix.cu [new file with mode: 0644]

diff --git a/include/deal.II/lac/cuda_sparse_matrix.h b/include/deal.II/lac/cuda_sparse_matrix.h
new file mode 100644 (file)
index 0000000..2d28bd6
--- /dev/null
@@ -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 <deal.II/base/config.h>
+#include <deal.II/base/subscriptor.h>
+
+#ifdef DEAL_II_WITH_CUDA
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/cuda_vector.h>
+#include <cusparse.h>
+
+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 <tt>@<float@> and
+   * @<double@></tt>.
+   *
+   * @ingroup Matrix1
+   * @author Bruno Turcksin
+   * @date 2018
+   */
+  template <typename Number>
+  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<Number> &sparse_matrix_host);
+
+    /**
+     * Move constructor. Create a new SparseMatrix by stealing the internal
+     * data.
+     */
+    SparseMatrix(CUDAWrappers::SparseMatrix<Number> &&);
+
+    /**
+     * Copy constructor is deleted.
+     */
+    SparseMatrix(const CUDAWrappers::SparseMatrix<Number> &) = 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<Number> &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<Number> &dst,
+               const LinearAlgebra::CUDAWrappers::Vector<Number> &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<Number> &dst,
+                const LinearAlgebra::CUDAWrappers::Vector<Number> &src) const;
+
+    /**
+     * Adding matrix-vector multiplication. Add $M \cdot src$ on $dst$
+     * with $M$ being this matrix.
+     */
+    void vmult_add(LinearAlgebra::CUDAWrappers::Vector<Number> &dst,
+                   const LinearAlgebra::CUDAWrappers::Vector<Number> &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<Number> &dst,
+                    const LinearAlgebra::CUDAWrappers::Vector<Number> &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<Number> &v) const;
+
+    /**
+     * Compute the matrix scalar product $\left(u,Mv\right)$.
+     */
+    Number matrix_scalar_product(const LinearAlgebra::CUDAWrappers::Vector<Number> &u,
+                                 const LinearAlgebra::CUDAWrappers::Vector<Number> &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<Number>       &dst,
+                    const LinearAlgebra::CUDAWrappers::Vector<Number> &x,
+                    const LinearAlgebra::CUDAWrappers::Vector<Number> &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<Number *, int *, int *, cusparseMatDescr_t>
+    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 <typename Number>
+  inline
+  unsigned int SparseMatrix<Number>::m() const
+  {
+    return n_rows;
+  }
+
+
+
+  template <typename Number>
+  inline
+  unsigned int SparseMatrix<Number>::n() const
+  {
+    return n_cols;
+  }
+
+
+
+  template <typename Number>
+  inline
+  std::size_t SparseMatrix<Number>::n_nonzero_elements() const
+  {
+    return nnz;
+  }
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
+#endif
index 8464b2426b3b1d2a1ea5d9d4bb81c620ecf019e8..dea7cf02f8a4a26ff91d1d15654df25a237b457d 100644 (file)
@@ -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 (file)
index 0000000..5008d6f
--- /dev/null
@@ -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 <deal.II/lac/cuda_sparse_matrix.h>
+#include <deal.II/base/cuda_size.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/lac/cuda_atomic.h>
+
+#ifdef DEAL_II_WITH_CUDA
+
+#include <cusparse.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace CUDAWrappers
+{
+  using ::dealii::CUDAWrappers::block_size;
+
+  namespace internal
+  {
+    template <typename Number>
+    __global__ void scale(Number                                         *val,
+                          const Number                                    a,
+                          const typename SparseMatrix<Number>::size_type  N)
+    {
+      const typename SparseMatrix<Number>::size_type idx = threadIdx.x +
+                                                           blockIdx.x * blockDim.x;
+      if (idx<N)
+        val[idx] *= a;
+    }
+
+
+
+    void csrmv(cusparseHandle_t handle, bool transpose, int m, int n, int nnz,
+               const cusparseMatDescr_t descr, const float *A_val_dev,
+               const int *A_row_ptr_dev, const int *A_column_index_dev,
+               const float *x, bool add,  float *y)
+    {
+      float alpha = 1.;
+      float beta = add ? 1. : 0.;
+      cusparseOperation_t cusparse_operation = transpose ?
+                                               CUSPARSE_OPERATION_TRANSPOSE :
+                                               CUSPARSE_OPERATION_NON_TRANSPOSE;
+
+      cusparseStatus_t error_code;
+      // This function performs y = alpha*op(A)*x + beta*y
+      error_code = cusparseScsrmv(handle, cusparse_operation, m, n, nnz,
+                                  &alpha, descr, A_val_dev, A_row_ptr_dev,
+                                  A_column_index_dev, x, &beta, y);
+      AssertCusparse(error_code);
+    }
+
+
+
+    void csrmv(cusparseHandle_t handle, bool transpose, int m, int n, int nnz,
+               const cusparseMatDescr_t descr, const double *A_val_dev,
+               const int *A_row_ptr_dev, const int *A_column_index_dev,
+               const double *x, bool add, double *y)
+    {
+      double alpha = 1.;
+      double beta = add ? 1. : 0.;
+      cusparseOperation_t cusparse_operation = transpose ?
+                                               CUSPARSE_OPERATION_TRANSPOSE :
+                                               CUSPARSE_OPERATION_NON_TRANSPOSE;
+
+      cusparseStatus_t error_code;
+      // This function performs y = alpha*op(A)*x + beta*y
+      error_code = cusparseDcsrmv(handle, cusparse_operation, m, n, nnz,
+                                  &alpha, descr, A_val_dev, A_row_ptr_dev,
+                                  A_column_index_dev, x, &beta, y);
+      AssertCusparse(error_code);
+    }
+
+
+
+    template <typename Number>
+    __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)
+    {
+      const typename SparseMatrix<Number>::size_type row = threadIdx.x +
+                                                           blockIdx.x * blockDim.x;
+
+      if (row<n_rows)
+        {
+          for (int j=row_ptr_dev[row]; j<row_ptr_dev[row+1] ; ++j)
+            dealii::LinearAlgebra::CUDAWrappers::atomicAdd_wrapper(
+              &sums[column_index_dev[j]], abs(val_dev[j]));
+        }
+    }
+
+
+
+    template <typename Number>
+    __global__ void linfty_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)
+    {
+      const typename SparseMatrix<Number>::size_type row = threadIdx.x +
+                                                           blockIdx.x * blockDim.x;
+
+      if (row<n_rows)
+        {
+          sums[row] = (Number) 0.;
+          for (int j=row_ptr_dev[row]; j<row_ptr_dev[row+1] ; ++j)
+            sums[row] += abs(val_dev[j]);
+        }
+    }
+  }
+
+
+
+  template <typename Number>
+  SparseMatrix<Number>::SparseMatrix()
+    :
+    nnz(0),
+    n_rows(0),
+    val_dev(nullptr),
+    column_index_dev(nullptr),
+    row_ptr_dev(nullptr),
+    descr(nullptr)
+  {}
+
+
+
+  template <typename Number>
+  SparseMatrix<Number>::SparseMatrix(cusparseHandle_t handle,
+                                     const ::dealii::SparseMatrix<Number> &sparse_matrix_host)
+    :
+    val_dev(nullptr),
+    column_index_dev(nullptr),
+    row_ptr_dev(nullptr),
+    descr(nullptr)
+  {
+    reinit(handle, sparse_matrix_host);
+  }
+
+
+
+  template <typename Number>
+  SparseMatrix<Number>::SparseMatrix(CUDAWrappers::SparseMatrix<Number> &&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 <typename Number>
+  SparseMatrix<Number>::~SparseMatrix<Number>()
+  {
+    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 <typename Number>
+  void SparseMatrix<Number>::reinit(cusparseHandle_t handle,
+                                    const ::dealii::SparseMatrix<Number> &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<Number> val;
+    val.reserve(nnz);
+    std::vector<int> column_index;
+    column_index.reserve(nnz);
+    std::vector<int> 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 <typename Number>
+  SparseMatrix<Number> &SparseMatrix<Number>::operator*= (const Number factor)
+  {
+    AssertIsFinite(factor);
+    const int n_blocks = 1 + (nnz-1)/block_size;
+    internal::scale<Number> <<<n_blocks,block_size>>>(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 <typename Number>
+  SparseMatrix<Number> &SparseMatrix<Number>::operator/= (const Number factor)
+  {
+    AssertIsFinite(factor);
+    Assert(factor!=Number(0.), ExcZero());
+    const int n_blocks = 1 + (nnz-1)/block_size;
+    internal::scale<Number> <<<n_blocks,block_size>>>(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 <typename Number>
+  void SparseMatrix<Number>::vmult(LinearAlgebra::CUDAWrappers::Vector<Number> &dst,
+                                   const
+                                   LinearAlgebra::CUDAWrappers::Vector<Number> &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 <typename Number>
+  void SparseMatrix<Number>::Tvmult(LinearAlgebra::CUDAWrappers::Vector<Number> &dst,
+                                    const LinearAlgebra::CUDAWrappers::Vector<Number> &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 <typename Number>
+  void SparseMatrix<Number>::vmult_add(LinearAlgebra::CUDAWrappers::Vector<Number> &dst,
+                                       const LinearAlgebra::CUDAWrappers::Vector<Number> &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 <typename Number>
+  void SparseMatrix<Number>::Tvmult_add(LinearAlgebra::CUDAWrappers::Vector<Number> &dst,
+                                        const LinearAlgebra::CUDAWrappers::Vector<Number> &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 <typename Number>
+  Number SparseMatrix<Number>::matrix_norm_square(const LinearAlgebra::CUDAWrappers::Vector<Number> &v) const
+  {
+    LinearAlgebra::CUDAWrappers::Vector<Number> tmp = v;
+    vmult(tmp, v);
+
+    return v*tmp;
+  }
+
+
+
+  template <typename Number>
+  Number SparseMatrix<Number>::matrix_scalar_product(const LinearAlgebra::CUDAWrappers::Vector<Number> &u,
+                                                     const LinearAlgebra::CUDAWrappers::Vector<Number> &v) const
+  {
+    LinearAlgebra::CUDAWrappers::Vector<Number> tmp = v;
+    vmult(tmp, v);
+
+    return u*tmp;
+  }
+
+
+
+  template <typename Number>
+  Number SparseMatrix<Number>::residual(LinearAlgebra::CUDAWrappers::Vector<Number>       &dst,
+                                        const LinearAlgebra::CUDAWrappers::Vector<Number> &x,
+                                        const LinearAlgebra::CUDAWrappers::Vector<Number> &b) const
+  {
+    vmult(dst, x);
+    dst.sadd(-1., 1., b);
+
+    return dst.l2_norm();
+  }
+
+
+
+  template <typename Number>
+  Number SparseMatrix<Number>::l1_norm() const
+  {
+    LinearAlgebra::CUDAWrappers::Vector<real_type> column_sums(n_cols);
+    const int n_blocks = 1 + (nnz-1)/block_size;
+    internal::l1_norm<Number> <<<n_blocks,block_size>>>(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 <typename Number>
+  Number SparseMatrix<Number>::linfty_norm() const
+  {
+    LinearAlgebra::CUDAWrappers::Vector<real_type> row_sums(n_rows);
+    const int n_blocks = 1 + (nnz-1)/block_size;
+    internal::linfty_norm<Number> <<<n_blocks,block_size>>>(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 <typename Number>
+  Number SparseMatrix<Number>::frobenius_norm() const
+  {
+    LinearAlgebra::CUDAWrappers::Vector<real_type> matrix_values(nnz);
+    cudaError_t cuda_error = cudaMemcpy(matrix_values.get_values(), val_dev,
+                                        nnz*sizeof(Number),
+                                        cudaMemcpyDeviceToDevice);
+
+    return matrix_values.l2_norm();
+  }
+
+
+
+  template <typename Number>
+  std::tuple<Number *, int *, int *, cusparseMatDescr_t>
+  SparseMatrix<Number>::get_cusparse_matrix()
+  {
+    return std::make_tuple(val_dev, column_index_dev, row_ptr_dev, descr);
+  }
+
+
+
+  template class SparseMatrix<float>;
+  template class SparseMatrix<double>;
+}
+DEAL_II_NAMESPACE_CLOSE
+
+#endif

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.