GpuSparseContext.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2026 Rasmus Munk Larsen <rmlarsen@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
// SPDX-License-Identifier: MPL-2.0
 
// GPU sparse matrix-vector (SpMV) and sparse matrix-dense matrix (SpMM) multiply
// via cuSPARSE.
//
// SparseContext owns the cuSPARSE descriptors and device buffers. It takes
// SparseMatrix<Scalar, ColMajor> (CSC), implicitly converting RowMajor input, and
// can borrow a gpu::Context so that sparse products share a stream with BLAS-1
// operations — which removes the cross-stream event waits in solvers like CG.
// It also takes BlockSparseMatrix: as BSR (internal::BsrBinding) where
// internal::use_cusparse_bsr holds, as the scalar-level toSparse() copy otherwise.
//
// Caching: host-input calls re-upload the values *and* index arrays on every
// call. Host pointer identity cannot detect a sparsity pattern rewritten in
// place or assigned into the same allocations (SparseMatrix reuses them for
// same-shape assignments), so the structure is never assumed unchanged. What
// is cached are the cuSPARSE descriptors (sparse descriptor keyed on
// dimensions + nonzero count + block size, dense descriptors keyed on shape)
// and the SpMV/SpMM workspace-size queries. For repeated products with no
// re-upload at all, use deviceView().
//
// Not thread-safe: concurrent multiply* calls on one instance race on the
// cuSPARSE handle, the bound stream, and the cached buffers. Use one per thread.
 
#ifndef EIGEN_GPU_SPARSE_CONTEXT_H
#define EIGEN_GPU_SPARSE_CONTEXT_H
 
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
 
#include <cstdint>
 
#include "./CuSparseSupport.h"
#include "./FwdDecl.h"
 
namespace Eigen {
namespace gpu {
 
namespace internal {
template <int BlockRows, int BlockCols>
struct use_cusparse_bsr
    : Eigen::internal::bool_constant<(EIGEN_HAS_CUSPARSE_BSR != 0 && BlockRows == BlockCols && BlockRows >= 2)> {};
}  // namespace internal
 
#if EIGEN_HAS_CUSPARSE_BSR
namespace internal {

template <typename Scalar, typename StorageIndex>
struct BsrArrays {
  Index brows;
  Index bcols;
  Index bnnz;
  Index blockSize;
  const StorageIndex* outer;
  const StorageIndex* inner;
  const Scalar* values;
 
  Index rows() const { return brows * blockSize; }
  Index cols() const { return bcols * blockSize; }
  Index nonZeros() const { return bnnz * blockSize * blockSize; }
};

template <typename Bsm>
BsrArrays<typename Bsm::Scalar, typename Bsm::StorageIndex> bsr_of(const Bsm& M) {
  return BsrArrays<typename Bsm::Scalar, typename Bsm::StorageIndex>{
      /*brows=*/M.blockOuterSize(), /*bcols=*/M.blockInnerSize(),
      /*bnnz=*/M.nonZeroBlocks(),   /*blockSize=*/Index(Bsm::BlockRows),
      /*outer=*/M.outerIndexPtr(),  /*inner=*/M.innerIndexPtr(),
      /*values=*/M.valuePtr()};
}

template <typename Bsm>
class BsrBinding {
  EIGEN_STATIC_ASSERT((use_cusparse_bsr<int(Bsm::BlockRows), int(Bsm::BlockCols)>::value),
                      CUSPARSE_BSR_REQUIRES_SQUARE_BLOCKS_OF_SIZE_AT_LEAST_2)
 
 public:
  using Scalar = typename Bsm::Scalar;
  using Arrays = BsrArrays<Scalar, typename Bsm::StorageIndex>;
 
  BsrBinding(const Bsm& A, GpuOp op) {
    if (op == GpuOp::ConjTrans && !NumTraits<Scalar>::IsComplex) op = GpuOp::Trans;
    const bool direct = Bsm::IsRowMajor ? (op == GpuOp::NoTrans) : (op == GpuOp::Trans);
    if (direct) {
      arrays_ = bsr_of(A);
      return;
    }
    if (op == GpuOp::ConjTrans && !Bsm::IsRowMajor) {
      // ColMajor conj(A) reads as BSR of conj(A)^T = A^H.
      copy_ = A.unaryExpr(Eigen::internal::scalar_conjugate_op<Scalar>());
    } else if (op == GpuOp::ConjTrans) {
      copy_ = A.adjoint();
    } else {
      copy_ = A.transpose();
    }
    arrays_ = bsr_of(copy_);
  }
 
  BsrBinding(const BsrBinding&) = delete;
  BsrBinding& operator=(const BsrBinding&) = delete;
 
  const Arrays& arrays() const { return arrays_; }
 
 private:
  // transpose()/adjoint() swap the block dimensions, which are equal here, so
  // the copy has A's own type in either storage order.
  Bsm copy_;
  Arrays arrays_{};
};
 
}  // namespace internal
#endif  // EIGEN_HAS_CUSPARSE_BSR

template <typename Scalar_>
class SpMVExpr {
 public:
  using Scalar = Scalar_;
  SpMVExpr(const DeviceSparseView<Scalar>& view, const DeviceMatrix<Scalar>& x) : view_(view), x_(x) {}
  const DeviceSparseView<Scalar>& view() const { return view_; }
  const DeviceMatrix<Scalar>& x() const { return x_; }
 
 private:
  const DeviceSparseView<Scalar>& view_;
  const DeviceMatrix<Scalar>& x_;
};

template <typename Scalar_>
class SpMVAffineExpr {
 public:
  using Scalar = Scalar_;
  SpMVAffineExpr(const SpMVExpr<Scalar>& product, Scalar alpha, Scalar beta, const DeviceMatrix<Scalar>& addend)
      : view_(product.view()), x_(product.x()), alpha_(alpha), beta_(beta), addend_(addend) {}
  const DeviceSparseView<Scalar>& view() const { return view_; }
  const DeviceMatrix<Scalar>& x() const { return x_; }
  Scalar alpha() const { return alpha_; }
  Scalar beta() const { return beta_; }
  const DeviceMatrix<Scalar>& addend() const { return addend_; }
 
