oxicuda_sparse/ops/

spmv.rs

//! Sparse matrix-vector multiplication (SpMV).
//!
//! Computes `y = alpha * A * x + beta * y` where `A` is a sparse CSR matrix
//! and `x`, `y` are dense vectors stored as raw device pointers.
//!
//! Three kernel strategies are available:
//! - **Scalar**: one thread per row (best for very sparse rows, < 4 nnz/row)
//! - **Vector**: one warp per row with shuffle reduction (best for moderate sparsity)
//! - **Adaptive**: auto-selects based on average nnz per row

use std::sync::Arc;

use oxicuda_blas::GpuFloat;
use oxicuda_driver::Module;
use oxicuda_driver::ffi::CUdeviceptr;
use oxicuda_launch::{Kernel, LaunchParams, grid_size_for};
use oxicuda_ptx::prelude::*;

use crate::error::{SparseError, SparseResult};
use crate::format::CsrMatrix;
use crate::handle::SparseHandle;
use crate::ptx_helpers::{
    add_float, emit_warp_reduce_sum, fma_float, load_float_imm, load_global_float, mul_float,
    ptx_suffix, reinterpret_bits_to_float, store_global_float,
};

/// Algorithm selection for SpMV.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
pub enum SpMVAlgo {
    /// One thread per row. Best for very sparse matrices (< 4 nnz/row).
    Scalar,
    /// One warp (32 threads) per row with warp shuffle reduction.
    /// Best for moderate sparsity (4-64 nnz/row).
    Vector,
    /// Automatically selects Scalar or Vector based on the matrix structure.
    Adaptive,
}

/// Default block size for scalar SpMV.
const SPMV_SCALAR_BLOCK: u32 = 256;

/// Default block size for vector SpMV (must be a multiple of warp size 32).
const SPMV_VECTOR_BLOCK: u32 = 256;

/// Threshold for auto-selecting Vector over Scalar (nnz per row).
const VECTOR_THRESHOLD: f64 = 4.0;

/// Resolve the [`SpMVAlgo::Adaptive`] selection to a concrete kernel given the
/// matrix's average non-zeros per row.
///
/// This pure function contains the kernel-selection heuristic so it can be
/// tested independently of GPU device memory.
///
/// * `avg_nnz_per_row < VECTOR_THRESHOLD` → [`SpMVAlgo::Scalar`]
/// * `avg_nnz_per_row >= VECTOR_THRESHOLD` → [`SpMVAlgo::Vector`]
#[inline]
pub(crate) fn resolve_adaptive(avg_nnz_per_row: f64) -> SpMVAlgo {
    if avg_nnz_per_row >= VECTOR_THRESHOLD {
        SpMVAlgo::Vector
    } else {
        SpMVAlgo::Scalar
    }
}

/// Sparse matrix-vector multiplication: `y = alpha * A * x + beta * y`.
///
/// # Arguments
///
/// * `handle` -- Sparse handle providing stream and device context.
/// * `algo` -- Algorithm selection strategy.
/// * `alpha` -- Scalar multiplier for `A * x`.
/// * `a` -- Sparse CSR matrix `A`.
/// * `x_ptr` -- Device pointer to dense vector `x` of length `A.cols()`.
/// * `beta` -- Scalar multiplier for existing `y`.
/// * `y_ptr` -- Device pointer to dense vector `y` of length `A.rows()`.
///
/// # Errors
///
/// Returns [`SparseError::PtxGeneration`] if kernel generation fails.
/// Returns [`SparseError::Cuda`] on kernel launch failure.
#[allow(clippy::too_many_arguments)]
pub fn spmv<T: GpuFloat>(
    handle: &SparseHandle,
    algo: SpMVAlgo,
    alpha: T,
    a: &CsrMatrix<T>,
    x_ptr: CUdeviceptr,
    beta: T,
    y_ptr: CUdeviceptr,
) -> SparseResult<()> {
    if a.rows() == 0 || a.cols() == 0 {
        return Ok(());
    }

    // Resolve Adaptive algorithm
    let effective_algo = match algo {
        SpMVAlgo::Adaptive => resolve_adaptive(a.avg_nnz_per_row()),
        other => other,
    };

    match effective_algo {
        SpMVAlgo::Scalar => spmv_scalar(handle, alpha, a, x_ptr, beta, y_ptr),
        SpMVAlgo::Vector => spmv_vector(handle, alpha, a, x_ptr, beta, y_ptr),
        SpMVAlgo::Adaptive => {
            // Already resolved above; unreachable
            spmv_scalar(handle, alpha, a, x_ptr, beta, y_ptr)
        }
    }
}

/// Scalar SpMV: one thread per row.
fn spmv_scalar<T: GpuFloat>(
    handle: &SparseHandle,
    alpha: T,
    a: &CsrMatrix<T>,
    x_ptr: CUdeviceptr,
    beta: T,
    y_ptr: CUdeviceptr,
) -> SparseResult<()> {
    let ptx = emit_spmv_scalar::<T>(handle.sm_version())?;
    let module = Arc::new(Module::from_ptx(&ptx)?);
    let kernel = Kernel::from_module(module, "spmv_scalar")?;

    let block_size = SPMV_SCALAR_BLOCK;
    let grid_size = grid_size_for(a.rows(), block_size);
    let params = LaunchParams::new(grid_size, block_size);

    kernel.launch(
        &params,
        handle.stream(),
        &(
            a.row_ptr().as_device_ptr(),
            a.col_idx().as_device_ptr(),
            a.values().as_device_ptr(),
            x_ptr,
            y_ptr,
            alpha.to_bits_u64(),
            beta.to_bits_u64(),
            a.rows(),
        ),
    )?;

    Ok(())
}

/// Vector SpMV: one warp (32 threads) per row.
fn spmv_vector<T: GpuFloat>(
    handle: &SparseHandle,
    alpha: T,
    a: &CsrMatrix<T>,
    x_ptr: CUdeviceptr,
    beta: T,
    y_ptr: CUdeviceptr,
) -> SparseResult<()> {
    let ptx = emit_spmv_vector::<T>(handle.sm_version())?;
    let module = Arc::new(Module::from_ptx(&ptx)?);
    let kernel = Kernel::from_module(module, "spmv_vector")?;

    let block_size = SPMV_VECTOR_BLOCK;
    // Each warp handles one row; warps_per_block = block_size / 32
    let warps_per_block = block_size / 32;
    let grid_size = grid_size_for(a.rows(), warps_per_block);
    let params = LaunchParams::new(grid_size, block_size);

    kernel.launch(
        &params,
        handle.stream(),
        &(
            a.row_ptr().as_device_ptr(),
            a.col_idx().as_device_ptr(),
            a.values().as_device_ptr(),
            x_ptr,
            y_ptr,
            alpha.to_bits_u64(),
            beta.to_bits_u64(),
            a.rows(),
        ),
    )?;

    Ok(())
}

