experimental: cuda oxide

CUDA_BENCHMARK.md

@@ -0,0 +1,91 @@
+# fGN sampling: CPU vs CUDA backends
+
+Head-to-head throughput of the fractional-Gaussian-noise (fGN) sampling backends,
+generated by `benches/fgn_cuda_compare.rs` on an **NVIDIA RTX 4070 SUPER** (Ada,
+sm_89, 12 GB).
+
+## Setup
+
+- **cpu** — `Cpu` backend: SIMD + rayon, `ndrustfft` circulant-embedding FFT.
+- **cuda_native** — `CudaNative`: cudarc + cuFFT + NVRTC Philox.
+- **gpu_cuda** — `CubeCl`: cubecl, hand-written shared-memory radix FFT.
+- Harness: criterion, `f32`. **Throughput unit: `Melem/s` = 10⁶ generated fGN
+  values per second** (`Gelem/s` = 10⁹). Higher is faster. `n` = path length,
+  `m` = number of paths; a batch produces `n·m` values.
+- Reported value = criterion's median (middle of the [lo, mid, hi] estimate).
+
+```bash
+cargo bench --bench fgn_cuda_compare --features "cuda-native,gpu-cuda"
+```
+
+## Single path (`sample`, m = 1) — Melem/s
+
+| n | cpu | cuda_native | gpu_cuda |
+|--:|--:|--:|--:|
+| 256 | **129.9** | 5.6 | 2.3 |
+| 512 | **129.4** | 9.1 | 3.7 |
+| 1024 | **129.6** | 12.8 | 6.9 |
+| 2048 | **122.4** | 25.4 | 11.5 |
+| 4096 | **119.7** | 49.0 | 19.3 |
+| 8192 | **116.5** | 94.5 | 27.5 |
+| 16384 | 113.9 | **149.1** | 31.9 |
+| 32768 | 108.2 | **312.3** | 39.7 |
+| 65536 | 77.4 | **346.9** | 40.4 |
+| 131072 | 74.7 | **218.8** | 23.2 |
+| 262144 | 69.6 | **148.0** | 25.2 |
+
+`cuda_native` (cuFFT) overtakes the CPU at **n ≈ 16384**; it peaks at n = 65536
+(~**4.5× CPU**) and stays ahead to n = 262144 (~2.1× CPU). `gpu_cuda` is slowest
+throughout.
+
+## Batch (`sample_par`, n × m) — Melem/s
+
+| n | m | cpu | cuda_native | gpu_cuda |
+|--:|--:|--:|--:|--:|
+| 1024 | 1024 | **1726** | 627 | 47 |
+| 1024 | 16384 | **1751** | 557 | 49 |
+| 1024 | 131072 | **1198** | 471 | 48 |
+| 4096 | 16384 | 269 | **581** | 48 |
+| 4096 | 65536 | 118 | **383** | 42 |
+| 16384 | 4096 | **228** | 225 | 30 |
+| 16384 | 16384 | 59 | **183** | 32 |
+| 65536 | 1024 | **632** | 301 | 22 |
+| 65536 | 4096 | **504** | 236 | 28 |
+| 131072 | 512 | **382** | 288 | 22 |
+| 131072 | 2048 | **369** | 221 | 25 |
+| 262144 | 256 | **278** | 260 | 21 |
+| 262144 | 1024 | 279 | **288** | 28 |
+
+## Regimes (who wins, and when)
+
+- **Short paths (n ≤ 1024), any m:** CPU dominates (1.2–1.75 **Gelem/s**) — the
+  FFTs fit cache and rayon saturates the cores; GPU launch/transfer overhead is
+  pure loss.
+- **Many *and* long paths (n ≥ 4096 with large m):** **cuda_native wins ~3×** —
+  e.g. 4096×16384 (581 vs 269), 4096×65536 (383 vs 118), 16384×16384 (183 vs 59).
+  This is the batched-cuFFT sweet spot.
+- **Very long paths, modest m (n ≥ 65536):** CPU and cuda_native are close;
+  the per-call host↔device copy + plan setup erodes the GPU's lead, so the winner
+  flips case by case (e.g. CPU wins 65536×1024; they tie at 262144×1024).
+- **gpu_cuda (cubecl)** never wins — its hand-written FFT runs at a near-constant
+  ~20–49 Melem/s regardless of size, ~3–30× behind the others. Replacing it with
+  a cuFFT-class kernel is the obvious optimisation.
+
+## Takeaways
+
+- Use the **CPU** for short paths and for latency-sensitive single draws below
+  n≈16k.
+- Use **cuda_native (cuFFT)** for large single paths (n ≥ 16k) and for
+  generating **many long paths at once** (the 3× regime above).
+- **gpu_cuda** is not competitive today; it is kept for the cross-platform
+  (wgpu/Metal) path, not for raw CUDA speed.
+
+## cuda-oxide backend
+
+The experimental `cuda-oxide` backend is **not in this benchmark**: its kernels
+use libdevice transcendentals (`ln`/`sqrt`/`cos`/`sin`), and `nvvmCompileProgram`
+rejects the codegen's opaque-pointer NVVM IR (`parse expected type`). This
+reproduces with cuda-oxide v0.1.0's **own** `math_atan` example inside its **exact
+pinned nix environment** (CUDA 13.2 + LLVM 22 + nightly-2026-04-03), so it is an
+upstream limitation. The backend is maintained on a separate experimental branch
+until that is fixed.