 private:
  const DeviceSparseView<Scalar>& view_;
  const DeviceMatrix<Scalar>& x_;
  Scalar alpha_;
  Scalar beta_;
  const DeviceMatrix<Scalar>& addend_;
};
 
template <typename S>
SpMVAffineExpr<S> operator-(const DeviceMatrix<S>& b, const SpMVExpr<S>& p) {
  return SpMVAffineExpr<S>(p, S(-1), S(1), b);
}
 
template <typename S>
SpMVAffineExpr<S> operator+(const DeviceMatrix<S>& b, const SpMVExpr<S>& p) {
  return SpMVAffineExpr<S>(p, S(1), S(1), b);
}
 
template <typename S>
SpMVAffineExpr<S> operator+(const SpMVExpr<S>& p, const DeviceMatrix<S>& b) {
  return b + p;
}
 
template <typename S>
SpMVAffineExpr<S> operator-(const SpMVExpr<S>& p, const DeviceMatrix<S>& b) {
  return SpMVAffineExpr<S>(p, S(1), S(-1), b);
}
 
}  // namespace gpu
 
namespace internal {
// DeviceSparseView is a matrix-free matrix type for Eigen's iterative solvers: it looks like a
// SparseMatrix to the traits machinery, and IterativeSolverBase stores it by pointer because
// Ref<> cannot bind to it (the same mechanism as the matrix-free example in the documentation).
template <typename Scalar_>
struct traits<gpu::DeviceSparseView<Scalar_>> : traits<SparseMatrix<Scalar_, ColMajor, int>> {};
}  // namespace internal
 
namespace gpu {

template <typename Scalar_>
class DeviceSparseView : public EigenBase<DeviceSparseView<Scalar_>> {
 public:
  using Scalar = Scalar_;
  using RealScalar = typename NumTraits<Scalar>::Real;
  using StorageIndex = int;
  using SpMat = SparseMatrix<Scalar, ColMajor, int>;
  // What IterativeSolverBase and ConjugateGradient read from a matrix type.
  static constexpr int ColsAtCompileTime = Dynamic;
  static constexpr int MaxColsAtCompileTime = Dynamic;
  static constexpr bool IsRowMajor = false;
 
  DeviceSparseView(SparseContext<Scalar>& ctx, Index rows, Index cols, uint64_t generation)
      : ctx_(ctx), rows_(rows), cols_(cols), generation_(generation) {}

  SpMVExpr<Scalar> operator*(const DeviceMatrix<Scalar>& x) const { return SpMVExpr<Scalar>(*this, x); }
 
  Index rows() const { return rows_; }
  Index cols() const { return cols_; }
  Index nonZeros() const { return ctx_.nonZeros(); }
  const SparseContext<Scalar>& context() const { return ctx_; }

  uint64_t generation() const { return generation_; }
 
 private:
  SparseContext<Scalar>& ctx_;
  Index rows_;
  Index cols_;
  uint64_t generation_;
};
 
template <typename Scalar_>
class SparseContext {
 public:
  using Scalar = Scalar_;
  using RealScalar = typename NumTraits<Scalar>::Real;
  using StorageIndex = int;
  using SpMat = SparseMatrix<Scalar, ColMajor, StorageIndex>;
  using DenseVector = Matrix<Scalar, Dynamic, 1>;
  using DenseMatrix = Matrix<Scalar, Dynamic, Dynamic, ColMajor>;
  template <int Options, int BlockRows, int BlockCols>
  using BlockSpMat = BlockSparseMatrix<Scalar, Options, BlockRows, BlockCols, StorageIndex>;

  SparseContext() : owns_handle_(true) {
    EIGEN_CUDA_RUNTIME_CHECK(cudaStreamCreate(&stream_));
    owns_stream_ = true;
    EIGEN_CUSPARSE_CHECK(cusparseCreate(&handle_));
    EIGEN_CUSPARSE_CHECK(cusparseSetStream(handle_, stream_));
  }

  explicit SparseContext(Context& ctx)
      : stream_(ctx.stream()), handle_(ctx.cusparseHandle()), owns_stream_(false), owns_handle_(false) {}
 
  ~SparseContext() {
    destroy_spmat_descriptor(/*checked=*/false);
    destroy_dense_descriptors();
    if (owns_handle_ && handle_) (void)cusparseDestroy(handle_);
    if (owns_stream_ && stream_) (void)cudaStreamDestroy(stream_);
  }
 
  SparseContext(const SparseContext&) = delete;
  SparseContext& operator=(const SparseContext&) = delete;

  DeviceSparseView<Scalar> deviceView(const SpMat& A) {
    eigen_assert(A.isCompressed());
    upload_sparse(A);
    return DeviceSparseView<Scalar>(*this, A.rows(), A.cols(), generation_);
  }

  uint64_t uploadGeneration() const { return generation_; }

  Index nonZeros() const { return cached_nnz_; }

  template <typename InputType, typename Rhs>
  DenseVector multiply(const SparseMatrixBase<InputType>& A, const MatrixBase<Rhs>& x) {
    return multiply_host_return(A, x, GpuOp::NoTrans);
  }

  template <typename InputType, typename Rhs, typename Dest>
  void multiply(const SparseMatrixBase<InputType>& A, const MatrixBase<Rhs>& x, MatrixBase<Dest>& y,
                Scalar alpha = Scalar(1), Scalar beta = Scalar(0), GpuOp op = GpuOp::NoTrans) {
    const InputType& input = A.derived();
    internal::check_storage_index_bounds<StorageIndex>(input.rows(), input.cols(), input.nonZeros());
    SpMat storage;
    const SpMat& mat = internal::bind_sparse<SpMat>(input, storage);
    multiply_host_impl(mat, x.derived(), y.derived(), alpha, beta, internal::to_cusparse_op<Scalar>(op));
  }

  template <typename InputType>
  void multiply(const SparseMatrixBase<InputType>& A, const DeviceMatrix<Scalar>& d_x, DeviceMatrix<Scalar>& d_y) {
    multiply(A, d_x, d_y, Scalar(1), Scalar(0), GpuOp::NoTrans);
  }

  template <typename InputType>
  void multiply(const SparseMatrixBase<InputType>& A, const DeviceMatrix<Scalar>& d_x, DeviceMatrix<Scalar>& d_y,
                Scalar alpha, Scalar beta, GpuOp op = GpuOp::NoTrans) {
    const InputType& input = A.derived();
    internal::check_storage_index_bounds<StorageIndex>(input.rows(), input.cols(), input.nonZeros());
    SpMat storage;
    const SpMat& mat = internal::bind_sparse<SpMat>(input, storage);
    upload_sparse(mat);
    spmv_device_exec(d_x, d_y, alpha, beta, op);
  }

  template <typename InputType, typename Rhs>
  DenseVector multiplyT(const SparseMatrixBase<InputType>& A, const MatrixBase<Rhs>& x) {
    return multiply_host_return(A, x, GpuOp::Trans);
  }

  template <typename InputType, typename Rhs>
  DenseVector multiplyAdjoint(const SparseMatrixBase<InputType>& A, const MatrixBase<Rhs>& x) {
    return multiply_host_return(A, x, GpuOp::ConjTrans);
  }