/// Generates PTX for scalar SpMV (one thread per row).
fn emit_spmv_scalar<T: GpuFloat>(sm: SmVersion) -> SparseResult<String> {
    let elem_bytes = T::size_u32();

    KernelBuilder::new("spmv_scalar")
        .target(sm)
        .param("row_ptr", PtxType::U64)
        .param("col_idx", PtxType::U64)
        .param("values", PtxType::U64)
        .param("x_ptr", PtxType::U64)
        .param("y_ptr", PtxType::U64)
        .param("alpha_bits", PtxType::U64)
        .param("beta_bits", PtxType::U64)
        .param("num_rows", PtxType::U32)
        .body(move |b| {
            let gid = b.global_thread_id_x();
            let num_rows = b.load_param_u32("num_rows");

            let gid_inner = gid.clone();
            b.if_lt_u32(gid, num_rows, move |b| {
                let row = gid_inner;
                let row_ptr_base = b.load_param_u64("row_ptr");
                let col_idx_base = b.load_param_u64("col_idx");
                let values_base = b.load_param_u64("values");
                let x_ptr = b.load_param_u64("x_ptr");
                let y_ptr = b.load_param_u64("y_ptr");
                let alpha_bits = b.load_param_u64("alpha_bits");
                let beta_bits = b.load_param_u64("beta_bits");

                let alpha = reinterpret_bits_to_float::<T>(b, alpha_bits);
                let beta = reinterpret_bits_to_float::<T>(b, beta_bits);

                // Load row_ptr[row] and row_ptr[row+1] (i32 = 4 bytes)
                let rp_addr = b.byte_offset_addr(row_ptr_base.clone(), row.clone(), 4);
                let row_start = b.load_global_i32(rp_addr);

                let row_plus_1 = b.alloc_reg(PtxType::U32);
                b.raw_ptx(&format!("add.u32 {row_plus_1}, {row}, 1;"));
                let rp_addr_next = b.byte_offset_addr(row_ptr_base, row_plus_1, 4);
                let row_end = b.load_global_i32(rp_addr_next);

                // Initialize accumulator
                let acc = load_float_imm::<T>(b, 0.0);

                // Loop over non-zeros in this row
                let loop_label = b.fresh_label("spmv_loop");
                let done_label = b.fresh_label("spmv_done");

                let k = b.alloc_reg(PtxType::U32);
                // Convert row_start (i32) to u32
                let rs_u32 = b.alloc_reg(PtxType::U32);
                b.raw_ptx(&format!("mov.b32 {rs_u32}, {row_start};"));
                b.raw_ptx(&format!("mov.u32 {k}, {rs_u32};"));

                let re_u32 = b.alloc_reg(PtxType::U32);
                b.raw_ptx(&format!("mov.b32 {re_u32}, {row_end};"));

                b.label(&loop_label);
                // Exit the loop when k >= row_end. Use the structured `branch_if`
                // so the branch target is emitted with the same `$`-prefix
                // convention as `b.label`/`b.branch`; a raw `bra L__...` would not
                // match the `$L__...:` label definition and `ptxas` would reject it
                // ("Unknown symbol").
                let pred = b.alloc_reg(PtxType::Pred);
                b.raw_ptx(&format!("setp.hs.u32 {pred}, {k}, {re_u32};"));
                b.branch_if(pred, &done_label);

                // Load col_idx[k] (i32 = 4 bytes)
                let ci_addr = b.byte_offset_addr(col_idx_base.clone(), k.clone(), 4);
                let col = b.load_global_i32(ci_addr);
                let col_u32 = b.alloc_reg(PtxType::U32);
                b.raw_ptx(&format!("mov.b32 {col_u32}, {col};"));

                // Load values[k]
                let v_addr = b.byte_offset_addr(values_base.clone(), k.clone(), elem_bytes);
                let val = load_global_float::<T>(b, v_addr);

                // Load x[col]
                let x_addr = b.byte_offset_addr(x_ptr.clone(), col_u32, elem_bytes);
                let x_val = load_global_float::<T>(b, x_addr);

                // acc += val * x_val
                let new_acc = fma_float::<T>(b, val, x_val, acc.clone());
                let mov_suffix = ptx_suffix::<T>();
                b.raw_ptx(&format!("mov.{mov_suffix} {acc}, {new_acc};"));

                // k++
                b.raw_ptx(&format!("add.u32 {k}, {k}, 1;"));
                b.branch(&loop_label);
                b.label(&done_label);

                // Compute y = alpha * acc + beta * y_old
                let y_addr = b.byte_offset_addr(y_ptr, row, elem_bytes);
                let y_old = load_global_float::<T>(b, y_addr.clone());

                let alpha_acc = mul_float::<T>(b, alpha, acc);
                let beta_y = mul_float::<T>(b, beta, y_old);
                let result = add_float::<T>(b, alpha_acc, beta_y);

                store_global_float::<T>(b, y_addr, result);
            });

            b.ret();
        })
        .build()
        .map_err(|e| SparseError::PtxGeneration(e.to_string()))
}