Cargo.toml

@@ -10,6 +10,10 @@ members = [
     "stochastic-rs-viz",
     "stochastic-rs-py",
 ]
+# Device-only cuda-oxide kernel crate: its own `[workspace]` root (built
+# standalone with the cuda-oxide codegen backend to emit NVVM IR), excluded so
+# the umbrella workspace does not claim it as a member.
+exclude = ["stochastic-rs-stochastic/fgn-oxide-kernels"]
 
 [workspace.package]
 version = "2.3.0"
@@ -61,7 +65,6 @@ owens-t = "0.1.5"
 
 # General utility
 anyhow = "1.0.89"
-either = "1.15.0"
 parking_lot = "0.12.5"
 ordered-float = "5.0.0"
 chrono = "0.4.38"
@@ -91,6 +94,9 @@ cubecl = { version = "0.9.0", default-features = false }
 cubecl-cuda = "0.9.0"
 cubecl-wgpu = "0.9.0"
 cudarc = { version = "0.19.2", features = ["cuda-12080", "cuda-version-from-build-system"] }
+cuda-core = { git = "https://github.com/NVlabs/cuda-oxide.git", rev = "4a56e4220aab8ce5d085a411e7f806cebb647d14" }
+cuda-device = { git = "https://github.com/NVlabs/cuda-oxide.git", rev = "4a56e4220aab8ce5d085a411e7f806cebb647d14" }
+cuda-host = { git = "https://github.com/NVlabs/cuda-oxide.git", rev = "4a56e4220aab8ce5d085a411e7f806cebb647d14" }
 gpu-fft = "1.1.1"
 metal = "0.33.0"
 
@@ -141,7 +147,9 @@ cubecl = { workspace = true, optional = true, default-features = false }
 cubecl-cuda = { workspace = true, optional = true }
 cubecl-wgpu = { workspace = true, optional = true }
 cudarc = { workspace = true, optional = true, features = [ "cuda-12080", "cuda-version-from-build-system", ] }
-either = { workspace = true }
+cuda-core = { workspace = true, optional = true }
+cuda-device = { workspace = true, optional = true }
+cuda-host = { workspace = true, optional = true }
 flate2 = { workspace = true }
 gpu-fft = { workspace = true, optional = true }
 gauss-quad = { workspace = true }
@@ -223,6 +231,16 @@ name = "fgn_cuda_native"
 harness = false
 required-features = ["cuda-native"]
 
+[[bench]]
+name = "fgn_cuda_oxide"
+harness = false
+required-features = ["cuda-oxide-experimental"]
+
+[[bench]]
+name = "fgn_cuda_compare"
+harness = false
+required-features = ["cuda-native", "gpu-cuda"]
+
 [[bench]]
 name = "fgn_all_backends"
 harness = false