  template <typename InputType, typename Rhs>
  DenseMatrix multiplyMat(const SparseMatrixBase<InputType>& A, const MatrixBase<Rhs>& X, GpuOp op = GpuOp::NoTrans) {
    const InputType& input = A.derived();
    internal::check_storage_index_bounds<StorageIndex>(input.rows(), input.cols(), input.nonZeros());
    SpMat storage;
    const SpMat& mat = internal::bind_sparse<SpMat>(input, storage);
    const DenseMatrix rhs(X.derived());
 
    const cusparseOperation_t cu_op = internal::to_cusparse_op<Scalar>(op);
    const Index m = (op == GpuOp::NoTrans) ? mat.rows() : mat.cols();
    const Index k = (op == GpuOp::NoTrans) ? mat.cols() : mat.rows();
    eigen_assert(k == rhs.rows());
 
    const Index n = rhs.cols();
    if (m == 0 || n == 0 || mat.nonZeros() == 0) return DenseMatrix::Zero(m, n);
 
    DenseMatrix Y = DenseMatrix::Zero(m, n);
    spmm_impl(mat, rhs, Y, Scalar(1), Scalar(0), cu_op);
    return Y;
  }
 
  // BlockSparseMatrix overloads: same contracts as the SparseMatrix ones,
  // dispatched on internal::use_cusparse_bsr<BlockRows, BlockCols> to a BSR
  // upload of op(A) or to the scalar-level CSC path.

  template <int Options, int BlockRows, int BlockCols>
  DeviceSparseView<Scalar> deviceView(const BlockSpMat<Options, BlockRows, BlockCols>& A) {
    return device_view_block(A, internal::use_cusparse_bsr<BlockRows, BlockCols>());
  }

  template <int Options, int BlockRows, int BlockCols, typename Rhs>
  DenseVector multiply(const BlockSpMat<Options, BlockRows, BlockCols>& A, const MatrixBase<Rhs>& x) {
    return multiply_host_return_block(A, x, GpuOp::NoTrans);
  }

  template <int Options, int BlockRows, int BlockCols, typename Rhs, typename Dest>
  void multiply(const BlockSpMat<Options, BlockRows, BlockCols>& A, const MatrixBase<Rhs>& x, MatrixBase<Dest>& y,
                Scalar alpha = Scalar(1), Scalar beta = Scalar(0), GpuOp op = GpuOp::NoTrans) {
    multiply_host_block(A, x.derived(), y.derived(), alpha, beta, op,
                        internal::use_cusparse_bsr<BlockRows, BlockCols>());
  }

  template <int Options, int BlockRows, int BlockCols>
  void multiply(const BlockSpMat<Options, BlockRows, BlockCols>& A, const DeviceMatrix<Scalar>& d_x,
                DeviceMatrix<Scalar>& d_y) {
    multiply(A, d_x, d_y, Scalar(1), Scalar(0), GpuOp::NoTrans);
  }

  template <int Options, int BlockRows, int BlockCols>
  void multiply(const BlockSpMat<Options, BlockRows, BlockCols>& A, const DeviceMatrix<Scalar>& d_x,
                DeviceMatrix<Scalar>& d_y, Scalar alpha, Scalar beta, GpuOp op = GpuOp::NoTrans) {
    multiply_device_block(A, d_x, d_y, alpha, beta, op, internal::use_cusparse_bsr<BlockRows, BlockCols>());
  }

  template <int Options, int BlockRows, int BlockCols, typename Rhs>
  DenseVector multiplyT(const BlockSpMat<Options, BlockRows, BlockCols>& A, const MatrixBase<Rhs>& x) {
    return multiply_host_return_block(A, x, GpuOp::Trans);
  }

  template <int Options, int BlockRows, int BlockCols, typename Rhs>
  DenseVector multiplyAdjoint(const BlockSpMat<Options, BlockRows, BlockCols>& A, const MatrixBase<Rhs>& x) {
    return multiply_host_return_block(A, x, GpuOp::ConjTrans);
  }

  template <int Options, int BlockRows, int BlockCols, typename Rhs>
  DenseMatrix multiplyMat(const BlockSpMat<Options, BlockRows, BlockCols>& A, const MatrixBase<Rhs>& X,
                          GpuOp op = GpuOp::NoTrans) {
    return multiply_mat_block(A, X, op, internal::use_cusparse_bsr<BlockRows, BlockCols>());
  }
 
  cudaStream_t stream() const { return stream_; }
 
 private:
  cudaStream_t stream_ = nullptr;
  cusparseHandle_t handle_ = nullptr;
  bool owns_stream_ = false;
  bool owns_handle_ = false;
 
  // Cached device buffers for sparse matrix (grow-only).
  internal::DeviceBuffer d_outerPtr_;
  internal::DeviceBuffer d_innerIdx_;
  internal::DeviceBuffer d_values_;
 
  // Cached device buffers for host-API dense vectors (grow-only).
  internal::DeviceBuffer d_x_;
  internal::DeviceBuffer d_y_;
 
  mutable internal::DeviceBuffer d_workspace_;
 
  // Cached cuSPARSE sparse matrix descriptor, keyed on the scalar-level shape
  // and nonzero count plus the block size (0 for CSC/CSR, B for BSR).
  cusparseSpMatDescr_t spmat_desc_ = nullptr;
  Index cached_rows_ = -1;
  Index cached_cols_ = -1;
  Index cached_nnz_ = -1;
  Index cached_block_size_ = 0;
 
  // Bumped on every sparse upload; DeviceSparseViews record it at creation so
  // a stale view (its data replaced by a later upload) asserts at evaluation.
  uint64_t generation_ = 0;
 
  // Cached dense-vector/matrix descriptors, re-pointed per call and recreated
  // only when the shape changes.
  mutable cusparseDnVecDescr_t x_vec_desc_ = nullptr;
  mutable cusparseDnVecDescr_t y_vec_desc_ = nullptr;
  mutable int64_t x_vec_size_ = -1;
  mutable int64_t y_vec_size_ = -1;
  mutable cusparseDnMatDescr_t x_mat_desc_ = nullptr;
  mutable cusparseDnMatDescr_t y_mat_desc_ = nullptr;
  mutable int64_t x_mat_rows_ = -1, x_mat_cols_ = -1;
  mutable int64_t y_mat_rows_ = -1, y_mat_cols_ = -1;
 