/// Generates PTX for vector SpMV (one warp per row).
fn emit_spmv_vector<T: GpuFloat>(sm: SmVersion) -> SparseResult<String> {
    let elem_bytes = T::size_u32();

    KernelBuilder::new("spmv_vector")
        .target(sm)
        .param("row_ptr", PtxType::U64)
        .param("col_idx", PtxType::U64)
        .param("values", PtxType::U64)
        .param("x_ptr", PtxType::U64)
        .param("y_ptr", PtxType::U64)
        .param("alpha_bits", PtxType::U64)
        .param("beta_bits", PtxType::U64)
        .param("num_rows", PtxType::U32)
        .body(move |b| {
            // Each warp handles one row. Warp ID = global_thread_id / 32
            let tid_global = b.global_thread_id_x();
            let num_rows = b.load_param_u32("num_rows");

            // Lane within warp (0..31)
            let lane = b.alloc_reg(PtxType::U32);
            b.raw_ptx(&format!("and.b32 {lane}, {tid_global}, 31;"));

            // Warp ID = tid_global >> 5
            let warp_id = b.alloc_reg(PtxType::U32);
            b.raw_ptx(&format!("shr.u32 {warp_id}, {tid_global}, 5;"));

            let warp_id_inner = warp_id.clone();
            let lane_inner = lane.clone();
            b.if_lt_u32(warp_id, num_rows, move |b| {
                let row = warp_id_inner;
                let lane = lane_inner;

                let row_ptr_base = b.load_param_u64("row_ptr");
                let col_idx_base = b.load_param_u64("col_idx");
                let values_base = b.load_param_u64("values");
                let x_ptr = b.load_param_u64("x_ptr");
                let y_ptr = b.load_param_u64("y_ptr");
                let alpha_bits = b.load_param_u64("alpha_bits");
                let beta_bits = b.load_param_u64("beta_bits");

                let alpha = reinterpret_bits_to_float::<T>(b, alpha_bits);
                let beta = reinterpret_bits_to_float::<T>(b, beta_bits);

                // Load row bounds
                let rp_addr = b.byte_offset_addr(row_ptr_base.clone(), row.clone(), 4);
                let row_start_i32 = b.load_global_i32(rp_addr);
                let row_start = b.alloc_reg(PtxType::U32);
                b.raw_ptx(&format!("mov.b32 {row_start}, {row_start_i32};"));

                let row_plus_1 = b.alloc_reg(PtxType::U32);
                b.raw_ptx(&format!("add.u32 {row_plus_1}, {row}, 1;"));
                let rp_addr_next = b.byte_offset_addr(row_ptr_base, row_plus_1, 4);
                let row_end_i32 = b.load_global_i32(rp_addr_next);
                let row_end = b.alloc_reg(PtxType::U32);
                b.raw_ptx(&format!("mov.b32 {row_end}, {row_end_i32};"));

                // Each lane starts at row_start + lane, stride 32
                let acc = load_float_imm::<T>(b, 0.0);

                let k = b.alloc_reg(PtxType::U32);
                b.raw_ptx(&format!("add.u32 {k}, {row_start}, {lane};"));

                let loop_label = b.fresh_label("spmv_vloop");
                let done_label = b.fresh_label("spmv_vdone");

                b.label(&loop_label);
                // Exit the loop when k >= row_end via the structured `branch_if`
                // (see the scalar kernel for why a raw `bra` would be rejected).
                let pred = b.alloc_reg(PtxType::Pred);
                b.raw_ptx(&format!("setp.hs.u32 {pred}, {k}, {row_end};"));
                b.branch_if(pred, &done_label);

                // Load col and value
                let ci_addr = b.byte_offset_addr(col_idx_base.clone(), k.clone(), 4);
                let col_i32 = b.load_global_i32(ci_addr);
                let col_u32 = b.alloc_reg(PtxType::U32);
                b.raw_ptx(&format!("mov.b32 {col_u32}, {col_i32};"));

                let v_addr = b.byte_offset_addr(values_base.clone(), k.clone(), elem_bytes);
                let val = load_global_float::<T>(b, v_addr);

                let x_addr = b.byte_offset_addr(x_ptr.clone(), col_u32, elem_bytes);
                let x_val = load_global_float::<T>(b, x_addr);

                let new_acc = fma_float::<T>(b, val, x_val, acc.clone());
                let mov_suffix = ptx_suffix::<T>();
                b.raw_ptx(&format!("mov.{mov_suffix} {acc}, {new_acc};"));

                // k += 32 (warp width)
                b.raw_ptx(&format!("add.u32 {k}, {k}, 32;"));
                b.branch(&loop_label);
                b.label(&done_label);

                // Warp shuffle reduction
                let reduced = emit_warp_reduce_sum::<T>(b, acc);

                // Only lane 0 writes the row result; every other lane skips ahead.
                // Branch when lane != 0 using the structured `branch_if` so the
                // target matches the `$`-prefixed label definition.
                let not_lane_0 = b.alloc_reg(PtxType::Pred);
                b.raw_ptx(&format!("setp.ne.u32 {not_lane_0}, {lane}, 0;"));
                let skip_label = b.fresh_label("spmv_skip");
                b.branch_if(not_lane_0, &skip_label);

                let y_addr = b.byte_offset_addr(y_ptr, row, elem_bytes);
                let y_old = load_global_float::<T>(b, y_addr.clone());

                let alpha_acc = mul_float::<T>(b, alpha, reduced);
                let beta_y = mul_float::<T>(b, beta, y_old);
                let result = add_float::<T>(b, alpha_acc, beta_y);
                store_global_float::<T>(b, y_addr, result);

                b.label(&skip_label);
            });

            b.ret();
        })
        .build()
        .map_err(|e| SparseError::PtxGeneration(e.to_string()))
}

#[cfg(test)]
mod tests {
    use super::*;

    /// Runs `ptxas -arch=sm_86` on the given PTX, writing it through
    /// [`std::env::temp_dir`]. Returns `Some(())` on success, `Some` error text
    /// on assembler rejection, or `None` if `ptxas` is not on `PATH` (so callers
    /// can skip gracefully on machines without the CUDA toolkit).
    fn try_ptxas(name: &str, ptx: &str) -> Option<Result<(), String>> {
        use std::io::Write;
        let path = std::env::temp_dir().join(format!("oxicuda_spmv_{name}.ptx"));
        {
            let mut f = std::fs::File::create(&path).expect("test: create temp PTX file");
            f.write_all(ptx.as_bytes())
                .expect("test: write temp PTX file");
        }
        match std::process::Command::new("ptxas")
            .arg("-arch=sm_86")
            .arg(&path)
            .arg("-o")
            .arg("/dev/null")
            .output()
        {
            Ok(out) if out.status.success() => Some(Ok(())),
            Ok(out) => Some(Err(String::from_utf8_lossy(&out.stderr).into_owned())),
            // ptxas missing (no CUDA toolkit): skip gracefully.
            Err(_) => None,
        }
    }

