Eigen/GPU [5/5]: BLAS-1 ops, DeviceScalar, device-resident SpMV, and CG interop

libeigen/eigen!2415

Co-authored-by: Rasmus Munk Larsen <rlarsen@nvidia.com>
Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>

.gitignore

@@ -40,3 +40,4 @@ Makefile
 !Eigen/Core
 !Eigen/src/Core
 CLAUDE.md
+.DS_Store

Eigen/src/IterativeLinearSolvers/ConjugateGradient.h

@@ -32,7 +32,10 @@ EIGEN_DONT_INLINE void conjugate_gradient(const MatrixType& mat, const Rhs& rhs,
                                           Index& iters, typename Dest::RealScalar& tol_error) {
   typedef typename Dest::RealScalar RealScalar;
   typedef typename Dest::Scalar Scalar;
-  typedef Matrix<Scalar, Dynamic, 1> VectorType;
+  // Use Dest's plain (owning) type as VectorType. For CPU Matrix/Map this
+  // resolves to Matrix<Scalar,Dynamic,1>. For GPU DeviceMatrix, PlainObject
+  // is DeviceMatrix itself (already owning).
+  typedef typename Dest::PlainObject VectorType;
 
   RealScalar tol = tol_error;
   Index maxIters = iters;

unsupported/Eigen/GPU

@@ -42,6 +42,7 @@
 
 #ifdef EIGEN_USE_GPU
 // IWYU pragma: begin_exports
+#include "src/GPU/DeviceScalar.h"
 #include "src/GPU/DeviceMatrix.h"
 #include "src/GPU/GpuContext.h"
 #include "src/GPU/DeviceExpr.h"
@@ -57,8 +58,10 @@
 #include "src/GPU/CuFftSupport.h"
 #include "src/GPU/GpuFFT.h"
 #include "src/GPU/CuSparseSupport.h"
+#ifdef EIGEN_SPARSECORE_MODULE_H
 #include "src/GPU/GpuSparseContext.h"
-#ifdef EIGEN_CUDSS
+#endif
+#if defined(EIGEN_CUDSS) && defined(EIGEN_SPARSECORE_MODULE_H)
 #include "src/GPU/CuDssSupport.h"
 #include "src/GPU/GpuSparseSolverBase.h"
 #include "src/GPU/GpuSparseLLT.h"

unsupported/Eigen/src/GPU/CuBlasSupport.h

@@ -22,6 +22,7 @@
 
 #include "./GpuSupport.h"
 #include <cublas_v2.h>
+#include <cublasLt.h>
 #include <cstring>
 
 namespace Eigen {
@@ -265,6 +266,101 @@ inline cublasStatus_t cublasXdgmm(cublasHandle_t h, cublasSideMode_t side, int m
                      reinterpret_cast<const cuDoubleComplex*>(x), incx, reinterpret_cast<cuDoubleComplex*>(C), ldc);
 }
 
+// ---- cuBLAS Level-1 wrappers ------------------------------------------------
+// Type-dispatched wrappers for BLAS-1 vector operations: dot, axpy, nrm2, scal, copy.
+// These work with CUBLAS_POINTER_MODE_HOST or CUBLAS_POINTER_MODE_DEVICE depending
+// on the caller's configuration. For device pointer mode, scalar result pointers
+// (dot, nrm2) must point to device memory.
+
+// dot: result = x^T * y (real) or x^H * y (complex conjugate dot)
+inline cublasStatus_t cublasXdot(cublasHandle_t h, int n, const float* x, int incx, const float* y, int incy,
+                                 float* result) {
+  return cublasSdot(h, n, x, incx, y, incy, result);
+}
+inline cublasStatus_t cublasXdot(cublasHandle_t h, int n, const double* x, int incx, const double* y, int incy,
+                                 double* result) {
+  return cublasDdot(h, n, x, incx, y, incy, result);
+}
+inline cublasStatus_t cublasXdot(cublasHandle_t h, int n, const std::complex<float>* x, int incx,
+                                 const std::complex<float>* y, int incy, std::complex<float>* result) {
+  return cublasCdotc(h, n, reinterpret_cast<const cuComplex*>(x), incx, reinterpret_cast<const cuComplex*>(y), incy,
+                     reinterpret_cast<cuComplex*>(result));
+}
+inline cublasStatus_t cublasXdot(cublasHandle_t h, int n, const std::complex<double>* x, int incx,
+                                 const std::complex<double>* y, int incy, std::complex<double>* result) {
+  return cublasZdotc(h, n, reinterpret_cast<const cuDoubleComplex*>(x), incx,
+                     reinterpret_cast<const cuDoubleComplex*>(y), incy, reinterpret_cast<cuDoubleComplex*>(result));
+}
+
+// nrm2: result = ||x||_2 (always returns real)
+inline cublasStatus_t cublasXnrm2(cublasHandle_t h, int n, const float* x, int incx, float* result) {
+  return cublasSnrm2(h, n, x, incx, result);
+}
+inline cublasStatus_t cublasXnrm2(cublasHandle_t h, int n, const double* x, int incx, double* result) {
+  return cublasDnrm2(h, n, x, incx, result);
+}
+inline cublasStatus_t cublasXnrm2(cublasHandle_t h, int n, const std::complex<float>* x, int incx, float* result) {
+  return cublasScnrm2(h, n, reinterpret_cast<const cuComplex*>(x), incx, result);
+}
+inline cublasStatus_t cublasXnrm2(cublasHandle_t h, int n, const std::complex<double>* x, int incx, double* result) {
+  return cublasDznrm2(h, n, reinterpret_cast<const cuDoubleComplex*>(x), incx, result);
+}
+
+// axpy: y += alpha * x
+inline cublasStatus_t cublasXaxpy(cublasHandle_t h, int n, const float* alpha, const float* x, int incx, float* y,
+                                  int incy) {
+  return cublasSaxpy(h, n, alpha, x, incx, y, incy);
+}
+inline cublasStatus_t cublasXaxpy(cublasHandle_t h, int n, const double* alpha, const double* x, int incx, double* y,
+                                  int incy) {
+  return cublasDaxpy(h, n, alpha, x, incx, y, incy);
+}
+inline cublasStatus_t cublasXaxpy(cublasHandle_t h, int n, const std::complex<float>* alpha,
+                                  const std::complex<float>* x, int incx, std::complex<float>* y, int incy) {
+  cuComplex a;
+  std::memcpy(&a, alpha, sizeof(a));
+  return cublasCaxpy(h, n, &a, reinterpret_cast<const cuComplex*>(x), incx, reinterpret_cast<cuComplex*>(y), incy);
+}
+inline cublasStatus_t cublasXaxpy(cublasHandle_t h, int n, const std::complex<double>* alpha,
+                                  const std::complex<double>* x, int incx, std::complex<double>* y, int incy) {
+  cuDoubleComplex a;
+  std::memcpy(&a, alpha, sizeof(a));
+  return cublasZaxpy(h, n, &a, reinterpret_cast<const cuDoubleComplex*>(x), incx, reinterpret_cast<cuDoubleComplex*>(y),
+                     incy);
+}
+
+// SCAL with complex alpha on complex vectors (Cscal/Zscal). The real-alpha
+// overloads (Sscal/Dscal/Csscal/Zdscal) live above with the FFT-scaling forms.
+inline cublasStatus_t cublasXscal(cublasHandle_t h, int n, const std::complex<float>* alpha, std::complex<float>* x,
+                                  int incx) {
+  cuComplex a;
+  std::memcpy(&a, alpha, sizeof(a));
+  return cublasCscal(h, n, &a, reinterpret_cast<cuComplex*>(x), incx);
+}
+inline cublasStatus_t cublasXscal(cublasHandle_t h, int n, const std::complex<double>* alpha, std::complex<double>* x,
+                                  int incx) {
+  cuDoubleComplex a;
+  std::memcpy(&a, alpha, sizeof(a));
+  return cublasZscal(h, n, &a, reinterpret_cast<cuDoubleComplex*>(x), incx);
+}
+
+// copy: y = x
+inline cublasStatus_t cublasXcopy(cublasHandle_t h, int n, const float* x, int incx, float* y, int incy) {
+  return cublasScopy(h, n, x, incx, y, incy);
+}
+inline cublasStatus_t cublasXcopy(cublasHandle_t h, int n, const double* x, int incx, double* y, int incy) {
+  return cublasDcopy(h, n, x, incx, y, incy);
+}
+inline cublasStatus_t cublasXcopy(cublasHandle_t h, int n, const std::complex<float>* x, int incx,
+                                  std::complex<float>* y, int incy) {
+  return cublasCcopy(h, n, reinterpret_cast<const cuComplex*>(x), incx, reinterpret_cast<cuComplex*>(y), incy);
+}
+inline cublasStatus_t cublasXcopy(cublasHandle_t h, int n, const std::complex<double>* x, int incx,
+                                  std::complex<double>* y, int incy) {
+  return cublasZcopy(h, n, reinterpret_cast<const cuDoubleComplex*>(x), incx, reinterpret_cast<cuDoubleComplex*>(y),
+                     incy);
+}
+
 }  // namespace internal
 }  // namespace gpu
 }  // namespace Eigen

unsupported/Eigen/src/GPU/DeviceDispatch.h

@@ -15,6 +15,8 @@
 #include "./InternalHeaderCheck.h"
 
 #include <climits>
+#include <cstdint>
+#include <limits>
 
 #include "./DeviceExpr.h"
 #include "./DeviceBlasExpr.h"
@@ -43,6 +45,10 @@ void dispatch_gemm(
   const DeviceMatrix<Scalar>& A = traits_lhs::matrix(expr.lhs());
   const DeviceMatrix<Scalar>& B = traits_rhs::matrix(expr.rhs());
 
+  // cuBLAS GEMM: C must not alias A or B (undefined behavior).
+  eigen_assert(dst.data() != A.data() && "GEMM: output aliases left operand (use a temporary)");
+  eigen_assert(dst.data() != B.data() && "GEMM: output aliases right operand (use a temporary)");
+
   constexpr cublasOperation_t transA = to_cublas_op(traits_lhs::op);
   constexpr cublasOperation_t transB = to_cublas_op(traits_rhs::op);
 
@@ -464,6 +470,306 @@ void SelfAdjointView<Scalar_, UpLo_>::rankUpdate(const DeviceMatrix<Scalar_>& A,
   internal::dispatch_syrk(Context::threadLocal(), matrix(), expr, alpha, beta);
 }
 
+// ---- Helper: scoped CUBLAS_POINTER_MODE_DEVICE ------------------------------
+// Saves the current pointer mode, switches to device, runs the callable,
+// then restores the original mode. Used by dot, norm, operator*=(DeviceScalar),
+// and operator+=(DeviceScaledDevice).
+
+namespace internal {
+template <typename F>
+void with_device_pointer_mode(cublasHandle_t h, F&& f) {
+  cublasPointerMode_t prev;
+  EIGEN_CUBLAS_CHECK(cublasGetPointerMode(h, &prev));
+  EIGEN_CUBLAS_CHECK(cublasSetPointerMode(h, CUBLAS_POINTER_MODE_DEVICE));
+  f();
+  EIGEN_CUBLAS_CHECK(cublasSetPointerMode(h, prev));
+}
+}  // namespace internal
+
+// ---- DeviceMatrix BLAS-1 out-of-line definitions ----------------------------
+// Defined here because they need the full Context definition.
+// All methods take an explicit Context& so callers can ensure same-stream
+// execution (zero event overhead when all operations share one context).
+//
+// Reduction methods (dot, norm, squaredNorm) use CUBLAS_POINTER_MODE_DEVICE:
+// the scalar result is written to device memory and stays there until read
+// via DeviceScalar's implicit conversion to Scalar (which syncs).
+
+namespace internal {
+// BLAS-1 cuBLAS wrappers take int counts. Index is ptrdiff_t on 64-bit
+// systems, so guard against silent truncation for matrices with > INT_MAX
+// elements.
+inline int blas1_int_size(Index rows, Index cols) {
+  const int64_t total = static_cast<int64_t>(rows) * static_cast<int64_t>(cols);
+  eigen_assert(total <= static_cast<int64_t>((std::numeric_limits<int>::max)()) &&
+               "cuBLAS BLAS-1 length exceeds int range");
+  return static_cast<int>(total);
+}
+}  // namespace internal
+
+template <typename Scalar_>
+DeviceScalar<typename DeviceMatrix<Scalar_>::Scalar> DeviceMatrix<Scalar_>::dot(Context& ctx,
+                                                                                const DeviceMatrix& other) const {
+  const int n = internal::blas1_int_size(rows_, cols_);
+  eigen_assert(n == internal::blas1_int_size(other.rows_, other.cols_));
+  DeviceScalar<Scalar> result(Scalar(0), ctx.stream());
+  if (n > 0) {
+    waitReady(ctx.stream());
+    other.waitReady(ctx.stream());
+    internal::with_device_pointer_mode(ctx.cublasHandle(), [&] {
+      EIGEN_CUBLAS_CHECK(
+          internal::cublasXdot(ctx.cublasHandle(), n, data_.get(), 1, other.data_.get(), 1, result.devicePtr()));
+    });
+  }
+  return result;
+}
+
+namespace internal {
+// Real: dot(x,x) returns DeviceScalar<Scalar> which IS DeviceScalar<RealScalar>.
+// Move-construct without any sync.
+template <typename Scalar, typename RealScalar>
+typename std::enable_if<std::is_same<Scalar, RealScalar>::value, DeviceScalar<RealScalar>>::type squaredNorm_from_dot(
+    DeviceScalar<Scalar>&& d, cudaStream_t) {
+  return std::move(d);
+}
+// Complex: must sync to extract the real part (DeviceScalar arithmetic is real-only).
+template <typename Scalar, typename RealScalar>
+typename std::enable_if<!std::is_same<Scalar, RealScalar>::value, DeviceScalar<RealScalar>>::type squaredNorm_from_dot(
+    DeviceScalar<Scalar>&& d, cudaStream_t stream) {
+  return DeviceScalar<RealScalar>(numext::real(Scalar(d)), stream);
+}
+}  // namespace internal
+
+template <typename Scalar_>
+DeviceScalar<typename NumTraits<Scalar_>::Real> DeviceMatrix<Scalar_>::squaredNorm(Context& ctx) const {
+  // Use dot(x,x) instead of nrm2()^2: dot kernel is ~4.5x faster than nrm2
+  // (nrm2 uses a numerically careful scaled-sum-of-squares algorithm that is
+  // unnecessary for CG convergence checks).
+  using RealScalar = typename NumTraits<Scalar_>::Real;
+  return internal::squaredNorm_from_dot<Scalar_, RealScalar>(dot(ctx, *this), ctx.stream());
+}
+
+template <typename Scalar_>
+DeviceScalar<typename NumTraits<Scalar_>::Real> DeviceMatrix<Scalar_>::norm(Context& ctx) const {
+  using RealScalar = typename NumTraits<Scalar>::Real;
+  const int n = internal::blas1_int_size(rows_, cols_);
+  DeviceScalar<RealScalar> result(RealScalar(0), ctx.stream());
+  if (n > 0) {
+    waitReady(ctx.stream());
+    internal::with_device_pointer_mode(ctx.cublasHandle(), [&] {
+      EIGEN_CUBLAS_CHECK(internal::cublasXnrm2(ctx.cublasHandle(), n, data_.get(), 1, result.devicePtr()));
+    });
+  }
+  return result;
+}
+
+template <typename Scalar_>
+void DeviceMatrix<Scalar_>::setZero(cudaStream_t stream) {
+  if (sizeInBytes() > 0) {
+    waitReady(stream);
+    EIGEN_CUDA_RUNTIME_CHECK(cudaMemsetAsync(data_.get(), 0, sizeInBytes(), stream));
+    recordReady(stream);
+  }
+}
+
+template <typename Scalar_>
+void DeviceMatrix<Scalar_>::setZero(Context& ctx) {
+  setZero(ctx.stream());
+}
+
+template <typename Scalar_>
+void DeviceMatrix<Scalar_>::addScaled(Context& ctx, Scalar alpha, const DeviceMatrix& x) {
+  const int n = internal::blas1_int_size(rows_, cols_);
+  eigen_assert(n == internal::blas1_int_size(x.rows_, x.cols_));
+  if (n > 0) {
+    waitReady(ctx.stream());
+    x.waitReady(ctx.stream());
+    EIGEN_CUBLAS_CHECK(internal::cublasXaxpy(ctx.cublasHandle(), n, &alpha, x.data_.get(), 1, data_.get(), 1));
+    recordReady(ctx.stream());
+  }
+}
+
+template <typename Scalar_>
+void DeviceMatrix<Scalar_>::scale(Context& ctx, Scalar alpha) {
+  const int n = internal::blas1_int_size(rows_, cols_);
+  if (n > 0) {
+    waitReady(ctx.stream());
+    EIGEN_CUBLAS_CHECK(internal::cublasXscal(ctx.cublasHandle(), n, &alpha, data_.get(), 1));
+    recordReady(ctx.stream());
+  }
+}
+
+template <typename Scalar_>
+void DeviceMatrix<Scalar_>::copyFrom(Context& ctx, const DeviceMatrix& other) {
+  // Wait on *this before resize — resize may free the old buffer while another
+  // stream is still reading it.
+  if (!empty()) waitReady(ctx.stream());
+  resize(other.rows_, other.cols_);
+  const int n = internal::blas1_int_size(rows_, cols_);
+  if (n > 0) {
+    other.waitReady(ctx.stream());
+    EIGEN_CUBLAS_CHECK(internal::cublasXcopy(ctx.cublasHandle(), n, other.data_.get(), 1, data_.get(), 1));
+    recordReady(ctx.stream());
+  }
+}
+
+// ---- BLAS-1 operator overloads for CG compatibility -------------------------
+
+// this += alpha * x  (axpy)
+template <typename Scalar_>
+DeviceMatrix<Scalar_>& DeviceMatrix<Scalar_>::operator+=(const Scaled<DeviceMatrix>& expr) {
+  addScaled(Context::threadLocal(), expr.scalar(), internal::device_expr_traits<DeviceMatrix>::matrix(expr.inner()));
+  return *this;
+}
+
+// this -= alpha * x  (axpy with negated alpha)
+template <typename Scalar_>
+DeviceMatrix<Scalar_>& DeviceMatrix<Scalar_>::operator-=(const Scaled<DeviceMatrix>& expr) {
+  addScaled(Context::threadLocal(), -expr.scalar(), internal::device_expr_traits<DeviceMatrix>::matrix(expr.inner()));
+  return *this;
+}
+
+// this += x  (axpy with alpha=1)
+template <typename Scalar_>
+DeviceMatrix<Scalar_>& DeviceMatrix<Scalar_>::operator+=(const DeviceMatrix& other) {
+  Scalar one(1);
+  addScaled(Context::threadLocal(), one, other);
+  return *this;
+}
+
+// this -= x  (axpy with alpha=-1)
+template <typename Scalar_>
+DeviceMatrix<Scalar_>& DeviceMatrix<Scalar_>::operator-=(const DeviceMatrix& other) {
+  Scalar neg_one(-1);
+  addScaled(Context::threadLocal(), neg_one, other);
+  return *this;
+}
+
+// this *= alpha  (scal, host pointer)
+template <typename Scalar_>
+DeviceMatrix<Scalar_>& DeviceMatrix<Scalar_>::operator*=(Scalar alpha) {
+  scale(Context::threadLocal(), alpha);
+  return *this;
+}
+
+// this *= alpha  (scal, device pointer — avoids host sync)
+template <typename Scalar_>
+DeviceMatrix<Scalar_>& DeviceMatrix<Scalar_>::operator*=(const DeviceScalar<Scalar>& alpha) {
+  const int n = internal::blas1_int_size(rows_, cols_);
+  if (n > 0) {
+    auto& ctx = Context::threadLocal();
+    waitReady(ctx.stream());
+    internal::with_device_pointer_mode(ctx.cublasHandle(), [&] {
+      EIGEN_CUBLAS_CHECK(internal::cublasXscal(ctx.cublasHandle(), n, alpha.devicePtr(), data_.get(), 1));
+    });
+    recordReady(ctx.stream());
+  }
+  return *this;
+}
+
+// this += DeviceScalar * x  (axpy with CUBLAS_POINTER_MODE_DEVICE)
+template <typename Scalar_>
+DeviceMatrix<Scalar_>& DeviceMatrix<Scalar_>::operator+=(const DeviceScaledDevice<Scalar_>& expr) {
+  const int n = internal::blas1_int_size(rows_, cols_);
+  const auto& x = expr.matrix();
+  eigen_assert(n == internal::blas1_int_size(x.rows_, x.cols_));
+  if (n > 0) {
+    auto& ctx = Context::threadLocal();
+    waitReady(ctx.stream());
+    x.waitReady(ctx.stream());
+    internal::with_device_pointer_mode(ctx.cublasHandle(), [&] {
+      EIGEN_CUBLAS_CHECK(
+          internal::cublasXaxpy(ctx.cublasHandle(), n, expr.alpha().devicePtr(), x.data_.get(), 1, data_.get(), 1));
+    });
+    recordReady(ctx.stream());
+  }
+  return *this;
+}
+
+// this -= DeviceScalar * x  (axpy with negated device scalar)
+template <typename Scalar_>
+DeviceMatrix<Scalar_>& DeviceMatrix<Scalar_>::operator-=(const DeviceScaledDevice<Scalar_>& expr) {
+  auto neg_alpha = -expr.alpha();
+  DeviceScaledDevice<Scalar_> neg_expr(neg_alpha, expr.matrix());
+  return operator+=(neg_expr);
+}
+
+// this = alpha * A + beta * B  (cuBLAS geam)
+template <typename Scalar_>
+DeviceMatrix<Scalar_>& DeviceMatrix<Scalar_>::operator=(const DeviceAddExpr<Scalar_>& expr) {
+  auto& ctx = Context::threadLocal();
+  const auto& A = expr.A();
+  const auto& B = expr.B();
+  eigen_assert(A.rows() == B.rows() && A.cols() == B.cols());
+  const int m = static_cast<int>(A.rows());
+  const int n = static_cast<int>(A.cols());
+  // Wait on *this before resize — resize may free the old buffer while another
+  // stream is still reading it.
+  if (!empty()) waitReady(ctx.stream());
+  resize(A.rows(), A.cols());
+  if (m > 0 && n > 0) {
+    A.waitReady(ctx.stream());
+    B.waitReady(ctx.stream());
+    Scalar_ alpha = expr.alpha();
+    Scalar_ beta = expr.beta();
+    EIGEN_CUBLAS_CHECK(internal::cublasXgeam(ctx.cublasHandle(), CUBLAS_OP_N, CUBLAS_OP_N, m, n, &alpha, A.data(), m,
+                                             &beta, B.data(), m, data_.get(), m));
+    recordReady(ctx.stream());
+  }
+  return *this;
+}
+
+// cwiseProduct (allocating).
+template <typename Scalar_>
+DeviceMatrix<Scalar_> DeviceMatrix<Scalar_>::cwiseProduct(Context& ctx, const DeviceMatrix& other) const {
+  const int n = internal::blas1_int_size(rows_, cols_);
+  eigen_assert(n == internal::blas1_int_size(other.rows_, other.cols_));
+  DeviceMatrix result(rows_, cols_);
+  if (n > 0) {
+    waitReady(ctx.stream());
+    other.waitReady(ctx.stream());
+    internal::device_cwiseProduct(data_.get(), other.data_.get(), result.data_.get(), n, ctx.stream());
+    result.recordReady(ctx.stream());
+  }
+  return result;
+}
+
+// In-place cwiseProduct: this = a .* b (reuses this buffer, no allocation).
+template <typename Scalar_>
+void DeviceMatrix<Scalar_>::cwiseProduct(Context& ctx, const DeviceMatrix& a, const DeviceMatrix& b) {
+  const int n = internal::blas1_int_size(a.rows_, a.cols_);
+  eigen_assert(n == internal::blas1_int_size(b.rows_, b.cols_));
+  if (!empty()) waitReady(ctx.stream());
+  resize(a.rows_, a.cols_);
+  if (n > 0) {
+    a.waitReady(ctx.stream());
+    b.waitReady(ctx.stream());
+    internal::device_cwiseProduct(a.data_.get(), b.data_.get(), data_.get(), n, ctx.stream());
+    recordReady(ctx.stream());
+  }
+}
+
+// Convenience overloads using thread-local default Context.
+template <typename Scalar_>
+DeviceScalar<typename DeviceMatrix<Scalar_>::Scalar> DeviceMatrix<Scalar_>::dot(const DeviceMatrix& other) const {
+  return dot(Context::threadLocal(), other);
+}
+
+template <typename Scalar_>
+DeviceScalar<typename NumTraits<Scalar_>::Real> DeviceMatrix<Scalar_>::squaredNorm() const {
+  return squaredNorm(Context::threadLocal());
+}
+
+template <typename Scalar_>
+DeviceScalar<typename NumTraits<Scalar_>::Real> DeviceMatrix<Scalar_>::norm() const {
+  return norm(Context::threadLocal());
+}
+
+template <typename Scalar_>
+void DeviceMatrix<Scalar_>::setZero() {
+  setZero(Context::threadLocal());
+}
+
 }  // namespace gpu
 }  // namespace Eigen
 

unsupported/Eigen/src/GPU/DeviceExpr.h

@@ -185,6 +185,87 @@ GemmExpr<Lhs, Rhs> operator*(const Lhs& a, const Rhs& b) {
   return {a, b};
 }
 
