cuda_sparse_matrix.h

// ---------------------------------------------------------------------
//
// 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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef dealii_cuda_sparse_matrix_h
#define dealii_cuda_sparse_matrix_h

#include <deal.II/base/config.h>

#include <deal.II/base/subscriptor.h>

#include <iomanip>

#ifdef DEAL_II_COMPILER_CUDA_AWARE
#  include <deal.II/base/cuda.h>

#  include <deal.II/lac/cuda_vector.h>
#  include <deal.II/lac/sparse_matrix.h>

#  include <cusparse.h>

DEAL_II_NAMESPACE_OPEN

namespace CUDAWrappers
{
  template <typename Number>
  class SparseMatrix : public virtual Subscriptor
  {
  public:
    using size_type = int;

    using value_type = Number;

    using real_type = Number;

    SparseMatrix();

    SparseMatrix(Utilities::CUDA::Handle &             handle,
                 const ::SparseMatrix<Number> &sparse_matrix_host);

    SparseMatrix(CUDAWrappers::SparseMatrix<Number> &&);

    SparseMatrix(const CUDAWrappers::SparseMatrix<Number> &) = delete;

    ~SparseMatrix();

    SparseMatrix &
    operator=(CUDAWrappers::SparseMatrix<Number> &&);

    SparseMatrix &
    operator=(const CUDAWrappers::SparseMatrix<Number> &) = delete;

    void
    reinit(Utilities::CUDA::Handle &             handle,
           const ::SparseMatrix<Number> &sparse_matrix_host);

    size_type
    m() const;

    size_type
    n() const;

    std::size_t
    n_nonzero_elements() const;

    template <class StreamType>
    void
    print(StreamType &out,
          const bool  across         = false,
          const bool  diagonal_first = true) const;

    void
    print_formatted(std::ostream &     out,
                    const unsigned int precision   = 3,
                    const bool         scientific  = true,
                    const unsigned int width       = 0,
                    const char *       zero_string = " ",
                    const double       denominator = 1.) const;

    SparseMatrix &
    operator*=(const Number factor);

    SparseMatrix &
    operator/=(const Number factor);

    void
    vmult(LinearAlgebra::CUDAWrappers::Vector<Number> &      dst,
          const LinearAlgebra::CUDAWrappers::Vector<Number> &src) const;

    void
    Tvmult(LinearAlgebra::CUDAWrappers::Vector<Number> &      dst,
           const LinearAlgebra::CUDAWrappers::Vector<Number> &src) const;

    void
    vmult_add(LinearAlgebra::CUDAWrappers::Vector<Number> &      dst,
              const LinearAlgebra::CUDAWrappers::Vector<Number> &src) const;

    void
    Tvmult_add(LinearAlgebra::CUDAWrappers::Vector<Number> &      dst,
               const LinearAlgebra::CUDAWrappers::Vector<Number> &src) const;

    Number
    matrix_norm_square(
      const LinearAlgebra::CUDAWrappers::Vector<Number> &v) const;

    Number
    matrix_scalar_product(
      const LinearAlgebra::CUDAWrappers::Vector<Number> &u,
      const LinearAlgebra::CUDAWrappers::Vector<Number> &v) const;

    Number
    residual(LinearAlgebra::CUDAWrappers::Vector<Number> &      dst,
             const LinearAlgebra::CUDAWrappers::Vector<Number> &x,
             const LinearAlgebra::CUDAWrappers::Vector<Number> &b) const;

    Number
    l1_norm() const;

    Number
    linfty_norm() const;

    Number
    frobenius_norm() const;

    std::tuple<Number *, int *, int *, cusparseMatDescr_t>
    get_cusparse_matrix() const;

  private:
    cusparseHandle_t cusparse_handle;

    int nnz;

    int n_rows;

    int n_cols;

    std::unique_ptr<Number[], void (*)(Number *)> val_dev;

    std::unique_ptr<int[], void (*)(int *)> column_index_dev;

    std::unique_ptr<int[], void (*)(int *)> row_ptr_dev;

    cusparseMatDescr_t descr;
  };