    /// The f64 scalar and vector SpMV kernels must be well-formed double-precision
    /// PTX: 64-bit value registers, `0D`-prefixed f64 immediates (never an f32
    /// `0F00000000` zero), no illegal `shfl.sync.down.b64`, and — when `ptxas` is
    /// available — they must assemble for `sm_86` (RTX A4000). Regression guard for
    /// the "CUDA: invalid PTX" failure of `cuda_spmv_csr`.
    #[test]
    fn spmv_f64_ptx_well_formed_and_assembles() {
        let scalar = emit_spmv_scalar::<f64>(SmVersion::Sm86).expect("f64 scalar PTX");
        let vector = emit_spmv_vector::<f64>(SmVersion::Sm86).expect("f64 vector PTX");

        for (name, ptx) in [("scalar", &scalar), ("vector", &vector)] {
            // f64 value registers are declared 64-bit (`.b64` is the reg-class of f64).
            assert!(
                ptx.contains(".reg .b64 %f"),
                "f64 {name} kernel must declare 64-bit %f value registers:\n{ptx}"
            );
            // No f32 zero immediate and no f32-typed float math may leak into the
            // f64 value path (the GEMM-class antipattern this guards against).
            assert!(
                !ptx.contains("0F00000000"),
                "f64 {name} kernel must not materialize an f32 0.0 immediate:\n{ptx}"
            );
            assert!(
                !ptx.contains(".f32"),
                "f64 {name} kernel must not contain any .f32-typed instruction:\n{ptx}"
            );
            // The f64 zero immediate uses the 16-hex-digit 0D form.
            assert!(
                ptx.contains("0D0000000000000000"),
                "f64 {name} kernel must materialize the f64 0.0 immediate (0D…):\n{ptx}"
            );
            // shfl only supports b32; a b64 shuffle is rejected by ptxas.
            assert!(
                !ptx.contains("shfl.sync.down.b64"),
                "f64 {name} kernel must not emit shfl.sync.down.b64:\n{ptx}"
            );
            // Genuine double-precision arithmetic must be present.
            assert!(
                ptx.contains("fma.rn.f64") && ptx.contains("ld.global.f64"),
                "f64 {name} kernel must use f64 fma/load instructions:\n{ptx}"
            );
        }

        // The vector kernel performs the warp reduction via paired b32 shuffles.
        assert!(
            vector.contains("shfl.sync.down.b32"),
            "f64 vector kernel must reduce via b32 shuffles:\n{vector}"
        );

        // Assemble with ptxas when present; both kernels must be accepted by sm_86.
        for (name, ptx) in [("scalar", &scalar), ("vector", &vector)] {
            match try_ptxas(&format!("f64_{name}"), ptx) {
                Some(Ok(())) => {}
                Some(Err(stderr)) => {
                    panic!("ptxas rejected the f64 {name} SpMV kernel:\n{stderr}\nPTX:\n{ptx}")
                }
                None => {
                    // ptxas unavailable: textual well-formedness checks above stand in.
                }
            }
        }
    }

    /// The f32 SpMV kernels must keep assembling too (the label fix is shared).
    #[test]
    fn spmv_f32_ptx_assembles() {
        let scalar = emit_spmv_scalar::<f32>(SmVersion::Sm86).expect("f32 scalar PTX");
        let vector = emit_spmv_vector::<f32>(SmVersion::Sm86).expect("f32 vector PTX");
        for (name, ptx) in [("scalar", &scalar), ("vector", &vector)] {
            assert!(
                ptx.contains(".reg .b32 %f"),
                "f32 {name} kernel must declare 32-bit %f value registers:\n{ptx}"
            );
            if let Some(Err(stderr)) = try_ptxas(&format!("f32_{name}"), ptx) {
                panic!("ptxas rejected the f32 {name} SpMV kernel:\n{stderr}\nPTX:\n{ptx}");
            }
        }
    }

    #[test]
    fn spmv_algo_auto_select() {
        // avg_nnz < threshold => Scalar
        // Verify VECTOR_THRESHOLD is set to a reasonable value for algorithm selection.
        let threshold = VECTOR_THRESHOLD;
        assert!(threshold > 3.0);
    }

    #[test]
    fn spmv_scalar_ptx_generates() {
        let ptx = emit_spmv_scalar::<f32>(SmVersion::Sm80);
        assert!(ptx.is_ok());
        let ptx = ptx.expect("test: PTX gen should succeed");
        assert!(ptx.contains(".entry spmv_scalar"));
        assert!(ptx.contains(".target sm_80"));
    }

    #[test]
    fn spmv_vector_ptx_generates() {
        let ptx = emit_spmv_vector::<f32>(SmVersion::Sm80);
        assert!(ptx.is_ok());
        let ptx = ptx.expect("test: PTX gen should succeed");
        assert!(ptx.contains(".entry spmv_vector"));
    }

    #[test]
    fn spmv_scalar_ptx_f64() {
        let ptx = emit_spmv_scalar::<f64>(SmVersion::Sm80);
        assert!(ptx.is_ok());
    }

    #[test]
    fn spmv_vector_ptx_f64() {
        let ptx = emit_spmv_vector::<f64>(SmVersion::Sm80);
        assert!(ptx.is_ok());
    }

    // -----------------------------------------------------------------------
    // Task 5a: Auto-selection heuristic tests (CPU-only, no GPU required)
    // -----------------------------------------------------------------------

    /// Very sparse rows (avg_nnz ≈ 1.5, well below threshold 4.0) → Scalar.
    #[test]
    fn test_spmv_selects_scalar_for_very_sparse() {
        // 100 rows, 150 nnz → avg = 1.5
        let avg = 150.0_f64 / 100.0;
        assert!(avg < VECTOR_THRESHOLD);
        assert_eq!(resolve_adaptive(avg), SpMVAlgo::Scalar);
    }

    /// Moderate density (avg_nnz = 32, above threshold 4.0) → Vector.
    #[test]
    fn test_spmv_selects_vector_for_moderate_density() {
        let avg = 32.0_f64;
        assert!(avg >= VECTOR_THRESHOLD);
        assert_eq!(resolve_adaptive(avg), SpMVAlgo::Vector);
    }

    /// Dense rows (avg_nnz = 128, well above threshold) → Vector.
    #[test]
    fn test_spmv_selects_vector_for_dense() {
        let avg = 128.0_f64;
        assert!(avg >= VECTOR_THRESHOLD);
        assert_eq!(resolve_adaptive(avg), SpMVAlgo::Vector);
    }

    /// Boundary: just below threshold → Scalar; at threshold → Vector.
    #[test]
    fn test_spmv_selection_boundary_conditions() {
        // Just below threshold (3.9999…)
        let just_below = VECTOR_THRESHOLD - f64::EPSILON * VECTOR_THRESHOLD;
        assert_eq!(resolve_adaptive(just_below), SpMVAlgo::Scalar);

        // Exactly at threshold
        assert_eq!(resolve_adaptive(VECTOR_THRESHOLD), SpMVAlgo::Vector);

        // Slightly above threshold
        let just_above = VECTOR_THRESHOLD + f64::EPSILON * VECTOR_THRESHOLD;
        assert_eq!(resolve_adaptive(just_above), SpMVAlgo::Vector);
    }