  // Cached workspace-size query results, indexed by cusparseOperation_t.
  // Invalidated when the sparse descriptor or a dense-descriptor shape changes.
  static constexpr size_t kWsUnknown = static_cast<size_t>(-1);
  mutable size_t spmv_ws_size_[3] = {kWsUnknown, kWsUnknown, kWsUnknown};
  mutable size_t spmm_ws_size_[3] = {kWsUnknown, kWsUnknown, kWsUnknown};
 
  static constexpr cusparseIndexType_t kIndexType =
      (sizeof(StorageIndex) == 4) ? CUSPARSE_INDEX_32I : CUSPARSE_INDEX_64I;
  static constexpr cudaDataType_t kValueType = internal::cuda_data_type<Scalar>::value;
 
  // Empty-operand result on the host: y <- beta * y, no device work.
  template <typename Dest>
  static void scale_host(Dest& y, Scalar beta) {
    if (beta == Scalar(0))
      y.setZero();
    else
      y *= beta;
  }
 
  // Shared host-input SpMV entry: y = op(A) * x into a fresh vector.
  template <typename InputType, typename Rhs>
  DenseVector multiply_host_return(const SparseMatrixBase<InputType>& A, const MatrixBase<Rhs>& x, GpuOp op) {
    const InputType& input = A.derived();
    internal::check_storage_index_bounds<StorageIndex>(input.rows(), input.cols(), input.nonZeros());
    SpMat storage;
    const SpMat& mat = internal::bind_sparse<SpMat>(input, storage);
    DenseVector y((op == GpuOp::NoTrans) ? mat.rows() : mat.cols());
    y.setZero();
    multiply_host_impl(mat, x.derived(), y, Scalar(1), Scalar(0), internal::to_cusparse_op<Scalar>(op));
    return y;
  }
 
  template <typename RhsDerived, typename DestDerived>
  void multiply_host_impl(const SpMat& A, const RhsDerived& x, DestDerived& y, Scalar alpha, Scalar beta,
                          cusparseOperation_t op) {
    eigen_assert(A.isCompressed());
 
    const Index m = A.rows();
    const Index n = A.cols();
    const Index nnz = A.nonZeros();
    const Index x_size = (op == CUSPARSE_OPERATION_NON_TRANSPOSE) ? n : m;
    const Index y_size = (op == CUSPARSE_OPERATION_NON_TRANSPOSE) ? m : n;
 
    eigen_assert(x.size() == x_size);
    eigen_assert(y.size() == y_size);
 
    if (m == 0 || n == 0 || nnz == 0) {
      scale_host(y, beta);
      return;
    }
 
    upload_sparse(A);
    spmv_host_cached(x, y, x_size, y_size, alpha, beta, op);
  }
 
  // BlockSparseMatrix primitives, dispatched on use_cusparse_bsr: the true_type
  // overloads bind op(A) as BSR (internal::BsrBinding) and run cuSPARSE with
  // NON_TRANSPOSE, the false_type ones forward the scalar-level toSparse() copy
  // to the CSC overloads.
  template <typename Bsm, typename Rhs>
  DenseVector multiply_host_return_block(const Bsm& A, const MatrixBase<Rhs>& x, GpuOp op) {
    DenseVector y((op == GpuOp::NoTrans) ? A.rows() : A.cols());
    y.setZero();
    multiply_host_block(A, x.derived(), y, Scalar(1), Scalar(0), op,
                        internal::use_cusparse_bsr<int(Bsm::BlockRows), int(Bsm::BlockCols)>());
    return y;
  }
 
  template <typename Bsm>
  DeviceSparseView<Scalar> device_view_block(const Bsm& A, std::false_type) {
    return deviceView(SpMat(A.toSparse()));
  }
 
  template <typename Bsm, typename RhsDerived, typename DestDerived>
  void multiply_host_block(const Bsm& A, const RhsDerived& x, DestDerived& y, Scalar alpha, Scalar beta, GpuOp op,
                           std::false_type) {
    multiply(A.toSparse(), x, y, alpha, beta, op);
  }
 
  template <typename Bsm>
  void multiply_device_block(const Bsm& A, const DeviceMatrix<Scalar>& d_x, DeviceMatrix<Scalar>& d_y, Scalar alpha,
                             Scalar beta, GpuOp op, std::false_type) {
    multiply(A.toSparse(), d_x, d_y, alpha, beta, op);
  }
 
  template <typename Bsm, typename Rhs>
  DenseMatrix multiply_mat_block(const Bsm& A, const MatrixBase<Rhs>& X, GpuOp op, std::false_type) {
    return multiplyMat(A.toSparse(), X, op);
  }
 
#if EIGEN_HAS_CUSPARSE_BSR
  template <typename Bsm>
  DeviceSparseView<Scalar> device_view_block(const Bsm& A, std::true_type) {
    internal::check_storage_index_bounds<StorageIndex>(A.rows(), A.cols(), A.nonZeros());
    const internal::BsrBinding<Bsm> bound(A, GpuOp::NoTrans);
    upload_bsr(bound.arrays());
    return DeviceSparseView<Scalar>(*this, A.rows(), A.cols(), generation_);
  }
 
  template <typename Bsm, typename RhsDerived, typename DestDerived>
  void multiply_host_block(const Bsm& A, const RhsDerived& x, DestDerived& y, Scalar alpha, Scalar beta, GpuOp op,
                           std::true_type) {
    internal::check_storage_index_bounds<StorageIndex>(A.rows(), A.cols(), A.nonZeros());
    const internal::BsrBinding<Bsm> bound(A, op);
    const internal::BsrArrays<Scalar, StorageIndex>& opA = bound.arrays();
 
    eigen_assert(x.size() == opA.cols());
    eigen_assert(y.size() == opA.rows());
 
    if (opA.rows() == 0 || opA.cols() == 0 || opA.bnnz == 0) {
      scale_host(y, beta);
      return;
    }
 
    upload_bsr(opA);
    spmv_host_cached(x, y, opA.cols(), opA.rows(), alpha, beta, CUSPARSE_OPERATION_NON_TRANSPOSE);
  }
 
  template <typename Bsm>
  void multiply_device_block(const Bsm& A, const DeviceMatrix<Scalar>& d_x, DeviceMatrix<Scalar>& d_y, Scalar alpha,
                             Scalar beta, GpuOp op, std::true_type) {
    internal::check_storage_index_bounds<StorageIndex>(A.rows(), A.cols(), A.nonZeros());
    const internal::BsrBinding<Bsm> bound(A, op);
    upload_bsr(bound.arrays());
    spmv_device_exec(d_x, d_y, alpha, beta, GpuOp::NoTrans);
  }
 