+// ---- DeviceScaledDevice: DeviceScalar * DeviceMatrix → device-pointer axpy ---
+// Like Scaled but carries a DeviceScalar (device pointer) instead of
+// a host scalar. operator+= dispatches to cuBLAS axpy with POINTER_MODE_DEVICE.
+
+template <typename Scalar_>
+class DeviceScaledDevice {
+ public:
+  using Scalar = Scalar_;
+  DeviceScaledDevice(const DeviceScalar<Scalar>& alpha, const DeviceMatrix<Scalar>& mat) : alpha_(alpha), mat_(mat) {}
+  const DeviceScalar<Scalar>& alpha() const { return alpha_; }
+  const DeviceMatrix<Scalar>& matrix() const { return mat_; }
+
+ private:
+  const DeviceScalar<Scalar>& alpha_;
+  const DeviceMatrix<Scalar>& mat_;
+};
+
+// DeviceScalar * DeviceMatrix → DeviceScaledDevice
+template <typename S>
+DeviceScaledDevice<S> operator*(const DeviceScalar<S>& alpha, const DeviceMatrix<S>& m) {
+  return {alpha, m};
+}
+
+// ---- DeviceAddExpr: a + b → cublasXgeam -------------------------------------
+// Captures `DeviceMatrix + Scaled<DeviceMatrix>` (and reverse).
+// Dispatched to geam: C = alpha * A + beta * B.
+//
+// Note: These operator+/- overloads are intentionally free functions on
+// DeviceMatrix, not Eigen expression templates. DeviceMatrix does not inherit
+// from MatrixBase, so there is no ambiguity with Eigen's own operator+/-.
+// If DeviceMatrix is ever made an Eigen expression type, these would need to
+// be revisited.
+
+template <typename Scalar_>
+class DeviceAddExpr {
+ public:
+  using Scalar = Scalar_;
+  DeviceAddExpr(Scalar alpha, const DeviceMatrix<Scalar>& A, Scalar beta, const DeviceMatrix<Scalar>& B)
+      : alpha_(alpha), A_(A), beta_(beta), B_(B) {}
+  Scalar alpha() const { return alpha_; }
+  Scalar beta() const { return beta_; }
+  const DeviceMatrix<Scalar>& A() const { return A_; }
+  const DeviceMatrix<Scalar>& B() const { return B_; }
+
+ private:
+  Scalar alpha_;
+  const DeviceMatrix<Scalar>& A_;
+  Scalar beta_;
+  const DeviceMatrix<Scalar>& B_;
+};
+
+// DeviceMatrix + DeviceMatrix → DeviceAddExpr (alpha=1, beta=1)
+template <typename S>
+DeviceAddExpr<S> operator+(const DeviceMatrix<S>& a, const DeviceMatrix<S>& b) {
+  return {S(1), a, S(1), b};
+}
+
+// DeviceMatrix + Scaled<DeviceMatrix> → DeviceAddExpr (alpha=1, beta=scaled)
+template <typename S>
+DeviceAddExpr<S> operator+(const DeviceMatrix<S>& a, const Scaled<DeviceMatrix<S>>& b) {
+  return {S(1), a, b.scalar(), b.inner()};
+}
+
+// Scaled<DeviceMatrix> + DeviceMatrix → DeviceAddExpr (alpha=scaled, beta=1)
+template <typename S>
+DeviceAddExpr<S> operator+(const Scaled<DeviceMatrix<S>>& a, const DeviceMatrix<S>& b) {
+  return {a.scalar(), a.inner(), S(1), b};
+}
+
+// DeviceMatrix - DeviceMatrix → DeviceAddExpr (alpha=1, beta=-1)
+template <typename S>
+DeviceAddExpr<S> operator-(const DeviceMatrix<S>& a, const DeviceMatrix<S>& b) {
+  return {S(1), a, S(-1), b};
+}
+
+// DeviceMatrix - Scaled<DeviceMatrix> → DeviceAddExpr (alpha=1, beta=-scaled)
+template <typename S>
+DeviceAddExpr<S> operator-(const DeviceMatrix<S>& a, const Scaled<DeviceMatrix<S>>& b) {
+  return {S(1), a, -b.scalar(), b.inner()};
+}
+
 }  // namespace gpu
 }  // namespace Eigen
 

unsupported/Eigen/src/GPU/DeviceMatrix.h

@@ -57,6 +57,16 @@ template <typename>
 class Assignment;
 template <typename, typename>
 class GemmExpr;
+template <typename>
+class Scaled;
+template <typename>
+class SpMVExpr;
+template <typename>
+class DeviceAddExpr;
+template <typename>
+class DeviceScaledDevice;
+template <typename>
+class DeviceScalar;
 template <typename, int>
 class LltSolveExpr;
 template <typename>