@@ -318,9 +336,19 @@ name = "dual_stream_compare"
 harness = false
 required-features = ["dual-stream-rng"]
 
+[[example]]
+name = "cuda_fgn_analysis"
+required-features = ["cuda-native"]
+
 [features]
 ai = ["dep:stochastic-rs-ai", "stochastic-rs-ai/quant"]
 cuda-native = ["dep:cudarc", "cudarc/cufft", "stochastic-rs-stochastic/cuda-native"]
+cuda-oxide-experimental = [
+    "dep:cuda-core",
+    "dep:cuda-device",
+    "dep:cuda-host",
+    "stochastic-rs-stochastic/cuda-oxide-experimental",
+]
 default = []
 # Experimental: opt in to the dual-stream RNG (`SimdRngDual` /
 # `SimdNormalDual`). ~5–11% speedup on Ziggurat-based Normal/Exp bulk

README.md

@@ -149,7 +149,7 @@ cargo bench --features cuda-native --bench fgn_cuda_native
 
 Single path:
 
-| n      | CPU `sample` | CUDA `sample_cuda_native(1)` | Speedup    |
+| n      | CPU `sample` | CUDA `.on(Device::CudaNative).sample()` | Speedup    |
 |-------:|-------------:|------------------------------:|-----------:|
 |  1,024 |       8.1 µs |                          46 µs|       0.18× |
 |  4,096 |        35 µs |                          84 µs|       0.42× |
@@ -158,7 +158,7 @@ Single path:
 
 Batch:
 
-| n, m         | CPU `sample_par` | CUDA `sample_cuda_native` | Speedup  |
+| n, m         | CPU `sample_par` | CUDA `.on(Device::CudaNative).sample_par` | Speedup  |
 |--------------|------------------:|---------------------------:|---------:|
 |   4,096, 32  |          147 µs   |                     117 µs |  **1.3×** |
 |   4,096, 512 |         1.78 ms   |                    2.37 ms |   0.75×   |

benches/fgn_accelerate.rs

@@ -6,6 +6,7 @@ use criterion::Criterion;
 use criterion::criterion_group;
 use criterion::criterion_main;
 use stochastic_rs::simd_rng::Unseeded;
+use stochastic_rs::stochastic::device::Accelerate;
 use stochastic_rs::stochastic::noise::fgn::Fgn;
 use stochastic_rs::traits::ProcessExt;
 