  template <typename Bsm, typename Rhs>
  DenseMatrix multiply_mat_block(const Bsm& A, const MatrixBase<Rhs>& X, GpuOp op, std::true_type) {
    internal::check_storage_index_bounds<StorageIndex>(A.rows(), A.cols(), A.nonZeros());
    const internal::BsrBinding<Bsm> bound(A, op);
    const internal::BsrArrays<Scalar, StorageIndex>& opA = bound.arrays();
    const DenseMatrix rhs(X.derived());
    eigen_assert(opA.cols() == rhs.rows());
 
    const Index n = rhs.cols();
    if (opA.rows() == 0 || opA.cols() == 0 || n == 0 || opA.bnnz == 0) return DenseMatrix::Zero(opA.rows(), n);
 
    DenseMatrix Y = DenseMatrix::Zero(opA.rows(), n);
    upload_bsr(opA);
    spmm_host_cached(rhs, Y, opA.rows(), opA.cols(), Scalar(1), Scalar(0), CUSPARSE_OPERATION_NON_TRANSPOSE);
    return Y;
  }
#endif  // EIGEN_HAS_CUSPARSE_BSR
 
  // y = alpha * op(D) * x + beta * y against the cached upload D, staging the
  // host vectors through d_x_ / d_y_. x_size and y_size are the column and row
  // counts of op(D).
  template <typename RhsDerived, typename DestDerived>
  void spmv_host_cached(const RhsDerived& x, DestDerived& y, Index x_size, Index y_size, Scalar alpha, Scalar beta,
                        cusparseOperation_t op) {
    ensure_buffer(d_x_, static_cast<size_t>(x_size) * sizeof(Scalar));
    // Ref binds in place when x is already a contiguous vector; only genuine
    // expressions are evaluated into the Ref's internal temporary.
    const Ref<const DenseVector> x_ref(x);
    EIGEN_CUDA_RUNTIME_CHECK(
        cudaMemcpyAsync(d_x_.get(), x_ref.data(), x_size * sizeof(Scalar), cudaMemcpyHostToDevice, stream_));
 
    ensure_buffer(d_y_, static_cast<size_t>(y_size) * sizeof(Scalar));
    if (beta != Scalar(0)) {
      const Ref<const DenseVector> y_ref(y);
      EIGEN_CUDA_RUNTIME_CHECK(
          cudaMemcpyAsync(d_y_.get(), y_ref.data(), y_size * sizeof(Scalar), cudaMemcpyHostToDevice, stream_));
    }
 
    exec_spmv(x_size, y_size, d_x_.get(), d_y_.get(), alpha, beta, op);
 
    EIGEN_CUDA_RUNTIME_CHECK(
        cudaMemcpyAsync(y.data(), d_y_.get(), y_size * sizeof(Scalar), cudaMemcpyDeviceToHost, stream_));
    EIGEN_CUDA_RUNTIME_CHECK(cudaStreamSynchronize(stream_));
  }
 
 public:
  void spmv_device_exec(const DeviceMatrix<Scalar>& d_x, DeviceMatrix<Scalar>& d_y, Scalar alpha = Scalar(1),
                        Scalar beta = Scalar(0), GpuOp op = GpuOp::NoTrans) const {
    eigen_assert(spmat_desc_ && "sparse matrix not uploaded — call deviceView() or multiply() first");
    check_op_against_upload(op);
    // cuSPARSE SpMV: y must not alias x (undefined behavior).
    eigen_assert(d_x.data() != d_y.data() && "SpMV: output aliases input vector");
 
    const cusparseOperation_t cu_op = internal::to_cusparse_op<Scalar>(op);
    const Index m = cached_rows_;
    const Index n = cached_cols_;
    const Index x_size = (cu_op == CUSPARSE_OPERATION_NON_TRANSPOSE) ? n : m;
    const Index y_size = (cu_op == CUSPARSE_OPERATION_NON_TRANSPOSE) ? m : n;
 
    eigen_assert(d_x.rows() * d_x.cols() == x_size);
 
    if (m == 0 || n == 0 || cached_nnz_ == 0) {
      // Empty A reduces SpMV to y <- beta*y; SparseContext owns no cuBLAS
      // handle for the scale, so the beta != 0 case must be handled by the caller.
      eigen_assert(beta == Scalar(0) && "SpMV with empty A and beta != 0 is unsupported; scale d_y externally");
      if (d_y.rows() * d_y.cols() != y_size) d_y.resize(y_size, 1);
      d_y.setZero(stream_);
      return;
    }
 
    // Ensure d_y is allocated.
    if (d_y.rows() * d_y.cols() != y_size) {
      d_y.resize(y_size, 1);
    }
 
    // Wait for input data to be ready on this stream.
    d_x.waitReady(stream_);
    d_y.waitReady(stream_);
 
    exec_spmv(x_size, y_size, const_cast<void*>(static_cast<const void*>(d_x.data())), static_cast<void*>(d_y.data()),
              alpha, beta, cu_op);
 
    d_y.recordReady(stream_);
  }

  void spmm_device_exec(const DeviceMatrix<Scalar>& d_X, DeviceMatrix<Scalar>& d_Y, Scalar alpha = Scalar(1),
                        Scalar beta = Scalar(0), GpuOp op = GpuOp::NoTrans) const {
    eigen_assert(spmat_desc_ && "sparse matrix not uploaded — call deviceView() or multiply() first");
    check_op_against_upload(op);
    eigen_assert(d_X.data() != d_Y.data() && "SpMM: output aliases input matrix");
 
    const cusparseOperation_t cu_op = internal::to_cusparse_op<Scalar>(op);
    const bool transposed = (cu_op != CUSPARSE_OPERATION_NON_TRANSPOSE);
    const Index m_op = transposed ? cached_cols_ : cached_rows_;
    const Index k_op = transposed ? cached_rows_ : cached_cols_;
    const Index n = d_X.cols();
 
    eigen_assert(d_X.rows() == k_op);
 
    if (m_op == 0 || n == 0 || cached_nnz_ == 0) {
      eigen_assert(beta == Scalar(0) && "SpMM with empty A and beta != 0 is unsupported; scale d_Y externally");
      if (d_Y.rows() != m_op || d_Y.cols() != n) d_Y.resize(m_op, n);
      d_Y.setZero(stream_);
      return;
    }
 
    if (d_Y.rows() != m_op || d_Y.cols() != n) {
      d_Y.resize(m_op, n);
    }
 
    d_X.waitReady(stream_);
    d_Y.waitReady(stream_);
 
    exec_spmm(m_op, k_op, n, const_cast<void*>(static_cast<const void*>(d_X.data())), static_cast<void*>(d_Y.data()),
              alpha, beta, cu_op);
 
    d_Y.recordReady(stream_);
  }
 