@@ -183,6 +193,8 @@ template <typename Scalar_>
 class DeviceMatrix {
  public:
   using Scalar = Scalar_;
+  using RealScalar = typename NumTraits<Scalar>::Real;
+  using PlainObject = DeviceMatrix;  // owning type (for CG template compatibility)
   using PlainMatrix = Eigen::Matrix<Scalar, Dynamic, Dynamic, ColMajor>;
 
   // ---- Construction / destruction ------------------------------------------
@@ -190,6 +202,18 @@ class DeviceMatrix {
   /** Default: empty (0x0, no allocation). */
   DeviceMatrix() = default;
 
+  /** Allocate uninitialized column vector of given size.
+   * Matches Matrix<Scalar,Dynamic,1>(n) for CG template compatibility. */
+  explicit DeviceMatrix(Index n) : rows_(n), cols_(1) {
+    eigen_assert(n >= 0);
+    size_t bytes = sizeInBytes();
+    if (bytes > 0) {
+      void* p = nullptr;
+      EIGEN_CUDA_RUNTIME_CHECK(cudaMalloc(&p, bytes));
+      data_.reset(static_cast<Scalar*>(p));
+    }
+  }
+
   /** Allocate uninitialized device memory for a rows x cols matrix. */
   DeviceMatrix(Index rows, Index cols) : rows_(rows), cols_(cols) {
     eigen_assert(rows >= 0 && cols >= 0);
@@ -468,6 +492,105 @@ class DeviceMatrix {
   template <int UpLo>
   DeviceMatrix& operator=(const SymmExpr<Scalar, UpLo>& expr);
 
+  // ---- BLAS Level-1 operations ----------------------------------------------
+  // DeviceMatrix is always dense (lda == rows), so a vector is simply a
+  // DeviceMatrix with cols == 1. These BLAS-1 methods operate on the flat
+  // rows*cols element array, making them work for both vectors and matrices.
+  //
+  // All methods take an explicit Context& for stream/handle control.
+  // When everything uses the same context, event waits are skipped (same-stream).
+  // Defined out-of-line in DeviceDispatch.h (needs Context).
+
+  /** Dot product: this^H * other. Returns DeviceScalar — the result stays
+   * on device until read via implicit conversion to Scalar (which syncs).
+   * When used with `auto`, no sync occurs until the value is needed. */
+  DeviceScalar<Scalar> dot(Context& ctx, const DeviceMatrix& other) const;
+
+  /** Squared L2 norm via dot(x, x). Returns DeviceScalar (no sync until read).
+   * For real types, the result stays on device. For complex types, falls back
+   * to host sync (DeviceScalar arithmetic is real-only). */
+  DeviceScalar<typename NumTraits<Scalar>::Real> squaredNorm(Context& ctx) const;
+
+  /** L2 norm. Returns DeviceScalar (no host sync). */
+  DeviceScalar<typename NumTraits<Scalar>::Real> norm(Context& ctx) const;
+
+  /** Set all elements to zero. */
+  void setZero(Context& ctx);
+  void setZero(cudaStream_t stream);
+
+  /** this += alpha * x (cuBLAS axpy). Requires same total size. */
+  void addScaled(Context& ctx, Scalar alpha, const DeviceMatrix& x);
+
+  /** this *= alpha (cuBLAS scal). */
+  void scale(Context& ctx, Scalar alpha);
+
+  /** Deep copy: this = other (cuBLAS copy). Resizes if needed. */
+  void copyFrom(Context& ctx, const DeviceMatrix& other);
+
+  // Convenience overloads using the thread-local default Context.
+  DeviceScalar<Scalar> dot(const DeviceMatrix& other) const;
+  DeviceScalar<typename NumTraits<Scalar>::Real> squaredNorm() const;
+  DeviceScalar<typename NumTraits<Scalar>::Real> norm() const;
+  void setZero();
+
+  // ---- BLAS-1 operator overloads for CG/iterative solver compatibility ------
+  // These allow CG code like `x += alpha * p` to work with DeviceMatrix.
+  // `alpha * DeviceMatrix` already returns `Scaled<DeviceMatrix<Scalar>>`
+  // (defined in DeviceExpr.h). These operators dispatch to cuBLAS axpy/scal.
+  // Defined out-of-line in DeviceDispatch.h.
+
+  /** this += alpha * x (cuBLAS axpy). For `x += alpha * p`. */
+  DeviceMatrix& operator+=(const Scaled<DeviceMatrix>& expr);
+
+  /** this -= alpha * x (cuBLAS axpy with negated alpha). For `r -= alpha * tmp`. */
+  DeviceMatrix& operator-=(const Scaled<DeviceMatrix>& expr);
+
+  /** this += x (cuBLAS axpy with alpha=1). */
+  DeviceMatrix& operator+=(const DeviceMatrix& other);
+
+  /** this -= x (cuBLAS axpy with alpha=-1). */
+  DeviceMatrix& operator-=(const DeviceMatrix& other);
+
+  /** this *= alpha (cuBLAS scal, host pointer mode). For `p *= beta`. */
+  DeviceMatrix& operator*=(Scalar alpha);
+
+  /** this *= alpha (cuBLAS scal, device pointer mode). Avoids host sync. */
+  DeviceMatrix& operator*=(const DeviceScalar<Scalar>& alpha);
+
+  /** Element-wise product: result[i] = this[i] * other[i].
+   * Returns a new DeviceMatrix. Defined out-of-line in DeviceDispatch.h. */
+  DeviceMatrix cwiseProduct(Context& ctx, const DeviceMatrix& other) const;
+
+  /** In-place element-wise product: this[i] = a[i] * b[i].
+   * Reuses this matrix's buffer when sizes match, avoiding cudaMalloc. */
+  void cwiseProduct(Context& ctx, const DeviceMatrix& a, const DeviceMatrix& b);
+
+  /** this += DeviceScalar * x (cuBLAS axpy with POINTER_MODE_DEVICE). */
+  DeviceMatrix& operator+=(const DeviceScaledDevice<Scalar>& expr);
+
+  /** this -= DeviceScalar * x (cuBLAS axpy with negated device scalar). */
+  DeviceMatrix& operator-=(const DeviceScaledDevice<Scalar>& expr);
+
+  /** Assign from an SpMV expression: d_y = d_A * d_x. */
+  DeviceMatrix& operator=(const SpMVExpr<Scalar>& expr);
+
+  /** Assign from an add expression: d_C = alpha * d_A + beta * d_B (cuBLAS geam). */
+  DeviceMatrix& operator=(const DeviceAddExpr<Scalar>& expr);
+
+  /** No-op — all DeviceMatrix operations are implicitly noalias.
+   *
+   * Unlike Eigen's Matrix, where omitting .noalias() triggers a copy to a
+   * temporary for safety, DeviceMatrix dispatches directly to NVIDIA library
+   * calls which have no built-in aliasing protection. Every assignment
+   * (`d_C = d_A * d_B`, `d_y = d_A * d_x`, etc.) behaves as if .noalias()
+   * were specified. The caller must ensure operands don't alias the
+   * destination for GEMM and SpMV. geam (`d_C = d_A + alpha * d_B`) is
+   * safe with aliasing. Debug asserts catch violations.
+   *
+   * This method exists so that `tmp.noalias() = mat * p` compiles for both
+   * Matrix and DeviceMatrix. */
+  DeviceMatrix& noalias() { return *this; }
+
   // ---- Ownership transfer ---------------------------------------------------
 
   /** Adopt an existing device pointer. Caller relinquishes ownership. */

unsupported/Eigen/src/GPU/DeviceScalar.h

@@ -0,0 +1,125 @@
+// 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
+
+// Device-resident scalar for deferred host synchronization.
+//
+// DeviceScalar<Scalar> wraps a single value in device memory. Reductions
+// (dot, nrm2) write results directly to device memory via
+// CUBLAS_POINTER_MODE_DEVICE, deferring host sync until the value is read.
+//
+// Implicit conversion to Scalar triggers cudaStreamSynchronize + download.
+// In CG, this reduces 3 syncs/iter to effectively 1: the first conversion
+// syncs the stream, subsequent conversions in the same expression just
+// download (the stream is already flushed).
+//
+// Usage:
+//   auto dot_val = d_x.dot(d_y);       // DeviceScalar, no sync
+//   auto norm_val = d_r.squaredNorm();  // DeviceScalar, no sync
+//   Scalar alpha = absNew / dot_val;    // sync here (both values downloaded)
+//   d_x += alpha * d_p;                 // host-scalar axpy (as before)
+
+#ifndef EIGEN_GPU_DEVICE_SCALAR_H
+#define EIGEN_GPU_DEVICE_SCALAR_H
+
+// IWYU pragma: private
+#include "./InternalHeaderCheck.h"
+
+#include "./GpuSupport.h"
+#include "./DeviceScalarOps.h"
+
+namespace Eigen {
+namespace gpu {
+
+template <typename Scalar_>
+class DeviceScalar {
+ public:
+  using Scalar = Scalar_;
+
+  /** Allocate uninitialized device scalar. Contents are undefined until written
+   * (e.g., by cuBLAS dot/nrm2 with POINTER_MODE_DEVICE). Consistent with
+   * DeviceMatrix(rows, cols) which also does not zero-initialize. */
+  explicit DeviceScalar(cudaStream_t stream = nullptr) : d_val_(sizeof(Scalar)), stream_(stream) {}
+
+  DeviceScalar(Scalar host_val, cudaStream_t stream) : d_val_(sizeof(Scalar)), stream_(stream) {
+    EIGEN_CUDA_RUNTIME_CHECK(cudaMemcpyAsync(d_val_.get(), &host_val, sizeof(Scalar), cudaMemcpyHostToDevice, stream_));
+  }
+
+  DeviceScalar(DeviceScalar&& o) noexcept : d_val_(std::move(o.d_val_)), stream_(o.stream_) { o.stream_ = nullptr; }
+
+  DeviceScalar& operator=(DeviceScalar&& o) noexcept {
+    if (this != &o) {
+      d_val_ = std::move(o.d_val_);
+      stream_ = o.stream_;
+      o.stream_ = nullptr;
+    }
+    return *this;
+  }
+
+  DeviceScalar(const DeviceScalar&) = delete;
+  DeviceScalar& operator=(const DeviceScalar&) = delete;
+
+  /** Download from device. Synchronizes the stream on first call;
+   * subsequent calls in the same expression are cheap (stream already flushed). */
+  Scalar get() const {
+    Scalar result;
+    EIGEN_CUDA_RUNTIME_CHECK(cudaMemcpyAsync(&result, d_val_.get(), sizeof(Scalar), cudaMemcpyDeviceToHost, stream_));
+    EIGEN_CUDA_RUNTIME_CHECK(cudaStreamSynchronize(stream_));
+    return result;
+  }
+
+  /** Implicit conversion — allows `Scalar alpha = deviceScalar` and
+   * `if (deviceScalar < threshold)`. Triggers sync. */
+  operator Scalar() const { return get(); }
+
+  Scalar* devicePtr() { return static_cast<Scalar*>(d_val_.get()); }
+  const Scalar* devicePtr() const { return static_cast<const Scalar*>(d_val_.get()); }
+  cudaStream_t stream() const { return stream_; }
+
+  // ---- Device-side arithmetic (no host sync) ---------------------------------
+  // Uses NPP from DeviceScalarOps.h. All results stay on device.
+  // Currently supports real types only (float, double). Complex types
+  // fall back to implicit conversion (host sync) for division.
+  //
+  // Note: DeviceScalar has no cross-stream readiness tracking. All
+  // operations must be on the same CUDA stream. This is the natural
+  // pattern in iterative solvers where one GpuContext owns all work.
+
+  friend DeviceScalar operator/(const DeviceScalar& a, const DeviceScalar& b) {
+    eigen_assert(a.stream_ == b.stream_ && "DeviceScalar operator/: operands must share the same stream");
+    DeviceScalar result(a.stream_);
+    gpu::internal::device_scalar_div(a.devicePtr(), b.devicePtr(), result.devicePtr(), a.stream_);
+    return result;
+  }
+
+  friend DeviceScalar operator/(Scalar a, const DeviceScalar& b) {
+    DeviceScalar d_a(a, b.stream_);
+    return d_a / b;
+  }
+
+  friend DeviceScalar operator/(const DeviceScalar& a, Scalar b) {
+    DeviceScalar d_b(b, a.stream_);
+    return a / d_b;
+  }
+
+  DeviceScalar operator-() const {
+    DeviceScalar result(stream_);
+    gpu::internal::device_scalar_neg(devicePtr(), result.devicePtr(), stream_);
+    return result;
+  }
+
+ private:
+  internal::DeviceBuffer d_val_;
+  cudaStream_t stream_ = nullptr;
+};
+
+}  // namespace gpu
+}  // namespace Eigen
+
+#endif  // EIGEN_GPU_DEVICE_SCALAR_H

unsupported/Eigen/src/GPU/DeviceScalarOps.h

@@ -0,0 +1,115 @@
+// 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
+
+// Device-resident scalar and element-wise operations via NPP signals.
+// Header-only — no custom CUDA kernels needed. Uses nppsDiv, nppsMul,
+// nppsMulC from the NPP library (CUDA::npps, part of the CUDA toolkit).
+
+#ifndef EIGEN_GPU_DEVICE_SCALAR_OPS_H
+#define EIGEN_GPU_DEVICE_SCALAR_OPS_H
+
+#include <cuda_runtime.h>
+#include <npps_arithmetic_and_logical_operations.h>
+
+#include "./GpuSupport.h"
+
+namespace Eigen {
+namespace gpu {
+namespace internal {
+
+// ---- NppStreamContext helper ------------------------------------------------
+
+inline NppStreamContext make_npp_stream_ctx(cudaStream_t stream) {
+  // NPP requires nCudaDeviceId / device attributes to match the device that
+  // owns 'stream' at this call. We query each time (cheap relative to the NPP
+  // launch itself) so multi-device or borrowed-stream callers stay correct.
+  NppStreamContext ctx = {};
+  ctx.hStream = stream;
+#if CUDART_VERSION >= 12080
+  // cudaStreamGetDevice (added in CUDA 12.8) returns the device that owns the
+  // stream regardless of the calling thread's current device — safe for
+  // borrowed streams.
+  EIGEN_CUDA_RUNTIME_CHECK(cudaStreamGetDevice(stream, &ctx.nCudaDeviceId));
+#else
+  // Older CUDA runtimes lack cudaStreamGetDevice. Callers using borrowed
+  // streams from a different device must cudaSetDevice() first.
+  EIGEN_CUDA_RUNTIME_CHECK(cudaGetDevice(&ctx.nCudaDeviceId));
+#endif
+  EIGEN_CUDA_RUNTIME_CHECK(cudaDeviceGetAttribute(&ctx.nCudaDevAttrComputeCapabilityMajor,
+                                                  cudaDevAttrComputeCapabilityMajor, ctx.nCudaDeviceId));
+  EIGEN_CUDA_RUNTIME_CHECK(cudaDeviceGetAttribute(&ctx.nCudaDevAttrComputeCapabilityMinor,
+                                                  cudaDevAttrComputeCapabilityMinor, ctx.nCudaDeviceId));
+  EIGEN_CUDA_RUNTIME_CHECK(
+      cudaDeviceGetAttribute(&ctx.nMultiProcessorCount, cudaDevAttrMultiProcessorCount, ctx.nCudaDeviceId));
+  EIGEN_CUDA_RUNTIME_CHECK(cudaDeviceGetAttribute(&ctx.nMaxThreadsPerMultiProcessor,
+                                                  cudaDevAttrMaxThreadsPerMultiProcessor, ctx.nCudaDeviceId));
+  EIGEN_CUDA_RUNTIME_CHECK(
+      cudaDeviceGetAttribute(&ctx.nMaxThreadsPerBlock, cudaDevAttrMaxThreadsPerBlock, ctx.nCudaDeviceId));
+  int shared_mem_per_block = 0;
+  EIGEN_CUDA_RUNTIME_CHECK(
+      cudaDeviceGetAttribute(&shared_mem_per_block, cudaDevAttrMaxSharedMemoryPerBlock, ctx.nCudaDeviceId));
+  ctx.nSharedMemPerBlock = static_cast<size_t>(shared_mem_per_block);
+  EIGEN_CUDA_RUNTIME_CHECK(cudaStreamGetFlags(stream, &ctx.nStreamFlags));
+  return ctx;
+}
+
+// ---- Scalar division: c = a / b (device-resident, async) --------------------
+
+inline void device_scalar_div(const float* a, const float* b, float* c, cudaStream_t stream) {
+  NppStreamContext npp_ctx = make_npp_stream_ctx(stream);
+  nppsDiv_32f_Ctx(b, a, c, 1, npp_ctx);  // NPP: pDst[i] = pSrc2[i] / pSrc1[i]
+}
+
+inline void device_scalar_div(const double* a, const double* b, double* c, cudaStream_t stream) {
+  NppStreamContext npp_ctx = make_npp_stream_ctx(stream);
+  nppsDiv_64f_Ctx(b, a, c, 1, npp_ctx);  // NPP: pDst[i] = pSrc2[i] / pSrc1[i]
+}
+
+// ---- Scalar negation: c = -a (device-resident, async) -----------------------
+
+inline void device_scalar_neg(const float* a, float* c, cudaStream_t stream) {
+  NppStreamContext npp_ctx = make_npp_stream_ctx(stream);
+  nppsMulC_32f_Ctx(a, -1.0f, c, 1, npp_ctx);
+}
+
+inline void device_scalar_neg(const double* a, double* c, cudaStream_t stream) {
+  NppStreamContext npp_ctx = make_npp_stream_ctx(stream);
+  nppsMulC_64f_Ctx(a, -1.0, c, 1, npp_ctx);
+}
+
+// ---- Element-wise vector multiply: c[i] = a[i] * b[i] ----------------------
+
+inline void device_cwiseProduct(const float* a, const float* b, float* c, int n, cudaStream_t stream) {
+  NppStreamContext npp_ctx = make_npp_stream_ctx(stream);
+  nppsMul_32f_Ctx(a, b, c, static_cast<size_t>(n), npp_ctx);
+}
+
+inline void device_cwiseProduct(const double* a, const double* b, double* c, int n, cudaStream_t stream) {
+  NppStreamContext npp_ctx = make_npp_stream_ctx(stream);
+  nppsMul_64f_Ctx(a, b, c, static_cast<size_t>(n), npp_ctx);
+}
+
+// ---- Element-wise vector division: c[i] = a[i] / b[i] ----------------------
+
+inline void device_cwiseQuotient(const float* a, const float* b, float* c, int n, cudaStream_t stream) {
+  NppStreamContext npp_ctx = make_npp_stream_ctx(stream);
+  nppsDiv_32f_Ctx(b, a, c, static_cast<size_t>(n), npp_ctx);  // NPP: dst = src2 / src1
+}
+
+inline void device_cwiseQuotient(const double* a, const double* b, double* c, int n, cudaStream_t stream) {
+  NppStreamContext npp_ctx = make_npp_stream_ctx(stream);
+  nppsDiv_64f_Ctx(b, a, c, static_cast<size_t>(n), npp_ctx);
+}
+
+}  // namespace internal
+}  // namespace gpu
+}  // namespace Eigen
+
+#endif  // EIGEN_GPU_DEVICE_SCALAR_OPS_H

unsupported/Eigen/src/GPU/GpuContext.h

@@ -11,8 +11,8 @@
 // Unified GPU execution context.
 //
 // gpu::Context owns a CUDA stream and NVIDIA library handles (cuBLAS
-// eagerly, cuSOLVER lazily on first use). It is the entry point for all
-// GPU operations on gpu::DeviceMatrix.
+// eagerly, cuSOLVER / cuBLASLt / cuSPARSE lazily on first use). It is the
+// entry point for all GPU operations on gpu::DeviceMatrix.
 //
 // The cuSOLVER handle is created on the first call to cusolverHandle()
 // so that translation units which only use cuFFT or cuBLAS paths (e.g.
@@ -33,6 +33,7 @@
 
 #include "./CuBlasSupport.h"
 #include "./CuSolverSupport.h"
+#include <cusparse.h>
 
 namespace Eigen {
 namespace gpu {
@@ -44,27 +45,34 @@ namespace gpu {
  * Each Context instance creates a dedicated CUDA stream and a cuBLAS handle
  * bound to that stream. The cuSOLVER handle is created on first use via
  * cusolverHandle(); translation units that never call it do not require
- * cuSOLVER at link time. Multiple contexts enable concurrent execution on
- * independent streams.
+ * cuSOLVER at link time. cuBLASLt and cuSPARSE handles are similarly lazy.
+ * Multiple contexts enable concurrent execution on independent streams.
  *
  * A lazily-created thread-local default is available via threadLocal() for
  * simple single-stream usage. A single Context is not thread-safe — use one
  * per thread, or external synchronization, since cuBLAS / cuSOLVER handles
- * are not thread-safe per handle and cusolverHandle() lazy-init is racy.
+ * are not thread-safe per handle and lazy-init of secondary handles is racy.
  */
 class Context {
  public:
+  /** Create a new context with a dedicated CUDA stream. */
   Context() {
     EIGEN_CUDA_RUNTIME_CHECK(cudaStreamCreate(&stream_));
-    EIGEN_CUBLAS_CHECK(cublasCreate(&cublas_));
-    EIGEN_CUBLAS_CHECK(cublasSetStream(cublas_, stream_));
+    init_cublas();
   }
 
+  /** Create a context on an existing stream (e.g., stream 0 = nullptr).
+   * The caller retains ownership of the stream — this context will not destroy it. */
+  explicit Context(cudaStream_t stream) : stream_(stream), owns_stream_(false) { init_cublas(); }
+
   ~Context() {
-    // Indirect call keeps cusolverDnDestroy out of TUs that never call cusolverHandle().
+    // Indirect calls keep cusolverDnDestroy / cusparseDestroy out of TUs that
+    // never call cusolverHandle() / cusparseHandle() (e.g. the cufft test).
+    if (cusparse_destroyer_) (void)cusparse_destroyer_(cusparse_);
     if (cusolver_destroyer_) (void)cusolver_destroyer_(cusolver_);
+    if (cublas_lt_) (void)cublasLtDestroy(cublas_lt_);
     if (cublas_) (void)cublasDestroy(cublas_);
-    if (stream_) (void)cudaStreamDestroy(stream_);
+    if (owns_stream_ && stream_) (void)cudaStreamDestroy(stream_);
   }
 
   // Non-copyable, non-movable (owns library handles).
@@ -73,7 +81,10 @@ class Context {
   Context(Context&&) = delete;
   Context& operator=(Context&&) = delete;
 
-  /** Lazily-created thread-local default context.
+  /** Get the thread-local default context.
+   *
+   * If setThreadLocal() has been called, returns that context.
+   * Otherwise lazily creates a new context with a dedicated stream.
    *
    * \note The thread-local instance is destroyed when the thread exits (or at
    * static destruction time for the main thread). On some CUDA driver
@@ -83,10 +94,17 @@ class Context {
    * To avoid this, call cudaDeviceReset() only after all Context instances
    * (including thread-local ones) have been destroyed. */
   static Context& threadLocal() {
+    Context* override = tl_override_ptr();
+    if (override) return *override;
     thread_local Context ctx;
     return ctx;
   }
 
+  /** Override the thread-local default context for this thread.
+   * The caller retains ownership of \p ctx — it must outlive all uses.
+   * Pass nullptr to restore the lazily-created default. */
+  static void setThreadLocal(Context* ctx) { tl_override_ptr() = ctx; }
+
   cudaStream_t stream() const { return stream_; }
   cublasHandle_t cublasHandle() const { return cublas_; }
 
@@ -100,13 +118,55 @@ class Context {
     return cusolver_;
   }
 
+  /** cuBLASLt handle (lazy-initialized on first GEMM call). */
+  cublasLtHandle_t cublasLtHandle() const {
+    if (!cublas_lt_) {
+      EIGEN_CUBLAS_CHECK(cublasLtCreate(&cublas_lt_));
+    }
+    return cublas_lt_;
+  }
+
+  /** Workspace buffer for cublasLtMatmul (grown lazily by cublaslt_gemm).
+   * Not thread-safe — all GEMM calls must be on this context's stream. */
+  internal::DeviceBuffer* gemmWorkspace() const { return &gemm_workspace_; }
+
+  /** cuSPARSE handle (lazy-initialized on first call). */
+  cusparseHandle_t cusparseHandle() const {
+    if (!cusparse_) {
+      cusparseStatus_t s1 = cusparseCreate(&cusparse_);
+      eigen_assert(s1 == CUSPARSE_STATUS_SUCCESS && "cusparseCreate failed");
+      EIGEN_UNUSED_VARIABLE(s1);
+      cusparseStatus_t s2 = cusparseSetStream(cusparse_, stream_);
+      eigen_assert(s2 == CUSPARSE_STATUS_SUCCESS && "cusparseSetStream failed");
+      EIGEN_UNUSED_VARIABLE(s2);
+      cusparse_destroyer_ = &destroyCusparse;
+    }
+    return cusparse_;
+  }
+
  private:
   static cusolverStatus_t destroyCusolver(cusolverDnHandle_t h) { return cusolverDnDestroy(h); }
+  static cusparseStatus_t destroyCusparse(cusparseHandle_t h) { return cusparseDestroy(h); }
 
   cudaStream_t stream_ = nullptr;
   cublasHandle_t cublas_ = nullptr;
   cusolverDnHandle_t cusolver_ = nullptr;
   cusolverStatus_t (*cusolver_destroyer_)(cusolverDnHandle_t) = nullptr;
+  mutable cublasLtHandle_t cublas_lt_ = nullptr;  // lazy
+  mutable cusparseHandle_t cusparse_ = nullptr;   // lazy
+  mutable cusparseStatus_t (*cusparse_destroyer_)(cusparseHandle_t) = nullptr;
+  mutable internal::DeviceBuffer gemm_workspace_;  // lazy
+  bool owns_stream_ = true;
+
+  static Context*& tl_override_ptr() {
+    thread_local Context* ptr = nullptr;
+    return ptr;
+  }
+
+  void init_cublas() {
+    EIGEN_CUBLAS_CHECK(cublasCreate(&cublas_));
+    EIGEN_CUBLAS_CHECK(cublasSetStream(cublas_, stream_));
+  }
 };
 
 }  // namespace gpu

unsupported/Eigen/src/GPU/GpuEigenSolver.h

@@ -120,6 +120,10 @@ class SelfAdjointEigenSolver {
   Index cols() const { return n_; }
   Index rows() const { return n_; }
 
+  // TODO: Add device-side accessors (deviceEigenvalues(), deviceEigenvectors())
+  // returning DeviceMatrix views of the internal buffers, so users can chain
+  // GPU operations without round-tripping through host memory.
+
   /** Eigenvalues in ascending order. Downloads from device. */
   RealVector eigenvalues() const {
     solver_ctx_.sync_info();

unsupported/Eigen/src/GPU/GpuSVD.h

@@ -19,6 +19,7 @@
 //   SVD<double> svd(A, ComputeThinU | ComputeThinV);
 //   VectorXd S = svd.singularValues();
 //   MatrixXd U = svd.matrixU();       // m×k or m×m
+//   MatrixXd V = svd.matrixV();         // n×k or n×n (matches JacobiSVD)
 //   MatrixXd VT = svd.matrixVT();      // k×n or n×n (this is V^T)
 //   MatrixXd X = svd.solve(B);        // pseudoinverse
 //   MatrixXd X = svd.solve(B, k);     // truncated (top k triplets)
@@ -159,6 +160,10 @@ class SVD {
   Index rows() const { return transposed_ ? n_ : m_; }
   Index cols() const { return transposed_ ? m_ : n_; }
 
+  // TODO: Add device-side accessors (deviceU(), deviceVT(), deviceSingularValues())
+  // returning DeviceMatrix views of the internal buffers, so users can chain
+  // GPU operations without round-tripping through host memory.
+
   /** Singular values (always available). Downloads from device on each call. */
   RealVector singularValues() const {
     solver_ctx_.sync_info();

unsupported/Eigen/src/GPU/GpuSparseContext.h

@@ -11,22 +11,34 @@
 // GPU sparse matrix-vector multiply (SpMV) and sparse matrix-dense matrix
 // multiply (SpMM) via cuSPARSE.
 //
-// SparseContext manages a cuSPARSE handle and device buffers. It accepts
-// Eigen SparseMatrix<Scalar, ColMajor> (CSC) and performs SpMV/SpMM on the
-// GPU. RowMajor input is implicitly converted to ColMajor.
+// SparseContext manages cuSPARSE descriptors and device buffers. It accepts
+// Eigen SparseMatrix<Scalar, ColMajor> (CSC) and performs SpMV/SpMM on the GPU.
+// RowMajor input is implicitly converted to ColMajor.
+//
+// Can borrow a Context for same-stream execution with BLAS-1 ops (zero
+// event overhead in iterative solvers like CG).
 //
 // Thread safety: not thread-safe. Concurrent multiply* calls on a single
 // SparseContext race on the cuSPARSE handle, the bound stream, and the
 // cached device buffers. Use one SparseContext per thread.
 //
 // Usage:
-//   SparseContext<double> ctx;
+//   // Standalone (own stream):
+//   gpu::SparseContext<double> ctx;
 //   VectorXd y = ctx.multiply(A, x);                  // y = A * x
 //   ctx.multiply(A, x, y, 2.0, 1.0);                  // y = 2*A*x + y
 //   ctx.multiply(A, x, y, 1.0, 0.0, gpu::GpuOp::ConjTrans);  // y = A^H * x
 //   VectorXd z = ctx.multiplyT(A, x);                 // z = A^T * x
 //   VectorXcd w = ctx.multiplyAdjoint(A, x);          // w = A^H * x (complex)
 //   MatrixXd Y = ctx.multiplyMat(A, X);               // Y = A * X (multiple RHS)
+//
+//   // Shared context (same stream as BLAS-1 ops):
+//   gpu::Context gpu_ctx;
+//   gpu::SparseContext<double> sparse_ctx(gpu_ctx);
+//   VectorXd y = sparse_ctx.multiply(A, x);
+//
+//   // Device-resident (no host roundtrip):
+//   sparse_ctx.multiply(A, d_x, d_y);                 // DeviceMatrix in/out
 
 #ifndef EIGEN_GPU_SPARSE_CONTEXT_H
 #define EIGEN_GPU_SPARSE_CONTEXT_H
@@ -39,6 +51,57 @@
 namespace Eigen {
 namespace gpu {
 
+// Forward declarations.
+template <typename Scalar_>
+class SparseContext;
+template <typename Scalar_>
+class DeviceSparseView;
+
+/** SpMV expression: DeviceSparseView * DeviceMatrix → SpMVExpr.
+ * Evaluated by DeviceMatrix::operator=(SpMVExpr). */
+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_;
+};
+
+/** Device-resident sparse matrix view. Returned by SparseContext::deviceView().
+ * Lightweight handle referencing the context's cached device data.
+ *
+ * \warning One SparseContext caches one sparse matrix at a time.
+ * Creating a second deviceView on the same context overwrites the first.
+ * For multiple simultaneous sparse matrices, use separate SparseContext
+ * instances (they can share a Context for same-stream execution).
+ *
+ * Supports `d_y = d_A * d_x` via SpMVExpr. */
+template <typename Scalar_>
+class DeviceSparseView {
+ public:
+  using Scalar = Scalar_;
+  using SpMat = SparseMatrix<Scalar, ColMajor, int>;
+
+  DeviceSparseView(SparseContext<Scalar>& ctx, Index rows, Index cols) : ctx_(ctx), rows_(rows), cols_(cols) {}
+
+  /** SpMV expression: d_A * d_x. Evaluated by DeviceMatrix::operator=. */
+  SpMVExpr<Scalar> operator*(const DeviceMatrix<Scalar>& x) const { return SpMVExpr<Scalar>(*this, x); }
+
+  Index rows() const { return rows_; }
+  Index cols() const { return cols_; }
+  const SparseContext<Scalar>& context() const { return ctx_; }
+
+ private:
+  SparseContext<Scalar>& ctx_;
+  Index rows_;
+  Index cols_;
+};
+
 template <typename Scalar_>
 class SparseContext {
  public:
@@ -49,22 +112,47 @@ class SparseContext {
   using DenseVector = Matrix<Scalar, Dynamic, 1>;
   using DenseMatrix = Matrix<Scalar, Dynamic, Dynamic, ColMajor>;
 
-  SparseContext() {
+  /** Standalone: creates own stream and cuSPARSE handle. */
+  SparseContext() : owns_handle_(true) {
     EIGEN_CUDA_RUNTIME_CHECK(cudaStreamCreate(&stream_));
+    owns_stream_ = true;
     EIGEN_CUSPARSE_CHECK(cusparseCreate(&handle_));
     EIGEN_CUSPARSE_CHECK(cusparseSetStream(handle_, stream_));
   }
 
+  /** Borrow a Context: shares stream and cuSPARSE handle.
+   * The Context must outlive this SparseContext. */
+  explicit SparseContext(Context& ctx)
+      : stream_(ctx.stream()), handle_(ctx.cusparseHandle()), owns_stream_(false), owns_handle_(false) {}
+
   ~SparseContext() {
     destroy_descriptors_unchecked();
-    if (handle_) (void)cusparseDestroy(handle_);
-    if (stream_) (void)cudaStreamDestroy(stream_);
+    if (owns_handle_ && handle_) (void)cusparseDestroy(handle_);
+    if (owns_stream_ && stream_) (void)cudaStreamDestroy(stream_);
   }
 
   SparseContext(const SparseContext&) = delete;
   SparseContext& operator=(const SparseContext&) = delete;
 
-  // ---- SpMV: y = A * x -----------------------------------------------------
+  // ---- Device sparse view (for expression syntax: d_y = d_A * d_x) ----------
+
+  /** Upload a sparse matrix to device and return a lightweight view.
+   * The sparse data is uploaded immediately and cached in this context.
+   * The returned view can be used for repeated SpMV without re-uploading.
+   * If the matrix values change, call deviceView() again to re-upload.
+   *
+   * \warning One context caches one matrix. Calling deviceView() again
+   * overwrites the previous upload. For multiple simultaneous matrices,
+   * use separate SparseContext instances sharing the same Context.
+   *
+   * Supports `d_y = d_A * d_x` expression syntax. */
+  DeviceSparseView<Scalar> deviceView(const SpMat& A) {
+    eigen_assert(A.isCompressed());
+    upload_sparse(A);
+    return DeviceSparseView<Scalar>(*this, A.rows(), A.cols());
+  }
+
+  // ---- SpMV: y = A * x (host vectors) --------------------------------------
 
   /** Compute y = A * x. Returns y as a new dense vector. */
   template <typename InputType, typename Rhs>
@@ -74,23 +162,41 @@ class SparseContext {
     const SpMat mat(input);
     DenseVector y(mat.rows());
     y.setZero();
-    multiply_impl(mat, x.derived(), y, Scalar(1), Scalar(0), CUSPARSE_OPERATION_NON_TRANSPOSE);
+    multiply_host_impl(mat, x.derived(), y, Scalar(1), Scalar(0), CUSPARSE_OPERATION_NON_TRANSPOSE);
     return y;
   }
 
-  /** Compute y = alpha * op(A) * x + beta * y (in-place). */
+  /** Compute y = alpha * op(A) * x + beta * y (in-place, host vectors). */
   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();
     check_storage_index_bounds(input.rows(), input.cols(), input.nonZeros());
     const SpMat mat(input);
-    multiply_impl(mat, x.derived(), y.derived(), alpha, beta, internal::to_cusparse_op_for_scalar<Scalar>(op));
+    multiply_host_impl(mat, x.derived(), y.derived(), alpha, beta, internal::to_cusparse_op_for_scalar<Scalar>(op));
   }
 
-  // ---- SpMV transpose: y = A^T * x -----------------------------------------
+  // ---- SpMV: y = A * x (DeviceMatrix, no host roundtrip) -------------------
 
-  /** Compute y = A^T * x. Returns y as a new dense vector. */
+  /** Compute d_y = A * d_x. Device-resident, no host transfer.
+   * Sparse matrix A is uploaded to device (cached). Dense vectors stay on device. */
+  template <typename InputType>
+  void multiply(const SparseMatrixBase<InputType>& A, const DeviceMatrix<Scalar>& d_x, DeviceMatrix<Scalar>& d_y) {
+    const SpMat mat(A.derived());
+    multiply_device_impl(mat, d_x, d_y, Scalar(1), Scalar(0), CUSPARSE_OPERATION_NON_TRANSPOSE);
+  }
+
+  /** Compute d_y = alpha * op(A) * d_x + beta * d_y (DeviceMatrix, in-place). */
+  template <typename InputType>
+  void multiply(const SparseMatrixBase<InputType>& A, const DeviceMatrix<Scalar>& d_x, DeviceMatrix<Scalar>& d_y,
+                Scalar alpha, Scalar beta, cusparseOperation_t op = CUSPARSE_OPERATION_NON_TRANSPOSE) {
+    const SpMat mat(A.derived());
+    multiply_device_impl(mat, d_x, d_y, alpha, beta, op);
+  }
+
+  // ---- SpMV transpose -------------------------------------------------------
+
+  /** Compute y = A^T * x (host vectors). */
   template <typename InputType, typename Rhs>
   DenseVector multiplyT(const SparseMatrixBase<InputType>& A, const MatrixBase<Rhs>& x) {
     const InputType& input = A.derived();
@@ -98,7 +204,7 @@ class SparseContext {
     const SpMat mat(input);
     DenseVector y(mat.cols());
     y.setZero();
-    multiply_impl(mat, x.derived(), y, Scalar(1), Scalar(0), CUSPARSE_OPERATION_TRANSPOSE);
+    multiply_host_impl(mat, x.derived(), y, Scalar(1), Scalar(0), CUSPARSE_OPERATION_TRANSPOSE);
     return y;
   }
 
@@ -112,8 +218,8 @@ class SparseContext {
     const SpMat mat(input);
     DenseVector y(mat.cols());
     y.setZero();
-    multiply_impl(mat, x.derived(), y, Scalar(1), Scalar(0),
-                  internal::to_cusparse_op_for_scalar<Scalar>(GpuOp::ConjTrans));
+    multiply_host_impl(mat, x.derived(), y, Scalar(1), Scalar(0),
+                       internal::to_cusparse_op_for_scalar<Scalar>(GpuOp::ConjTrans));
     return y;
   }
 
@@ -147,32 +253,37 @@ class SparseContext {
  private:
   cudaStream_t stream_ = nullptr;
   cusparseHandle_t handle_ = nullptr;
+  bool owns_stream_ = false;
+  bool owns_handle_ = false;
 
-  // Cached device buffers (grow-only).
+  // Cached device buffers for sparse matrix (grow-only).
   internal::DeviceBuffer d_outerPtr_;
   internal::DeviceBuffer d_innerIdx_;
   internal::DeviceBuffer d_values_;
-  internal::DeviceBuffer d_x_;
-  internal::DeviceBuffer d_y_;
-  internal::DeviceBuffer d_workspace_;
   size_t d_outerPtr_size_ = 0;
   size_t d_innerIdx_size_ = 0;
   size_t d_values_size_ = 0;
+
+  // Cached device buffers for host-API dense vectors (grow-only).
+  internal::DeviceBuffer d_x_;
+  internal::DeviceBuffer d_y_;
   size_t d_x_size_ = 0;
   size_t d_y_size_ = 0;
-  size_t d_workspace_size_ = 0;
 
-  // Cached cuSPARSE descriptors.
+  mutable internal::DeviceBuffer d_workspace_;
+  mutable size_t d_workspace_size_ = 0;
+
+  // Cached cuSPARSE sparse matrix descriptor.
   cusparseSpMatDescr_t spmat_desc_ = nullptr;
   Index cached_rows_ = -1;
   Index cached_cols_ = -1;
   Index cached_nnz_ = -1;
 
-  // ---- SpMV implementation --------------------------------------------------
+  // ---- SpMV with host vectors (upload/download per call) --------------------
 
   template <typename RhsDerived, typename DestDerived>
-  void multiply_impl(const SpMat& A, const RhsDerived& x, DestDerived& y, Scalar alpha, Scalar beta,
-                     cusparseOperation_t op) {
+  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();
@@ -192,16 +303,13 @@ class SparseContext {
       return;
     }
 
-    // Upload sparse matrix to device.
     upload_sparse(A);
 
-    // Upload x to device.
     ensure_buffer(d_x_, d_x_size_, static_cast<size_t>(x_size) * sizeof(Scalar));
     const DenseVector x_tmp(x);
     EIGEN_CUDA_RUNTIME_CHECK(
         cudaMemcpyAsync(d_x_.get(), x_tmp.data(), x_size * sizeof(Scalar), cudaMemcpyHostToDevice, stream_));
 
-    // Upload y to device (for beta != 0).
     ensure_buffer(d_y_, d_y_size_, static_cast<size_t>(y_size) * sizeof(Scalar));
     if (beta != Scalar(0)) {
       const DenseVector y_tmp(y);
@@ -209,26 +317,111 @@ class SparseContext {
           cudaMemcpyAsync(d_y_.get(), y_tmp.data(), y_size * sizeof(Scalar), cudaMemcpyHostToDevice, stream_));
     }
 
-    // Create dense vector descriptors.
-    constexpr cudaDataType_t dtype = internal::cuda_data_type<Scalar>::value;
-    cusparseDnVecDescr_t x_desc = nullptr, y_desc = nullptr;
-    EIGEN_CUSPARSE_CHECK(cusparseCreateDnVec(&x_desc, x_size, d_x_.get(), dtype));
-    EIGEN_CUSPARSE_CHECK(cusparseCreateDnVec(&y_desc, y_size, d_y_.get(), dtype));
+    exec_spmv(x_size, y_size, d_x_.get(), d_y_.get(), alpha, beta, op);
 
-    // Query workspace size.
-    size_t ws_size = 0;
-    EIGEN_CUSPARSE_CHECK(cusparseSpMV_bufferSize(handle_, op, &alpha, spmat_desc_, x_desc, &beta, y_desc, dtype,
-                                                 CUSPARSE_SPMV_ALG_DEFAULT, &ws_size));
-    ensure_buffer(d_workspace_, d_workspace_size_, ws_size);
-
-    // Execute SpMV.
-    EIGEN_CUSPARSE_CHECK(cusparseSpMV(handle_, op, &alpha, spmat_desc_, x_desc, &beta, y_desc, dtype,
-                                      CUSPARSE_SPMV_ALG_DEFAULT, d_workspace_.get()));
-
-    // Download result.
     EIGEN_CUDA_RUNTIME_CHECK(
         cudaMemcpyAsync(y.data(), d_y_.get(), y_size * sizeof(Scalar), cudaMemcpyDeviceToHost, stream_));
     EIGEN_CUDA_RUNTIME_CHECK(cudaStreamSynchronize(stream_));
+  }
+
+  // ---- SpMV with DeviceMatrix (no host transfer) ----------------------------
+
+  // Called by public multiply(A, d_x, d_y) — always re-uploads A.
+  void multiply_device_impl(const SpMat& A, const DeviceMatrix<Scalar>& d_x, DeviceMatrix<Scalar>& d_y, Scalar alpha,
+                            Scalar beta, cusparseOperation_t op) {
+    upload_sparse(A);
+    spmv_device_exec(d_x, d_y, alpha, beta, op);
+  }
+
+ public:
+  /** Execute SpMV using the already-uploaded sparse matrix (no re-upload).
+   * Used by SpMVExpr (d_y = d_A * d_x) for cached deviceView() paths.
+   * The sparse matrix must have been uploaded via deviceView() or multiply(). */
+  void spmv_device_exec(const DeviceMatrix<Scalar>& d_x, DeviceMatrix<Scalar>& d_y, Scalar alpha = Scalar(1),
+                        Scalar beta = Scalar(0), cusparseOperation_t op = CUSPARSE_OPERATION_NON_TRANSPOSE) const {
+    eigen_assert(spmat_desc_ && "sparse matrix not uploaded — call deviceView() or multiply() first");
+    // cuSPARSE SpMV: y must not alias x (undefined behavior).
+    eigen_assert(d_x.data() != d_y.data() && "SpMV: output aliases input vector");
+
+    const Index m = cached_rows_;
+    const Index n = cached_cols_;
+    const Index x_size = (op == CUSPARSE_OPERATION_NON_TRANSPOSE) ? n : m;
+    const Index y_size = (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 beta != 0 must be done 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, 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); inverted on 11.x
+  // (descriptor is CSR of A^T).
+  static cusparseOperation_t descriptor_op(cusparseOperation_t user_op) {
+    if (!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;
+    }
+  }
+
+  // ---- Shared SpMV execution ------------------------------------------------
+
+  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 {
+    constexpr cudaDataType_t dtype = internal::cuda_data_type<Scalar>::value;
+    const cusparseOperation_t cu_op = descriptor_op(op);
+    cusparseDnVecDescr_t x_desc = nullptr, y_desc = nullptr;
+    EIGEN_CUSPARSE_CHECK(cusparseCreateDnVec(&x_desc, x_size, d_x_ptr, dtype));
+    EIGEN_CUSPARSE_CHECK(cusparseCreateDnVec(&y_desc, y_size, d_y_ptr, dtype));
+
+    size_t ws_size = 0;
+    EIGEN_CUSPARSE_CHECK(cusparseSpMV_bufferSize(handle_, cu_op, &alpha, spmat_desc_, x_desc, &beta, y_desc, dtype,
+                                                 CUSPARSE_SPMV_ALG_DEFAULT, &ws_size));
+    ensure_buffer(d_workspace_, d_workspace_size_, ws_size);
+
+    EIGEN_CUSPARSE_CHECK(cusparseSpMV(handle_, cu_op, &alpha, spmat_desc_, x_desc, &beta, y_desc, dtype,
+                                      CUSPARSE_SPMV_ALG_DEFAULT, d_workspace_.get()));
 
     EIGEN_CUSPARSE_CHECK(cusparseDestroyDnVec(x_desc));
     EIGEN_CUSPARSE_CHECK(cusparseDestroyDnVec(y_desc));
@@ -269,24 +462,21 @@ class SparseContext {
       EIGEN_CUDA_RUNTIME_CHECK(cudaMemcpyAsync(d_y_.get(), Y.data(), y_bytes, cudaMemcpyHostToDevice, stream_));
     }
 
-    // Create dense matrix descriptors.
     constexpr cudaDataType_t dtype = internal::cuda_data_type<Scalar>::value;
+    const cusparseOperation_t cu_op = descriptor_op(op);
     cusparseDnMatDescr_t x_desc = nullptr, y_desc = nullptr;
     // Eigen is column-major, so ld = rows.
     EIGEN_CUSPARSE_CHECK(cusparseCreateDnMat(&x_desc, k_op, n, k_op, d_x_.get(), dtype, CUSPARSE_ORDER_COL));
     EIGEN_CUSPARSE_CHECK(cusparseCreateDnMat(&y_desc, m_op, n, m_op, d_y_.get(), dtype, CUSPARSE_ORDER_COL));
 
-    // Query workspace.
     size_t ws_size = 0;
-    EIGEN_CUSPARSE_CHECK(cusparseSpMM_bufferSize(handle_, op, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, spmat_desc_,
-                                                 x_desc, &beta, y_desc, dtype, CUSPARSE_SPMM_ALG_DEFAULT, &ws_size));
+    EIGEN_CUSPARSE_CHECK(cusparseSpMM_bufferSize(handle_, cu_op, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, spmat_desc_,
+                                                 x_desc, &beta, y_desc, dtype, kSpMMAlg, &ws_size));
     ensure_buffer(d_workspace_, d_workspace_size_, ws_size);
 
-    // Execute SpMM.
-    EIGEN_CUSPARSE_CHECK(cusparseSpMM(handle_, op, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, spmat_desc_, x_desc, &beta,
-                                      y_desc, dtype, CUSPARSE_SPMM_ALG_DEFAULT, d_workspace_.get()));
+    EIGEN_CUSPARSE_CHECK(cusparseSpMM(handle_, cu_op, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, spmat_desc_, x_desc,
+                                      &beta, y_desc, dtype, kSpMMAlg, d_workspace_.get()));
 
-    // Download result.
     EIGEN_CUDA_RUNTIME_CHECK(cudaMemcpyAsync(Y.data(), d_y_.get(), y_bytes, cudaMemcpyDeviceToHost, stream_));
     EIGEN_CUDA_RUNTIME_CHECK(cudaStreamSynchronize(stream_));
 
@@ -307,28 +497,17 @@ class SparseContext {
   }
 
   void upload_sparse(const SpMat& A) {
-    // cuSPARSE 12.0+ accepts CSC directly for both SpMV and SpMM. cuSPARSE
-    // 11.x SpMM rejects CSC ("unsupported matrix format for matA (CSC)") and
-    // SpMV with CONJUGATE_TRANSPOSE on CSC+complex silently demotes to
-    // TRANSPOSE. CSR works on every version, so on 11.x we transpose-copy
-    // the user's ColMajor input into a RowMajor (CSR) representation.
-#if CUSPARSE_VERSION >= 12000
+    // 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.
     upload_compressed_arrays(A.rows(), A.cols(), A.nonZeros(),
-                             /*outer_count=*/A.cols() + 1, A.outerIndexPtr(), A.innerIndexPtr(), A.valuePtr(),
-                             /*is_csr=*/false);
-#else
-    using CsrMat = SparseMatrix<Scalar, RowMajor, StorageIndex>;
-    const CsrMat csr(A);
-    upload_compressed_arrays(csr.rows(), csr.cols(), csr.nonZeros(),
-                             /*outer_count=*/csr.rows() + 1, csr.outerIndexPtr(), csr.innerIndexPtr(), csr.valuePtr(),
-                             /*is_csr=*/true);
-    // cudaMemcpyAsync from pageable host memory blocks the host until the
-    // source is consumed, so the CsrMat temporary's lifetime is sufficient.
-#endif
+                             /*outer_count=*/A.cols() + 1, A.outerIndexPtr(), A.innerIndexPtr(), A.valuePtr());
   }
 
   void upload_compressed_arrays(Index m, Index n, Index nnz, Index outer_count, const StorageIndex* host_outer,
-                                const StorageIndex* host_inner, const Scalar* host_values, bool is_csr) {
+                                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>(nnz) * sizeof(StorageIndex);
     const size_t val_bytes = static_cast<size_t>(nnz) * sizeof(Scalar);
@@ -349,8 +528,11 @@ class SparseContext {
       constexpr cusparseIndexType_t idx_type = (sizeof(StorageIndex) == 4) ? CUSPARSE_INDEX_32I : CUSPARSE_INDEX_64I;
       constexpr cudaDataType_t val_type = internal::cuda_data_type<Scalar>::value;
 
-      if (is_csr) {
-        EIGEN_CUSPARSE_CHECK(cusparseCreateCsr(&spmat_desc_, m, n, nnz, d_outerPtr_.get(), d_innerIdx_.get(),
+      if (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(), idx_type, idx_type, CUSPARSE_INDEX_BASE_ZERO,
                                                val_type));
       } else {
@@ -361,7 +543,7 @@ class SparseContext {
       cached_rows_ = m;
       cached_cols_ = n;
       cached_nnz_ = nnz;
-    } else if (is_csr) {
+    } else if (kUseCsrOfTranspose) {
       EIGEN_CUSPARSE_CHECK(cusparseCsrSetPointers(spmat_desc_, d_outerPtr_.get(), d_innerIdx_.get(), d_values_.get()));
     } else {
       EIGEN_CUSPARSE_CHECK(cusparseCscSetPointers(spmat_desc_, d_outerPtr_.get(), d_innerIdx_.get(), d_values_.get()));
@@ -389,7 +571,7 @@ class SparseContext {
     cached_nnz_ = -1;
   }
 
-  void ensure_buffer(internal::DeviceBuffer& buf, size_t& current_size, size_t needed) {
+  void ensure_buffer(internal::DeviceBuffer& buf, size_t& current_size, size_t needed) const {
     if (needed > current_size) {
       if (buf) EIGEN_CUDA_RUNTIME_CHECK(cudaStreamSynchronize(stream_));
       buf = internal::DeviceBuffer(needed);
@@ -398,6 +580,17 @@ class SparseContext {
   }
 };
 
+// ---- DeviceMatrix::operator=(SpMVExpr) out-of-line definition ----------------
+// Defined here because it needs the full SparseContext definition.
+
+template <typename Scalar_>
+DeviceMatrix<Scalar_>& DeviceMatrix<Scalar_>::operator=(const SpMVExpr<Scalar_>& expr) {
+  // Use spmv_device_exec — the sparse matrix was already uploaded by deviceView().
+  // No re-upload on repeated SpMV with the same view.
+  expr.view().context().spmv_device_exec(expr.x(), *this, Scalar_(1), Scalar_(0), CUSPARSE_OPERATION_NON_TRANSPOSE);
+  return *this;
+}
+
 }  // namespace gpu
 }  // namespace Eigen
 

unsupported/Eigen/src/GPU/GpuSupport.h

@@ -22,6 +22,7 @@
 #include "./InternalHeaderCheck.h"
 
 #include <cuda_runtime.h>
+#include <vector>
 
 #include <limits>
 #include <memory>
@@ -80,6 +81,70 @@ struct CudaFreeHostDeleter {
   }
 };
 
+// ---- Thread-local pool of small device buffers ------------------------------
+// Recycles allocations up to kSmallBufferThreshold bytes (e.g., DeviceScalar)
+// to avoid cudaMalloc/cudaFree overhead. Larger allocations bypass the pool.
+
+template <size_t SmallBufferThreshold = 256, size_t MaxPoolSize = 64>
+struct DeviceBufferPool {
+  static constexpr size_t kSmallBufferThreshold = SmallBufferThreshold;
+  static constexpr size_t kMaxPoolSize = MaxPoolSize;
+
+  struct Entry {
+    void* ptr;
+    size_t bytes;
+  };
+
+  ~DeviceBufferPool() {
+    for (auto& e : free_list_) (void)cudaFree(e.ptr);
+  }
+
+  void* allocate(size_t bytes) {
+    for (size_t i = 0; i < free_list_.size(); ++i) {
+      if (free_list_[i].bytes >= bytes) {
+        void* p = free_list_[i].ptr;
+        free_list_[i] = free_list_.back();
+        free_list_.pop_back();
+        return p;
+      }
+    }
+    void* p = nullptr;
+    EIGEN_CUDA_RUNTIME_CHECK(cudaMalloc(&p, bytes));
+    return p;
+  }
+
+  void deallocate(void* p, size_t bytes) {
+    if (free_list_.size() < kMaxPoolSize) {
+      free_list_.push_back({p, bytes});
+    } else {
+      (void)cudaFree(p);
+    }
+  }
+
+  static DeviceBufferPool& threadLocal() {
+    thread_local DeviceBufferPool pool;
+    return pool;
+  }
+
+ private:
+  std::vector<Entry> free_list_;
+};
+
+// Stateful deleter that returns small buffers to the thread-local pool and
+// cudaFree's larger ones. size==0 means "always cudaFree" (for adopted ptrs).
+struct PooledCudaFreeDeleter {
+  size_t size = 0;
+
+  void operator()(void* p) const noexcept {
+    if (!p) return;
+    if (size > 0 && size <= DeviceBufferPool<>::kSmallBufferThreshold) {
+      DeviceBufferPool<>::threadLocal().deallocate(p, size);
+    } else {
+      (void)cudaFree(p);
+    }
+  }
+};
+
 // ---- RAII: device buffer ----------------------------------------------------
 
 class DeviceBuffer {
@@ -89,8 +154,12 @@ class DeviceBuffer {
   explicit DeviceBuffer(size_t bytes) {
     if (bytes > 0) {
       void* p = nullptr;
-      EIGEN_CUDA_RUNTIME_CHECK(cudaMalloc(&p, bytes));
-      ptr_.reset(p);
+      if (bytes <= DeviceBufferPool<>::kSmallBufferThreshold) {
+        p = DeviceBufferPool<>::threadLocal().allocate(bytes);
+      } else {
+        EIGEN_CUDA_RUNTIME_CHECK(cudaMalloc(&p, bytes));
+      }
+      ptr_ = std::unique_ptr<void, PooledCudaFreeDeleter>(p, PooledCudaFreeDeleter{bytes});
     }
   }
 
@@ -98,15 +167,18 @@ class DeviceBuffer {
   void* release() noexcept { return ptr_.release(); }
   explicit operator bool() const noexcept { return static_cast<bool>(ptr_); }
 
+  size_t size() const noexcept { return ptr_.get_deleter().size; }
+
   // Adopt an existing device pointer. Caller relinquishes ownership.
+  // Adopted buffers bypass the pool on destruction (deleter size == 0).
   static DeviceBuffer adopt(void* p) noexcept {
     DeviceBuffer b;
-    b.ptr_.reset(p);
+    b.ptr_ = std::unique_ptr<void, PooledCudaFreeDeleter>(p, PooledCudaFreeDeleter{});
     return b;
   }
 
  private:
-  std::unique_ptr<void, CudaFreeDeleter> ptr_;
+  std::unique_ptr<void, PooledCudaFreeDeleter> ptr_;
 };
 
 // ---- RAII: pinned host buffer -----------------------------------------------

unsupported/Eigen/src/GPU/README.md

@@ -1,9 +1,8 @@
 # Eigen GPU Module (`unsupported/Eigen/GPU`)
 
 GPU-accelerated linear algebra for Eigen users, dispatching to NVIDIA CUDA
-Math Libraries (cuBLAS, cuSOLVER, cuFFT, cuSPARSE, cuDSS). Requires CUDA
-11.4+; cuDSS features require CUDA 12.0+ and a separate cuDSS install.
-Header-only.
+Math Libraries (cuBLAS, cuSOLVER, cuFFT, cuSPARSE, cuDSS). Requires CUDA 11.4+;
+cuDSS features require CUDA 12.0+ and a separate cuDSS install. Header-only.
 
 ## Why this module
 
@@ -32,8 +31,9 @@ learning CUDA library APIs directly.
 
 **CPU and GPU coexist.** There is no global compile-time switch that replaces
 CPU implementations (unlike `EIGEN_USE_LAPACKE`). Users choose GPU solvers
-explicitly -- `gpu::LLT<double>` vs `LLT<MatrixXd>`, `gpu::SparseLLT<double>` vs
-`SimplicialLLT<SparseMatrix<double>>` -- and both coexist in the same binary.
+explicitly -- `gpu::LLT<double>` vs `Eigen::LLT<MatrixXd>`,
+`gpu::SparseLLT<double>` vs `SimplicialLLT<SparseMatrix<double>>` -- and both
+coexist in the same binary.
 This also lets users keep the factored matrix on device across multiple solves,
 something impossible with compile-time replacement.
 
@@ -70,6 +70,16 @@ expression template system (`MatrixBase`, `CwiseBinaryOp`, etc.).
 supported expression maps to a single NVIDIA library call. There is no
 coefficient-level evaluation, lazy fusion, or packet operations.
 
+**Interoperability where useful.** `DeviceMatrix` provides the same operator
+signatures as `Matrix` for common vector operations: `+=`, `-=`, `*=`,
+`dot()`, `squaredNorm()`, `norm()`, `setZero()`, and `noalias()`. This makes
+`DeviceMatrix` usable as a drop-in `VectorType` in Eigen algorithm templates
+that rely on these operations. For example, Eigen's `conjugate_gradient()`
+template works with `DeviceMatrix` with a single typedef change -- no
+modifications to the algorithm or the expression template system. Conjugate
+gradient is just the motivating example; we are open to expanding operator
+coverage as needed to support other high-level Eigen algorithms on the GPU.
+
 **Explicit over implicit.** Host-device transfers, stream management, and
 library handle lifetimes are visible in the API. There are no hidden
 allocations or synchronizations except where documented (e.g., `toHost()` must
@@ -77,11 +87,12 @@ synchronize to deliver data to the host).
 
 ## Key concepts
 
-### `DeviceMatrix<Scalar>`
+### `gpu::DeviceMatrix<Scalar>`
 
 A typed RAII wrapper for a dense column-major matrix in GPU device memory.
 This is the GPU counterpart of Eigen's `MatrixX<Scalar>`. A vector is simply
-a `DeviceMatrix` with one column.
+a `DeviceMatrix` with one column. All public GPU classes live in `namespace
+Eigen::gpu`.
 
 ```cpp
 // Upload from host
@@ -104,13 +115,40 @@ MatrixXd C = transfer.get();
 `selfadjointView<UpLo>()`, `llt()`, `lu()`. These return lightweight
 expression objects that are evaluated when assigned.
 
+For BLAS Level-1 operations, `DeviceMatrix` also provides `dot()`, `norm()`,
+`squaredNorm()`, `setZero()`, `noalias()`, and arithmetic operators
+(`+=`, `-=`, `*=`) that dispatch to cuBLAS `axpy`, `nrm2`, `dot`, `scal`,
+and `geam`. These are the operations needed by iterative solvers.
+
+### `gpu::DeviceScalar<Scalar>`
+
+A device-resident scalar value. Reductions like `dot()`, `norm()`, and
+`squaredNorm()` return `DeviceScalar` instead of a host scalar, deferring
+the host synchronization until the value is actually needed:
+
+```cpp
+auto dot_val = d_x.dot(d_y);          // DeviceScalar -- no sync
+auto norm_sq = d_r.squaredNorm();      // DeviceScalar -- no sync
+Scalar alpha = dot_val / norm_sq;      // sync here (implicit conversion)
+d_x += alpha * d_p;                    // host scalar * DeviceMatrix (axpy)
+```
+
+Division between `DeviceScalar` values (real types only) is performed on
+device via NPP, avoiding extra synchronizations. Small device allocations
+(including `DeviceScalar`) are recycled through a thread-local
+`DeviceBufferPool` to avoid `cudaMalloc`/`cudaFree` overhead in tight loops.
+Pool contract: recycled pointers are safe for same-thread, same-stream reuse
+(the typical iterative-solver pattern, where one `gpu::Context` drives all
+work). Mixing pooled buffers across threads or CUDA streams without external
+synchronization is not supported.
+
 ### `gpu::Context`
 
 Every GPU operation needs a CUDA stream and library handles (cuBLAS eagerly,
-cuSOLVER lazily on first use). `gpu::Context` bundles these together. A
-single `Context` is not thread-safe -- use one per thread (or external
-synchronization), since the underlying cuBLAS and cuSOLVER handles are not
-thread-safe per handle.
+cuSOLVER / cuBLASLt / cuSPARSE lazily on first use). `gpu::Context` bundles
+these together. A single `Context` is not thread-safe -- use one per thread
+(or external synchronization), since the underlying NVIDIA library handles
+are not thread-safe per handle.
 
 For simple usage, you don't need to create one -- a per-thread default context
 is created lazily on first use:
@@ -129,16 +167,32 @@ d_C1.device(ctx1) = d_A1 * d_B1;   // runs on stream 1
 d_C2.device(ctx2) = d_A2 * d_B2;   // runs on stream 2 (concurrently)
 ```
 
+To integrate with existing CUDA code, borrow an existing stream:
+
+```cpp
+gpu::Context ctx(my_existing_stream);  // wraps stream, does not take ownership
+```
+
+To override the thread-local default (e.g., in CG where all ops share one
+context):
+
+```cpp
+gpu::Context ctx;
+gpu::Context::setThreadLocal(&ctx);    // all threadLocal() calls return ctx
+// ... GPU operations ...
+gpu::Context::setThreadLocal(nullptr); // restore lazy-created default
+```
+
 ### Linking
 
 The module is header-only, but each feature pulls in the corresponding NVIDIA
-library at link time. cuSOLVER is created lazily on first use, so a
-translation unit that only uses cuBLAS, cuFFT, cuSPARSE, or cuDSS does not
-need to link cuSOLVER:
+library at link time. cuSOLVER, cuBLASLt, and cuSPARSE are created lazily on
+first use, so a translation unit that only uses cuBLAS or cuFFT does not need
+to link the others:
 
 | Feature                                 | Link flags                |
 |-----------------------------------------|---------------------------|
-| `DeviceMatrix`, GEMM, TRSM, SYMM, SYRK  | `-lcublas`                |
+| `DeviceMatrix`, GEMM, TRSM, SYMM, SYRK  | `-lcublas -lcublasLt`     |
 | Dense solvers (LLT, LU, QR, SVD, EVD)   | `-lcusolver -lcublas`     |
 | FFT (`gpu::FFT`)                        | `-lcufft -lcublas`        |
 | SpMV / SpMM (`gpu::SparseContext`)      | `-lcusparse -lcublas`     |
@@ -164,7 +218,7 @@ d_C = d_A * d_B.transpose();
 
 // Scaled and accumulated
 d_C += 2.0 * d_A * d_B;             // alpha=2, beta=1
-d_C.device(ctx) -= d_A * d_B;       // alpha=-1, beta=1 (requires explicit context)
+d_C.device(ctx) -= d_A * d_B;       // alpha=-1, beta=1 (GEMM requires explicit context for -=)
 
 // Triangular solve (TRSM)
 d_X = d_A.triangularView<Lower>().solve(d_B);
@@ -176,6 +230,30 @@ d_C = d_A.selfadjointView<Lower>() * d_B;
 d_C.selfadjointView<Lower>().rankUpdate(d_A);  // C += A * A^H
 ```
 
+### BLAS Level-1 operations
+
+```cpp
+// Dot product and norms (return DeviceScalar -- no sync until read)
+auto dot_val = d_x.dot(d_y);          // cublasDdot / cublasCdotc
+auto norm_val = d_r.norm();            // cublasDnrm2
+double n = norm_val;                   // implicit conversion triggers sync
+
+// Vector arithmetic (cuBLAS axpy / geam)
+d_x += alpha * d_p;                    // axpy: x = x + alpha * p
+d_x -= alpha * d_p;                    // axpy: x = x - alpha * p
+d_x *= alpha;                          // scal: x = alpha * x
+d_r.setZero();                         // cudaMemsetAsync
+
+// DeviceScalar arithmetic (stays on device, real types only)
+auto alpha = absNew / dot_val;         // device-side division via NPP
+d_x += alpha * d_p;                    // DeviceScalar * DeviceMatrix (axpy with device pointer)
+
+// Matrix add/subtract (cuBLAS geam)
+gpu::DeviceMatrix<double> d_C = d_A + d_B;       // C = A + B
+d_C = d_A + 2.0 * d_B;                          // C = A + 2*B
+d_C = d_A - d_B;                                 // C = A - B
+```
+
 ### Dense solvers (cuSOLVER)
 
 **One-shot expression syntax** -- Convenient, re-factorizes each time:
@@ -214,6 +292,7 @@ gpu::SVD<double> svd;
 svd.compute(d_A, ComputeThinU | ComputeThinV);
 VectorXd S = svd.singularValues();   // downloads to host
 MatrixXd U = svd.matrixU();          // downloads to host
+MatrixXd V = svd.matrixV();          // V (matches JacobiSVD)
 MatrixXd VT = svd.matrixVT();        // V^T (matches cuSOLVER)
 
 // SVD: device-side views (no D2H transfer; svd must outlive the views)
@@ -299,16 +378,71 @@ MatrixXcf C = fft.inv2(B);         // 2D inverse (scaled by 1/(rows*cols))
 SparseMatrix<double> A = ...;
 VectorXd x = ...;
 
-gpu::SparseContext<double> ctx;
-VectorXd y = ctx.multiply(A, x);            // y = A * x
-VectorXd z = ctx.multiplyT(A, x);           // z = A^T * x
-ctx.multiply(A, x, y, 2.0, 1.0);            // y = 2*A*x + y
-ctx.multiply(A, x, y, 1.0, 0.0,             // y = A^H * x (Hermitian SpMV)
-             gpu::GpuOp::ConjTrans);
+// Host vectors (upload/download handled internally)
+gpu::SparseContext<double> spmv;
+VectorXd y = spmv.multiply(A, x);           // y = A * x
+VectorXd z = spmv.multiplyT(A, x);          // z = A^T * x
+spmv.multiply(A, x, y, 2.0, 1.0);           // y = 2*A*x + y
+spmv.multiply(A, x, y, 1.0, 0.0,            // y = A^H * x (Hermitian SpMV)
+              gpu::GpuOp::ConjTrans);
 
 // Multiple RHS (SpMM)
-MatrixXd Y = ctx.multiplyMat(A, X);                       // Y = A * X
-MatrixXd Z = ctx.multiplyMat(A, X, gpu::GpuOp::Trans);    // Z = A^T * X
+MatrixXd Y = spmv.multiplyMat(A, X);                      // Y = A * X
+MatrixXd Z = spmv.multiplyMat(A, X, gpu::GpuOp::Trans);   // Z = A^T * X
+
+// Device-resident SpMV (sparse matrix cached on device)
+gpu::Context ctx;
+gpu::SparseContext<double> spmv_dev(ctx);   // share gpu::Context for same-stream
+auto d_A = spmv_dev.deviceView(A);          // upload sparse matrix once
+d_y = d_A * d_x;                            // operator syntax, stays on device
+```
+
+### Eigen algorithm interop (example: Conjugate gradient)
+
+The BLAS-1 operators and `DeviceSparseView` make `DeviceMatrix` usable as a
+vector type in GPU implementations of algorithms like conjugate gradient.
+Conjugate gradient is the motivating example -- the GPU CG mirrors Eigen's
+`conjugate_gradient()` line for line, with only one host sync per iteration
+(the convergence check). All scalar intermediates (`alpha`, `beta`, `absNew`)
+stay on device as `DeviceScalar` values:
+
+```cpp
+gpu::Context ctx;
+gpu::Context::setThreadLocal(&ctx);
+gpu::SparseContext<double> spmv(ctx);
+auto mat = spmv.deviceView(A);              // upload sparse matrix once
+
+auto rhs = gpu::DeviceMatrix<double>::fromHost(b, ctx.stream());
+gpu::DeviceMatrix<double> x(n, 1);
+x.setZero();
+gpu::DeviceMatrix<double> residual(n, 1);
+residual.copyFrom(ctx, rhs);                // r = b (x=0)
+gpu::DeviceMatrix<double> p(n, 1);
+p.copyFrom(ctx, residual);                  // p = r
+gpu::DeviceMatrix<double> z(n, 1), tmp(n, 1);
+
+auto absNew = residual.dot(p);              // DeviceScalar -- no sync
+
+while (i < maxIters) {
+  tmp.noalias() = mat * p;                   // SpMV, device-resident
+
+  auto alpha = absNew / p.dot(tmp);          // DeviceScalar / DeviceScalar -- no sync
+
+  x += alpha * p;                            // DeviceScalar * DeviceMatrix axpy -- no sync
+  residual -= alpha * tmp;                   // DeviceScalar * DeviceMatrix axpy -- no sync
+
+  residualNorm2 = residual.squaredNorm();    // THE one sync per iteration
+  if (residualNorm2 < threshold) break;
+
+  z.copyFrom(ctx, residual);                 // no preconditioner: z = r
+  auto absOld = std::move(absNew);           // no sync, no alloc
+  absNew = residual.dot(z);                  // DeviceScalar -- no sync
+  auto beta = absNew / absOld;               // DeviceScalar / DeviceScalar -- no sync
+
+  p *= beta;                                 // DeviceScalar scal -- no sync
+  p += z;                                    // axpy -- no sync
+}
+MatrixXd result = x.toHost();
 ```
 
 ### Precision control
@@ -334,6 +468,7 @@ Mandatory sync points:
 - `fromHost()` -- Synchronizes to complete the upload before returning
 - `toHost()` / `HostTransfer::get()` -- Must deliver data to host
 - `info()` -- Must read the factorization status
+- `DeviceScalar` implicit conversion -- Downloads scalar from device
 
 **Cross-stream safety** is automatic. `DeviceMatrix` tracks write completion
 via CUDA events. When a matrix written on stream A is read on stream B, the
@@ -383,6 +518,18 @@ noted otherwise).
 | `X = A.triangularView<L>().solve(B)` | `cublasXtrsm` | side=L, uplo, diag=NonUnit |
 | `C = A.selfadjointView<L>() * B` | `cublasXsymm` / `cublasXhemm` | side=L, uplo |
 | `C.selfadjointView<L>().rankUpdate(A)` | `cublasXsyrk` / `cublasXherk` | uplo, trans=N |
+| `C = A + B` | `cublasXgeam` | alpha=1, beta=1 |
+| `C = A + alpha * B` | `cublasXgeam` | alpha=1, beta from scaled |
+| `C = A - B` | `cublasXgeam` | alpha=1, beta=-1 |
+| `C = A - alpha * B` | `cublasXgeam` | alpha=1, beta=-scaled |
+| `x += alpha * y` | `cublasXaxpy` | alpha (host scalar) |
+| `x += dAlpha * y` | `cublasXaxpy` | alpha (DeviceScalar, device pointer mode) |
+| `x -= alpha * y` | `cublasXaxpy` | alpha negated |
+| `x *= alpha` | `cublasXscal` | alpha (host or DeviceScalar) |
+| `x.dot(y)` | `cublasXdot` / `cublasXdotc` | returns `DeviceScalar` |
+| `x.norm()` | `cublasXnrm2` | returns `DeviceScalar<RealScalar>` |
+| `x.squaredNorm()` | `cublasXdot(x, x)` | returns `DeviceScalar<RealScalar>` |
+| `d_y = view * d_x` | `cusparseSpMV` | device-resident SpMV |
 
 ### `DeviceMatrix<Scalar>`
 
@@ -393,13 +540,14 @@ one column.
 ```cpp
 // Construction
 DeviceMatrix<Scalar>()                                   // Empty (0x0)
+DeviceMatrix<Scalar>(Index n)                            // Allocate column vector (n x 1)
 DeviceMatrix<Scalar>(rows, cols)                         // Allocate uninitialized
 
 // Upload / download / pointer adoption
 static DeviceMatrix fromHost(matrix, stream=nullptr)           // -> DeviceMatrix (syncs)
-static DeviceMatrix fromHostAsync(ptr, rows, cols, outerStride, s)  // -> DeviceMatrix (no sync, caller manages ptr lifetime)
-static DeviceMatrix adopt(Scalar* device_ptr, rows, cols)      // Owning wrapper over a raw device pointer
-static DeviceMatrix view(Scalar* device_ptr, rows, cols)       // Non-owning view (does not free on destruction)
+static DeviceMatrix fromHostAsync(ptr, rows, cols, stream)         // -> DeviceMatrix (no sync, caller manages ptr lifetime)
+static DeviceMatrix adopt(Scalar* device_ptr, rows, cols)          // Owning wrapper over a raw device pointer
+static DeviceMatrix view(Scalar* device_ptr, rows, cols)           // Non-owning view (does not free on destruction)
 PlainMatrix        toHost(stream=nullptr)                      // -> host Matrix (syncs)
 HostTransfer       toHostAsync(stream=nullptr)                 // -> HostTransfer future (no sync)
 DeviceMatrix       clone(stream=nullptr)                       // -> DeviceMatrix (D2D copy, async)
@@ -420,6 +568,50 @@ LuExpr             lu()                                  // -> .solve(d_B) -> De
 TriangularView     triangularView<UpLo>()                // -> .solve(d_B) -> DeviceMatrix (TRSM)
 SelfAdjointView    selfadjointView<UpLo>()               // -> * d_B (SYMM), .rankUpdate(d_A) (SYRK)
 Assignment   device(gpu::Context& ctx)                // Bind assignment to explicit stream
+DeviceMatrix&      noalias()                             // No-op (all ops are implicitly noalias)
+
+// BLAS Level-1 (all have overloads with explicit gpu::Context& parameter)
+DeviceScalar<Scalar>     dot(const DeviceMatrix& other)  // cuBLAS dot/dotc -> DeviceScalar
+DeviceScalar<RealScalar> norm()                          // cuBLAS nrm2 -> DeviceScalar
+DeviceScalar<RealScalar>  squaredNorm()                    // dot(self, self) -> DeviceScalar (no sync)
+void                     setZero()                       // cudaMemsetAsync
+void                     addScaled(gpu::Context&, Scalar alpha, const DeviceMatrix& x)  // this += alpha * x (axpy)
+void                     scale(gpu::Context&, Scalar alpha)                              // this *= alpha (scal)
+void                     copyFrom(gpu::Context&, const DeviceMatrix& other)              // this = other (D2D copy)
+DeviceMatrix& operator+=(const Scaled<DeviceMatrix>&)    // axpy; spelled `mat += alpha * other`
+DeviceMatrix& operator-=(const Scaled<DeviceMatrix>&)    // axpy negated; spelled `mat -= alpha * other`
+DeviceMatrix& operator+=(const DeviceMatrix&)            // cuBLAS axpy (alpha=1)
+DeviceMatrix& operator-=(const DeviceMatrix&)            // cuBLAS axpy (alpha=-1)
+DeviceMatrix& operator+=(const DeviceScaledDevice<Scalar>&)  // axpy with device scalar; spelled `mat += d_alpha * other`
+DeviceMatrix& operator-=(const DeviceScaledDevice<Scalar>&)  // negated; spelled `mat -= d_alpha * other`
+DeviceMatrix& operator*=(Scalar)                         // cuBLAS scal (host pointer mode)
+DeviceMatrix& operator*=(const DeviceScalar<Scalar>&)    // cuBLAS scal (device pointer mode, no host sync)
+DeviceMatrix  cwiseProduct(gpu::Context&, const DeviceMatrix&)            // NPP nppsMul (float/double only)
+void          cwiseProduct(gpu::Context&, const DeviceMatrix&, const DeviceMatrix&)  // in-place: this = a .* b
+
+// geam expressions (evaluated on assignment)
+DeviceMatrix& operator=(const DeviceAddExpr&)            // C = A + B, C = A + alpha*B, C = A - B, etc.
+```
+
+### `DeviceScalar<Scalar>`
+
+Device-resident scalar. Returned by `dot()`, `norm()`, and `squaredNorm()`.
+Implicit conversion to `Scalar` triggers `cudaStreamSynchronize` + download.
+
+```cpp
+DeviceScalar(cudaStream_t stream = nullptr)              // Allocate uninitialized
+DeviceScalar(Scalar host_val, cudaStream_t stream)       // Upload host value
+
+Scalar         get()                                     // Download (syncs stream)
+               operator Scalar()                         // Implicit conversion (syncs)
+Scalar*        devicePtr()                               // Raw device pointer
+cudaStream_t   stream()
+
+// Device-side arithmetic (no host sync, real types only)
+DeviceScalar   operator/(DeviceScalar, DeviceScalar)     // NPP nppsDiv
+DeviceScalar   operator/(Scalar, DeviceScalar)           // upload + div
+DeviceScalar   operator/(DeviceScalar, Scalar)           // upload + div
+DeviceScalar   operator-()                               // NPP nppsMulC(-1)
 ```
 
 ### `gpu::Context`
@@ -430,13 +622,17 @@ across threads.
 
 ```cpp
 gpu::Context()                                             // Creates dedicated stream + cuBLAS handle
-                                                           // (cuSOLVER handle created lazily on first
-                                                           // call to cusolverHandle())
+                                                           // (cuSOLVER / cuBLASLt / cuSPARSE handles
+                                                           // are created lazily on first use)
+gpu::Context(cudaStream_t stream)                          // Borrow existing stream (not owned)
 static gpu::Context& threadLocal()                         // Per-thread default (lazy-created)
+static void        setThreadLocal(gpu::Context* ctx)       // Override thread-local default (nullptr restores)
 
 cudaStream_t       stream()
 cublasHandle_t     cublasHandle()
 cusolverDnHandle_t cusolverHandle()                        // Lazy: creates the handle on first call
+cublasLtHandle_t   cublasLtHandle()                        // Lazy-initialized
+cusparseHandle_t   cusparseHandle()                        // Lazy-initialized
 ```
 
 Non-copyable, non-movable (owns library handles). Translation units that
@@ -510,6 +706,7 @@ gpu::SVD&            compute(const DeviceMatrix& d_A, unsigned options = Compute
 
 RealVector         singularValues()                      // -> host vector (syncs, downloads)
 PlainMatrix        matrixU()                             // -> host Matrix (syncs, downloads)
+PlainMatrix        matrixV()                             // -> host Matrix (V = VT^H, matches JacobiSVD)
 PlainMatrix        matrixVT()                            // -> host Matrix (syncs, downloads V^T)
 
 DeviceMatrix       d_singularValues()                    // -> DeviceMatrix view (zero-copy)
@@ -526,11 +723,11 @@ Index              rows() / cols()
 cudaStream_t       stream()
 ```
 
-**Note:** `singularValues()`, `matrixU()`, and `matrixVT()` download to host
-on each call. The `d_*` accessors return non-owning `DeviceMatrix` views into
-the solver's internal buffers; the `gpu::SVD` object must outlive any view
-derived from it. For wide matrices (m < n) the U/V^T views are owning (one
-`cublasXgeam` adjoint pass).
+**Note:** `singularValues()`, `matrixU()`, `matrixV()`, and `matrixVT()`
+download to host on each call. The `d_*` accessors return non-owning
+`DeviceMatrix` views into the solver's internal buffers; the `gpu::SVD` object
+must outlive any view derived from it. For wide matrices (m < n) the U/V^T
+views are owning (one `cublasXgeam` adjoint pass).
 
 ### `gpu::SelfAdjointEigenSolver<Scalar>` -- Eigendecomposition (cuSOLVER)
 
@@ -632,13 +829,16 @@ the input scalar type (complex vs real).
 
 ### `gpu::SparseContext<Scalar>` -- SpMV/SpMM (cuSPARSE)
 
-Accepts `SparseMatrix<Scalar, ColMajor>`. All methods accept host data and
-return host data. Matrix dimensions and nonzero count must fit in `int`
-(cuSPARSE limitation; debug builds assert).
+Accepts `SparseMatrix<Scalar, ColMajor>`. Host-input methods accept host data
+and return host data; device-input methods (`deviceView()`, `multiply(A, d_x,
+d_y)`) operate on `DeviceMatrix`. Matrix dimensions and nonzero count must fit
+in `int` (cuSPARSE limitation; debug builds assert).
 
 ```cpp
 gpu::SparseContext()                                       // Creates own stream + cuSPARSE handle
+gpu::SparseContext(gpu::Context& ctx)                        // Borrow gpu::Context for same-stream execution
 
+// Host data in/out
 DenseVector        multiply(A, x)                                       // y = A * x
 void               multiply(A, x, y, alpha=1, beta=0,                   // y = alpha*op(A)*x + beta*y
                      op=GpuOp::NoTrans)
@@ -646,9 +846,25 @@ DenseVector        multiplyT(A, x)                                      // y = A
 DenseVector        multiplyAdjoint(A, x)                                // y = A^H * x
 DenseMatrix        multiplyMat(A, X, op=GpuOp::NoTrans)                 // Y = op(A) * X (SpMM)
 
+// DeviceMatrix in/out (sparse matrix re-uploaded each call)
+void               multiply(A, d_x, d_y)                                // SpMV with device vectors
+void               multiply(A, d_x, d_y, alpha, beta, op)
+
+// Device-resident sparse matrix (upload once, reuse)
+DeviceSparseView   deviceView(A)                                        // Upload sparse matrix, return view
+
 cudaStream_t       stream()
 ```
 
+### `DeviceSparseView<Scalar>` -- Device-resident sparse matrix
+
+Returned by `gpu::SparseContext::deviceView()`. Holds a sparse matrix on device
+for repeated SpMV without re-uploading.
+
+```cpp
+SpMVExpr           operator*(const DeviceMatrix& d_x)    // d_y = view * d_x (evaluated on assignment)
+```
+
 ### Aliasing
 
 Unlike Eigen's `Matrix`, where omitting `.noalias()` triggers a copy to a
@@ -656,7 +872,9 @@ temporary, DeviceMatrix dispatches directly to NVIDIA library calls which have
 no built-in aliasing protection. All operations are implicitly noalias.
 The caller must ensure operands don't alias the destination for GEMM, TRSM,
 SYMM/HEMM, and SYRK/HERK. Debug builds assert on these violations before
-dispatching to cuBLAS.
+dispatching to cuBLAS. `geam` expressions (`d_C = d_A + alpha * d_B`) are
+safe with aliasing. The `.noalias()` method exists as a no-op for Eigen
+template compatibility.
 
 ## Future work
 
@@ -695,30 +913,32 @@ dispatching to cuBLAS.
 
 | File | Depends on | Contents |
 |------|-----------|----------|
-| `GpuSupport.h` | `<cuda_runtime.h>` | Error macro, `DeviceBuffer`, `cuda_data_type<>` |
-| `DeviceMatrix.h` | `GpuSupport.h` | `DeviceMatrix<>`, `HostTransfer<>` |
-| `DeviceExpr.h` | `DeviceMatrix.h` | GEMM expression wrappers |
+| `GpuSupport.h` | `<cuda_runtime.h>` | Error macro, `DeviceBuffer`, `DeviceBufferPool`, `cuda_data_type<>` |
+| `DeviceMatrix.h` | `GpuSupport.h` | `gpu::DeviceMatrix<>`, `gpu::HostTransfer<>` |
+| `DeviceExpr.h` | `DeviceMatrix.h` | GEMM, geam, and device-scalar expression wrappers |
 | `DeviceBlasExpr.h` | `DeviceMatrix.h` | TRSM, SYMM, SYRK expression wrappers |
 | `DeviceSolverExpr.h` | `DeviceMatrix.h` | Solver expression wrappers (LLT, LU) |
-| `DeviceDispatch.h` | all above | All dispatch functions + `Assignment` |
+| `DeviceScalar.h` | `GpuSupport.h`, `DeviceScalarOps.h` | `gpu::DeviceScalar<>` (device-resident scalar) |
+| `DeviceScalarOps.h` | `<npps_*.h>` | Scalar div/neg/cwiseProduct via NPP |
+| `DeviceDispatch.h` | all above | All dispatch functions, BLAS-1 out-of-line defs, `gpu::Assignment` |
 | `GpuContext.h` | `CuBlasSupport.h`, `CuSolverSupport.h` | `gpu::Context` |
-| `CuBlasSupport.h` | `GpuSupport.h`, `<cublas_v2.h>` | cuBLAS error macro, op/compute type maps |
+| `CuBlasSupport.h` | `GpuSupport.h`, `<cublas_v2.h>`, `<cublasLt.h>` | cuBLAS error macro, type-specific wrappers |
 | `CuSolverSupport.h` | `GpuSupport.h`, `<cusolverDn.h>` | cuSOLVER params, fill-mode mapping |
 | `GpuSolverContext.h` | `CuSolverSupport.h`, `CuBlasSupport.h` | Shared solver context (stream, handles, scratch) |
-| `GpuLLT.h` | `GpuSolverContext.h` | Cached dense Cholesky factorization |
-| `GpuLU.h` | `GpuSolverContext.h` | Cached dense LU factorization |
-| `GpuQR.h` | `GpuSolverContext.h` | Dense QR decomposition |
-| `GpuSVD.h` | `GpuSolverContext.h` | Dense SVD decomposition |
-| `GpuEigenSolver.h` | `GpuSolverContext.h` | Self-adjoint eigenvalue decomposition |
+| `GpuLLT.h` | `GpuSolverContext.h` | `gpu::LLT<>` -- Cached dense Cholesky factorization |
+| `GpuLU.h` | `GpuSolverContext.h` | `gpu::LU<>` -- Cached dense LU factorization |
+| `GpuQR.h` | `GpuSolverContext.h` | `gpu::QR<>` -- Dense QR decomposition |
+| `GpuSVD.h` | `GpuSolverContext.h` | `gpu::SVD<>` -- Dense SVD decomposition |
+| `GpuEigenSolver.h` | `GpuSolverContext.h` | `gpu::SelfAdjointEigenSolver<>` |
 | `CuFftSupport.h` | `GpuSupport.h`, `<cufft.h>` | cuFFT error macro, type-dispatch wrappers |
-| `GpuFFT.h` | `CuFftSupport.h`, `CuBlasSupport.h`, `GpuContext.h` | 1D/2D FFT with plan caching |
+| `GpuFFT.h` | `CuFftSupport.h`, `CuBlasSupport.h`, `GpuContext.h` | `gpu::FFT<>` -- 1D/2D FFT with plan caching |
 | `CuSparseSupport.h` | `GpuSupport.h`, `<cusparse.h>` | cuSPARSE error macro |
-| `GpuSparseContext.h` | `CuSparseSupport.h` | SpMV/SpMM via cuSPARSE |
+| `GpuSparseContext.h` | `CuSparseSupport.h` | `gpu::SparseContext<>`, `gpu::DeviceSparseView<>` |
 | `CuDssSupport.h` | `GpuSupport.h`, `<cudss.h>` | cuDSS error macro, type traits (optional) |
 | `GpuSparseSolverBase.h` | `CuDssSupport.h` | CRTP base for sparse solvers (optional) |
-| `GpuSparseLLT.h` | `GpuSparseSolverBase.h` | Sparse Cholesky via cuDSS (optional) |
-| `GpuSparseLDLT.h` | `GpuSparseSolverBase.h` | Sparse LDL^T via cuDSS (optional) |
-| `GpuSparseLU.h` | `GpuSparseSolverBase.h` | Sparse LU via cuDSS (optional) |
+| `GpuSparseLLT.h` | `GpuSparseSolverBase.h` | `gpu::SparseLLT<>` -- Sparse Cholesky via cuDSS (optional) |
+| `GpuSparseLDLT.h` | `GpuSparseSolverBase.h` | `gpu::SparseLDLT<>` -- Sparse LDL^T via cuDSS (optional) |
+| `GpuSparseLU.h` | `GpuSparseSolverBase.h` | `gpu::SparseLU<>` -- Sparse LU via cuDSS (optional) |
 
 ## Building and testing
 
@@ -731,7 +951,7 @@ cmake -G Ninja -B build -S . \
 
 cmake --build build --target cublas cusolver_llt cusolver_lu \
   cusolver_qr cusolver_svd cusolver_eigen \
-  device_matrix cufft cusparse_spmv
+  device_matrix cufft cusparse_spmv cg
 ctest --test-dir build -L gpu --output-on-failure
 
 # Sparse solvers (cuDSS -- separate install required)
@@ -746,10 +966,24 @@ ctest --test-dir build -R '^cudss_' --output-on-failure
 
 ## Future enhancements
 
-- **Device-resident sparse matrix-vector products.** `gpu::SparseContext`
-  currently operates on host vectors and matrices, uploading and downloading
-  on each call. The key missing piece is a `DeviceSparseView` that holds a
-  sparse matrix on device and supports operator syntax (`d_y = d_A * d_x`)
-  with `DeviceMatrix` operands -- keeping the entire SpMV/SpMM pipeline on
-  device. This is essential for iterative solvers and any workflow that chains
-  sparse and dense operations without returning to the host.
+- **Batched API (`DeviceBatchMatrix`).** A strided batch of N identical-size
+  matrices dispatching to cuBLAS/cuSOLVER batched APIs (`cublasDgemmBatched`,
+  `cusolverDnXpotrfBatched`, etc.). This enables robotics and model-predictive
+  control workloads where many small independent systems are solved in
+  parallel.
+- **cuTENSOR for Tensor module.** Replace the hand-written GPU tensor
+  contraction and reduction kernels (~2300 lines in
+  `TensorContractionGpu.h` / `TensorReductionGpu.h`) with cuTENSOR dispatch,
+  following the same library-dispatch pattern used by `unsupported/Eigen/GPU`.
+- **Unified/zero-copy memory for Jetson.** Use `cudaMallocManaged` or
+  `cudaHostAllocMapped` to eliminate `fromHost()` / `toHost()` copies on
+  integrated GPUs (Jetson) where CPU and GPU share DRAM.
+- **Device-side Eigen interop.** Bridge between host-side `DeviceMatrix`
+  dispatch and device-side Eigen expression templates (Core + Tensor) running
+  inside CUDA kernels. Raw-pointer + `Map` / `TensorMap` as the zero-copy
+  interop surface.
+- **Per-stream CUDA memory pools.** Currently all streams in a `GpuContext`
+  share the default device memory pool. Attaching a dedicated
+  `cudaMemPool_t` per stream (`cudaDeviceSetMempool` /
+  `cudaMallocFromPoolAsync`) can reduce cross-stream allocator contention for
+  workloads that fan out many concurrent solves.

unsupported/benchmarks/GPU/CMakeLists.txt

@@ -12,8 +12,8 @@
 # SPDX-FileCopyrightText: The Eigen Authors
 # SPDX-License-Identifier: MPL-2.0
 
-cmake_minimum_required(VERSION 3.17)
-project(EigenGpuBenchmarks CXX)
+cmake_minimum_required(VERSION 3.18)
+project(EigenGpuBenchmarks CXX CUDA)
 
 find_package(benchmark REQUIRED)
 find_package(CUDAToolkit REQUIRED)
@@ -58,3 +58,37 @@ eigen_add_gpu_benchmark(bench_batching_float bench_batching.cpp DEFINITIONS SCAL
 # FFT benchmarks: 1D/2D C2C, R2C, C2R throughput and plan reuse.
 eigen_add_gpu_benchmark(bench_fft bench_fft.cpp LIBRARIES CUDA::cufft)
 eigen_add_gpu_benchmark(bench_fft_double bench_fft.cpp LIBRARIES CUDA::cufft DEFINITIONS SCALAR=double)
+
+# CG sync overhead benchmark: host vs device pointer mode for reductions.
+# Uses CUDA kernels for device scalar arithmetic.
+add_executable(bench_cg_sync bench_cg_sync.cu)
+target_include_directories(bench_cg_sync PRIVATE
+  ${EIGEN_SOURCE_DIR}
+  ${CUDAToolkit_INCLUDE_DIRS})
+target_link_libraries(bench_cg_sync PRIVATE
+  benchmark::benchmark benchmark::benchmark_main
+  CUDA::cudart CUDA::cusolver CUDA::cublas CUDA::cusparse CUDA::npps CUDA::nppc)
+target_compile_options(bench_cg_sync PRIVATE $<$<COMPILE_LANGUAGE:CUDA>:-O3 --expt-relaxed-constexpr>)
+target_compile_definitions(bench_cg_sync PRIVATE EIGEN_USE_GPU)
+
+# GPU CG vs CPU CG comparison benchmark.
+add_executable(bench_cg_vs_cpu bench_cg_vs_cpu.cu)
+target_include_directories(bench_cg_vs_cpu PRIVATE
+  ${EIGEN_SOURCE_DIR}
+  ${CUDAToolkit_INCLUDE_DIRS})
+target_link_libraries(bench_cg_vs_cpu PRIVATE
+  benchmark::benchmark benchmark::benchmark_main
+  CUDA::cudart CUDA::cusolver CUDA::cublas CUDA::cusparse CUDA::npps CUDA::nppc)
+target_compile_options(bench_cg_vs_cpu PRIVATE $<$<COMPILE_LANGUAGE:CUDA>:-O3 --expt-relaxed-constexpr>)
+target_compile_definitions(bench_cg_vs_cpu PRIVATE EIGEN_USE_GPU)
+
+# Bundle Adjustment benchmark: GPU CG vs CPU CG on real BAL datasets.
+add_executable(bench_ba bench_ba.cu)
+target_include_directories(bench_ba PRIVATE
+  ${EIGEN_SOURCE_DIR}
+  ${CUDAToolkit_INCLUDE_DIRS})
+target_link_libraries(bench_ba PRIVATE
+  benchmark::benchmark
+  CUDA::cudart CUDA::cusolver CUDA::cublas CUDA::cusparse CUDA::npps CUDA::nppc)
+target_compile_options(bench_ba PRIVATE $<$<COMPILE_LANGUAGE:CUDA>:-O3 --expt-relaxed-constexpr>)
+target_compile_definitions(bench_ba PRIVATE EIGEN_USE_GPU)

unsupported/benchmarks/GPU/ba_results.md

@@ -0,0 +1,154 @@
+<!--
+SPDX-FileCopyrightText: The Eigen Authors
+SPDX-License-Identifier: MPL-2.0
+-->
+
+# Bundle Adjustment: GPU CG vs CPU CG Results
+
+Benchmark of Eigen's GPU CG pipeline on normal equations arising from bundle
+adjustment (BAL datasets). Compares CPU `ConjugateGradient` (Jacobi preconditioner)
+against GPU CG using `DeviceMatrix` + `GpuSparseContext` + `DeviceScalar`.
+
+## Hardware
+
+- **CPU**: Intel Core i7-13700HX (Raptor Lake, 12 cores / 24 threads, single thread for Eigen CG)
+- **GPU**: NVIDIA GeForce RTX 4070 Laptop GPU (Ada Lovelace, 4608 CUDA cores, 8 GB GDDR6)
+- **CUDA**: 13.2 / Driver 595.79
+- **OS**: Ubuntu 24.04 (WSL2, kernel 6.6.87)
+
+## Software
+
+- Eigen: `eigen-gpu-cg` branch
+- Google Benchmark 1.9.1
+- Compiler: nvcc 13.2 + g++ 13.3
+- Normal equations: H = J^T*J + I (Levenberg-Marquardt damping lambda=1.0)
+- CG tolerance: 1e-8, max iterations: 10000
+
+## Method
+
+For each BAL problem file:
+1. Parse the BAL file (cameras, 3D points, 2D observations)
+2. Compute the full Jacobian J using the BAL camera model (Rodrigues rotation +
+   perspective projection + radial distortion) with central finite differences
+3. Form the normal equations H = J^T*J + lambda*I (sparse, symmetric positive definite)
+4. Solve H*dx = -J^T*r using CG with Jacobi preconditioner on CPU and GPU
+5. Report wall-clock time (mean of 3 repetitions)
+
+GPU CG uses: `GpuSparseContext` for SpMV, `DeviceMatrix` for vectors,
+`DeviceScalar` with `CUBLAS_POINTER_MODE_DEVICE` for dot/norm reductions,
+in-place `cwiseProduct` via NPP for Jacobi preconditioner application,
+device-pointer-mode `scal` to avoid host sync on the beta update.
+
+## Results
+
+### Summary table
+
+| Dataset | Cameras | Points | Obs | H size | H nnz | CG iters | CPU CG (ms) | GPU CG (ms) | Speedup |
+|---------|---------|--------|-----|--------|-------|----------|-------------|-------------|---------|
+| Ladybug-49 | 49 | 7,776 | 31,843 | 23,769 | 1.8M | 4,421 | 4,006 | 1,152 | **3.5x** |
+| Ladybug-138 | 138 | 19,878 | 85,217 | 60,876 | 4.8M | 7,008 | 21,498 | 3,553 | **6.1x** |
+| Ladybug-646 | 646 | 73,584 | 327,297 | 226,566 | 18.4M | 10,000* | 123,727 | 14,268 | **8.7x** |
+| Dubrovnik-356 | 356 | 226,730 | 1,255,268 | 683,394 | 69.8M | 4,308 | 216,149 | 24,493 | **8.8x** |
+
+\* Hit 10,000 iteration cap (poorly conditioned problem). Both CPU and GPU
+hit the same cap, so timing comparison remains valid.
+
+### Profile breakdown (Ladybug-138, nsys)
+
+GPU kernel time is dominated by SpMV (91%). The remaining 9% is BLAS-1
+operations (dot, axpy, scal) and NPP element-wise ops (cwiseProduct).
+
+| Kernel | Time (ms) | % | Calls |
+|--------|-----------|---|-------|
+| cuSPARSE csrmv (SpMV) | 2507 | 91.3% | 7,006 |
+| cuBLAS dot | 92 | 3.4% | 21,020 |
+| cuBLAS axpy (device ptr) | 27 | 1.0% | 14,012 |
+| cuSPARSE partition | 19 | 0.7% | 7,006 |
+| NPP cwiseProduct | 16 + 13 | 1.1% | 14,011 + 7,006 |
+| cuBLAS axpy (host ptr) | 12 | 0.5% | 7,005 |
+| cuBLAS scal (device ptr) | 11 | 0.4% | 7,005 |
+| NPP scalar ops | 7 | 0.2% | 7,006 |
+
+### Optimizations applied
+
+Three profiling-driven optimizations reduced GPU CG time by **1.8x**
+(6.5s → 3.6s on Ladybug-138):
+
+1. **In-place `cwiseProduct`**: The Jacobi preconditioner apply
+   (`z = invdiag .* residual`) was allocating a new DeviceMatrix every
+   iteration. Added `z.cwiseProduct(ctx, a, b)` that reuses `z`'s buffer.
+   Reduced `cudaMalloc` calls from 7,053 to 23 (saving 2.3s).
+
+2. **`squaredNorm` via `dot(x,x)`**: cuBLAS `nrm2` uses a numerically
+   careful scaled-sum-of-squares algorithm (29µs/call). Replaced with
+   `dot(x,x)` (6.4µs/call) — 4.5x faster per call, saving ~320ms.
+
+3. **Device-pointer `scal`**: `p *= beta` was converting `DeviceScalar`
+   beta to host (triggering a stream sync), then calling host-pointer-mode
+   scal. Added `operator*=(DeviceScalar)` that uses device-pointer-mode
+   scal, eliminating one sync per iteration. Halved `cudaStreamSynchronize`
+   calls from 14K to 7K.
+
+### Observations
+
+1. **GPU speedup scales with problem size**: from 3.5x on small problems
+   (24K variables) to 8.8x on large problems (683K variables). This is
+   expected — larger problems have more parallelism for the GPU to exploit.
+
+2. **Iteration counts match**: CPU and GPU CG converge in the same number
+   of iterations (within 1%), confirming numerical equivalence.
+
+3. **Bottleneck is SpMV**: CG iteration time is dominated (91%) by the
+   sparse matrix-vector product on H. Further speedup requires either
+   faster SpMV (e.g., block-sparse formats) or algorithmic improvements
+   (Schur complement, better preconditioners).
+
+4. **Remaining overhead**: CUDA API calls (cudaMemcpyAsync for 8-byte
+   DeviceScalar transfers) account for ~50% of non-kernel time. Batching
+   multiple scalar reductions into a single transfer would help.
+
+5. **Jacobi preconditioner is weak for BA**: The Ladybug-646 problem does
+   not converge in 10K iterations. Ceres uses block Jacobi or Schur
+   complement preconditioners that would also benefit from GPU acceleration.
+
+### Scaling plot data
+
+```text
+# n        nnz_H       cpu_ms    gpu_ms    speedup
+23769      1793475     4006      1152      3.48
+60876      4791762     21498     3553      6.05
+226566     18387948    123727    14268     8.67
+683394     69827066    216149    24493     8.82
+```
+
+## BAL datasets
+
+Downloaded from http://grail.cs.washington.edu/projects/bal/
+
+| File | Source |
+|------|--------|
+| problem-49-7776-pre.txt | Ladybug sequence |
+| problem-138-19878-pre.txt | Ladybug sequence |
+| problem-646-73584-pre.txt | Ladybug sequence |
+| problem-356-226730-pre.txt | Dubrovnik reconstruction |
+
+## Reproducing
+
+```bash
+# Build
+cmake -G Ninja -B build-bench-gpu -S unsupported/benchmarks/GPU -DCMAKE_CUDA_ARCHITECTURES=89
+cmake --build build-bench-gpu --target bench_ba
+
+# Download BAL datasets
+wget http://grail.cs.washington.edu/projects/bal/data/ladybug/problem-49-7776-pre.txt.bz2
+wget http://grail.cs.washington.edu/projects/bal/data/ladybug/problem-138-19878-pre.txt.bz2
+wget http://grail.cs.washington.edu/projects/bal/data/ladybug/problem-646-73584-pre.txt.bz2
+wget http://grail.cs.washington.edu/projects/bal/data/dubrovnik/problem-356-226730-pre.txt.bz2
+bunzip2 *.bz2
+
+# Run (one at a time)
+BAL_FILE=problem-49-7776-pre.txt ./build-bench-gpu/bench_ba --benchmark_repetitions=3
+BAL_FILE=problem-138-19878-pre.txt ./build-bench-gpu/bench_ba --benchmark_repetitions=3
+BAL_FILE=problem-646-73584-pre.txt ./build-bench-gpu/bench_ba --benchmark_repetitions=3
+BAL_FILE=problem-356-226730-pre.txt ./build-bench-gpu/bench_ba --benchmark_repetitions=3
+```

unsupported/benchmarks/GPU/bench_ba.cu

@@ -0,0 +1,535 @@
+// Bundle Adjustment benchmark: GPU CG vs CPU CG on real BAL datasets.
+//
+// Tests Eigen's GPU CG pipeline (gpu::DeviceMatrix + gpu::SparseContext + DeviceScalar)
+// on the normal equations (J^T*J) arising from bundle adjustment problems.
+//
+// Reads a BAL (Bundle Adjustment in the Large) format file, computes the
+// Jacobian and residual, forms the normal equations H = J^T*J + lambda*I,
+// then solves H*dx = -J^T*r with both CPU and GPU conjugate gradients.
+//
+// BAL format: http://grail.cs.washington.edu/projects/bal/
+//
+// Usage:
+//   cmake --build build-bench-gpu --target bench_gpu_ba
+//
+//   # Download a BAL dataset (bz2-compressed):
+//   wget http://grail.cs.washington.edu/projects/bal/data/ladybug/problem-49-7776-pre.txt.bz2
+//   bunzip2 problem-49-7776-pre.txt.bz2
+//
+//   # Run on a specific problem:
+//   BAL_FILE=problem-49-7776-pre.txt ./build-bench-gpu/bench_gpu_ba
+//
+//   # Append results to the log:
+//   BAL_FILE=problem-49-7776-pre.txt ./build-bench-gpu/bench_gpu_ba \
+//     --benchmark_format=console 2>&1 | tee -a benchmarks/GPU/ba_results.log
+// SPDX-FileCopyrightText: The Eigen Authors
+// SPDX-License-Identifier: MPL-2.0
+
+#include <benchmark/benchmark.h>
+
+#include <Eigen/Sparse>
+#include <Eigen/IterativeLinearSolvers>
+#include <unsupported/Eigen/GPU>
+
+#include <cmath>
+#include <cstdio>
+#include <fstream>
+#include <string>
+#include <vector>
+
+using namespace Eigen;
+
+// ============================================================================
+// BAL problem data
+// ============================================================================
+
+struct BALProblem {
+  int num_cameras = 0;
+  int num_points = 0;
+  int num_observations = 0;
+
+  // Observations: (camera_idx, point_idx, observed_x, observed_y).
+  std::vector<int> camera_index;
+  std::vector<int> point_index;
+  std::vector<double> observations_x;
+  std::vector<double> observations_y;
+
+  // Camera parameters: 9 per camera (Rodrigues r[3], translation t[3], f, k1, k2).
+  std::vector<double> cameras;  // [num_cameras * 9]
+
+  // 3D points: 3 per point.
+  std::vector<double> points;  // [num_points * 3]
+
+  const double* camera(int i) const { return &cameras[i * 9]; }
+  const double* point(int i) const { return &points[i * 3]; }
+
+  bool load(const std::string& filename) {
+    std::ifstream in(filename);
+    if (!in) {
+      fprintf(stderr, "ERROR: Cannot open BAL file: %s\n", filename.c_str());
+      return false;
+    }
+
+    in >> num_cameras >> num_points >> num_observations;
+    if (!in || num_cameras <= 0 || num_points <= 0 || num_observations <= 0) {
+      fprintf(stderr, "ERROR: Invalid BAL header in %s\n", filename.c_str());
+      return false;
+    }
+
+    camera_index.resize(num_observations);
+    point_index.resize(num_observations);
+    observations_x.resize(num_observations);
+    observations_y.resize(num_observations);
+    for (int i = 0; i < num_observations; ++i) {
+      in >> camera_index[i] >> point_index[i] >> observations_x[i] >> observations_y[i];
+    }
+
+    cameras.resize(num_cameras * 9);
+    for (int i = 0; i < num_cameras * 9; ++i) {
+      in >> cameras[i];
+    }
+
+    points.resize(num_points * 3);
+    for (int i = 0; i < num_points * 3; ++i) {
+      in >> points[i];
+    }
+
+    if (!in) {
+      fprintf(stderr, "ERROR: Truncated BAL file: %s\n", filename.c_str());
+      return false;
+    }
+
+    fprintf(stderr, "Loaded BAL: %d cameras, %d points, %d observations\n", num_cameras, num_points, num_observations);
+    return true;
+  }
+};
+
+// ============================================================================
+// Camera projection model (BAL convention)
+// ============================================================================
+
+// Rodrigues rotation: rotate point X by axis-angle vector omega.
+static void rodrigues_rotate(const double* omega, const double* X, double* result) {
+  double theta2 = omega[0] * omega[0] + omega[1] * omega[1] + omega[2] * omega[2];
+  if (theta2 > 1e-30) {
+    double theta = std::sqrt(theta2);
+    double costh = std::cos(theta);
+    double sinth = std::sin(theta);
+    double k = (1.0 - costh) / theta2;
+
+    // Cross product omega x X.
+    double wx = omega[1] * X[2] - omega[2] * X[1];
+    double wy = omega[2] * X[0] - omega[0] * X[2];
+    double wz = omega[0] * X[1] - omega[1] * X[0];
+
+    // Dot product omega . X.
+    double dot = omega[0] * X[0] + omega[1] * X[1] + omega[2] * X[2];
+
+    result[0] = X[0] * costh + wx * (sinth / theta) + omega[0] * dot * k;
+    result[1] = X[1] * costh + wy * (sinth / theta) + omega[1] * dot * k;
+    result[2] = X[2] * costh + wz * (sinth / theta) + omega[2] * dot * k;
+  } else {
+    // Small angle: R ≈ I + [omega]×.
+    result[0] = X[0] + omega[1] * X[2] - omega[2] * X[1];
+    result[1] = X[1] + omega[2] * X[0] - omega[0] * X[2];
+    result[2] = X[2] + omega[0] * X[1] - omega[1] * X[0];
+  }
+}
+
+// Project a 3D point through a camera, returning the 2D residual.
+// camera: [r0,r1,r2, t0,t1,t2, f, k1, k2]
+// point:  [X, Y, Z]
+// observed: [ox, oy]
+// residual: [rx, ry] = projected - observed
+static void project(const double* camera, const double* point, const double* observed, double* residual) {
+  // Rotate.
+  double P[3];
+  rodrigues_rotate(camera, point, P);
+
+  // Translate.
+  P[0] += camera[3];
+  P[1] += camera[4];
+  P[2] += camera[5];
+
+  // Normalize (BAL convention: negative z).
+  double xp = -P[0] / P[2];
+  double yp = -P[1] / P[2];
+
+  // Radial distortion.
+  double r2 = xp * xp + yp * yp;
+  double distortion = 1.0 + camera[7] * r2 + camera[8] * r2 * r2;
+
+  // Apply focal length.
+  double predicted_x = camera[6] * distortion * xp;
+  double predicted_y = camera[6] * distortion * yp;
+
+  residual[0] = predicted_x - observed[0];
+  residual[1] = predicted_y - observed[1];
+}
+
+// ============================================================================
+// Jacobian computation (numerical differentiation)
+// ============================================================================
+
+// Compute the 2x9 Jacobian block w.r.t. camera params and 2x3 block w.r.t.
+// point coords for a single observation, using central finite differences.
+static void compute_jacobian_block(const double* camera, const double* point, const double* observed,
+                                   double* J_cam,    // 2x9, row-major
+                                   double* J_point)  // 2x3, row-major
+{
+  constexpr double eps = 1e-8;
+
+  // Camera parameters (9).
+  double cam_pert[9];
+  std::copy(camera, camera + 9, cam_pert);
+  for (int j = 0; j < 9; ++j) {
+    double orig = cam_pert[j];
+    double rp[2], rm[2];
+
+    cam_pert[j] = orig + eps;
+    project(cam_pert, point, observed, rp);
+    cam_pert[j] = orig - eps;
+    project(cam_pert, point, observed, rm);
+    cam_pert[j] = orig;
+
+    J_cam[0 * 9 + j] = (rp[0] - rm[0]) / (2.0 * eps);
+    J_cam[1 * 9 + j] = (rp[1] - rm[1]) / (2.0 * eps);
+  }
+
+  // Point coordinates (3).
+  double pt_pert[3];
+  std::copy(point, point + 3, pt_pert);
+  for (int j = 0; j < 3; ++j) {
+    double orig = pt_pert[j];
+    double rp[2], rm[2];
+
+    pt_pert[j] = orig + eps;
+    project(camera, pt_pert, observed, rp);
+    pt_pert[j] = orig - eps;
+    project(camera, pt_pert, observed, rm);
+    pt_pert[j] = orig;
+
+    J_point[0 * 3 + j] = (rp[0] - rm[0]) / (2.0 * eps);
+    J_point[1 * 3 + j] = (rp[1] - rm[1]) / (2.0 * eps);
+  }
+}
+
+// ============================================================================
+// Build normal equations: H = J^T*J + lambda*I, g = -J^T*r
+// ============================================================================
+
+struct NormalEquations {
+  SparseMatrix<double, ColMajor, int> H;
+  VectorXd g;
+  VectorXd residual;
+  double residual_norm;
+  int jacobian_rows;
+  int jacobian_cols;
+  long jacobian_nnz;
+};
+
+static NormalEquations build_normal_equations(const BALProblem& problem, double lambda = 1.0) {
+  const int num_cam_params = problem.num_cameras * 9;
+  const int num_pt_params = problem.num_points * 3;
+  const int num_params = num_cam_params + num_pt_params;
+  const int num_residuals = problem.num_observations * 2;
+
+  fprintf(stderr, "Building Jacobian: %d x %d, %ld nonzeros\n", num_residuals, num_params,
+          (long)problem.num_observations * 24);
+
+  // Build J as a triplet list.
+  using Triplet = Eigen::Triplet<double>;
+  std::vector<Triplet> triplets;
+  triplets.reserve(problem.num_observations * 24);  // 2 rows × 12 nonzeros = 24 entries per obs
+
+  VectorXd residual(num_residuals);
+
+  for (int obs = 0; obs < problem.num_observations; ++obs) {
+    int ci = problem.camera_index[obs];
+    int pi = problem.point_index[obs];
+    double observed[2] = {problem.observations_x[obs], problem.observations_y[obs]};
+
+    // Compute residual.
+    double r[2];
+    project(problem.camera(ci), problem.point(pi), observed, r);
+    residual[obs * 2 + 0] = r[0];
+    residual[obs * 2 + 1] = r[1];
+
+    // Compute Jacobian blocks.
+    double J_cam[18], J_pt[6];  // 2x9 and 2x3
+    compute_jacobian_block(problem.camera(ci), problem.point(pi), observed, J_cam, J_pt);
+
+    // Insert camera block: rows [2*obs, 2*obs+1], cols [9*ci, 9*ci+8].
+    for (int row = 0; row < 2; ++row) {
+      for (int col = 0; col < 9; ++col) {
+        double val = J_cam[row * 9 + col];
+        if (val != 0.0) {
+          triplets.emplace_back(obs * 2 + row, ci * 9 + col, val);
+        }
+      }
+    }
+
+    // Insert point block: rows [2*obs, 2*obs+1], cols [num_cam_params + 3*pi, ...].
+    for (int row = 0; row < 2; ++row) {
+      for (int col = 0; col < 3; ++col) {
+        double val = J_pt[row * 3 + col];
+        if (val != 0.0) {
+          triplets.emplace_back(obs * 2 + row, num_cam_params + pi * 3 + col, val);
+        }
+      }
+    }
+  }
+
+  // Build sparse Jacobian.
+  SparseMatrix<double, ColMajor, int> J(num_residuals, num_params);
+  J.setFromTriplets(triplets.begin(), triplets.end());
+
+  fprintf(stderr, "Jacobian: %dx%d, nnz=%ld\n", (int)J.rows(), (int)J.cols(), (long)J.nonZeros());
+
+  // Form normal equations: H = J^T*J + lambda*I.
+  SparseMatrix<double, ColMajor, int> H = (J.transpose() * J).pruned();
+
+  // Add Levenberg-Marquardt damping.
+  for (int i = 0; i < num_params; ++i) {
+    H.coeffRef(i, i) += lambda;
+  }
+  H.makeCompressed();
+
+  // Gradient: g = -J^T * r.
+  VectorXd g = -(J.transpose() * residual);
+
+  double rnorm = residual.norm();
+  fprintf(stderr, "Normal equations: H is %dx%d, nnz=%ld, |r|=%.6e\n", (int)H.rows(), (int)H.cols(), (long)H.nonZeros(),
+          rnorm);
+
+  return {std::move(H), std::move(g), std::move(residual), rnorm, num_residuals, num_params, (long)J.nonZeros()};
+}
+
+// ============================================================================
+// Global problem state (loaded once before benchmarks run)
+// ============================================================================
+
+static BALProblem g_problem;
+static NormalEquations g_neq;
+static bool g_loaded = false;
+
+static void ensure_loaded() {
+  if (g_loaded) return;
+
+  const char* bal_file = std::getenv("BAL_FILE");
+  if (!bal_file) {
+    fprintf(stderr,
+            "ERROR: Set BAL_FILE environment variable to a BAL problem file.\n"
+            "  Download from: http://grail.cs.washington.edu/projects/bal/\n"
+            "  Example:\n"
+            "    wget http://grail.cs.washington.edu/projects/bal/data/ladybug/"
+            "problem-49-7776-pre.txt.bz2\n"
+            "    bunzip2 problem-49-7776-pre.txt.bz2\n"
+            "    BAL_FILE=problem-49-7776-pre.txt ./build-bench-gpu/bench_gpu_ba\n");
+    std::exit(1);
+  }
+
+  if (!g_problem.load(bal_file)) {
+    std::exit(1);
+  }
+
+  g_neq = build_normal_equations(g_problem);
+  g_loaded = true;
+}
+
+// ============================================================================
+// CPU CG benchmark
+// ============================================================================
+
+static void BM_BA_CPU_CG(benchmark::State& state) {
+  ensure_loaded();
+  const auto& H = g_neq.H;
+  const auto& g = g_neq.g;
+
+  ConjugateGradient<SparseMatrix<double, ColMajor, int>, Lower | Upper> cg;
+  cg.setMaxIterations(10000);
+  cg.setTolerance(1e-8);
+  cg.compute(H);
+
+  int last_iters = 0;
+  double last_error = 0;
+  for (auto _ : state) {
+    VectorXd dx = cg.solve(g);
+    benchmark::DoNotOptimize(dx.data());
+    last_iters = cg.iterations();
+    last_error = cg.error();
+  }
+
+  state.counters["n"] = H.rows();
+  state.counters["nnz"] = H.nonZeros();
+  state.counters["iters"] = last_iters;
+  state.counters["error"] = last_error;
+  state.counters["cameras"] = g_problem.num_cameras;
+  state.counters["points"] = g_problem.num_points;
+  state.counters["observations"] = g_problem.num_observations;
+}
+
+// ============================================================================
+// GPU CG benchmark (with Jacobi preconditioner)
+// ============================================================================
+
+static void cuda_warmup() {
+  static bool done = false;
+  if (!done) {
+    void* p;
+    cudaMalloc(&p, 1);
+    cudaFree(p);
+    done = true;
+  }
+}
+
+static void BM_BA_GPU_CG(benchmark::State& state) {
+  ensure_loaded();
+  cuda_warmup();
+
+  const auto& H = g_neq.H;
+  const auto& g = g_neq.g;
+  const Index n = H.rows();
+
+  // Extract inverse diagonal (Jacobi preconditioner).
+  using SpMat = SparseMatrix<double, ColMajor, int>;
+  VectorXd invdiag(n);
+  for (Index j = 0; j < H.outerSize(); ++j) {
+    SpMat::InnerIterator it(H, j);
+    while (it && it.index() != j) ++it;
+    if (it && it.index() == j && it.value() != 0.0)
+      invdiag(j) = 1.0 / it.value();
+    else
+      invdiag(j) = 1.0;
+  }
+
+  // Set up GPU context and upload data.
+  gpu::Context ctx;
+  gpu::Context::setThreadLocal(&ctx);
+  gpu::SparseContext<double> spmv_ctx(ctx);
+  auto mat = spmv_ctx.deviceView(H);
+  auto d_invdiag = gpu::DeviceMatrix<double>::fromHost(invdiag, ctx.stream());
+  auto d_g = gpu::DeviceMatrix<double>::fromHost(g, ctx.stream());
+
+  int last_iters = 0;
+  double last_error = 0;
+
+  for (auto _ : state) {
+    gpu::DeviceMatrix<double> d_x(n, 1);
+    d_x.setZero(ctx);
+    gpu::DeviceMatrix<double> residual(n, 1);
+    residual.copyFrom(ctx, d_g);
+
+    double rhsNorm2 = d_g.squaredNorm(ctx);
+    double threshold = 1e-8 * 1e-8 * rhsNorm2;
+    double residualNorm2 = residual.squaredNorm(ctx);
+
+    gpu::DeviceMatrix<double> p = d_invdiag.cwiseProduct(ctx, residual);
+    gpu::DeviceMatrix<double> z(n, 1), tmp(n, 1);
+
+    auto absNew = residual.dot(ctx, p);
+    Index i = 0;
+    Index maxIters = 10000;
+    while (i < maxIters) {
+      tmp.noalias() = mat * p;
+      auto alpha = absNew / p.dot(ctx, tmp);
+      d_x += alpha * p;
+      residual -= alpha * tmp;
+
+      residualNorm2 = residual.squaredNorm(ctx);
+      if (residualNorm2 < threshold) break;
+
+      z.cwiseProduct(ctx, d_invdiag, residual);  // in-place, no allocation
+      auto absOld = std::move(absNew);
+      absNew = residual.dot(ctx, z);
+      auto beta = absNew / absOld;
+
+      p *= beta;  // device-pointer scal, no host sync
+      p += z;
+      i++;
+    }
+    benchmark::DoNotOptimize(d_x.data());
+    last_iters = i;
+    last_error = std::sqrt(residualNorm2 / rhsNorm2);
+  }
+
+  gpu::Context::setThreadLocal(nullptr);
+
+  state.counters["n"] = n;
+  state.counters["nnz"] = H.nonZeros();
+  state.counters["iters"] = last_iters;
+  state.counters["error"] = last_error;
+  state.counters["cameras"] = g_problem.num_cameras;
+  state.counters["points"] = g_problem.num_points;
+  state.counters["observations"] = g_problem.num_observations;
+}
+
+// ============================================================================
+// CPU CG with Jacobi preconditioner (apples-to-apples comparison)
+// ============================================================================
+
+static void BM_BA_CPU_CG_Jacobi(benchmark::State& state) {
+  ensure_loaded();
+  const auto& H = g_neq.H;
+  const auto& g = g_neq.g;
+
+  // Eigen's DiagonalPreconditioner is effectively Jacobi.
+  ConjugateGradient<SparseMatrix<double, ColMajor, int>, Lower | Upper> cg;
+  cg.setMaxIterations(10000);
+  cg.setTolerance(1e-8);
+  cg.compute(H);
+
+  int last_iters = 0;
+  double last_error = 0;
+  for (auto _ : state) {
+    VectorXd dx = cg.solve(g);
+    benchmark::DoNotOptimize(dx.data());
+    last_iters = cg.iterations();
+    last_error = cg.error();
+  }
+
+  state.counters["n"] = H.rows();
+  state.counters["nnz"] = H.nonZeros();
+  state.counters["iters"] = last_iters;
+  state.counters["error"] = last_error;
+}
+
+// ============================================================================
+// Register benchmarks
+// ============================================================================
+
+BENCHMARK(BM_BA_CPU_CG)->Unit(benchmark::kMillisecond);
+BENCHMARK(BM_BA_CPU_CG_Jacobi)->Unit(benchmark::kMillisecond);
+BENCHMARK(BM_BA_GPU_CG)->Unit(benchmark::kMillisecond);
+
+// ============================================================================
+// Custom main: print summary after benchmarks
+// ============================================================================
+
+int main(int argc, char** argv) {
+  benchmark::Initialize(&argc, argv);
+
+  // Print problem info before benchmarks.
+  const char* bal_file = std::getenv("BAL_FILE");
+  if (bal_file) {
+    ensure_loaded();
+    fprintf(stderr,
+            "\n"
+            "=== Bundle Adjustment GPU CG Benchmark ===\n"
+            "BAL file:      %s\n"
+            "Cameras:       %d\n"
+            "Points:        %d\n"
+            "Observations:  %d\n"
+            "J size:        %d x %d, nnz=%ld\n"
+            "H size:        %d x %d, nnz=%ld\n"
+            "|residual|:    %.6e\n"
+            "==========================================\n\n",
+            bal_file, g_problem.num_cameras, g_problem.num_points, g_problem.num_observations, g_neq.jacobian_rows,
+            g_neq.jacobian_cols, g_neq.jacobian_nnz, (int)g_neq.H.rows(), (int)g_neq.H.cols(), (long)g_neq.H.nonZeros(),
+            g_neq.residual_norm);
+  }
+
+  benchmark::RunSpecifiedBenchmarks();
+  benchmark::Shutdown();
+  return 0;
+}

unsupported/benchmarks/GPU/bench_cg_sync.cu

@@ -0,0 +1,293 @@
+// Benchmark: GPU Conjugate Gradient via gpu::DeviceMatrix operators.
+//
+// Shows the path to running Eigen's CG on GPU with minimal code changes.
+// The gpu::DeviceMatrix benchmark mirrors Eigen's conjugate_gradient() line-by-line.
+// A raw cuBLAS device-pointer-mode implementation is included as a lower bound.
+//
+// The only change needed in Eigen's CG template to support gpu::DeviceMatrix:
+//   Line 34:  typedef Dest VectorType;  (instead of Matrix<Scalar, Dynamic, 1>)
+//
+// Usage:
+//   cmake --build build-bench-gpu --target bench_gpu_cg_sync
+//   ./build-bench-gpu/bench_gpu_cg_sync
+// SPDX-FileCopyrightText: The Eigen Authors
+// SPDX-License-Identifier: MPL-2.0
+
+#include <benchmark/benchmark.h>
+
+#include <Eigen/Sparse>
+#include <unsupported/Eigen/GPU>
+#include <cusparse.h>
+
+using namespace Eigen;
+
+using Scalar = double;
+using RealScalar = double;
+using Vec = Matrix<Scalar, Dynamic, 1>;
+using SpMat = SparseMatrix<Scalar, ColMajor, int>;
+
+static SpMat make_spd(Index n) {
+  SpMat A(n, n);
+  A.reserve(VectorXi::Constant(n, 3));
+  for (Index i = 0; i < n; ++i) {
+    A.insert(i, i) = 4.0;
+    if (i > 0) A.insert(i, i - 1) = -1.0;
+    if (i < n - 1) A.insert(i, i + 1) = -1.0;
+  }
+  A.makeCompressed();
+  return A;
+}
+
+static void cuda_warmup() {
+  static bool done = false;
+  if (!done) {
+    void* p;
+    cudaMalloc(&p, 1);
+    cudaFree(p);
+    done = true;
+  }
+}
+
+// ==========================================================================
+// GPU CG using gpu::DeviceMatrix operators — mirrors Eigen's conjugate_gradient()
+// ==========================================================================
+//
+// Compare with Eigen/src/IterativeLinearSolvers/ConjugateGradient.h lines 29-84.
+// Left column: Eigen CG code.  Right column: this benchmark.
+//
+//   Eigen CG                              GPU CG (this benchmark)
+//   --------                              -----------------------
+//   VectorType residual = rhs - mat * x;  residual.copyFrom(ctx, rhs);  [x=0 so r=b]
+//   RealScalar rhsNorm2 = rhs.sqNorm();   RealScalar rhsNorm2 = rhs.squaredNorm();
+//   ...
+//   tmp.noalias() = mat * p;              tmp.noalias() = mat * p;  [identical]
+//   Scalar alpha = absNew / p.dot(tmp);   Scalar alpha = absNew / p.dot(tmp);  [identical]
+//   x += alpha * p;                       x += alpha * p;  [identical]
+//   residual -= alpha * tmp;              residual -= alpha * tmp;  [identical]
+//   residualNorm2 = residual.sqNorm();    residualNorm2 = residual.squaredNorm();  [identical]
+//   ...
+//   p = z + beta * p;                     p *= beta; p += z;  [equivalent, no alloc]
+
+static void BM_CG_DeviceMatrixOps(benchmark::State& state) {
+  cuda_warmup();
+  const Index n = state.range(0);
+
+  SpMat A = make_spd(n);
+  Vec b = Vec::Random(n);
+
+  // One shared context: SpMV + BLAS-1 on same stream, zero event overhead.
+  gpu::Context ctx;
+  gpu::Context::setThreadLocal(&ctx);
+  gpu::SparseContext<Scalar> spmv(ctx);
+  auto mat = spmv.deviceView(A);
+
+  // Upload RHS once.
+  auto rhs = gpu::DeviceMatrix<Scalar>::fromHost(b, ctx.stream());
+
+  for (auto _ : state) {
+    // --- Eigen CG lines 34-63: initialization ---
+    //   typedef Dest VectorType;                               // GPU CHANGE: was Matrix<Scalar,Dynamic,1>
+    //   VectorType residual = rhs - mat * x;                   // x=0, so residual = rhs
+    gpu::DeviceMatrix<Scalar> x(n, 1);
+    x.setZero();
+    gpu::DeviceMatrix<Scalar> residual(n, 1);
+    residual.copyFrom(ctx, rhs);
+
+    //   RealScalar rhsNorm2 = rhs.squaredNorm();
+    RealScalar rhsNorm2 = rhs.squaredNorm();
+    if (rhsNorm2 == 0) continue;
+
+    RealScalar tol = 1e-10;
+    const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
+    RealScalar threshold = numext::maxi(RealScalar(tol * tol * rhsNorm2), considerAsZero);
+
+    //   RealScalar residualNorm2 = residual.squaredNorm();
+    RealScalar residualNorm2 = residual.squaredNorm();
+    if (residualNorm2 < threshold) continue;
+
+    //   VectorType p(n);
+    //   p = precond.solve(residual);                           // no preconditioner: p = residual
+    gpu::DeviceMatrix<Scalar> p(n, 1);
+    p.copyFrom(ctx, residual);
+
+    //   VectorType z(n), tmp(n);
+    gpu::DeviceMatrix<Scalar> z(n, 1), tmp(n, 1);
+
+    //   auto absNew = numext::real(residual.dot(p));
+    //   gpu::DeviceScalar — stays on device, no sync.
+    auto absNew = residual.dot(p);  // gpu::DeviceScalar, no sync
+
+    //   while (i < maxIters) {
+    Index maxIters = 200;
+    Index i = 0;
+    while (i < maxIters) {
+      //     tmp.noalias() = mat * p;
+      tmp.noalias() = mat * p;  // SpMV, device-resident
+
+      //     auto alpha = absNew / p.dot(tmp);
+      //   gpu::DeviceScalar / gpu::DeviceScalar → device kernel, no sync!
+      auto alpha = absNew / p.dot(tmp);  // gpu::DeviceScalar, no sync
+
+      //     x += alpha * p;
+      //   gpu::DeviceScalar * gpu::DeviceMatrix → device-pointer axpy, no sync!
+      x += alpha * p;
+
+      //     residual -= alpha * tmp;
+      residual -= alpha * tmp;  // device-pointer axpy, no sync
+
+      //     residualNorm2 = residual.squaredNorm();
+      residualNorm2 = residual.squaredNorm();  // THE one sync per iteration
+
+      //     if (residualNorm2 < threshold) break;
+      if (residualNorm2 < threshold) break;
+
+      //     z = precond.solve(residual);
+      z.copyFrom(ctx, residual);  // no preconditioner
+
+      //     auto absOld = std::move(absNew);
+      auto absOld = std::move(absNew);  // no sync, no alloc
+
+      //     absNew = numext::real(residual.dot(z));
+      absNew = residual.dot(z);  // gpu::DeviceScalar, no sync
+
+      //     auto beta = absNew / absOld;
+      //   gpu::DeviceScalar / gpu::DeviceScalar → device kernel, no sync!
+      auto beta = absNew / absOld;  // gpu::DeviceScalar, no sync
+
+      //     p = z + beta * p;
+      p *= beta;  // device-pointer scal, no host sync
+      p += z;
+
+      i++;
+    }
+  }
+
+  gpu::Context::setThreadLocal(nullptr);
+  state.SetItemsProcessed(state.iterations() * 200);
+}
+
+BENCHMARK(BM_CG_DeviceMatrixOps)->RangeMultiplier(4)->Range(1 << 10, 1 << 20);
+
+// ==========================================================================
+// Raw cuBLAS device-pointer-mode CG (1 sync/iter) — performance lower bound
+// ==========================================================================
+
+__global__ void scalar_div_kernel(const Scalar* a, const Scalar* b, Scalar* out) { *out = *a / *b; }
+__global__ void scalar_neg_kernel(const Scalar* in, Scalar* out) { *out = -(*in); }
+
+static void BM_CG_DevicePointerMode(benchmark::State& state) {
+  cuda_warmup();
+  const Index n = state.range(0);
+  const int maxIters = 200;
+
+  SpMat A = make_spd(n);
+  Vec b = Vec::Random(n);
+
+  cudaStream_t stream;
+  cudaStreamCreate(&stream);
+  cublasHandle_t cublas;
+  cublasCreate(&cublas);
+  cublasSetStream(cublas, stream);
+
+  cusparseHandle_t cusparse;
+  cusparseCreate(&cusparse);
+  cusparseSetStream(cusparse, stream);
+
+  internal::DeviceBuffer d_outer((n + 1) * sizeof(int));
+  internal::DeviceBuffer d_inner(A.nonZeros() * sizeof(int));
+  internal::DeviceBuffer d_vals(A.nonZeros() * sizeof(Scalar));
+  cudaMemcpy(d_outer.ptr, A.outerIndexPtr(), (n + 1) * sizeof(int), cudaMemcpyHostToDevice);
+  cudaMemcpy(d_inner.ptr, A.innerIndexPtr(), A.nonZeros() * sizeof(int), cudaMemcpyHostToDevice);
+  cudaMemcpy(d_vals.ptr, A.valuePtr(), A.nonZeros() * sizeof(Scalar), cudaMemcpyHostToDevice);
+
+  cusparseSpMatDescr_t matA;
+  cusparseCreateCsc(&matA, n, n, A.nonZeros(), d_outer.ptr, d_inner.ptr, d_vals.ptr, CUSPARSE_INDEX_32I,
+                    CUSPARSE_INDEX_32I, CUSPARSE_INDEX_BASE_ZERO, CUDA_R_64F);
+
+  internal::DeviceBuffer d_tmp_buf(n * sizeof(Scalar));
+  cusparseDnVecDescr_t tmp_x, tmp_y;
+  cusparseCreateDnVec(&tmp_x, n, d_tmp_buf.ptr, CUDA_R_64F);
+  cusparseCreateDnVec(&tmp_y, n, d_tmp_buf.ptr, CUDA_R_64F);
+  Scalar spmv_alpha = 1.0, spmv_beta = 0.0;
+  size_t ws_size = 0;
+  cusparseSpMV_bufferSize(cusparse, CUSPARSE_OPERATION_NON_TRANSPOSE, &spmv_alpha, matA, tmp_x, &spmv_beta, tmp_y,
+                          CUDA_R_64F, CUSPARSE_SPMV_ALG_DEFAULT, &ws_size);
+  internal::DeviceBuffer d_workspace(ws_size);
+  cusparseDestroyDnVec(tmp_x);
+  cusparseDestroyDnVec(tmp_y);
+
+  internal::DeviceBuffer d_x(n * sizeof(Scalar)), d_r(n * sizeof(Scalar));
+  internal::DeviceBuffer d_p(n * sizeof(Scalar)), d_tmp(n * sizeof(Scalar));
+  internal::DeviceBuffer d_b(n * sizeof(Scalar));
+  internal::DeviceBuffer d_absNew(sizeof(Scalar)), d_absOld(sizeof(Scalar));
+  internal::DeviceBuffer d_pdot(sizeof(Scalar)), d_alpha(sizeof(Scalar));
+  internal::DeviceBuffer d_neg_alpha(sizeof(Scalar)), d_beta(sizeof(Scalar));
+  internal::DeviceBuffer d_rnorm(sizeof(RealScalar));
+
+  cudaMemcpy(d_b.ptr, b.data(), n * sizeof(Scalar), cudaMemcpyHostToDevice);
+
+  auto spmv = [&](Scalar* x_ptr, Scalar* y_ptr) {
+    cusparseDnVecDescr_t vx, vy;
+    cusparseCreateDnVec(&vx, n, x_ptr, CUDA_R_64F);
+    cusparseCreateDnVec(&vy, n, y_ptr, CUDA_R_64F);
+    cusparseSpMV(cusparse, CUSPARSE_OPERATION_NON_TRANSPOSE, &spmv_alpha, matA, vx, &spmv_beta, vy, CUDA_R_64F,
+                 CUSPARSE_SPMV_ALG_DEFAULT, d_workspace.ptr);
+    cusparseDestroyDnVec(vx);
+    cusparseDestroyDnVec(vy);
+  };
+
+  for (auto _ : state) {
+    cudaMemsetAsync(static_cast<Scalar*>(d_x.ptr), 0, n * sizeof(Scalar), stream);
+    cudaMemcpyAsync(d_r.ptr, d_b.ptr, n * sizeof(Scalar), cudaMemcpyDeviceToDevice, stream);
+    cudaMemcpyAsync(d_p.ptr, d_b.ptr, n * sizeof(Scalar), cudaMemcpyDeviceToDevice, stream);
+
+    cublasSetPointerMode(cublas, CUBLAS_POINTER_MODE_DEVICE);
+    cublasDdot(cublas, n, static_cast<Scalar*>(d_r.ptr), 1, static_cast<Scalar*>(d_p.ptr), 1,
+               static_cast<Scalar*>(d_absNew.ptr));
+
+    for (int i = 0; i < maxIters; ++i) {
+      spmv(static_cast<Scalar*>(d_p.ptr), static_cast<Scalar*>(d_tmp.ptr));
+
+      cublasDdot(cublas, n, static_cast<Scalar*>(d_p.ptr), 1, static_cast<Scalar*>(d_tmp.ptr), 1,
+                 static_cast<Scalar*>(d_pdot.ptr));
+
+      scalar_div_kernel<<<1, 1, 0, stream>>>(static_cast<Scalar*>(d_absNew.ptr), static_cast<Scalar*>(d_pdot.ptr),
+                                             static_cast<Scalar*>(d_alpha.ptr));
+      scalar_neg_kernel<<<1, 1, 0, stream>>>(static_cast<Scalar*>(d_alpha.ptr), static_cast<Scalar*>(d_neg_alpha.ptr));
+
+      cublasDaxpy(cublas, n, static_cast<Scalar*>(d_alpha.ptr), static_cast<Scalar*>(d_p.ptr), 1,
+                  static_cast<Scalar*>(d_x.ptr), 1);
+      cublasDaxpy(cublas, n, static_cast<Scalar*>(d_neg_alpha.ptr), static_cast<Scalar*>(d_tmp.ptr), 1,
+                  static_cast<Scalar*>(d_r.ptr), 1);
+
+      cublasDnrm2(cublas, n, static_cast<Scalar*>(d_r.ptr), 1, static_cast<RealScalar*>(d_rnorm.ptr));
+
+      RealScalar rnorm;
+      cudaMemcpyAsync(&rnorm, d_rnorm.ptr, sizeof(RealScalar), cudaMemcpyDeviceToHost, stream);
+      cudaStreamSynchronize(stream);
+      if (rnorm * rnorm < 1e-20) break;
+
+      cudaMemcpyAsync(d_absOld.ptr, d_absNew.ptr, sizeof(Scalar), cudaMemcpyDeviceToDevice, stream);
+      cublasDdot(cublas, n, static_cast<Scalar*>(d_r.ptr), 1, static_cast<Scalar*>(d_r.ptr), 1,
+                 static_cast<Scalar*>(d_absNew.ptr));
+
+      scalar_div_kernel<<<1, 1, 0, stream>>>(static_cast<Scalar*>(d_absNew.ptr), static_cast<Scalar*>(d_absOld.ptr),
+                                             static_cast<Scalar*>(d_beta.ptr));
+
+      cublasDscal(cublas, n, static_cast<Scalar*>(d_beta.ptr), static_cast<Scalar*>(d_p.ptr), 1);
+      cublasSetPointerMode(cublas, CUBLAS_POINTER_MODE_HOST);
+      Scalar one = 1.0;
+      cublasDaxpy(cublas, n, &one, static_cast<Scalar*>(d_r.ptr), 1, static_cast<Scalar*>(d_p.ptr), 1);
+      cublasSetPointerMode(cublas, CUBLAS_POINTER_MODE_DEVICE);
+    }
+    cudaStreamSynchronize(stream);
+  }
+
+  state.SetItemsProcessed(state.iterations() * maxIters);
+  cusparseDestroySpMat(matA);
+  cusparseDestroy(cusparse);
+  cublasDestroy(cublas);
+  cudaStreamDestroy(stream);
+}
+
+BENCHMARK(BM_CG_DevicePointerMode)->RangeMultiplier(4)->Range(1 << 10, 1 << 20);

unsupported/benchmarks/GPU/bench_cg_vs_cpu.cu

@@ -0,0 +1,218 @@
+// Benchmark: GPU CG vs CPU CG on realistic sparse systems.
+//
+// Tests 2D Laplacian (5-point stencil) and 3D Laplacian (7-point stencil)
+// in both float and double precision.
+//
+// Usage:
+//   cmake --build build-bench-gpu --target bench_gpu_cg_vs_cpu
+//   ./build-bench-gpu/bench_gpu_cg_vs_cpu
+// SPDX-FileCopyrightText: The Eigen Authors
+// SPDX-License-Identifier: MPL-2.0
+
+#include <benchmark/benchmark.h>
+
+#include <Eigen/Sparse>
+#include <Eigen/IterativeLinearSolvers>
+#include <unsupported/Eigen/GPU>
+
+using namespace Eigen;
+
+// ---- Sparse matrix generators -----------------------------------------------
+
+template <typename Scalar>
+SparseMatrix<Scalar, ColMajor, int> make_laplacian_2d(int grid_n) {
+  using SpMat = SparseMatrix<Scalar, ColMajor, int>;
+  const int n = grid_n * grid_n;
+  SpMat A(n, n);
+  A.reserve(VectorXi::Constant(n, 5));
+  for (int i = 0; i < grid_n; ++i) {
+    for (int j = 0; j < grid_n; ++j) {
+      int idx = i * grid_n + j;
+      A.insert(idx, idx) = Scalar(4);
+      if (i > 0) A.insert(idx, idx - grid_n) = Scalar(-1);
+      if (i < grid_n - 1) A.insert(idx, idx + grid_n) = Scalar(-1);
+      if (j > 0) A.insert(idx, idx - 1) = Scalar(-1);
+      if (j < grid_n - 1) A.insert(idx, idx + 1) = Scalar(-1);
+    }
+  }
+  A.makeCompressed();
+  return A;
+}
+
+template <typename Scalar>
+SparseMatrix<Scalar, ColMajor, int> make_laplacian_3d(int grid_n) {
+  using SpMat = SparseMatrix<Scalar, ColMajor, int>;
+  const int n = grid_n * grid_n * grid_n;
+  const int n2 = grid_n * grid_n;
+  SpMat A(n, n);
+  A.reserve(VectorXi::Constant(n, 7));
+  for (int i = 0; i < grid_n; ++i) {
+    for (int j = 0; j < grid_n; ++j) {
+      for (int k = 0; k < grid_n; ++k) {
+        int idx = i * n2 + j * grid_n + k;
+        A.insert(idx, idx) = Scalar(6);
+        if (i > 0) A.insert(idx, idx - n2) = Scalar(-1);
+        if (i < grid_n - 1) A.insert(idx, idx + n2) = Scalar(-1);
+        if (j > 0) A.insert(idx, idx - grid_n) = Scalar(-1);
+        if (j < grid_n - 1) A.insert(idx, idx + grid_n) = Scalar(-1);
+        if (k > 0) A.insert(idx, idx - 1) = Scalar(-1);
+        if (k < grid_n - 1) A.insert(idx, idx + 1) = Scalar(-1);
+      }
+    }
+  }
+  A.makeCompressed();
+  return A;
+}
+
+static void cuda_warmup() {
+  static bool done = false;
+  if (!done) {
+    void* p;
+    cudaMalloc(&p, 1);
+    cudaFree(p);
+    done = true;
+  }
+}
+
+// ---- CPU CG -----------------------------------------------------------------
+
+template <typename Scalar, typename MatGen>
+void run_cpu_cg(benchmark::State& state, MatGen make_matrix) {
+  using SpMat = SparseMatrix<Scalar, ColMajor, int>;
+  using Vec = Matrix<Scalar, Dynamic, 1>;
+  using RealScalar = typename NumTraits<Scalar>::Real;
+
+  const int grid_n = state.range(0);
+  SpMat A = make_matrix(grid_n);
+  Vec b = Vec::Random(A.rows());
+
+  ConjugateGradient<SpMat, Lower | Upper> cg;
+  cg.setMaxIterations(10000);
+  cg.setTolerance(RealScalar(1e-8));
+  cg.compute(A);
+
+  int last_iters = 0;
+  for (auto _ : state) {
+    Vec x = cg.solve(b);
+    benchmark::DoNotOptimize(x.data());
+    last_iters = cg.iterations();
+  }
+  state.counters["n"] = A.rows();
+  state.counters["nnz"] = A.nonZeros();
+  state.counters["iters"] = last_iters;
+  state.counters["error"] = cg.error();
+}
+
+// ---- GPU CG -----------------------------------------------------------------
+
+template <typename Scalar, typename MatGen>
+void run_gpu_cg(benchmark::State& state, MatGen make_matrix) {
+  using SpMat = SparseMatrix<Scalar, ColMajor, int>;
+  using Vec = Matrix<Scalar, Dynamic, 1>;
+  using RealScalar = typename NumTraits<Scalar>::Real;
+
+  cuda_warmup();
+  const int grid_n = state.range(0);
+  SpMat A = make_matrix(grid_n);
+  const Index n = A.rows();
+  Vec b = Vec::Random(n);
+
+  // Extract inverse diagonal.
+  Vec invdiag(n);
+  for (Index j = 0; j < A.outerSize(); ++j) {
+    typename SpMat::InnerIterator it(A, j);
+    while (it && it.index() != j) ++it;
+    if (it && it.index() == j && it.value() != Scalar(0))
+      invdiag(j) = Scalar(1) / it.value();
+    else
+      invdiag(j) = Scalar(1);
+  }
+
+  gpu::Context ctx;
+  gpu::Context::setThreadLocal(&ctx);
+  gpu::SparseContext<Scalar> spmv_ctx(ctx);
+  auto mat = spmv_ctx.deviceView(A);
+  auto d_invdiag = gpu::DeviceMatrix<Scalar>::fromHost(invdiag, ctx.stream());
+  auto d_b = gpu::DeviceMatrix<Scalar>::fromHost(b, ctx.stream());
+
+  int last_iters = 0;
+  RealScalar last_error = 0;
+
+  for (auto _ : state) {
+    gpu::DeviceMatrix<Scalar> d_x(n, 1);
+    d_x.setZero(ctx);
+    gpu::DeviceMatrix<Scalar> residual(n, 1);
+    residual.copyFrom(ctx, d_b);
+
+    RealScalar rhsNorm2 = d_b.squaredNorm(ctx);
+    RealScalar tol = RealScalar(1e-8);
+    RealScalar threshold = tol * tol * rhsNorm2;
+    RealScalar residualNorm2 = residual.squaredNorm(ctx);
+
+    gpu::DeviceMatrix<Scalar> p = d_invdiag.cwiseProduct(ctx, residual);
+    gpu::DeviceMatrix<Scalar> z(n, 1), tmp(n, 1);
+
+    auto absNew = residual.dot(ctx, p);
+    Index i = 0;
+    Index maxIters = 10000;
+    while (i < maxIters) {
+      tmp.noalias() = mat * p;
+      auto alpha = absNew / p.dot(ctx, tmp);
+      d_x += alpha * p;
+      residual -= alpha * tmp;
+
+      residualNorm2 = residual.squaredNorm(ctx);
+      if (residualNorm2 < threshold) break;
+
+      z.cwiseProduct(ctx, d_invdiag, residual);
+      auto absOld = std::move(absNew);
+      absNew = residual.dot(ctx, z);
+      auto beta = absNew / absOld;
+
+      p *= beta;
+      p += z;
+      i++;
+    }
+    benchmark::DoNotOptimize(d_x.data());
+    last_iters = i;
+    last_error = numext::sqrt(residualNorm2 / rhsNorm2);
+  }
+
+  gpu::Context::setThreadLocal(nullptr);
+  state.counters["n"] = n;
+  state.counters["nnz"] = A.nonZeros();
+  state.counters["iters"] = last_iters;
+  state.counters["error"] = last_error;
+}
+
+// ---- 2D Laplacian, double ---------------------------------------------------
+
+static void BM_CG_CPU_2D_double(benchmark::State& state) { run_cpu_cg<double>(state, make_laplacian_2d<double>); }
+static void BM_CG_GPU_2D_double(benchmark::State& state) { run_gpu_cg<double>(state, make_laplacian_2d<double>); }
+
+BENCHMARK(BM_CG_CPU_2D_double)->ArgsProduct({{32, 64, 128, 256, 512}});
+BENCHMARK(BM_CG_GPU_2D_double)->ArgsProduct({{32, 64, 128, 256, 512}});
+
+// ---- 2D Laplacian, float ----------------------------------------------------
+
+static void BM_CG_CPU_2D_float(benchmark::State& state) { run_cpu_cg<float>(state, make_laplacian_2d<float>); }
+static void BM_CG_GPU_2D_float(benchmark::State& state) { run_gpu_cg<float>(state, make_laplacian_2d<float>); }
+
+BENCHMARK(BM_CG_CPU_2D_float)->ArgsProduct({{32, 64, 128, 256, 512}});
+BENCHMARK(BM_CG_GPU_2D_float)->ArgsProduct({{32, 64, 128, 256, 512}});
+
+// ---- 3D Laplacian, double ---------------------------------------------------
+
+static void BM_CG_CPU_3D_double(benchmark::State& state) { run_cpu_cg<double>(state, make_laplacian_3d<double>); }
+static void BM_CG_GPU_3D_double(benchmark::State& state) { run_gpu_cg<double>(state, make_laplacian_3d<double>); }
+
+BENCHMARK(BM_CG_CPU_3D_double)->ArgsProduct({{16, 32, 48, 64}});
+BENCHMARK(BM_CG_GPU_3D_double)->ArgsProduct({{16, 32, 48, 64}});
+
+// ---- 3D Laplacian, float ----------------------------------------------------
+
+static void BM_CG_CPU_3D_float(benchmark::State& state) { run_cpu_cg<float>(state, make_laplacian_3d<float>); }
+static void BM_CG_GPU_3D_float(benchmark::State& state) { run_gpu_cg<float>(state, make_laplacian_3d<float>); }
+
+BENCHMARK(BM_CG_CPU_3D_float)->ArgsProduct({{16, 32, 48, 64}});
+BENCHMARK(BM_CG_GPU_3D_float)->ArgsProduct({{16, 32, 48, 64}});

unsupported/test/GPU/CMakeLists.txt

@@ -3,7 +3,7 @@
 # SPDX-FileCopyrightText: The Eigen Authors
 # SPDX-License-Identifier: MPL-2.0
 
-find_package(CUDAToolkit)
+find_package(CUDAToolkit QUIET)
 
 # ei_add_gpu_test(<name> [EXTRA_LIBS <lib>...] [EXTRA_INCLUDES <dir>...] [EXTRA_DEFINES <def>...])
 #
@@ -57,17 +57,18 @@ function(ei_add_gpu_test test_name)
   endforeach()
 endfunction()
 
-# DeviceMatrix core: CUDA runtime only.
-ei_add_gpu_test(device_matrix)
+# DeviceMatrix core: CUDA runtime + cuBLAS + cuSOLVER (for BLAS-1 ops via GpuContext).
+ei_add_gpu_test(device_matrix
+  EXTRA_LIBS CUDA::cublas CUDA::cusolver CUDA::npps CUDA::nppc)
 
-# cuBLAS integration (BLAS-3, triangular solves, rank-k updates, ...).
-option(EIGEN_TEST_CUBLAS "Test cuBLAS integration" ON)
+# cuBLAS integration. OFF by default; later MRs in the stack opt in.
+option(EIGEN_TEST_CUBLAS "Test cuBLAS integration" OFF)
 if(EIGEN_TEST_CUBLAS AND TARGET CUDA::cublas)
   ei_add_gpu_test(cublas EXTRA_LIBS CUDA::cublas CUDA::cusolver)
 endif()
 
-# cuSOLVER dense factorizations (LLT, LU).
-option(EIGEN_TEST_CUSOLVER "Test cuSOLVER integration" ON)
+# cuSOLVER dense factorizations (LLT, LU, QR, SVD, self-adjoint eigensolver).
+option(EIGEN_TEST_CUSOLVER "Test cuSOLVER integration" OFF)
 if(EIGEN_TEST_CUSOLVER AND TARGET CUDA::cusolver)
   foreach(_cusolver_test IN ITEMS cusolver_llt cusolver_lu cusolver_qr cusolver_svd cusolver_eigen)
     ei_add_gpu_test(${_cusolver_test} EXTRA_LIBS CUDA::cusolver CUDA::cublas)
@@ -79,10 +80,18 @@ if(TARGET CUDA::cufft)
   ei_add_gpu_test(cufft EXTRA_LIBS CUDA::cufft CUDA::cublas)
 endif()
 
-# cuSPARSE SpMV (part of the CUDA toolkit).
+# cuSPARSE SpMV (part of the CUDA toolkit) and the end-to-end CG test that
+# exercises Eigen::ConjugateGradient with DeviceMatrix + device-resident SpMV.
 if(TARGET CUDA::cusparse)
   ei_add_gpu_test(cusparse_spmv
     EXTRA_LIBS CUDA::cusparse CUDA::cublas CUDA::cusolver)
+  ei_add_gpu_test(cg
+    EXTRA_LIBS CUDA::cusparse CUDA::cublas CUDA::cusolver CUDA::npps CUDA::nppc)
+endif()
+
+option(EIGEN_TEST_CUSPARSE "Test cuSPARSE integration" OFF)
+if(EIGEN_TEST_CUSPARSE AND TARGET CUDA::cusparse)
+  ei_add_gpu_test(cusparse EXTRA_LIBS CUDA::cusparse)
 endif()
 
 # cuDSS sparse direct solvers -- distributed separately from the CUDA Toolkit.

unsupported/test/GPU/cg.cpp

@@ -0,0 +1,231 @@
+// 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
+
+// End-to-end test: CG algorithm running on GPU via gpu::DeviceMatrix.
+//
+// Uses DeviceSparseView for SpMV, gpu::DeviceMatrix for vectors, DeviceScalar
+// for deferred reductions. Verifies correctness against CPU ConjugateGradient.
+
+#define EIGEN_USE_GPU
+#include "main.h"
+#include <Eigen/Sparse>
+#include <Eigen/IterativeLinearSolvers>
+#include <unsupported/Eigen/GPU>
+#include "gpu_test_helpers.h"
+
+using namespace Eigen;
+
+// ---- Helper: build a sparse SPD matrix --------------------------------------
+
+template <typename Scalar>
+SparseMatrix<Scalar, ColMajor, int> make_spd(Index n, double density = 0.1) {
+  using SpMat = SparseMatrix<Scalar, ColMajor, int>;
+  using RealScalar = typename NumTraits<Scalar>::Real;
+
+  SpMat R(n, n);
+  R.reserve(VectorXi::Constant(n, static_cast<int>(n * density) + 1));
+  for (Index j = 0; j < n; ++j) {
+    for (Index i = 0; i < n; ++i) {
+      if (i == j || (std::rand() / double(RAND_MAX)) < density) {
+        R.insert(i, j) = Scalar(std::rand() / double(RAND_MAX) - 0.5);
+      }
+    }
+  }
+  R.makeCompressed();
+  SpMat A = R.adjoint() * R;
+  for (Index i = 0; i < n; ++i) A.coeffRef(i, i) += Scalar(RealScalar(n));
+  A.makeCompressed();
+  return A;
+}
+
+// ---- GPU CG without preconditioner ------------------------------------------
+
+template <typename Scalar>
+void test_gpu_cg(Index n) {
+  using SpMat = SparseMatrix<Scalar, ColMajor, int>;
+  using Vec = Matrix<Scalar, Dynamic, 1>;
+  using RealScalar = typename NumTraits<Scalar>::Real;
+
+  SpMat A = make_spd<Scalar>(n);
+  Vec b = Vec::Random(n);
+
+  // CPU reference (identity preconditioner to match GPU).
+  ConjugateGradient<SpMat, Lower | Upper, IdentityPreconditioner> cpu_cg;
+  cpu_cg.setMaxIterations(1000);
+  cpu_cg.setTolerance(RealScalar(1e-8));
+  cpu_cg.compute(A);
+  Vec x_cpu = cpu_cg.solve(b);
+  VERIFY_IS_EQUAL(cpu_cg.info(), Success);
+
+  // GPU CG: mirrors Eigen's conjugate_gradient() using gpu::DeviceMatrix ops.
+  gpu::Context ctx;
+  gpu::Context::setThreadLocal(&ctx);
+  gpu::SparseContext<Scalar> spmv_ctx(ctx);
+  auto mat = spmv_ctx.deviceView(A);
+
+  auto d_b = gpu::DeviceMatrix<Scalar>::fromHost(b, ctx.stream());
+  gpu::DeviceMatrix<Scalar> d_x(n, 1);
+  d_x.setZero(ctx);
+
+  // r = b (since x=0)
+  gpu::DeviceMatrix<Scalar> residual(n, 1);
+  residual.copyFrom(ctx, d_b);
+
+  RealScalar rhsNorm2 = d_b.squaredNorm(ctx);
+  RealScalar tol = RealScalar(1e-8);
+  RealScalar threshold = tol * tol * rhsNorm2;
+  RealScalar residualNorm2 = residual.squaredNorm(ctx);
+
+  // p = r (no preconditioner)
+  gpu::DeviceMatrix<Scalar> p(n, 1);
+  p.copyFrom(ctx, residual);
+  gpu::DeviceMatrix<Scalar> z(n, 1), tmp(n, 1);
+
+  auto absNew = residual.dot(ctx, p);
+  Index maxIters = 1000;
+  Index i = 0;
+  while (i < maxIters) {
+    tmp.noalias() = mat * p;
+
+    auto alpha = absNew / p.dot(ctx, tmp);
+    d_x += alpha * p;
+    residual -= alpha * tmp;
+
+    residualNorm2 = residual.squaredNorm(ctx);
+    if (residualNorm2 < threshold) break;
+
+    // z = r (no preconditioner)
+    z.copyFrom(ctx, residual);
+
+    auto absOld = std::move(absNew);
+    absNew = residual.dot(ctx, z);
+    auto beta = absNew / absOld;
+
+    p *= beta;
+    p += z;
+    i++;
+  }
+
+  gpu::Context::setThreadLocal(nullptr);
+
+  Vec x_gpu = d_x.toHost(ctx.stream());
+
+  // Verify residual.
+  Vec r = A * x_gpu - b;
+  RealScalar relres = r.norm() / b.norm();
+  VERIFY(relres < RealScalar(1e-6));
+
+  // Compare with CPU.
+  RealScalar sol_tol = RealScalar(100) * RealScalar(n) * NumTraits<Scalar>::epsilon();
+  VERIFY((x_gpu - x_cpu).norm() / (x_cpu.norm() + RealScalar(1)) < sol_tol);
+}
+
+// ---- GPU CG with Jacobi preconditioner --------------------------------------
+
+template <typename Scalar>
+void test_gpu_cg_jacobi(Index n) {
+  using SpMat = SparseMatrix<Scalar, ColMajor, int>;
+  using Vec = Matrix<Scalar, Dynamic, 1>;
+  using RealScalar = typename NumTraits<Scalar>::Real;
+
+  SpMat A = make_spd<Scalar>(n);
+  Vec b = Vec::Random(n);
+
+  // CPU reference.
+  ConjugateGradient<SpMat, Lower | Upper> cpu_cg;
+  cpu_cg.setMaxIterations(1000);
+  cpu_cg.setTolerance(RealScalar(1e-8));
+  cpu_cg.compute(A);
+  Vec x_cpu = cpu_cg.solve(b);
+
+  // Extract inverse diagonal.
+  Vec invdiag(n);
+  for (Index j = 0; j < A.outerSize(); ++j) {
+    typename SpMat::InnerIterator it(A, j);
+    while (it && it.index() != j) ++it;
+    if (it && it.index() == j && it.value() != Scalar(0))
+      invdiag(j) = Scalar(1) / it.value();
+    else
+      invdiag(j) = Scalar(1);
+  }
+
+  // GPU CG with Jacobi preconditioner.
+  gpu::Context ctx;
+  gpu::Context::setThreadLocal(&ctx);
+  gpu::SparseContext<Scalar> spmv_ctx(ctx);
+  auto mat = spmv_ctx.deviceView(A);
+  auto d_invdiag = gpu::DeviceMatrix<Scalar>::fromHost(invdiag, ctx.stream());
+
+  auto d_b = gpu::DeviceMatrix<Scalar>::fromHost(b, ctx.stream());
+  gpu::DeviceMatrix<Scalar> d_x(n, 1);
+  d_x.setZero(ctx);
+
+  gpu::DeviceMatrix<Scalar> residual(n, 1);
+  residual.copyFrom(ctx, d_b);
+
+  RealScalar rhsNorm2 = d_b.squaredNorm(ctx);
+  RealScalar tol = RealScalar(1e-8);
+  RealScalar threshold = tol * tol * rhsNorm2;
+  RealScalar residualNorm2 = residual.squaredNorm(ctx);
+
+  // p = precond.solve(r) = invdiag .* r
+  gpu::DeviceMatrix<Scalar> p = d_invdiag.cwiseProduct(ctx, residual);
+  gpu::DeviceMatrix<Scalar> z(n, 1), tmp(n, 1);
+
+  auto absNew = residual.dot(ctx, p);
+  Index maxIters = 1000;
+  Index i = 0;
+  while (i < maxIters) {
+    tmp.noalias() = mat * p;
+
+    auto alpha = absNew / p.dot(ctx, tmp);
+    d_x += alpha * p;
+    residual -= alpha * tmp;
+
+    residualNorm2 = residual.squaredNorm(ctx);
+    if (residualNorm2 < threshold) break;
+
+    // z = precond.solve(r) = invdiag .* r
+    z.cwiseProduct(ctx, d_invdiag, residual);
+
+    auto absOld = std::move(absNew);
+    absNew = residual.dot(ctx, z);
+    auto beta = absNew / absOld;
+
+    p *= beta;
+    p += z;
+    i++;
+  }
+
+  gpu::Context::setThreadLocal(nullptr);
+
+  Vec x_gpu = d_x.toHost(ctx.stream());
+
+  Vec r = A * x_gpu - b;
+  RealScalar relres = r.norm() / b.norm();
+  VERIFY(relres < RealScalar(1e-6));
+
+  RealScalar sol_tol = RealScalar(100) * RealScalar(n) * NumTraits<Scalar>::epsilon();
+  VERIFY((x_gpu - x_cpu).norm() / (x_cpu.norm() + RealScalar(1)) < sol_tol);
+}
+
+EIGEN_DECLARE_TEST(gpu_cg) {
+  gpu_test::require_cusparse_context();
+
+  // Split by scalar so each part compiles in parallel.
+  CALL_SUBTEST_1(test_gpu_cg<double>(64));
+  CALL_SUBTEST_1(test_gpu_cg<double>(256));
+  CALL_SUBTEST_1(test_gpu_cg_jacobi<double>(64));
+  CALL_SUBTEST_1(test_gpu_cg_jacobi<double>(256));
+  CALL_SUBTEST_2(test_gpu_cg<float>(64));
+  CALL_SUBTEST_2(test_gpu_cg<float>(256));
+  CALL_SUBTEST_2(test_gpu_cg_jacobi<float>(64));
+  CALL_SUBTEST_2(test_gpu_cg_jacobi<float>(256));
+}

unsupported/test/GPU/cusparse_spmv.cpp

@@ -225,6 +225,105 @@ void test_empty() {
   VERIFY_IS_EQUAL(y.size(), 0);
 }
 
+// ---- gpu::DeviceMatrix SpMV (no host roundtrip) ----------------------------------
+
+template <typename Scalar>
+void test_spmv_device(Index n) {
+  using SpMat = SparseMatrix<Scalar, ColMajor, int>;
+  using Vec = Matrix<Scalar, Dynamic, 1>;
+  using RealScalar = typename NumTraits<Scalar>::Real;
+
+  SpMat A = make_sparse<Scalar>(n, n);
+  Vec x = Vec::Random(n);
+
+  // Use shared gpu::Context for same-stream execution.
+  gpu::Context gpu_ctx;
+  gpu::SparseContext<Scalar> ctx(gpu_ctx);
+
+  auto d_x = gpu::DeviceMatrix<Scalar>::fromHost(x, gpu_ctx.stream());
+  gpu::DeviceMatrix<Scalar> d_y;
+
+  ctx.multiply(A, d_x, d_y);
+
+  Vec y_gpu = d_y.toHost(gpu_ctx.stream());
+  Vec y_cpu = A * x;
+
+  RealScalar tol = RealScalar(10) * RealScalar(n) * NumTraits<Scalar>::epsilon();
+  VERIFY((y_gpu - y_cpu).norm() / (y_cpu.norm() + RealScalar(1)) < tol);
+}
+
+// ---- Expression syntax: d_y = d_A * d_x ------------------------------------
+
+template <typename Scalar>
+void test_spmv_expr(Index n) {
+  using SpMat = SparseMatrix<Scalar, ColMajor, int>;
+  using Vec = Matrix<Scalar, Dynamic, 1>;
+  using RealScalar = typename NumTraits<Scalar>::Real;
+
+  SpMat A = make_sparse<Scalar>(n, n);
+  Vec x = Vec::Random(n);
+
+  gpu::Context gpu_ctx;
+  gpu::SparseContext<Scalar> ctx(gpu_ctx);
+
+  // Upload sparse matrix and create device view.
+  auto d_A = ctx.deviceView(A);
+
+  // Upload x.
+  auto d_x = gpu::DeviceMatrix<Scalar>::fromHost(x, gpu_ctx.stream());
+
+  // Expression syntax: d_y = d_A * d_x
+  gpu::DeviceMatrix<Scalar> d_y;
+  d_y = d_A * d_x;
+
+  // Also test with noalias():
+  gpu::DeviceMatrix<Scalar> d_tmp;
+  d_tmp.noalias() = d_A * d_x;
+
+  Vec y_gpu = d_y.toHost(gpu_ctx.stream());
+  Vec tmp_gpu = d_tmp.toHost(gpu_ctx.stream());
+  Vec y_cpu = A * x;
+
+  RealScalar tol = RealScalar(10) * RealScalar(n) * NumTraits<Scalar>::epsilon();
+  VERIFY((y_gpu - y_cpu).norm() / (y_cpu.norm() + RealScalar(1)) < tol);
+  VERIFY((tmp_gpu - y_cpu).norm() / (y_cpu.norm() + RealScalar(1)) < tol);
+}
+
+// ---- deviceView overwrite: second view replaces first -----------------------
+
+template <typename Scalar>
+void test_deviceview_overwrite(Index n) {
+  using SpMat = SparseMatrix<Scalar, ColMajor, int>;
+  using Vec = Matrix<Scalar, Dynamic, 1>;
+  using RealScalar = typename NumTraits<Scalar>::Real;
+
+  SpMat A1 = make_sparse<Scalar>(n, n);
+  SpMat A2 = make_sparse<Scalar>(n, n);  // different random matrix
+
+  Vec x = Vec::Random(n);
+
+  gpu::Context gpu_ctx;
+  gpu::SparseContext<Scalar> ctx(gpu_ctx);
+
+  // First view: A1.
+  auto d_A1 = ctx.deviceView(A1);
+  auto d_x = gpu::DeviceMatrix<Scalar>::fromHost(x, gpu_ctx.stream());
+  gpu::DeviceMatrix<Scalar> d_y1;
+  d_y1 = d_A1 * d_x;
+  Vec y1_gpu = d_y1.toHost(gpu_ctx.stream());
+  Vec y1_cpu = A1 * x;
+  RealScalar tol = RealScalar(10) * RealScalar(n) * NumTraits<Scalar>::epsilon();
+  VERIFY((y1_gpu - y1_cpu).norm() / (y1_cpu.norm() + RealScalar(1)) < tol);
+
+  // Second view overwrites first: now uses A2.
+  auto d_A2 = ctx.deviceView(A2);
+  gpu::DeviceMatrix<Scalar> d_y2;
+  d_y2 = d_A2 * d_x;
+  Vec y2_gpu = d_y2.toHost(gpu_ctx.stream());
+  Vec y2_cpu = A2 * x;
+  VERIFY((y2_gpu - y2_cpu).norm() / (y2_cpu.norm() + RealScalar(1)) < tol);
+}
+
 // ---- Per-scalar driver ------------------------------------------------------
 
 template <typename Scalar>
@@ -234,12 +333,22 @@ void test_scalar() {
   CALL_SUBTEST(test_spmv<Scalar>(64, 128));  // wide
   CALL_SUBTEST(test_spmv_alpha_beta<Scalar>(64));
   CALL_SUBTEST(test_spmv_transpose<Scalar>(128, 64));
+  // cuSPARSE < 12 cannot represent A^H * x for complex scalars with the
+  // CSR-of-A^T trick used by SparseContext; SparseContext asserts in that
+  // case. Real-scalar ConjTrans is demoted to Trans and works fine.
+#if !defined(CUSPARSE_VERSION) || CUSPARSE_VERSION >= 12000
   CALL_SUBTEST(test_spmv_adjoint<Scalar>(128, 64));
+#else
+  if (!NumTraits<Scalar>::IsComplex) CALL_SUBTEST(test_spmv_adjoint<Scalar>(128, 64));
+#endif
   CALL_SUBTEST(test_spmm<Scalar>(64, 64, 4));
   CALL_SUBTEST(test_spmm_transpose<Scalar>(128, 64, 4));
   CALL_SUBTEST(test_identity<Scalar>(64));
   CALL_SUBTEST(test_reuse<Scalar>(64));
   CALL_SUBTEST(test_empty<Scalar>());
+  CALL_SUBTEST(test_spmv_device<Scalar>(64));
+  CALL_SUBTEST(test_spmv_expr<Scalar>(64));
+  CALL_SUBTEST(test_deviceview_overwrite<Scalar>(64));
 }
 
 EIGEN_DECLARE_TEST(gpu_cusparse_spmv) {

unsupported/test/GPU/device_matrix.cpp

@@ -8,11 +8,12 @@
 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 // SPDX-License-Identifier: MPL-2.0
 
-// Tests for DeviceMatrix and HostTransfer: typed RAII GPU memory wrapper.
+// Tests for gpu::DeviceMatrix and HostTransfer: typed RAII GPU memory wrapper.
 // No cuSOLVER dependency — only CUDA runtime.
 
 #define EIGEN_USE_GPU
 #include "main.h"
+#include <Eigen/Sparse>
 #include <unsupported/Eigen/GPU>
 
 using namespace Eigen;
@@ -244,6 +245,217 @@ void test_scalar() {
   CALL_SUBTEST(test_move_assign<Scalar>(64, 64));
 }
 
+// ---- BLAS-1: dot product ----------------------------------------------------
+
+template <typename Scalar>
+void test_blas1(Index n) {
+  using Vec = Matrix<Scalar, Dynamic, 1>;
+  using RealScalar = typename NumTraits<Scalar>::Real;
+
+  // All BLAS-1 ops share one gpu::Context — same stream, zero event overhead.
+  gpu::Context ctx;
+
+  RealScalar tol = RealScalar(10) * RealScalar(n) * NumTraits<Scalar>::epsilon();
+
+  // dot
+  {
+    Vec a = Vec::Random(n);
+    Vec b = Vec::Random(n);
+    auto d_a = gpu::DeviceMatrix<Scalar>::fromHost(a, ctx.stream());
+    auto d_b = gpu::DeviceMatrix<Scalar>::fromHost(b, ctx.stream());
+    Scalar gpu_dot = d_a.dot(ctx, d_b);
+    Scalar cpu_dot = a.dot(b);
+    VERIFY(numext::abs(gpu_dot - cpu_dot) < tol * numext::abs(cpu_dot) + tol);
+  }
+
+  // norm / squaredNorm
+  {
+    Vec a = Vec::Random(n);
+    auto d_a = gpu::DeviceMatrix<Scalar>::fromHost(a, ctx.stream());
+    RealScalar gpu_norm = d_a.norm(ctx);
+    RealScalar cpu_norm = a.norm();
+    VERIFY(numext::abs(gpu_norm - cpu_norm) < tol * cpu_norm + tol);
+    RealScalar gpu_sqnorm = d_a.squaredNorm(ctx);
+    RealScalar cpu_sqnorm = a.squaredNorm();
+    VERIFY(numext::abs(gpu_sqnorm - cpu_sqnorm) < tol * cpu_sqnorm + tol);
+  }
+
+  // addScaled (axpy)
+  {
+    Vec x = Vec::Random(n);
+    Vec y = Vec::Random(n);
+    Scalar alpha(2.5);
+    Vec y_ref = y + alpha * x;
+    auto d_y = gpu::DeviceMatrix<Scalar>::fromHost(y, ctx.stream());
+    auto d_x = gpu::DeviceMatrix<Scalar>::fromHost(x, ctx.stream());
+    d_y.addScaled(ctx, alpha, d_x);
+    Vec y_gpu = d_y.toHost(ctx.stream());
+    VERIFY((y_gpu - y_ref).norm() < tol * y_ref.norm() + tol);
+  }
+
+  // scale (scal)
+  {
+    Vec x = Vec::Random(n);
+    Scalar alpha(3.0);
+    Vec x_ref = alpha * x;
+    auto d_x = gpu::DeviceMatrix<Scalar>::fromHost(x, ctx.stream());
+    d_x.scale(ctx, alpha);
+    Vec x_gpu = d_x.toHost(ctx.stream());
+    VERIFY((x_gpu - x_ref).norm() < tol * x_ref.norm() + tol);
+  }
+
+  // copyFrom
+  {
+    Vec x = Vec::Random(n);
+    auto d_x = gpu::DeviceMatrix<Scalar>::fromHost(x, ctx.stream());
+    gpu::DeviceMatrix<Scalar> d_y;
+    d_y.copyFrom(ctx, d_x);
+    Vec y = d_y.toHost(ctx.stream());
+    VERIFY_IS_APPROX(y, x);
+  }
+
+  // setZero
+  {
+    Vec x = Vec::Random(n);
+    auto d_x = gpu::DeviceMatrix<Scalar>::fromHost(x, ctx.stream());
+    d_x.setZero(ctx);
+    Vec result = d_x.toHost(ctx.stream());
+    VERIFY_IS_EQUAL(result, Vec::Zero(n));
+  }
+}
+
+// ---- BLAS-1 operator overloads (CG-style) -----------------------------------
+
+template <typename Scalar>
+void test_cg_operators(Index n) {
+  using Vec = Matrix<Scalar, Dynamic, 1>;
+  using RealScalar = typename NumTraits<Scalar>::Real;
+
+  RealScalar tol = RealScalar(10) * RealScalar(n) * NumTraits<Scalar>::epsilon();
+
+  Vec x = Vec::Random(n);
+  Vec p = Vec::Random(n);
+  Vec tmp = Vec::Random(n);
+  Vec z = Vec::Random(n);
+  Scalar alpha(2.5);
+  Scalar beta(0.7);
+
+  // Test: x += alpha * p
+  {
+    Vec x_ref = x + alpha * p;
+    auto d_x = gpu::DeviceMatrix<Scalar>::fromHost(x);
+    auto d_p = gpu::DeviceMatrix<Scalar>::fromHost(p);
+    d_x += alpha * d_p;
+    Vec x_gpu = d_x.toHost();
+    VERIFY((x_gpu - x_ref).norm() < tol * x_ref.norm() + tol);
+  }
+
+  // Test: r -= alpha * tmp
+  {
+    Vec r = Vec::Random(n);
+    Vec r_ref = r - alpha * tmp;
+    auto d_r = gpu::DeviceMatrix<Scalar>::fromHost(r);
+    auto d_tmp = gpu::DeviceMatrix<Scalar>::fromHost(tmp);
+    d_r -= alpha * d_tmp;
+    Vec r_gpu = d_r.toHost();
+    VERIFY((r_gpu - r_ref).norm() < tol * r_ref.norm() + tol);
+  }
+
+  // Test: p = z + beta * p  (cuBLAS geam)
+  {
+    Vec p_copy = p;
+    Vec p_ref = z + beta * p_copy;
+    auto d_p = gpu::DeviceMatrix<Scalar>::fromHost(p_copy);
+    auto d_z = gpu::DeviceMatrix<Scalar>::fromHost(z);
+    d_p = d_z + beta * d_p;
+    Vec p_gpu = d_p.toHost();
+    VERIFY((p_gpu - p_ref).norm() < tol * p_ref.norm() + tol);
+  }
+
+  // Test: operator+= and operator-= with gpu::DeviceMatrix (no scalar)
+  {
+    Vec a = Vec::Random(n);
+    Vec b = Vec::Random(n);
+    Vec a_ref = a + b;
+    auto d_a = gpu::DeviceMatrix<Scalar>::fromHost(a);
+    auto d_b = gpu::DeviceMatrix<Scalar>::fromHost(b);
+    d_a += d_b;
+    VERIFY((d_a.toHost() - a_ref).norm() < tol * a_ref.norm() + tol);
+  }
+}
+
+// ---- gpu::DeviceScalar: deferred sync -------------------------------------------
+
+template <typename Scalar>
+void test_device_scalar() {
+  using Vec = Matrix<Scalar, Dynamic, 1>;
+  using RealScalar = typename NumTraits<Scalar>::Real;
+
+  const Index n = 256;
+  Vec a = Vec::Random(n);
+  Vec b = Vec::Random(n);
+
+  gpu::Context ctx;
+  auto d_a = gpu::DeviceMatrix<Scalar>::fromHost(a, ctx.stream());
+  auto d_b = gpu::DeviceMatrix<Scalar>::fromHost(b, ctx.stream());
+
+  // dot() returns gpu::DeviceScalar — implicit conversion to Scalar syncs.
+  Scalar gpu_dot = d_a.dot(ctx, d_b);
+  Scalar cpu_dot = a.dot(b);
+  RealScalar tol = RealScalar(10) * RealScalar(n) * NumTraits<Scalar>::epsilon();
+  VERIFY(numext::abs(gpu_dot - cpu_dot) < tol * numext::abs(cpu_dot) + tol);
+
+  // squaredNorm() returns host RealScalar directly (syncs internally).
+  RealScalar gpu_sqnorm = d_a.squaredNorm(ctx);
+  RealScalar cpu_sqnorm = a.squaredNorm();
+  VERIFY(numext::abs(gpu_sqnorm - cpu_sqnorm) < tol * cpu_sqnorm + tol);
+
+  // norm() returns gpu::DeviceScalar<RealScalar> — implicit conversion syncs.
+  RealScalar gpu_norm = d_a.norm(ctx);
+  RealScalar cpu_norm = a.norm();
+  VERIFY(numext::abs(gpu_norm - cpu_norm) < tol * cpu_norm + tol);
+
+  // Convenience overloads (thread-local context).
+  gpu::Context::setThreadLocal(&ctx);
+  Scalar gpu_dot2 = d_a.dot(d_b);
+  VERIFY(numext::abs(gpu_dot2 - cpu_dot) < tol * numext::abs(cpu_dot) + tol);
+  gpu::Context::setThreadLocal(nullptr);
+
+  // Empty vectors: dot and norm must return zero.
+  {
+    gpu::DeviceMatrix<Scalar> d_empty(0, 1);
+    gpu::DeviceMatrix<Scalar> d_empty2(0, 1);
+    Scalar empty_dot = d_empty.dot(ctx, d_empty2);
+    VERIFY_IS_EQUAL(empty_dot, Scalar(0));
+    RealScalar empty_sqnorm = d_empty.squaredNorm(ctx);
+    VERIFY_IS_EQUAL(empty_sqnorm, RealScalar(0));
+    RealScalar empty_norm = d_empty.norm(ctx);
+    VERIFY_IS_EQUAL(empty_norm, RealScalar(0));
+  }
+}
+
+// ---- cwiseProduct -----------------------------------------------------------
+
+template <typename Scalar>
+void test_cwiseProduct() {
+  using Vec = Matrix<Scalar, Dynamic, 1>;
+  using RealScalar = typename NumTraits<Scalar>::Real;
+
+  const Index n = 256;
+  Vec a = Vec::Random(n);
+  Vec b = Vec::Random(n);
+  Vec ref = a.array() * b.array();
+
+  gpu::Context ctx;
+  auto d_a = gpu::DeviceMatrix<Scalar>::fromHost(a, ctx.stream());
+  auto d_b = gpu::DeviceMatrix<Scalar>::fromHost(b, ctx.stream());
+  auto d_c = d_a.cwiseProduct(ctx, d_b);
+  Vec result = d_c.toHost(ctx.stream());
+
+  RealScalar tol = RealScalar(10) * RealScalar(n) * NumTraits<Scalar>::epsilon();
+  VERIFY((result - ref).norm() < tol * ref.norm() + tol);
+}
+
 EIGEN_DECLARE_TEST(gpu_device_matrix) {
   CALL_SUBTEST(test_default_construct());
   CALL_SUBTEST(test_empty());
@@ -256,4 +468,18 @@ EIGEN_DECLARE_TEST(gpu_device_matrix) {
   CALL_SUBTEST(test_scalar<double>());
   CALL_SUBTEST(test_scalar<std::complex<float>>());
   CALL_SUBTEST(test_scalar<std::complex<double>>());
+  CALL_SUBTEST(test_blas1<float>(256));
+  CALL_SUBTEST(test_blas1<double>(256));
+  CALL_SUBTEST(test_blas1<std::complex<float>>(256));
+  CALL_SUBTEST(test_blas1<std::complex<double>>(256));
+  CALL_SUBTEST(test_cg_operators<float>(256));
+  CALL_SUBTEST(test_cg_operators<double>(256));
+  CALL_SUBTEST(test_cg_operators<std::complex<float>>(256));
+  CALL_SUBTEST(test_cg_operators<std::complex<double>>(256));
+  CALL_SUBTEST(test_device_scalar<float>());
+  CALL_SUBTEST(test_device_scalar<double>());
+  CALL_SUBTEST(test_device_scalar<std::complex<float>>());
+  CALL_SUBTEST(test_device_scalar<std::complex<double>>());
+  CALL_SUBTEST(test_cwiseProduct<float>());
+  CALL_SUBTEST(test_cwiseProduct<double>());
 }