    /// Empty matrix (0.0 avg_nnz) is handled gracefully → Scalar (no Vector wasted).
    #[test]
    fn test_spmv_selection_empty_matrix() {
        assert_eq!(resolve_adaptive(0.0), SpMVAlgo::Scalar);
    }

    /// VECTOR_THRESHOLD sanity: must equal 4.0 (the spec-defined boundary).
    #[test]
    fn test_vector_threshold_sanity() {
        assert_eq!(
            VECTOR_THRESHOLD, 4.0,
            "VECTOR_THRESHOLD must be 4.0 per spec"
        );
        assert!(VECTOR_THRESHOLD.is_finite());
    }

    // -----------------------------------------------------------------------
    // Deepening: explicit avg_nnz_per_row bracket tests matching sparse
    // matrix categories from estimation.md and architecture notes.
    // -----------------------------------------------------------------------

    /// avg_nnz_per_row ≤ 2 (diagonal / identity matrices) → Scalar kernel.
    ///
    /// Models a diagonal matrix (1 nnz/row) — the most sparse real-world case.
    #[test]
    fn test_spmv_scalar_for_diagonal_matrix() {
        // 1000-row diagonal → avg = 1.0
        let avg = 1000.0_f64 / 1000.0;
        assert!(avg <= 2.0, "avg={avg} should be ≤ 2");
        assert_eq!(
            resolve_adaptive(avg),
            SpMVAlgo::Scalar,
            "diagonal matrices (avg ≤ 2) should use Scalar SpMV"
        );
    }

    /// avg_nnz_per_row ≤ 2, fractional (near-diagonal) → Scalar kernel.
    ///
    /// Models a tridiagonal-like matrix with ~2 nnz/row.
    #[test]
    fn test_spmv_scalar_for_tridiagonal_matrix() {
        // 1000 rows, 2000 nnz → avg = 2.0 (tridiagonal boundary)
        let avg = 2000.0_f64 / 1000.0;
        assert!(avg <= 2.0, "avg={avg} should be ≤ 2");
        assert_eq!(
            resolve_adaptive(avg),
            SpMVAlgo::Scalar,
            "near-diagonal matrices (avg ≤ 2) should use Scalar SpMV"
        );
    }

    /// avg_nnz_per_row in (2, 32] (moderate stencil / FEM) → Vector kernel.
    ///
    /// Models a 5-point 2D finite-difference stencil (avg ≈ 5 nnz/row).
    #[test]
    fn test_spmv_vector_for_5pt_stencil() {
        // 1000×1000 grid → 5_000_000 rows with ~5 nnz each
        let avg = 5.0_f64;
        assert!(avg > 2.0 && avg <= 32.0, "avg={avg} should be in (2, 32]");
        assert_eq!(
            resolve_adaptive(avg),
            SpMVAlgo::Vector,
            "5-point stencil (avg ≈ 5) should use Vector SpMV"
        );
    }

    /// avg_nnz_per_row ≈ 16 (7-point 3D stencil) → Vector kernel.
    #[test]
    fn test_spmv_vector_for_7pt_3d_stencil() {
        let avg = 7.0_f64;
        assert!(avg <= 32.0, "avg={avg} should be ≤ 32");
        assert_eq!(
            resolve_adaptive(avg),
            SpMVAlgo::Vector,
            "7-point 3D stencil (avg ≈ 7) should use Vector SpMV"
        );
    }

    /// avg_nnz_per_row exactly at VECTOR_THRESHOLD boundary (4.0) → Vector.
    ///
    /// Tests that the boundary is inclusive: avg = VECTOR_THRESHOLD selects
    /// Vector, not Scalar (i.e., `>=` rather than `>`).
    #[test]
    fn test_spmv_vector_at_exact_threshold() {
        let avg = VECTOR_THRESHOLD; // 4.0
        assert_eq!(
            resolve_adaptive(avg),
            SpMVAlgo::Vector,
            "avg == VECTOR_THRESHOLD should select Vector (inclusive boundary)"
        );
        // One ULP below threshold → Scalar
        let below = VECTOR_THRESHOLD - f64::MIN_POSITIVE;
        // May still be 4.0 due to float precision, so only check if strictly below
        if below < VECTOR_THRESHOLD {
            assert_eq!(
                resolve_adaptive(below),
                SpMVAlgo::Scalar,
                "avg strictly below VECTOR_THRESHOLD should select Scalar"
            );
        }
    }

    /// avg_nnz_per_row > 32 (dense row, graph networks) → Vector kernel.
    ///
    /// In the current two-class model, any avg ≥ VECTOR_THRESHOLD selects
    /// Vector regardless of whether avg is 5 or 500. This confirms that the
    /// "Adaptive" algorithm resolves correctly for highly dense rows.
    #[test]
    fn test_spmv_vector_for_high_density_rows() {
        // avg = 64: above the ≤ 32 bracket, still selects Vector
        let avg_64 = 64.0_f64;
        assert_eq!(
            resolve_adaptive(avg_64),
            SpMVAlgo::Vector,
            "high-density rows (avg = 64) should use Vector SpMV via Adaptive"
        );

        // avg = 256: very dense (near-dense matrix)
        let avg_256 = 256.0_f64;
        assert_eq!(
            resolve_adaptive(avg_256),
            SpMVAlgo::Vector,
            "near-dense rows (avg = 256) should use Vector SpMV via Adaptive"
        );
    }

    /// Adaptive algo resolves to the same result as calling resolve_adaptive
    /// directly for various avg_nnz values. Confirms SpMVAlgo::Adaptive is
    /// not accidentally treated as a concrete kernel variant.
    #[test]
    fn test_spmv_adaptive_algo_is_not_concrete() {
        // SpMVAlgo::Adaptive is a selection hint, not a concrete kernel.
        // resolve_adaptive must return Scalar or Vector, never Adaptive.
        let test_avgs = [0.0, 0.5, 1.0, 2.0, 3.99, 4.0, 4.01, 32.0, 64.0, 128.0];
        for avg in test_avgs {
            let resolved = resolve_adaptive(avg);
            assert!(
                matches!(resolved, SpMVAlgo::Scalar | SpMVAlgo::Vector),
                "resolve_adaptive({avg}) returned {resolved:?}, expected Scalar or Vector"
            );
        }
    }

