GPU Backend

cargo build --release --features "parallel gpu"
cargo nextest run --features "parallel gpu" --test golden_gpu --test golden_gpu_density_matrix

Constructing a GpuContext compiles the CUDA kernels through NVRTC, one to three seconds on a GTX 1080 Ti. The PTX is shared by every context in the process and cached on disk in prism-q-ptx under the user cache directory (XDG_CACHE_HOME, else LOCALAPPDATA, else HOME/.cache, else the OS temp directory), keyed by device arch, crate version, and a hash of the kernel source, so later processes skip the compile. A missing, unreadable, or corrupt cache file only costs a recompile; delete the directory to force one.

CUDA acceleration covers statevector execution, stabilizer execution, density-matrix execution, and compiled BTS sampling. Seven entry points are available:

• BackendKind::AutoGpu { context } (simulate(circuit).gpu_auto(ctx)). Automatic backend selection with the device opted in. The shape-based decision tree runs unchanged; a selected statevector or stabilizer workload that clears the family's qubit crossover and fits in VRAM runs on the device. Everything else, including a device allocation that fails at init, takes the identical CPU path (the soft VRAM fallback).
• BackendKind::StatevectorGpu { context }. Public dispatch path for statevector GPU execution. It routes through simulate(circuit).backend(kind).seed(seed).run(), keeps fusion and subsystem decomposition, and uses crate::gpu::min_qubits() (default 14, PRISM_GPU_MIN_QUBITS override) to keep small sub-circuits on CPU.
• BackendKind::StabilizerGpu { context }. Public dispatch path for stabilizer GPU execution. Gate application uses a device tableau and one word-grouped batched Clifford kernel (stab_apply_word_grouped). Measurement and reset keep pivot search, row cascade, phase fixup, and deterministic outcomes on the device. The default crossover stays conservative (STABILIZER_MIN_QUBITS_DEFAULT = 100_000, PRISM_STABILIZER_GPU_MIN_QUBITS override) until benchmarks justify lowering it. Direct backend benchmarks should use StabilizerBackend::with_gpu(ctx) to exclude diagnostic readbacks from probabilities(), export_tableau(), and export_statevector(). Golden tests cover every kernel path, including 500q GHZ measure-all.
• BackendKind::DensityMatrixGpu { context }. The exact mixture held in device memory. Explicit only: neither Auto nor AutoGpu selects it, there is no crossover, and there is no host fallback. init budgets the 4^n buffer against the free VRAM and errors before allocating when it does not fit, so an 11 GiB card holds 13 qubits (1 GiB at 13, 4 GiB at 14 plus scratch). The unitary half reuses the dense statevector kernels on the embedded 2n-qubit buffer, and every channel, measurement, and readout sweep runs as a kernel of its own. The noisy Simulate terminals answer from the device mixture exactly as DensityMatrix does. DensityMatrixBackend::new(seed).with_gpu(ctx) is the direct form.
• StatevectorBackend::new(seed).with_gpu(ctx). Direct statevector GPU opt-in. Every instruction routes to CUDA after the context is attached. No crossover or subsystem decomposition applies.
• StabilizerBackend::new(seed).with_gpu(ctx). Direct stabilizer GPU opt-in for kernel benchmarks and targeted correctness tests.
• run_shots_compiled_with_gpu (or CompiledSampler::with_gpu(ctx)). GPU BTS sampling for flat sparse parity. The path launches one kernel per 65_536-shot chunk, uses random bits generated on the host, and preserves the CPU sample_bts_meas_major layout. The sampler caches sparse parity CSR arrays, packed reference bits, and reusable scratch on the device. It is active only when num_shots >= BTS_MIN_SHOTS_DEFAULT (131_072 by default, PRISM_GPU_BTS_MIN_SHOTS override). sample_bulk_packed_device returns a DevicePackedShots handle. Marginals reduce to one counter per measurement row on the device. Exact counts use a bounded device hash reduction for up to 8 packed measurement words when the compact result is cheaper to transfer than the full shot matrix. Otherwise the API uses a host copy for correctness.