 private:
  // cuSPARSE 11.x's cusparseSpMM rejects CSC for matA (CSC support landed in
  // CUDA 12.0). On 11.x we register the same buffers as CSR-of-A^T (dims
  // swapped) and invert the user-facing op before each cuSPARSE call. On 12+
  // we keep the natural CSC path so users pay no extra cost.
#if !defined(CUSPARSE_VERSION) || CUSPARSE_VERSION < 12000
  static constexpr bool kUseCsrOfTranspose = true;
  static constexpr cusparseSpMMAlg_t kSpMMAlg = CUSPARSE_SPMM_CSR_ALG2;
#else
  static constexpr bool kUseCsrOfTranspose = false;
  static constexpr cusparseSpMMAlg_t kSpMMAlg = CUSPARSE_SPMM_ALG_DEFAULT;
#endif
 
  // Map a user-facing op on A to the cuSPARSE op on the cached descriptor.
  // Identity on cuSPARSE 12+ (descriptor is CSC of A, or BSR of op(A));
  // inverted on 11.x (descriptor is CSR of A^T).
  static cusparseOperation_t descriptor_op(cusparseOperation_t user_op) {
    EIGEN_IF_CONSTEXPR (!kUseCsrOfTranspose) return user_op;
    switch (user_op) {
      case CUSPARSE_OPERATION_NON_TRANSPOSE:
        return CUSPARSE_OPERATION_TRANSPOSE;
      case CUSPARSE_OPERATION_TRANSPOSE:
        return CUSPARSE_OPERATION_NON_TRANSPOSE;
      default:
        // CONJUGATE_TRANSPOSE on the CSR-of-A^T descriptor would compute
        // conj(A) * x, not A^H * x — not supported via this representation.
        eigen_assert(false && "CUSPARSE_OPERATION_CONJUGATE_TRANSPOSE not supported on cuSPARSE < 12.0");
        return user_op;
    }
  }
 
  static int op_index(cusparseOperation_t op) {
    switch (op) {
      case CUSPARSE_OPERATION_TRANSPOSE:
        return 1;
      case CUSPARSE_OPERATION_CONJUGATE_TRANSPOSE:
        return 2;
      default:
        return 0;
    }
  }
 
  void invalidate_ws_caches() const {
    for (int i = 0; i < 3; ++i) {
      spmv_ws_size_[i] = kWsUnknown;
      spmm_ws_size_[i] = kWsUnknown;
    }
  }
 
  // Recreate the descriptor when the size changes; otherwise just re-point it
  // at the new device buffer (cusparseDnVecSetValues is a host-side pointer
  // update, no GPU work).
  void update_dnvec(cusparseDnVecDescr_t& desc, int64_t& cur_size, int64_t size, void* ptr) const {
    if (!desc || cur_size != size) {
      if (desc) EIGEN_CUSPARSE_CHECK(cusparseDestroyDnVec(desc));
      EIGEN_CUSPARSE_CHECK(cusparseCreateDnVec(&desc, size, ptr, kValueType));
      cur_size = size;
      invalidate_ws_caches();
    } else {
      EIGEN_CUSPARSE_CHECK(cusparseDnVecSetValues(desc, ptr));
    }
  }
 
  void update_dnmat(cusparseDnMatDescr_t& desc, int64_t& cur_rows, int64_t& cur_cols, int64_t rows, int64_t cols,
                    void* ptr) const {
    if (!desc || cur_rows != rows || cur_cols != cols) {
      if (desc) EIGEN_CUSPARSE_CHECK(cusparseDestroyDnMat(desc));
      // Column-major with ld = rows.
      EIGEN_CUSPARSE_CHECK(cusparseCreateDnMat(&desc, rows, cols, rows, ptr, kValueType, CUSPARSE_ORDER_COL));
      cur_rows = rows;
      cur_cols = cols;
      invalidate_ws_caches();
    } else {
      EIGEN_CUSPARSE_CHECK(cusparseDnMatSetValues(desc, ptr));
    }
  }
 
  void destroy_dense_descriptors() {
    if (x_vec_desc_) (void)cusparseDestroyDnVec(x_vec_desc_);
    if (y_vec_desc_) (void)cusparseDestroyDnVec(y_vec_desc_);
    if (x_mat_desc_) (void)cusparseDestroyDnMat(x_mat_desc_);
    if (y_mat_desc_) (void)cusparseDestroyDnMat(y_mat_desc_);
    x_vec_desc_ = y_vec_desc_ = nullptr;
    x_mat_desc_ = y_mat_desc_ = nullptr;
  }
 
  void exec_spmv(Index x_size, Index y_size, void* d_x_ptr, void* d_y_ptr, Scalar alpha, Scalar beta,
                 cusparseOperation_t op) const {
    const cusparseOperation_t cu_op = descriptor_op(op);
    update_dnvec(x_vec_desc_, x_vec_size_, x_size, d_x_ptr);
    update_dnvec(y_vec_desc_, y_vec_size_, y_size, d_y_ptr);
 
    size_t& ws_size = spmv_ws_size_[op_index(cu_op)];
    if (ws_size == kWsUnknown) {
      EIGEN_CUSPARSE_CHECK(cusparseSpMV_bufferSize(handle_, cu_op, &alpha, spmat_desc_, x_vec_desc_, &beta, y_vec_desc_,
                                                   kValueType, CUSPARSE_SPMV_ALG_DEFAULT, &ws_size));
    }
    ensure_buffer(d_workspace_, ws_size);
 
    EIGEN_CUSPARSE_CHECK(cusparseSpMV(handle_, cu_op, &alpha, spmat_desc_, x_vec_desc_, &beta, y_vec_desc_, kValueType,
                                      CUSPARSE_SPMV_ALG_DEFAULT, d_workspace_.get()));
  }
 
  void exec_spmm(Index m_op, Index k_op, Index n, void* d_x_ptr, void* d_y_ptr, Scalar alpha, Scalar beta,
                 cusparseOperation_t op) const {
    const cusparseOperation_t cu_op = descriptor_op(op);
    // X is k_op x n, Y is m_op x n (column-major, post-op shapes).
    update_dnmat(x_mat_desc_, x_mat_rows_, x_mat_cols_, k_op, n, d_x_ptr);
    update_dnmat(y_mat_desc_, y_mat_rows_, y_mat_cols_, m_op, n, d_y_ptr);
 
    size_t& ws_size = spmm_ws_size_[op_index(cu_op)];
    if (ws_size == kWsUnknown) {
      EIGEN_CUSPARSE_CHECK(cusparseSpMM_bufferSize(handle_, cu_op, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha,
                                                   spmat_desc_, x_mat_desc_, &beta, y_mat_desc_, kValueType, kSpMMAlg,
                                                   &ws_size));
    }
    ensure_buffer(d_workspace_, ws_size);
 
    EIGEN_CUSPARSE_CHECK(cusparseSpMM(handle_, cu_op, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, spmat_desc_,
                                      x_mat_desc_, &beta, y_mat_desc_, kValueType, kSpMMAlg, d_workspace_.get()));
  }
 