    // -----------------------------------------------------------------------
    // Quality gate: CSR-Vector warp shuffle reduction simulation (CPU)
    // -----------------------------------------------------------------------

    /// Simulate a single-warp (32 threads) tree reduction of partial dot-products.
    ///
    /// In the Vector SpMV kernel each warp computes partial sums for the row
    /// elements it handles, then performs a binary tree (warp-shuffle) reduction
    /// to sum all 32 partial sums into a single row result.
    ///
    /// This test verifies the correctness of that reduction algorithm on the CPU.
    #[test]
    fn spmv_warp_reduction_sim_32_threads() {
        // 32 partial sums (one per thread in a warp)
        let partial: Vec<f64> = (0..32_u32).map(|i| f64::from(i * i + 1)).collect();
        let naive_sum: f64 = partial.iter().sum();

        // Simulate binary tree reduction (warp shuffle pattern):
        // stride 16, 8, 4, 2, 1
        let mut sums = partial.clone();
        let mut active = 32_usize;
        while active > 1 {
            let half = active / 2;
            for lane in 0..half {
                sums[lane] += sums[lane + half];
            }
            active = half;
        }
        let tree_sum = sums[0];

        assert!(
            (tree_sum - naive_sum).abs() < 1e-9,
            "Warp tree reduction ({tree_sum}) must match naive sum ({naive_sum})"
        );
    }

    /// Simulate a half-warp (16 threads) tree reduction.
    ///
    /// Verifies reduction correctness for the half-warp code path used when
    /// the row is shorter than a full warp.
    #[test]
    fn spmv_half_warp_reduction_sim_16_threads() {
        let partial: Vec<f64> = (0..16_u32).map(|i| f64::from(2 * i + 3)).collect();
        let naive_sum: f64 = partial.iter().sum();

        let mut sums = partial.clone();
        let mut active = 16_usize;
        while active > 1 {
            let half = active / 2;
            for lane in 0..half {
                sums[lane] += sums[lane + half];
            }
            active = half;
        }
        let tree_sum = sums[0];

        assert!(
            (tree_sum - naive_sum).abs() < 1e-9,
            "Half-warp tree reduction ({tree_sum}) must match naive sum ({naive_sum})"
        );
    }

    // -----------------------------------------------------------------------
    // Quality gate: SpMV numerical accuracy vs dense reference (CPU simulation)
    // -----------------------------------------------------------------------

    /// Dense-reference SpMV: computes y = A * x for a general dense matrix.
    fn dense_spmv(a_rows: usize, a_cols: usize, a: &[f64], x: &[f64]) -> Vec<f64> {
        let mut y = vec![0.0_f64; a_rows];
        for i in 0..a_rows {
            for j in 0..a_cols {
                y[i] += a[i * a_cols + j] * x[j];
            }
        }
        y
    }

    /// CSR SpMV simulation: computes y = A_csr * x on the CPU.
    fn csr_spmv_sim(
        nrows: usize,
        row_ptr: &[usize],
        col_idx: &[usize],
        values: &[f64],
        x: &[f64],
    ) -> Vec<f64> {
        let mut y = vec![0.0_f64; nrows];
        for i in 0..nrows {
            for idx in row_ptr[i]..row_ptr[i + 1] {
                y[i] += values[idx] * x[col_idx[idx]];
            }
        }
        y
    }

    /// SpMV for 4×4 identity matrix: y = I * x must equal x.
    ///
    /// This is the simplest correctness test: the identity provides a known
    /// reference where every output equals the corresponding input.
    #[test]
    fn spmv_numerical_accuracy_identity_4x4() {
        let n = 4_usize;
        // Identity matrix in CSR format
        let row_ptr = vec![0, 1, 2, 3, 4];
        let col_idx = vec![0, 1, 2, 3];
        let values = vec![1.0_f64; n];
        let x = vec![1.0_f64, 2.0, 3.0, 4.0];

        let y_csr = csr_spmv_sim(n, &row_ptr, &col_idx, &values, &x);
        let y_dense = dense_spmv(
            n,
            n,
            &[
                1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0,
            ],
            &x,
        );
        for i in 0..n {
            assert!(
                (y_csr[i] - y_dense[i]).abs() < 1e-13,
                "SpMV I×x: y_csr[{i}]={} != y_dense[{i}]={}",
                y_csr[i],
                y_dense[i],
            );
        }
    }

    /// SpMV for a 0.1% sparse 1000×1000 matrix with a known diagonal pattern.
    ///
    /// Only diagonal entries are set (1000 out of 1_000_000 possible entries = 0.1%).
    /// Result must equal x (diagonal matrix with ones = identity).
    #[test]
    fn spmv_very_sparse_0_1_percent_1000x1000() {
        let n = 1000_usize;
        // Diagonal matrix (0.1% density)
        let row_ptr: Vec<usize> = (0..=n).collect();
        let col_idx: Vec<usize> = (0..n).collect();
        let values: Vec<f64> = vec![2.0; n]; // diagonal value = 2
        let x: Vec<f64> = (0..n).map(|i| i as f64 * 0.001 + 1.0).collect();

        let y = csr_spmv_sim(n, &row_ptr, &col_idx, &values, &x);

        for i in 0..n {
            let expected = 2.0 * x[i];
            assert!(
                (y[i] - expected).abs() < 1e-10,
                "0.1% sparse SpMV row {i}: got {}, expected {expected}",
                y[i],
            );
        }
    }

    /// SpMV for a 10% sparse 100×100 matrix: banded structure.
    ///
    /// Uses a banded matrix with bandwidth 5 (5 non-zeros per row on average),
    /// giving approximately 10% density for a 100×100 system.
    #[test]
    fn spmv_moderate_10_percent_100x100() {
        let n = 100_usize;
        let bandwidth = 5_usize; // ±2 off-diagonal + diagonal

        let mut row_ptr = vec![0_usize; n + 1];
        let mut col_idx = Vec::new();
        let mut values = Vec::new();

        for i in 0..n {
            let start = i.saturating_sub(2);
            let end = (i + 3).min(n);
            for j in start..end {
                col_idx.push(j);
                values.push(if i == j { 4.0_f64 } else { -1.0 });
            }
            row_ptr[i + 1] = col_idx.len();
        }
        let _ = bandwidth; // document variable used in comments

        // x = [1, 1, 1, ..., 1]
        let x = vec![1.0_f64; n];
        let y_csr = csr_spmv_sim(n, &row_ptr, &col_idx, &values, &x);

        // Build dense matrix and compute reference
        let mut a_dense = vec![0.0_f64; n * n];
        for i in 0..n {
            let start = i.saturating_sub(2);
            let end = (i + 3).min(n);
            for j in start..end {
                a_dense[i * n + j] = if i == j { 4.0 } else { -1.0 };
            }
        }
        let y_dense = dense_spmv(n, n, &a_dense, &x);

        for i in 0..n {
            assert!(
                (y_csr[i] - y_dense[i]).abs() < 1e-10,
                "10% sparse SpMV row {i}: got {}, expected {}",
                y_csr[i],
                y_dense[i],
            );
        }
    }