  template <typename Number>
  inline SparseMatrix<Number>::size_type
  SparseMatrix<Number>::m() const
  {
    return n_rows;
  }



  template <typename Number>
  inline SparseMatrix<Number>::size_type
  SparseMatrix<Number>::n() const
  {
    return n_cols;
  }



  template <typename Number>
  inline std::size_t
  SparseMatrix<Number>::n_nonzero_elements() const
  {
    return nnz;
  }



  template <typename Number>
  template <class StreamType>
  inline void
  SparseMatrix<Number>::print(StreamType &out,
                              const bool  across,
                              const bool  diagonal_first) const
  {
    Assert(column_index_dev != nullptr, ExcNotInitialized());
    Assert(val_dev != nullptr, ExcNotInitialized());
    Assert(row_ptr_dev != nullptr, ExcNotInitialized());

    std::vector<int>    rows(n_rows + 1);
    std::vector<int>    cols(nnz);
    std::vector<double> val(nnz);
    Utilities::CUDA::copy_to_host(row_ptr_dev.get(), rows);
    Utilities::CUDA::copy_to_host(column_index_dev.get(), cols);
    Utilities::CUDA::copy_to_host(val_dev.get(), val);

    bool   has_diagonal = false;
    Number diagonal     = Number();

    for (size_type i = 0; i < n_rows; ++i)
      {
        if (diagonal_first)
          {
            // find the diagonal and print if it exists
            for (size_type j = rows[i]; j < rows[i + 1] && cols[j] <= i; ++j)
              {
                if (i == cols[j])
                  {
                    diagonal     = val[j];
                    has_diagonal = true;
                    if (across)
                      out << ' ' << i << ',' << i << ':' << diagonal;
                    else
                      out << '(' << i << ',' << i << ") " << diagonal
                          << std::endl;
                    break;
                  }
              }
          }
        for (size_type j = rows[i]; j < rows[i + 1]; ++j)
          {
            if (has_diagonal && i == cols[j])
              continue;
            if (across)
              out << ' ' << i << ',' << cols[j] << ':' << val[j];
            else
              out << "(" << i << "," << cols[j] << ") " << val[j] << std::endl;
          }
      }
    if (across)
      out << std::endl;
  }



  template <typename Number>
  void
  SparseMatrix<Number>::print_formatted(std::ostream &     out,
                                        const unsigned int precision,
                                        const bool         scientific,
                                        const unsigned int width_,
                                        const char *       zero_string,
                                        const double       denominator) const
  {
    Assert(column_index_dev != nullptr, ExcNotInitialized());
    Assert(val_dev != nullptr, ExcNotInitialized());
    Assert(row_ptr_dev != nullptr, ExcNotInitialized());

    std::vector<int>    rows(n_rows + 1);
    std::vector<int>    cols(nnz);
    std::vector<Number> val(nnz);
    Utilities::CUDA::copy_to_host(row_ptr_dev.get(), rows);
    Utilities::CUDA::copy_to_host(column_index_dev.get(), cols);
    Utilities::CUDA::copy_to_host(val_dev.get(), val);

    unsigned int width = width_;

    std::ios::fmtflags old_flags     = out.flags();
    unsigned int       old_precision = out.precision(precision);

    if (scientific)
      {
        out.setf(std::ios::scientific, std::ios::floatfield);
        if (!width)
          width = precision + 7;
      }
    else
      {
        out.setf(std::ios::fixed, std::ios::floatfield);
        if (!width)
          width = precision + 2;
      }

    for (size_type i = 0; i < n_rows; ++i)
      {
        size_type j = rows[i];
        for (size_type k = 0; k < n_cols; ++k)
          {
            if (k == cols[j])
              {
                out << std::setw(width) << val[j] * Number(denominator) << ' ';
                ++j;
              }
            else
              out << std::setw(width) << zero_string << ' ';
          }
        out << std::endl;
      };
    AssertThrow(out, ExcIO());

    // reset output format
    out.precision(old_precision);
    out.flags(old_flags);
  }
} // namespace CUDAWrappers

DEAL_II_NAMESPACE_CLOSE

#endif
#endif