perf(prover): device-only preprocessed tables and GPU commits for mid-size tables

Makefile

@@ -565,6 +565,8 @@ test-cuda-integration:
 test-cuda-fallback:
 	cargo test -p lambda-vm-prover --release --features test-cuda-faults \
 	    --test cuda_fallback_tests -- --ignored --nocapture --test-threads=1
+	cargo test -p lambda-vm-prover --release --features lambda-vm-prover/cuda \
+	    --test gpu_force_downgrade -- --ignored --nocapture --test-threads=1
 
 # The prover/stark/crypto/ecsm test suite with the GPU (cuda) path enabled (requires NVIDIA
 # GPU + nvcc). The GPU CI counterpart of CPU CI's sharded prover tests. Single-threaded: the

crypto/math-cuda/build.rs

@@ -72,6 +72,13 @@ fn to_real_arch(arch: &str) -> String {
     }
 }
 
+/// Single source for the barycentric multi-kernel eval-point cap. The CUDA
+/// side sizes a per-thread accumulator array with it (`BARY_MAX_K`, passed via
+/// `-D` below) and the Rust dispatch asserts against it (generated into
+/// `bary_consts.rs`) — defining it twice invites stack corruption in the
+/// kernel the day one side moves without the other.
+const BARY_MAX_EVAL_POINTS: usize = 8;
+
 fn compile_kernel(src: &str, out_name: &str, have_nvcc: bool) {
     let manifest_dir = PathBuf::from(env::var("CARGO_MANIFEST_DIR").unwrap());
     let out_dir = PathBuf::from(env::var("OUT_DIR").unwrap());
@@ -118,6 +125,7 @@ fn compile_kernel(src: &str, out_name: &str, have_nvcc: bool) {
 
     let mut cmd = Command::new(nvcc_path());
     cmd.args(["--cubin", "-O3", "-std=c++17", "-arch", &arch]);
+    cmd.arg(format!("-DBARY_MAX_K={BARY_MAX_EVAL_POINTS}"));
     // SASS→source line mapping for Nsight Compute. Unlike -G this does not
     // change codegen, but keep it opt-in so production cubins stay byte-stable.
     if env::var("LAMBDA_VM_NVCC_LINEINFO").is_ok_and(|v| v != "0" && !v.is_empty()) {
@@ -136,6 +144,19 @@ fn compile_kernel(src: &str, out_name: &str, have_nvcc: bool) {
 }
 
 fn main() {
+    // Rust-side mirror of the kernel cap; see BARY_MAX_EVAL_POINTS above.
+    let out_dir = PathBuf::from(env::var("OUT_DIR").unwrap());
+    fs::write(
+        out_dir.join("bary_consts.rs"),
+        format!(
+            "/// Compile-time cap of the multi kernels' per-thread accumulator array\n\
+             /// (`BARY_MAX_K` in barycentric.cu — single-sourced from build.rs).\n\
+             /// Callers with more evaluation points fall back to the per-point kernels.\n\
+             pub const BARY_MAX_EVAL_POINTS: usize = {BARY_MAX_EVAL_POINTS};\n"
+        ),
+    )
+    .expect("failed to write bary_consts.rs");
+
     // Headers aren't compiled, so emit rerun-if-changed to rebuild on
     // header edits.
     println!("cargo:rerun-if-changed=kernels/goldilocks.cuh");

crypto/math-cuda/kernels/barycentric.cu

@@ -191,6 +191,137 @@ extern "C" __global__ void barycentric_ext3_batched_strided(
     }
 }
 
+// Multi-eval-point + row-chunked barycentric. Two fixes over the *_strided
+// kernels above: (1) the LDE column data is read ONCE for all K evaluation
+// points (K inv_denom blocks, K accumulators) instead of once per point, and
+// (2) each column is split into `num_chunks` row ranges so the grid is
+// `num_cols * num_chunks` blocks instead of `num_cols` — the single-block-per-
+// column grid left most SMs idle at typical column counts. Blocks emit partial
+// sums; `barycentric_combine_partials` folds the chunk axis.
+//
+// `inv_denoms` holds K contiguous blocks of 3N u64 (ext3 interleaved), one per
+// evaluation point — the layout `compute_and_invert_denoms_ext3_dev` already
+// produces. Partials layout: `[(k*num_cols + col)*num_chunks + chunk]` ext3
+// interleaved, so the combine pass reads each (k, col)'s chunks contiguously.
+#ifndef BARY_MAX_K
+#error "BARY_MAX_K must be passed by build.rs (-DBARY_MAX_K=...) — single-sourced there"
+#endif
+
+extern "C" __global__ void barycentric_base_strided_multi(
+    const uint64_t *columns,
+    uint64_t col_stride,
+    uint64_t row_stride,
+    const uint64_t *coset_points,
+    const uint64_t *inv_denoms,
+    uint64_t n,
+    uint64_t k_points,
+    uint64_t num_chunks,
+    uint64_t *partials
+) {
+    uint64_t col = blockIdx.x;
+    uint64_t chunk = blockIdx.y;
+    const uint64_t *col_data = columns + col * col_stride;
+    uint64_t chunk_len = (n + num_chunks - 1) / num_chunks;
+    uint64_t start = chunk * chunk_len;
+    uint64_t end = start + chunk_len < n ? start + chunk_len : n;
+
+    ext3::Fe3 acc[BARY_MAX_K];
+    for (uint32_t k = 0; k < k_points; ++k) acc[k] = ext3::zero();
+
+    for (uint64_t i = start + threadIdx.x; i < end; i += BARY_BLOCK_DIM) {
+        uint64_t eval  = col_data[i * row_stride];
+        uint64_t point = coset_points[i];
+        uint64_t pe    = goldilocks::mul(point, eval);
+        for (uint32_t k = 0; k < k_points; ++k) {
+            const uint64_t *inv = inv_denoms + (uint64_t)k * 3 * n + i * 3;
+            ext3::Fe3 inv_d = ext3::make(inv[0], inv[1], inv[2]);
+            acc[k] = ext3::add(acc[k], ext3::mul_base(inv_d, pe));
+        }
+    }
+
+    for (uint32_t k = 0; k < k_points; ++k) {
+        ext3::Fe3 sum = block_reduce_ext3(acc[k]);
+        if (threadIdx.x == 0) {
+            uint64_t o = ((k * gridDim.x + col) * num_chunks + chunk) * 3;
+            partials[o + 0] = sum.a;
+            partials[o + 1] = sum.b;
+            partials[o + 2] = sum.c;
+        }
+        // block_reduce_ext3 reuses its shared buffers: every thread must be
+        // done reading round k's result before round k+1 overwrites them.
+        __syncthreads();
+    }
+}
+
+extern "C" __global__ void barycentric_ext3_strided_multi(
+    const uint64_t *columns,
+    uint64_t col_stride,
+    uint64_t row_stride,
+    const uint64_t *coset_points,
+    const uint64_t *inv_denoms,
+    uint64_t n,
+    uint64_t k_points,
+    uint64_t num_chunks,
+    uint64_t *partials
+) {
+    uint64_t col = blockIdx.x;
+    uint64_t chunk = blockIdx.y;
+    const uint64_t *slab_a = columns + (col * 3 + 0) * col_stride;
+    const uint64_t *slab_b = columns + (col * 3 + 1) * col_stride;
+    const uint64_t *slab_c = columns + (col * 3 + 2) * col_stride;
+    uint64_t chunk_len = (n + num_chunks - 1) / num_chunks;
+    uint64_t start = chunk * chunk_len;
+    uint64_t end = start + chunk_len < n ? start + chunk_len : n;
+
+    ext3::Fe3 acc[BARY_MAX_K];
+    for (uint32_t k = 0; k < k_points; ++k) acc[k] = ext3::zero();
+
+    for (uint64_t i = start + threadIdx.x; i < end; i += BARY_BLOCK_DIM) {
+        uint64_t lde_i = i * row_stride;
+        ext3::Fe3 eval = ext3::make(slab_a[lde_i], slab_b[lde_i], slab_c[lde_i]);
+        uint64_t point = coset_points[i];
+        ext3::Fe3 pe = ext3::mul_base(eval, point);
+        for (uint32_t k = 0; k < k_points; ++k) {
+            const uint64_t *inv = inv_denoms + (uint64_t)k * 3 * n + i * 3;
+            ext3::Fe3 inv_d = ext3::make(inv[0], inv[1], inv[2]);
+            acc[k] = ext3::add(acc[k], ext3::mul(pe, inv_d));
+        }
+    }
+
+    for (uint32_t k = 0; k < k_points; ++k) {
+        ext3::Fe3 sum = block_reduce_ext3(acc[k]);
+        if (threadIdx.x == 0) {
+            uint64_t o = ((k * gridDim.x + col) * num_chunks + chunk) * 3;
+            partials[o + 0] = sum.a;
+            partials[o + 1] = sum.b;
+            partials[o + 2] = sum.c;
+        }
+        __syncthreads();
+    }
+}
+
+// Fold the chunk axis of the multi kernels' partials: one thread per
+// (k, col) pair sums its `num_chunks` ext3 partials sequentially (the whole
+// buffer is tiny — K * cols * chunks). Output `out_ext3_int[k*num_cols+col]`,
+// same per-column layout as the single-point kernels, K blocks concatenated.
+extern "C" __global__ void barycentric_combine_partials(
+    const uint64_t *partials,
+    uint64_t num_chunks,
+    uint64_t total,
+    uint64_t *out_ext3_int
+) {
+    uint64_t idx = (uint64_t)blockIdx.x * blockDim.x + threadIdx.x;
+    if (idx >= total) return;
+    const uint64_t *row = partials + idx * num_chunks * 3;
+    ext3::Fe3 acc = ext3::zero();
+    for (uint64_t c = 0; c < num_chunks; ++c) {
+        acc = ext3::add(acc, ext3::make(row[c * 3 + 0], row[c * 3 + 1], row[c * 3 + 2]));
+    }
+    out_ext3_int[idx * 3 + 0] = acc.a;
+    out_ext3_int[idx * 3 + 1] = acc.b;
+    out_ext3_int[idx * 3 + 2] = acc.c;
+}
+
 // Gather full rows from a device-resident base-field LDE (`buf[col*col_stride +
 // row]`). One block per gathered row, threads stride over columns. Output is
 // row-major `out[q*num_cols + col]` for gathered-row slot `q` — directly the

crypto/math-cuda/src/barycentric.rs

@@ -350,6 +350,163 @@ pub fn barycentric_ext3_on_device_with_dev_inv_denoms(
     Ok(out)
 }
 
+include!(concat!(env!("OUT_DIR"), "/bary_consts.rs"));
+
+/// Row-chunk count for the multi kernels: enough `cols * chunks` blocks to
+/// occupy the device, without shrinking a chunk's row range below the point
+/// where launch + combine overhead dominates.
+fn bary_num_chunks(num_cols: usize, n: usize) -> usize {
+    let by_occupancy = (2048 / num_cols.max(1)).max(1);
+    let by_rows = (n / 8192).max(1);
+    by_occupancy.min(by_rows).min(64)
+}
+
+/// Multi-eval-point counterpart of
+/// [`barycentric_base_on_device_with_dev_inv_denoms`]: one pass over the LDE
+/// column data computes the barycentric sums for ALL `k_points` evaluation
+/// points (their inv_denom blocks live contiguously in `inv_denoms_dev`, the
+/// layout `compute_and_invert_denoms_ext3_dev` produces). Returns
+/// `3 * k_points * num_cols` u64: `k_points` concatenated per-column blocks,
+/// each in the same layout as the single-point kernels.
+pub fn barycentric_base_multi_on_device(
+    stream: &Arc<CudaStream>,
+    main_handle: &GpuLdeBase,
+    row_stride: usize,
+    coset_points_dev: &CudaSlice<u64>,
+    inv_denoms_dev: &CudaSlice<u64>,
+    n: usize,
+    k_points: usize,
+) -> Result<Vec<u64>> {
+    main_handle.wait_ready_on(stream)?;
+    assert!((1..=BARY_MAX_EVAL_POINTS).contains(&k_points));
+    assert!(coset_points_dev.len() >= n);
+    assert!(inv_denoms_dev.len() >= k_points * 3 * n);
+    let num_cols = main_handle.m;
+    if num_cols == 0 || n == 0 {
+        return Ok(vec![0; 3 * k_points * num_cols]);
+    }
+    let be = backend()?;
+    let num_chunks = bary_num_chunks(num_cols, n);
+    let total = k_points * num_cols;
+    let mut partials = stream.alloc_zeros::<u64>(total * num_chunks * 3)?;
+    let mut out_dev = stream.alloc_zeros::<u64>(3 * total)?;
+    let points_view = coset_points_dev.slice(0..n);
+    let inv_view = inv_denoms_dev.slice(0..k_points * 3 * n);
+
+    let col_stride_u64 = main_handle.lde_size as u64;
+    let row_stride_u64 = row_stride as u64;
+    let n_u64 = n as u64;
+    let k_u64 = k_points as u64;
+    let chunks_u64 = num_chunks as u64;
+    let total_u64 = total as u64;
+    let cfg = LaunchConfig {
+        grid_dim: (num_cols as u32, num_chunks as u32, 1),
+        block_dim: (BLOCK_DIM, 1, 1),
+        shared_mem_bytes: 0,
+    };
+    unsafe {
+        stream
+            .launch_builder(&be.barycentric_base_strided_multi)
+            .arg(main_handle.buf.as_ref())
+            .arg(&col_stride_u64)
+            .arg(&row_stride_u64)
+            .arg(&points_view)
+            .arg(&inv_view)
+            .arg(&n_u64)
+            .arg(&k_u64)
+            .arg(&chunks_u64)
+            .arg(&mut partials)
+            .launch(cfg)?;
+    }
+    let combine_cfg = LaunchConfig {
+        grid_dim: (total.div_ceil(BLOCK_DIM as usize) as u32, 1, 1),
+        block_dim: (BLOCK_DIM, 1, 1),
+        shared_mem_bytes: 0,
+    };
+    unsafe {
+        stream
+            .launch_builder(&be.barycentric_combine_partials)
+            .arg(&partials)
+            .arg(&chunks_u64)
+            .arg(&total_u64)
+            .arg(&mut out_dev)
+            .launch(combine_cfg)?;
+    }
+    let out = stream.clone_dtoh(&out_dev)?;
+    stream.synchronize()?;
+    Ok(out)
+}
+
+/// Ext3 counterpart of [`barycentric_base_multi_on_device`].
+pub fn barycentric_ext3_multi_on_device(
+    stream: &Arc<CudaStream>,
+    aux_handle: &GpuLdeExt3,
+    row_stride: usize,
+    coset_points_dev: &CudaSlice<u64>,
+    inv_denoms_dev: &CudaSlice<u64>,
+    n: usize,
+    k_points: usize,
+) -> Result<Vec<u64>> {
+    aux_handle.wait_ready_on(stream)?;
+    assert!((1..=BARY_MAX_EVAL_POINTS).contains(&k_points));
+    assert!(coset_points_dev.len() >= n);
+    assert!(inv_denoms_dev.len() >= k_points * 3 * n);
+    let num_cols = aux_handle.m;
+    if num_cols == 0 || n == 0 {
+        return Ok(vec![0; 3 * k_points * num_cols]);
+    }
+    let be = backend()?;
+    let num_chunks = bary_num_chunks(num_cols, n);
+    let total = k_points * num_cols;
+    let mut partials = stream.alloc_zeros::<u64>(total * num_chunks * 3)?;
+    let mut out_dev = stream.alloc_zeros::<u64>(3 * total)?;
+    let points_view = coset_points_dev.slice(0..n);
+    let inv_view = inv_denoms_dev.slice(0..k_points * 3 * n);
+
+    let col_stride_u64 = aux_handle.lde_size as u64;
+    let row_stride_u64 = row_stride as u64;
+    let n_u64 = n as u64;
+    let k_u64 = k_points as u64;
+    let chunks_u64 = num_chunks as u64;
+    let total_u64 = total as u64;
+    let cfg = LaunchConfig {
+        grid_dim: (num_cols as u32, num_chunks as u32, 1),
+        block_dim: (BLOCK_DIM, 1, 1),
+        shared_mem_bytes: 0,
+    };
+    unsafe {
+        stream
+            .launch_builder(&be.barycentric_ext3_strided_multi)
+            .arg(aux_handle.buf.as_ref())
+            .arg(&col_stride_u64)
+            .arg(&row_stride_u64)
+            .arg(&points_view)
+            .arg(&inv_view)
+            .arg(&n_u64)
+            .arg(&k_u64)
+            .arg(&chunks_u64)
+            .arg(&mut partials)
+            .launch(cfg)?;
+    }
+    let combine_cfg = LaunchConfig {
+        grid_dim: (total.div_ceil(BLOCK_DIM as usize) as u32, 1, 1),
+        block_dim: (BLOCK_DIM, 1, 1),
+        shared_mem_bytes: 0,
+    };
+    unsafe {
+        stream
+            .launch_builder(&be.barycentric_combine_partials)
+            .arg(&partials)
+            .arg(&chunks_u64)
+            .arg(&total_u64)
+            .arg(&mut out_dev)
+            .launch(combine_cfg)?;
+    }
+    let out = stream.clone_dtoh(&out_dev)?;
+    stream.synchronize()?;
+    Ok(out)
+}
+
 /// Gather full rows from a device-resident base-field LDE handle. `rows` are LDE
 /// row indices; returns their column values row-major (`rows.len() * main.m`
 /// u64, `out[q*num_cols + col]`) — i.e. the concatenation of

crypto/math-cuda/src/device.rs

@@ -207,6 +207,9 @@ pub struct Backend {
     pub barycentric_ext3_batched: CudaFunction,
     pub barycentric_base_batched_strided: CudaFunction,
     pub barycentric_ext3_batched_strided: CudaFunction,
+    pub barycentric_base_strided_multi: CudaFunction,
+    pub barycentric_ext3_strided_multi: CudaFunction,
+    pub barycentric_combine_partials: CudaFunction,
     pub gather_rows_base: CudaFunction,
     pub gather_rows_ext3: CudaFunction,
 
@@ -438,6 +441,9 @@ impl Backend {
                 .load_function("barycentric_base_batched_strided")?,
             barycentric_ext3_batched_strided: bary
                 .load_function("barycentric_ext3_batched_strided")?,
+            barycentric_base_strided_multi: bary.load_function("barycentric_base_strided_multi")?,
+            barycentric_ext3_strided_multi: bary.load_function("barycentric_ext3_strided_multi")?,
+            barycentric_combine_partials: bary.load_function("barycentric_combine_partials")?,
             gather_rows_base: bary.load_function("gather_rows_base")?,
             gather_rows_ext3: bary.load_function("gather_rows_ext3")?,
             deep_composition_ext3_row: deep.load_function("deep_composition_ext3_row")?,