    // -----------------------------------------------------------------------
    // Quality gate: auto-format selection thresholds
    // -----------------------------------------------------------------------

    /// Verify format selection thresholds cover the three density regimes.
    ///
    /// - avg_nnz < VECTOR_THRESHOLD (4.0) → Scalar
    /// - avg_nnz >= VECTOR_THRESHOLD       → Vector
    ///
    /// The test explicitly checks the three named brackets from the spec.
    #[test]
    fn spmv_format_selection_three_brackets() {
        // Very sparse (≤ 2): diagonal-like
        assert_eq!(
            resolve_adaptive(1.0),
            SpMVAlgo::Scalar,
            "avg_nnz=1.0 (≤ 2 bracket) must select Scalar"
        );
        assert_eq!(
            resolve_adaptive(2.0),
            SpMVAlgo::Scalar,
            "avg_nnz=2.0 (≤ 2 bracket) must select Scalar"
        );
        // Moderate (≤ 64): stencil-like (above VECTOR_THRESHOLD)
        assert_eq!(
            resolve_adaptive(5.0),
            SpMVAlgo::Vector,
            "avg_nnz=5.0 (≤ 64 bracket) must select Vector"
        );
        assert_eq!(
            resolve_adaptive(32.0),
            SpMVAlgo::Vector,
            "avg_nnz=32.0 (≤ 64 bracket) must select Vector"
        );
        // Dense (> 64): near-dense graph
        assert_eq!(
            resolve_adaptive(65.0),
            SpMVAlgo::Vector,
            "avg_nnz=65.0 (> 64 bracket) must select Vector (binary model)"
        );
        assert_eq!(
            resolve_adaptive(256.0),
            SpMVAlgo::Vector,
            "avg_nnz=256.0 (> 64 bracket) must select Vector"
        );
    }

    /// CPU-proxy throughput benchmark: SpMV on a synthetic 10k×10k 5-point stencil matrix.
    ///
    /// Simulates the type of sparse matrix found in the SuiteSparse collection.
    /// Measures CPU reference throughput and reports GFLOPS as a structural
    /// lower-bound for the GPU target (cuSPARSE comparison requires real hardware).
    #[test]
    fn spmv_suitesparse_proxy_throughput_10k() {
        // 2D 5-point Laplacian stencil on a 100×100 grid → 10k×10k sparse matrix.
        // Each interior row has 5 non-zeros; boundary rows have 3–4.
        let grid = 100_usize;
        let n = grid * grid; // 10_000 rows

        let mut row_ptr: Vec<usize> = Vec::with_capacity(n + 1);
        let mut col_idx: Vec<usize> = Vec::new();
        let mut values: Vec<f64> = Vec::new();

        row_ptr.push(0);
        for row in 0..n {
            let r = row / grid;
            let c = row % grid;
            // North neighbour
            if r > 0 {
                col_idx.push(row - grid);
                values.push(-1.0);
            }
            // West neighbour
            if c > 0 {
                col_idx.push(row - 1);
                values.push(-1.0);
            }
            // Self (diagonal = 4)
            col_idx.push(row);
            values.push(4.0);
            // East neighbour
            if c + 1 < grid {
                col_idx.push(row + 1);
                values.push(-1.0);
            }
            // South neighbour
            if r + 1 < grid {
                col_idx.push(row + grid);
                values.push(-1.0);
            }
            row_ptr.push(col_idx.len());
        }

        let nnz = col_idx.len();
        let x: Vec<f64> = (0..n).map(|i| (i as f64) * 0.0001 + 1.0).collect();

        // Warm-up pass
        let _ = csr_spmv_sim(n, &row_ptr, &col_idx, &values, &x);

        const ITERS: usize = 10;
        let start = std::time::Instant::now();
        let mut y = vec![0.0_f64; n];
        for _ in 0..ITERS {
            y = csr_spmv_sim(n, &row_ptr, &col_idx, &values, &x);
        }
        let elapsed_ns = start.elapsed().as_nanos() as f64;

        // 2 flops per non-zero (one multiply + one add)
        let total_flops = 2.0 * nnz as f64 * ITERS as f64;
        let gflops = total_flops / elapsed_ns; // (flops) / (ns) = GFlops/s

        println!(
            "SpMV SuiteSparse proxy (10k×10k 5-pt stencil, {} nnz, {} iters): {:.3} GFLOPS (CPU reference)",
            nnz, ITERS, gflops
        );

        // Sanity: result must be non-zero and throughput must be measurable
        assert!(y[n / 2] != 0.0, "SpMV result must be non-zero");
        assert!(
            gflops > 0.001,
            "SpMV CPU reference throughput unrealistically low: {:.6} GFLOPS",
            gflops
        );
    }
}

// ---------------------------------------------------------------------------
// On-device numeric validation (feature = "gpu-tests")
// ---------------------------------------------------------------------------

#[cfg(all(test, feature = "gpu-tests"))]
mod gpu_device_tests {
    use super::*;
    use crate::gpu_test_support::{assert_close, gpu_handle};
    use crate::host_csr::{f64_to_gpu, gpu_to_f64};
    use oxicuda_memory::DeviceBuffer;

    /// CPU oracle for `y = alpha * A * x + beta * y0` over a CSR matrix
    /// (row count is derived from `row_ptr`).
    fn cpu_csr_spmv(
        row_ptr: &[i32],
        col_idx: &[i32],
        values: &[f64],
        x: &[f64],
        y0: &[f64],
        alpha: f64,
        beta: f64,
    ) -> Vec<f64> {
        let rows = row_ptr.len() - 1;
        let mut y = vec![0.0_f64; rows];
        for (i, slot) in y.iter_mut().enumerate() {
            let start = row_ptr[i] as usize;
            let end = row_ptr[i + 1] as usize;
            let mut acc = 0.0_f64;
            for k in start..end {
                acc += values[k] * x[col_idx[k] as usize];
            }
            *slot = alpha * acc + beta * y0[i];
        }
        y
    }

