GPU support

Every process carries its sampling backend as a type parameter and is re-typed with .on::<B2>() (see backends). This page is the inventory: which processes have a device kernel today, on which backends, and what the rest do.

The support matrix

Feature flags: metal (MetalNative), cuda-native (CudaNative), cubecl-cuda / cubecl-wgpu (CubeCl on CUDA or on Metal / Vulkan / WebGPU through wgpu), accelerate. A marker only exists when its feature is compiled, so selecting a backend a process does not support — or one the build does not carry — is a compile error, never a silent fallback. The same holds for precision: the capability traits take the scalar (FgnBackend<T>, EulerBackend<T>) and a device implements them only for the precision its kernels compute in, so Gbm<f64> or Fgn<f64> on MetalNative / CubeCl does not compile. Nothing is computed in f32 behind an f64 type.

What the device path does

• fGN family. The circulant-embedding FFT that produces fractional Gaussian noise runs on the device; the process's own drift / diffusion recursion then runs on the host as before. One launch per sample_par batch.
• Euler engine. The whole Euler–Maruyama recursion of a path runs in one kernel thread; sample_par(m) is one launch for all m paths, sample() launches a single path. CIR uses full truncation (Lord, Koekkoek & van Dijk 2010) in the kernel. The CUDA and Metal kernels are one C body rendered per language, and the CubeCL kernel repeats its integer hash, so MetalNative, CudaNative and CubeCl agree with each other seed for seed up to f32 libm rounding.
• Host backends. Cpu is the process's own sampler; Accelerate is a CPU device that only differs for the fGN family (vDSP FFT). For every other process both are the same code path.

Large batches

A device handle holds at most batch_budget bytes of path data per launch (1 GiB by default, STOCHASTIC_RS_DEVICE_BATCH_BYTES, or CudaNative::default().with_batch_budget(n)). A larger sample_par(m) runs as several launches whose union is bit-identical to one launch, because the kernels hash the global path index; sample_map(m, f) maps each chunk in parallel on the host before the next chunk is launched, so the batch never has to fit in memory at once. On native CUDA the chunks alternate between two streams: one copies back through pinned memory while the next computes. The fGN devices chunk the same way, one seed per batch and an element offset per chunk, and keep the FFT plans and buffers of the last four (n, m) shapes, so alternating two sizes does not re-plan on every call.

More than one GPU

The handle names the device: CudaNative::new(1) (CUDA ordinal), MetalNative::new(1) (index into Device::all()), CubeCl::new(1) (the runtime's device index); .on_device(handle) puts a process on it, and handle.probe() reports the device opened. The defaults read STOCHASTIC_RS_DEVICE. From Python the name carries the ordinal: device="cuda-native:1", probe_device("metal:0").

When the device is missing

handle.probe() on any backend handle returns Ok(DeviceInfo) with the device's name and precisions, or Err(DeviceError) saying why it cannot be used (Unavailable, Compile, Launch). Every process has try_sample() and try_sample_par(m) on ProcessExt, returning the same error instead of panicking (always Ok on the host devices, bit-identical to sample*); the plain sample* calls panic with the error's message and never fall back to the host. From Python a device failure is a RuntimeError.

Precision and reproducibility

The device seed is drawn from the process's own seed source: a Deterministic process reproduces its device paths call after call, an Unseeded one draws fresh entropy — exactly as on the host. Host and device paths agree in distribution, not bit for bit, because the device draws its own stream.

From Python

The device-capable classes take device= (PyGbm(..., dtype="f32", device="metal"); the fGN family and the Euler engine's three processes), checked and probed at construction. PyGbm(0.05, 0.2, 253, x0=100.0, t=1.0, seed=7, device="metal", dtype="f32").sample_par(m) runs the Euler engine with device one of "cpu", "gpu" (the first compiled device backend), "cuda-native", "metal", "cubecl"; a backend the build does not carry raises ValueError with a rebuild hint. The array is float64 from the CPU and native CUDA paths and float32 from Metal and CubeCL, whose kernels compute in single precision. The published wheels are CPU-only; the device backends come with a source build — maturin develop --features metal,cubecl-wgpu on a Mac, --features cuda-native on an NVIDIA machine.

Trying the CUDA back-end without a CUDA machine

notebooks/colab_cuda_check.ipynb opens in Google Colab: on a T4 runtime it installs Rust, clones the repository, runs the native CUDA tests (the fGN pipeline and the Euler engine), builds the Python wheel with device="cuda-native" and plots what the GPU produced. The Rust tests build in 10–20 minutes on Colab's two cores, the wheel in another ten; the tests and the plots take seconds.

A run on a T4

notebooks/runs/colab_cuda_check_T4_2026-09-03.ipynb is the executed copy of a run on 2026-09-03 (Tesla T4, driver 580.82 / CUDA 13.0, toolkit 12.8; GitHub's own viewer may not render Colab outputs, nbviewer does). All 12 cuda_native tests passed, including the two that assert a chunked batch equals one launch path for path, in f32 and f64; the wheel built for CPython 3.13 and its device tests passed; probe_device reported {'backend': 'CudaNative', 'name': 'Tesla T4', 'precisions': ['f32', 'f64']}.

The variance of fractional Brownian motion sampled on the GPU sits on the t^{2H} lines for H = 0.3, 0.5, 0.8 exactly where the CPU samples sit, and the Euler engine's GBM terminal mean over 100k paths was 105.15 against the 105.13 forward. Wall time per call, warm, against Colab's two CPU cores:

The Euler case is slower on the T4 in f64 because that card executes double precision at 1/32 of its single-precision rate, and the 200 MB result is copied back and re-laid as rows on the host; dtype="f32" is the setting for the Euler engine on inference-class GPUs. The fGN pipeline is FFT-bound and wins already at these sizes; the gap widens with m.

Adding a kernel

A process joins the Euler engine by implementing EulerCoefficients (its EulerSpec family, initial value, grid, horizon, device seed) and switching its backend_switch! line from via host to via euler. A new drift / diffusion family is one EulerSpec variant plus one family branch in the CubeCL kernel (euler/gpu.rs) and one in the C body the CUDA and Metal kernels render from a single text (euler/kernel.rs). Widening a process's bound from HostBackend to EulerBackend is not a breaking change for its callers.

Not on a device yet, and what it would take: the two-state models (Heston, SABR, Bates, Heston-SLV — a correlated pair and a two-component stepper), time-dependent coefficients (CIR++ φ(t), Hull–White θ(t) — a curve buffer on the device), jump processes (a Poisson draw in the kernel), and the Volterra / Markov-lift family (a separate O(n·N') kernel).