  void spmm_impl(const SpMat& A, const DenseMatrix& X, DenseMatrix& Y, Scalar alpha, Scalar beta,
                 cusparseOperation_t op) {
    eigen_assert(A.isCompressed());
 
    // For op != NON_TRANSPOSE, Y = op(A) * X. The dense X / Y descriptors must
    // describe the *post-op* shapes: X has k_op rows (= input dim of op(A)),
    // Y has m_op rows (= output dim of op(A)).
    const bool transposed = (op != CUSPARSE_OPERATION_NON_TRANSPOSE);
    const Index m_op = transposed ? A.cols() : A.rows();
    const Index k_op = transposed ? A.rows() : A.cols();
    const Index n = X.cols();
    const Index nnz = A.nonZeros();
 
    if (m_op == 0 || n == 0 || k_op == 0 || nnz == 0) {
      scale_host(Y, beta);
      return;
    }
 
    upload_sparse(A);
    spmm_host_cached(X, Y, m_op, k_op, alpha, beta, op);
  }
 
  // Y = alpha * op(D) * X + beta * Y against the cached upload D, staging the
  // host matrices through d_x_ / d_y_. X is k_op x n, Y is m_op x n.
  void spmm_host_cached(const DenseMatrix& X, DenseMatrix& Y, Index m_op, Index k_op, Scalar alpha, Scalar beta,
                        cusparseOperation_t op) {
    const Index n = X.cols();
    const size_t x_bytes = static_cast<size_t>(k_op) * static_cast<size_t>(n) * sizeof(Scalar);
    const size_t y_bytes = static_cast<size_t>(m_op) * static_cast<size_t>(n) * sizeof(Scalar);
    ensure_buffer(d_x_, x_bytes);
    ensure_buffer(d_y_, y_bytes);
    EIGEN_CUDA_RUNTIME_CHECK(cudaMemcpyAsync(d_x_.get(), X.data(), x_bytes, cudaMemcpyHostToDevice, stream_));
    if (beta != Scalar(0)) {
      EIGEN_CUDA_RUNTIME_CHECK(cudaMemcpyAsync(d_y_.get(), Y.data(), y_bytes, cudaMemcpyHostToDevice, stream_));
    }
 
    exec_spmm(m_op, k_op, n, d_x_.get(), d_y_.get(), alpha, beta, op);
 
    EIGEN_CUDA_RUNTIME_CHECK(cudaMemcpyAsync(Y.data(), d_y_.get(), y_bytes, cudaMemcpyDeviceToHost, stream_));
    EIGEN_CUDA_RUNTIME_CHECK(cudaStreamSynchronize(stream_));
  }
 
  void upload_sparse(const SpMat& A) {
    // cuSPARSE 12.0+ accepts CSC directly. On cuSPARSE 11.x, cusparseSpMM
    // rejects CSC and CONJUGATE_TRANSPOSE on CSC+complex SpMV silently
    // demotes to TRANSPOSE. We register the same CSC buffers as CSR-of-A^T
    // (dims swapped) on 11.x and invert the op at exec time via
    // descriptor_op() — no transpose-copy required.
    const Index m = A.rows();
    const Index n = A.cols();
    const Index nnz = A.nonZeros();
    upload_arrays(/*outer_count=*/n + 1, /*inner_count=*/nnz, /*value_count=*/nnz, A.outerIndexPtr(), A.innerIndexPtr(),
                  A.valuePtr());
    if (descriptor_key_matches(m, n, nnz, /*block_size=*/0)) return;
 
    destroy_spmat_descriptor(/*checked=*/true);
    EIGEN_IF_CONSTEXPR (kUseCsrOfTranspose) {
      // cuSPARSE 11.x: cusparseSpMM rejects CSC for matA. CSC of A and CSR of
      // A^T share the same buffers, so register the data as CSR-of-A^T (dims
      // swapped) and invert the op in exec_spmv / spmm_impl via descriptor_op.
      EIGEN_CUSPARSE_CHECK(cusparseCreateCsr(&spmat_desc_, n, m, nnz, d_outerPtr_.get(), d_innerIdx_.get(),
                                             d_values_.get(), kIndexType, kIndexType, CUSPARSE_INDEX_BASE_ZERO,
                                             kValueType));
    } else {
      EIGEN_CUSPARSE_CHECK(cusparseCreateCsc(&spmat_desc_, m, n, nnz, d_outerPtr_.get(), d_innerIdx_.get(),
                                             d_values_.get(), kIndexType, kIndexType, CUSPARSE_INDEX_BASE_ZERO,
                                             kValueType));
    }
    set_descriptor_key(m, n, nnz, /*block_size=*/0);
  }
 
#if EIGEN_HAS_CUSPARSE_BSR
  void upload_bsr(const internal::BsrArrays<Scalar, StorageIndex>& opA) {
    upload_arrays(/*outer_count=*/opA.brows + 1, /*inner_count=*/opA.bnnz, /*value_count=*/opA.nonZeros(), opA.outer,
                  opA.inner, opA.values);
    if (descriptor_key_matches(opA.rows(), opA.cols(), opA.nonZeros(), opA.blockSize)) return;
 
    destroy_spmat_descriptor(/*checked=*/true);
    // Row-major blocks: cusparseSpMM accepts no other block layout for BSR.
    EIGEN_CUSPARSE_CHECK(cusparseCreateBsr(&spmat_desc_, opA.brows, opA.bcols, opA.bnnz, opA.blockSize, opA.blockSize,
                                           d_outerPtr_.get(), d_innerIdx_.get(), d_values_.get(), kIndexType,
                                           kIndexType, CUSPARSE_INDEX_BASE_ZERO, kValueType, CUSPARSE_ORDER_ROW));
    set_descriptor_key(opA.rows(), opA.cols(), opA.nonZeros(), opA.blockSize);
  }
#endif  // EIGEN_HAS_CUSPARSE_BSR
 