    /// Drive the production `spmv` op for one element type and compare to the
    /// CPU oracle.
    #[allow(clippy::too_many_arguments)]
    fn run_spmv<T: GpuFloat>(
        algo: SpMVAlgo,
        rows: u32,
        cols: u32,
        row_ptr: &[i32],
        col_idx: &[i32],
        values: &[f64],
        x: &[f64],
        y0: &[f64],
        alpha: f64,
        beta: f64,
        tol: f64,
        tag: &str,
    ) {
        let Some(handle) = gpu_handle() else {
            return;
        };
        let dev_values: Vec<T> = values.iter().map(|&v| f64_to_gpu::<T>(v)).collect();
        let a = CsrMatrix::<T>::from_host(rows, cols, row_ptr, col_idx, &dev_values)
            .expect("test: build CSR");
        let dev_x: Vec<T> = x.iter().map(|&v| f64_to_gpu::<T>(v)).collect();
        let dev_y: Vec<T> = y0.iter().map(|&v| f64_to_gpu::<T>(v)).collect();
        let x_buf = DeviceBuffer::from_host(&dev_x).expect("test: upload x");
        let y_buf = DeviceBuffer::from_host(&dev_y).expect("test: upload y");

        spmv::<T>(
            &handle,
            algo,
            f64_to_gpu::<T>(alpha),
            &a,
            x_buf.as_device_ptr(),
            f64_to_gpu::<T>(beta),
            y_buf.as_device_ptr(),
        )
        .expect("test: spmv launch");
        handle.stream().synchronize().expect("test: sync");

        let mut out = vec![T::gpu_zero(); rows as usize];
        y_buf.copy_to_host(&mut out).expect("test: download y");
        let got: Vec<f64> = out.iter().map(|&v| gpu_to_f64(v)).collect();
        let want = cpu_csr_spmv(row_ptr, col_idx, values, x, y0, alpha, beta);
        assert_close(&got, &want, tol, tag);
    }

    /// 4x4 symmetric tridiagonal-ish matrix with a dense-ish last row.
    fn matrix_4x4() -> (u32, u32, Vec<i32>, Vec<i32>, Vec<f64>) {
        // [ 2 -1  0  0]
        // [-1  2 -1  0]
        // [ 0 -1  2 -1]
        // [ 3  0 -1  4]
        let row_ptr = vec![0, 2, 5, 8, 11];
        let col_idx = vec![0, 1, 0, 1, 2, 1, 2, 3, 0, 2, 3];
        let values = vec![2.0, -1.0, -1.0, 2.0, -1.0, -1.0, 2.0, -1.0, 3.0, -1.0, 4.0];
        (4, 4, row_ptr, col_idx, values)
    }

    /// Wider matrix (5x5, ~8 nnz/row average via banded structure) to exercise
    /// the warp/vector kernel path.
    fn matrix_6x6_banded() -> (u32, u32, Vec<i32>, Vec<i32>, Vec<f64>) {
        let n = 6usize;
        let mut row_ptr = vec![0i32];
        let mut col_idx = Vec::new();
        let mut values = Vec::new();
        for i in 0..n {
            let lo = i.saturating_sub(2);
            let hi = (i + 3).min(n);
            for j in lo..hi {
                col_idx.push(j as i32);
                values.push(if i == j { 5.0 } else { -1.0 + 0.1 * (i as f64) });
            }
            row_ptr.push(col_idx.len() as i32);
        }
        (n as u32, n as u32, row_ptr, col_idx, values)
    }

    #[test]
    fn spmv_scalar_f64_alpha_beta() {
        let (r, c, rp, ci, v) = matrix_4x4();
        let x = vec![1.0, 2.0, 3.0, 4.0];
        let y0 = vec![10.0, 20.0, 30.0, 40.0];
        run_spmv::<f64>(
            SpMVAlgo::Scalar,
            r,
            c,
            &rp,
            &ci,
            &v,
            &x,
            &y0,
            2.5,
            -0.75,
            1e-10,
            "spmv_scalar_f64",
        );
    }

    #[test]
    fn spmv_vector_f64_alpha_beta() {
        let (r, c, rp, ci, v) = matrix_6x6_banded();
        let x: Vec<f64> = (0..r as usize).map(|i| 0.5 + i as f64).collect();
        let y0: Vec<f64> = (0..r as usize).map(|i| 100.0 - i as f64).collect();
        run_spmv::<f64>(
            SpMVAlgo::Vector,
            r,
            c,
            &rp,
            &ci,
            &v,
            &x,
            &y0,
            1.5,
            0.25,
            1e-10,
            "spmv_vector_f64",
        );
    }

    #[test]
    fn spmv_scalar_f32_alpha_beta() {
        let (r, c, rp, ci, v) = matrix_4x4();
        let x = vec![1.0, 2.0, 3.0, 4.0];
        let y0 = vec![10.0, 20.0, 30.0, 40.0];
        run_spmv::<f32>(
            SpMVAlgo::Scalar,
            r,
            c,
            &rp,
            &ci,
            &v,
            &x,
            &y0,
            2.0,
            0.5,
            1e-4,
            "spmv_scalar_f32",
        );
    }

    #[test]
    fn spmv_vector_f32_alpha_beta() {
        let (r, c, rp, ci, v) = matrix_6x6_banded();
        let x: Vec<f64> = (0..r as usize).map(|i| 0.5 + i as f64).collect();
        let y0: Vec<f64> = (0..r as usize).map(|i| 7.0 + i as f64).collect();
        run_spmv::<f32>(
            SpMVAlgo::Vector,
            r,
            c,
            &rp,
            &ci,
            &v,
            &x,
            &y0,
            1.25,
            -0.5,
            1e-4,
            "spmv_vector_f32",
        );
    }

    #[test]
    fn spmv_beta_zero_overwrites_garbage() {
        // beta = 0 must fully overwrite the prior y (incl. any stale content).
        let (r, c, rp, ci, v) = matrix_4x4();
        let x = vec![1.0, 1.0, 1.0, 1.0];
        let y0 = vec![1e9, -1e9, 5e8, -5e8];
        run_spmv::<f64>(
            SpMVAlgo::Scalar,
            r,
            c,
            &rp,
            &ci,
            &v,
            &x,
            &y0,
            1.0,
            0.0,
            1e-10,
            "spmv_beta_zero",
        );
    }
}