@@ -17,13 +18,14 @@ fn bench_single(c: &mut Criterion) {
 
   for &n in &[1024usize, 4096, 16384, 65536] {
     let fgn = Fgn::new(0.7f32, n, None, Unseeded);
-    let _ = fgn.sample_accelerate(1);
+    let dev = Fgn::new(0.7f32, n, None, Unseeded).on::<Accelerate>();
+    let _ = dev.sample();
 
     g.bench_with_input(BenchmarkId::new("cpu", n), &n, |b, _| {
       b.iter(|| black_box(fgn.sample()));
     });
     g.bench_with_input(BenchmarkId::new("accelerate", n), &n, |b, _| {
-      b.iter(|| black_box(fgn.sample_accelerate(1).unwrap()));
+      b.iter(|| black_box(dev.sample()));
     });
   }
   g.finish();
@@ -44,7 +46,8 @@ fn bench_batch(c: &mut Criterion) {
   ];
   for &(n, m) in &cases {
     let fgn = Fgn::new(0.7f32, n, None, Unseeded);
-    let _ = fgn.sample_accelerate(m);
+    let dev = Fgn::new(0.7f32, n, None, Unseeded).on::<Accelerate>();
+    let _ = dev.sample_par(m);
     let label = format!("n={n},m={m}");
 
     g.bench_with_input(BenchmarkId::new("cpu", &label), &(n, m), |b, &(_, m)| {
@@ -54,7 +57,7 @@ fn bench_batch(c: &mut Criterion) {
       BenchmarkId::new("accelerate", &label),
       &(n, m),
       |b, &(_, m)| {
-        b.iter(|| black_box(fgn.sample_accelerate(m).unwrap()));
+        b.iter(|| black_box(dev.sample_par(m)));
       },
     );
   }

benches/fgn_all_backends.rs

@@ -6,6 +6,9 @@ use criterion::Criterion;
 use criterion::criterion_group;
 use criterion::criterion_main;
 use stochastic_rs::simd_rng::Unseeded;
+use stochastic_rs::stochastic::device::Accelerate;
+use stochastic_rs::stochastic::device::CubeCl;
+use stochastic_rs::stochastic::device::MetalNative;
 use stochastic_rs::stochastic::noise::fgn::Fgn;
 use stochastic_rs::traits::ProcessExt;
 
@@ -17,23 +20,26 @@ fn bench_single(c: &mut Criterion) {
 
   for &n in &[1024usize, 4096, 16384, 65536] {
     let fgn = Fgn::new(0.7f32, n, None, Unseeded);
+    let dev_gpu = Fgn::new(0.7f32, n, None, Unseeded).on::<CubeCl>();
+    let dev_metal = Fgn::new(0.7f32, n, None, Unseeded).on::<MetalNative>();
+    let dev_accel = Fgn::new(0.7f32, n, None, Unseeded).on::<Accelerate>();
 
     // warmup
-    let _ = fgn.sample_gpu(1);
-    let _ = fgn.sample_metal(1);
-    let _ = fgn.sample_accelerate(1);
+    let _ = dev_gpu.sample();
+    let _ = dev_metal.sample();
+    let _ = dev_accel.sample();
 
     g.bench_with_input(BenchmarkId::new("cpu", n), &n, |b, _| {
       b.iter(|| black_box(fgn.sample()));
     });
     g.bench_with_input(BenchmarkId::new("gpu_cubecl", n), &n, |b, _| {
-      b.iter(|| black_box(fgn.sample_gpu(1).unwrap()));
+      b.iter(|| black_box(dev_gpu.sample()));
     });
     g.bench_with_input(BenchmarkId::new("metal", n), &n, |b, _| {
-      b.iter(|| black_box(fgn.sample_metal(1).unwrap()));
+      b.iter(|| black_box(dev_metal.sample()));
     });
     g.bench_with_input(BenchmarkId::new("accelerate", n), &n, |b, _| {
-      b.iter(|| black_box(fgn.sample_accelerate(1).unwrap()));
+      b.iter(|| black_box(dev_accel.sample()));
     });
   }
   g.finish();
@@ -54,22 +60,25 @@ fn bench_batch(c: &mut Criterion) {
   ];
   for &(n, m) in &cases {
     let fgn = Fgn::new(0.7f32, n, None, Unseeded);
-    let _ = fgn.sample_gpu(m);
-    let _ = fgn.sample_metal(m);
-    let _ = fgn.sample_accelerate(m);
+    let dev_gpu = Fgn::new(0.7f32, n, None, Unseeded).on::<CubeCl>();
+    let dev_metal = Fgn::new(0.7f32, n, None, Unseeded).on::<MetalNative>();
+    let dev_accel = Fgn::new(0.7f32, n, None, Unseeded).on::<Accelerate>();
+    let _ = dev_gpu.sample_par(m);
+    let _ = dev_metal.sample_par(m);
+    let _ = dev_accel.sample_par(m);
     let label = format!("n={n},m={m}");
 
     g.bench_with_input(BenchmarkId::new("cpu", &label), &m, |b, &m| {
       b.iter(|| black_box(fgn.sample_par(m)));
     });
     g.bench_with_input(BenchmarkId::new("gpu_cubecl", &label), &m, |b, &m| {
-      b.iter(|| black_box(fgn.sample_gpu(m).unwrap()));
+      b.iter(|| black_box(dev_gpu.sample_par(m)));
     });
     g.bench_with_input(BenchmarkId::new("metal", &label), &m, |b, &m| {
-      b.iter(|| black_box(fgn.sample_metal(m).unwrap()));
+      b.iter(|| black_box(dev_metal.sample_par(m)));
     });
     g.bench_with_input(BenchmarkId::new("accelerate", &label), &m, |b, &m| {
-      b.iter(|| black_box(fgn.sample_accelerate(m).unwrap()));
+      b.iter(|| black_box(dev_accel.sample_par(m)));
     });
   }
   g.finish();