  void upload_arrays(Index outer_count, Index inner_count, Index value_count, const StorageIndex* host_outer,
                     const StorageIndex* host_inner, const Scalar* host_values) {
    const size_t outer_bytes = static_cast<size_t>(outer_count) * sizeof(StorageIndex);
    const size_t inner_bytes = static_cast<size_t>(inner_count) * sizeof(StorageIndex);
    const size_t val_bytes = static_cast<size_t>(value_count) * sizeof(Scalar);
 
    // Values *and* index arrays are re-uploaded unconditionally: host pointer
    // identity cannot detect a same-shape/same-nnz pattern rewritten in place
    // or assigned into the same allocations (SparseMatrix reuses them), so a
    // structure cache keyed on pointers would silently serve stale indices.
    // Only the cuSPARSE descriptor and the workspace-size queries are cached,
    // keyed on (rows, cols, nnz, block size). Every upload invalidates
    // outstanding DeviceSparseViews via the generation counter.
    ++generation_;
    ensure_buffer(d_values_, val_bytes);
    ensure_buffer(d_outerPtr_, outer_bytes);
    ensure_buffer(d_innerIdx_, inner_bytes);
    EIGEN_CUDA_RUNTIME_CHECK(cudaMemcpyAsync(d_values_.get(), host_values, val_bytes, cudaMemcpyHostToDevice, stream_));
    EIGEN_CUDA_RUNTIME_CHECK(
        cudaMemcpyAsync(d_outerPtr_.get(), host_outer, outer_bytes, cudaMemcpyHostToDevice, stream_));
    EIGEN_CUDA_RUNTIME_CHECK(
        cudaMemcpyAsync(d_innerIdx_.get(), host_inner, inner_bytes, cudaMemcpyHostToDevice, stream_));
  }
 
  // Same shape, nnz and format: the grow-only device buffers cannot have been
  // reallocated, so the existing descriptor still points at the freshly
  // written data.
  bool descriptor_key_matches(Index m, Index n, Index nnz, Index block_size) const {
    return m == cached_rows_ && n == cached_cols_ && nnz == cached_nnz_ && block_size == cached_block_size_;
  }
 
  void set_descriptor_key(Index m, Index n, Index nnz, Index block_size) {
    cached_rows_ = m;
    cached_cols_ = n;
    cached_nnz_ = nnz;
    cached_block_size_ = block_size;
  }
 
  // A BSR descriptor holds op(A) as bound at upload time, and cuSPARSE runs no
  // transposed BSR product.
  void check_op_against_upload(GpuOp op) const {
    eigen_assert((cached_block_size_ == 0 || op == GpuOp::NoTrans) &&
                 "cuSPARSE runs BSR products only with op == NoTrans; pass op to multiply(A, d_x, d_y, ...) instead");
    EIGEN_UNUSED_VARIABLE(op);
  }
 
  // Destroy the sparse-matrix descriptor and reset the cache identity.
  // `checked` selects checking (mid-lifetime rebuilds) vs swallowing failure
  // (noexcept destructor).
  void destroy_spmat_descriptor(bool checked) {
    if (spmat_desc_) {
      if (checked)
        EIGEN_CUSPARSE_CHECK(cusparseDestroySpMat(spmat_desc_));
      else
        (void)cusparseDestroySpMat(spmat_desc_);
      spmat_desc_ = nullptr;
    }
    set_descriptor_key(-1, -1, -1, 0);
    invalidate_ws_caches();
  }
 
  void ensure_buffer(internal::DeviceBuffer& buf, size_t needed) const {
    if (needed > buf.size()) {
      if (buf) EIGEN_CUDA_RUNTIME_CHECK(cudaStreamSynchronize(stream_));
      buf = internal::DeviceBuffer(needed);
    }
  }
};
 
// Defined here because it needs the full SparseContext definition.
 
template <typename Scalar_>
DeviceMatrix<Scalar_>& DeviceMatrix<Scalar_>::operator=(const SpMVExpr<Scalar_>& expr) {
  // Uses the sparse matrix already uploaded by deviceView() — no re-upload on
  // repeated products with the same view. A stale view (the context has since
  // uploaded again, replacing the cached data) is caught here.
  eigen_assert(expr.view().generation() == expr.view().context().uploadGeneration() &&
               "DeviceSparseView is stale: its SparseContext has since uploaded another sparse matrix");
  if (expr.x().cols() <= 1) {
    expr.view().context().spmv_device_exec(expr.x(), *this, Scalar_(1), Scalar_(0), GpuOp::NoTrans);
  } else {
    expr.view().context().spmm_device_exec(expr.x(), *this, Scalar_(1), Scalar_(0), GpuOp::NoTrans);
  }
  return *this;
}
 
template <typename Scalar_>
DeviceMatrix<Scalar_>& DeviceMatrix<Scalar_>::operator=(const SpMVAffineExpr<Scalar_>& expr) {
  const DeviceSparseView<Scalar_>& view = expr.view();
  const DeviceMatrix& addend = expr.addend();
  eigen_assert(view.generation() == view.context().uploadGeneration() &&
               "DeviceSparseView is stale: its SparseContext has since uploaded another sparse matrix");
  eigen_assert(addend.rows() == view.rows() && addend.cols() == expr.x().cols() &&
               "SpMVAffineExpr: the addend must have the shape of the product");
  eigen_assert(&expr.x() != this && "SpMVAffineExpr: the destination aliases the dense operand");
  // The product accumulates into a copy of the addend; when the addend is the destination
  // itself (d_r = d_r - d_A * d_x) cuSPARSE's beta accumulates in place.
  if (&addend != this) copyFrom(Context::threadLocal(), addend);
  // With no stored entries the expression is beta * addend, which spmv_device_exec cannot
  // form: SparseContext owns no cuBLAS handle to scale d_y by beta.
  if (view.nonZeros() == 0 || view.rows() == 0 || view.cols() == 0) {
    if (expr.beta() != Scalar_(1)) scale(Context::threadLocal(), expr.beta());
    return *this;
  }
  if (expr.x().cols() <= 1) {
    view.context().spmv_device_exec(expr.x(), *this, expr.alpha(), expr.beta(), GpuOp::NoTrans);
  } else {
    view.context().spmm_device_exec(expr.x(), *this, expr.alpha(), expr.beta(), GpuOp::NoTrans);
  }
  return *this;
}
 
template <typename Scalar_>
DeviceMatrix<Scalar_>::DeviceMatrix(const SpMVExpr<Scalar_>& expr) : DeviceMatrix() {
  *this = expr;
}
 
template <typename Scalar_>
DeviceMatrix<Scalar_>::DeviceMatrix(const SpMVAffineExpr<Scalar_>& expr) : DeviceMatrix() {
  *this = expr;
}
}  // namespace gpu
}  // namespace Eigen
 
#endif  // EIGEN_GPU_SPARSE_CONTEXT_H