When a GPU context is attached, Backend::init allocates state on the device instead of a host Vec<Complex64> and every instruction routes to a CUDA kernel. On the hard statevector path (StatevectorGpu, with_gpu), a state that does not fit the currently free VRAM is rejected at init with an error naming the requested and free device memory; GpuContext::max_qubits_for_statevector reports the advisory cap from free memory.

The four BackendKind entry points are also reachable from Python, from a build carrying the gpu feature. See Python Bindings.

Module layout (src/gpu/)

Kernel coverage

Every variant in the Gate enum has a dedicated kernel. Batched variants (BatchPhase, BatchRzz, DiagonalBatch, MultiFused { all_diagonal: true }) use LUT kernels that consume the same host table builders as the CPU path. Non-diagonal MultiFused uses a shared memory tiled kernel (apply_multi_fused_tiled, TILE_Q = 10, TILE_SIZE = 1024) over a chosen set of ten qubits per pass: the five lowest qubits, which keep a warp's loads contiguous, plus up to five of the sub-gates' higher targets. A MultiFused over n qubits therefore takes about (n - 5) / 5 passes, each applying its sub-gates in shared memory. A pass with fewer than three sub-gates falls back to per gate launches. Multi2q still launches once per sub-gate; rare in practice.

PTX template substitution: the CUDA C source lives in .cu files beside the Rust launchers and reaches KERNEL_SOURCE_TEMPLATE through include_str!. kernels/dense.cu carries placeholders such as {{BP_TABLE_SIZE}} and {{TILE_Q}}. The kernel_source() function substitutes them at device construction from the Rust constants in src/backend/statevector/kernels.rs, keeping CPU and GPU in sync.

Correctness

tests/golden_gpu.rs compares GPU amplitudes against the CPU statevector within 1e-12 for every gate variant, the fusion paths, and the BackendKind::StatevectorGpu public dispatch path at the crossover boundary; tests/golden_gpu_density_matrix.rs does the same for the full density matrix buffer. Both suites skip when no device opens, and no CI job opens one. scripts/test-gpu.ps1 runs them with PRISM_REQUIRE_GPU=1 so a missing or unusable device fails the run.

Shot reproducibility

Two limits bound what a seed guarantees, and neither is visible from the golden equality tests.

• CPU against GPU: agreement in distribution, not bit for bit. Both paths draw the same RNG stream for a given shot seed, but the device reduction that produces a measurement probability sums in tree order with FMA contraction, so it can differ from the host sum in the last ulp, and a uniform draw landing between the two flips that outcome and every outcome after it. Where the probability is a dyadic rational (0.5, 1.0, and the amplitudes reachable from Clifford gates) both sums are exact and the shots do match, which is what statevector_gpu_mid_measure_shots_match_cpu pins. statevector_gpu_shot_frequencies_match_cpu_off_dyadic pins the general case: equal frequencies within 5 sigma after Rx(0.3).
• GPU BTS sampling: reproducible at a fixed Rayon thread count. Above MIN_PAR_DRAWS random bits per chunk, fill_random_bits seeds one stream per worker and partitions the draws by rayon::current_num_threads(), so the same seed on a host with a different worker count produces different shots. Below that threshold the serial single-stream path runs and the seed reproduces outright. Pin RAYON_NUM_THREADS when byte-identical shot payloads matter across machines.

Current limits

• Device placement is silent. Circuits below the crossover run on the host, and the AutoGpu soft VRAM fallback degrades to host execution without a report; nothing user-visible says whether a run executed on the device.
• Stabilizer probabilities(), export_tableau(), and export_statevector() read back to the CPU.
• Every trajectory shot rebuilds the backend, reallocating the device buffer (measured at 0.1 ms per shot at 20 qubits, so not a practical cost). Custom Kraus branch probabilities come from an on-device reduced-density-matrix reduction (rdm_qubit), not a full-state readback.
• Kernel design and crossover analysis live in the module docstrings on src/gpu/kernels/dense.rs.