Introduction

qsv is a high-performance quantum statevector simulator written in Rust, built as a study in performance engineering. Off-the-shelf simulators already exist (qsim, Qiskit-Aer, QuEST, Yao.jl); qsv's purpose is not to be another one, but to demonstrate how one is optimized — cache-, SIMD-, and threading-aware design driven by profiling and benchmarked honestly against the established tools.

The thesis that organizes everything

Statevector simulation is memory-bandwidth-bound, not compute-bound.

Applying a 1-qubit gate streams the entire $2^N$-amplitude array while doing only ~2 complex multiplies per 16-byte amplitude — an arithmetic intensity of roughly 0.13 FLOP/byte, deep in the bandwidth-bound region of the roofline. Every optimization decision in qsv is justified by one question:

Does this reduce bytes moved per gate, or raise arithmetic intensity per byte moved?

This reframes "squeeze all the performance" as primarily a memory-traffic and cache problem. See How we optimize for the full argument and the evidence behind it.

What's here

• A tutorial for using qsv as a library.
• An architecture overview of the crate and its pluggable backends.
• How we optimize — the optimization strategy and milestone narrative (the centerpiece).
• The core kernel — the insert_zero_bit indexing trick every production simulator shares, explained from the ground up.
• Benchmarking & profiling — how we measure.
• Research notes — distilled findings from qsim, Qiskit-Aer, QuEST, Yao.jl, cuStateVec, spinoza, and Algorithms for Modern Hardware.

Status

Foundations and the first optimization milestones are in place: a Structure-of-Arrays statevector, the universal bit-shift kernel, a pluggable Backend trait with three implementations (a naive oracle plus two optimized backends), and a differential test suite that validates every kernel against the oracle. See the roadmap for the milestone-by-milestone plan.

Installation

Prerequisites

• Rust 1.86+ (stable). The project pins rust-version = "1.86".
• A C toolchain is not required — qsv-core is pure Rust with no external dependencies.

Check your toolchain:

rustc --version   # 1.86.0 or newer

Build & test

git clone <repo-url> statevec-sim
cd statevec-sim

cargo build --workspace            # debug build
cargo build --release --workspace  # optimized build
cargo test  --workspace            # unit + oracle + differential tests

A clean run should report all tests passing, with cargo clippy --workspace --all-targets -- -D warnings and cargo fmt --check both clean (these gate CI).

Run the CLI

The qsv binary is a smoke runner that prepares a GHZ state and prints its outcome probabilities:

cargo run --release --bin qsv -- 24   # GHZ on 24 qubits (256 MB statevector)

Optional tooling

These power the documentation and the benchmarking/profiling workflow:

# Documentation (this book). mdBook 0.5+ needs rustc 1.88; on 1.86 pin the 0.4 line:
cargo install mdbook --version "^0.4" --locked
mdbook serve docs        # live-reload at http://localhost:3000
mdbook build docs        # render to docs/book/

# Profiling (macOS / Apple Silicon)
cargo install samply --locked          # sampling profiler, Firefox-profiler UI
cargo install cargo-instruments        # Xcode Instruments integration

# Benchmarks use criterion (pulled automatically by `cargo bench`).

Memory ceiling

A statevector of $N$ qubits in f64 needs $2^N\times 16$ bytes (real + imaginary). On a 36 GB machine the practical in-place ceiling is ~30 qubits (f64, 16 GB) or ~31 (f32). qsv updates the state in place precisely because at that size there is no room for a second buffer.

Tutorial

This walks through using qsv-core as a library: building circuits, running them on a backend, and reading out results.

A first circuit: the Bell state

#![allow(unused)]
fn main() {
use qsv_core::prelude::*;

// 2-qubit register, then H(0); CX(0,1)  ->  (|00> + |11>)/√2
let mut circuit = Circuit::<f64>::new(2);
circuit.h(0).cx(0, 1);

let state = BitShiftBackend.execute(&circuit);

let inv_sqrt2 = std::f64::consts::FRAC_1_SQRT_2;
assert!((state.amplitude(0b00).re - inv_sqrt2).abs() < 1e-12);
assert!((state.amplitude(0b11).re - inv_sqrt2).abs() < 1e-12);
}

The builder is fluent and chainable. Qubit 0 is the least significant bit of the basis index, so |00⟩ is index 0 and |11⟩ is index 3.

Gates

Single-qubit: h, x, y, z, s, t, sx, and parametric rx(θ), ry(θ), rz(θ), phase(λ). Two-qubit (control first): cx(c,t), cz(a,b), cphase(c,t,λ), swap(a,b), rzz(a,b,θ).

#![allow(unused)]
fn main() {
let mut c = Circuit::<f64>::new(3);
c.h(0)
 .rx(1, std::f64::consts::FRAC_PI_2)
 .cx(0, 1)
 .cz(1, 2)
 .rzz(0, 2, 0.7);   // QAOA-style two-qubit rotation
}

Arbitrary unitaries can be pushed directly as a dense matrix on an ordered qubit list with circuit.push(dense_gate, &qubits) — this is how a Toffoli or a fused gate is applied.

Reading out results

#![allow(unused)]
fn main() {
let backend = BitShiftBackend;
let state = backend.execute(&circuit);

let amp = state.amplitude(5);          // Cplx<f64> at basis index 5
let p   = state.probabilities();       // Vec<f64> of |ψ_i|²
let norm = state.norm_sqr();           // ⟨ψ|ψ⟩, == 1 for a valid state
}

Built-in circuit generators

The circuits module provides standard circuits (also used by the tests and benchmarks):

#![allow(unused)]
fn main() {
use qsv_core::circuits::{ghz, qft, random_circuit};

let g = ghz(5);                        // (|0…0> + |1…1>)/√2
let f = qft(8);                        // Quantum Fourier Transform
let r = random_circuit(10, 200, 42);   // 200 random gates, reproducible (seed 42)

// QFT of |0…0> is the uniform superposition:
let s = BitShiftBackend.execute(&qft(4));
let expected = 1.0 / ((1usize << 4) as f64);          // |amplitude|² for every state
assert!((s.amplitude(0).norm_sqr() - expected).abs() < 1e-12);
}

Choosing a backend

qsv exposes several backends, all implementing the same Backend trait — so you can swap them with a one-line change. They exist side by side so each optimization milestone stays runnable and benchmarkable, and so every fast kernel can be validated against the slow one.

CpuBackend::parallel() (the Default) multithreads above a size threshold; CpuBackend::serial() forces single-threaded. Build with --no-default-features to drop the parallel/simd features (and the rayon/wide dependencies) entirely.

Gate fusion

fuse rewrites a circuit so adjacent gates spanning a few qubits become single composite gates — fewer passes over the state, the dominant optimization on structured circuits. It is a pure Circuit -> Circuit transform, so it composes with any backend:

#![allow(unused)]
fn main() {
use qsv_core::prelude::*;
use qsv_core::circuits::qft;

let circuit = qft(16);
let fused = fuse(&circuit, &FusionConfig::default()); // default max_qubits = 4
let state = CpuBackend::default().execute(&fused);     // same result, fewer passes
}

fuse never changes the computed state (verified against the oracle across max_qubits settings); it only reduces the number of gate applications.

#![allow(unused)]
fn main() {
// Identical circuit, different engine — the seam in action.
let a = RefBackend.execute(&circuit);
let b = BitShiftBackend.execute(&circuit);
// a and b agree amplitude-for-amplitude (this is exactly what the test suite checks).
}

Choosing f32 instead of f64 is a type parameter: Circuit::<f32>::new(n) — half the memory traffic, lower precision.

Precision and generics

Everything is generic over the Real trait (f64 by default, f32 available). Kernels are monomorphized per scalar type — there is no dynamic dispatch in the hot path.

Architecture overview

qsv is a Cargo workspace, which keeps the optimization-critical library isolated (tiny, auditable, dependency-free) from heavier bench/binding crates.

crates/qsv-core   the product: state, gates, circuit, backends, fusion
crates/qsv-cli    `qsv` binary (smoke runner; QASM3 + sampling later)
crates/qsv-bench  criterion benchmarks + profiling binaries
crates/qsv-cuda   optional GPU backend (cudarc + NVRTC; `--features cuda`)
docs/             this mdBook (design, research, tutorial)
_local/           shallow clones of reference simulators (git-ignored)

qsv-core modules

#![deny(unsafe_code)] is set crate-wide; only the hot-path modules opt back in with a localized #[allow(unsafe_code)] and a // SAFETY: justification, keeping the unsafe surface tiny.

Load-bearing abstractions

StateVector<R> — Structure of Arrays

Amplitudes are stored as two separate arrays, re: Vec<R> and im: Vec<R>, not as an interleaved Vec<Complex>. For the bandwidth-bound complex-multiply kernel this is the right layout: a SIMD load of re and a SIMD load of im each yield a register of like-typed values, so the multiply is a straight broadcast-FMA chain with no lane shuffles (no NEON ld2 / x86 unpck de-interleave). All updates are in place — at 30 qubits the vector is 16 GB and there is no room for a second buffer.

Backend<R> — the pluggable seam

Everything above this trait (circuit, gates, fusion) is backend-agnostic; a backend owns the amplitude storage and the gate kernels. This is the seam behind which a future CUDA/cuTile or Metal backend will live without touching the circuit layer. It is kept leak-proof by:

• an associated type State — CPU uses the host StateVector, a GPU backend would use an opaque device handle; no method ever hands out a &mut [R] to host memory;
• reductions (probabilities, and later sample/expectation) are backend methods, so a GPU computes them on-device rather than copying back;
• a single download as the only device→host crossing;
• a default execute(&Circuit) that a GPU overrides to batch a whole circuit.

The trait is validated by having three implementors today — RefBackend, ReshapeBackend, BitShiftBackend — which proves it encodes no CPU-only assumptions.

Gate — zero hot-path allocation

Gate matrices are small and fixed-size; parametric gates (rx(θ), …) materialize on the stack, and only fused/arbitrary unitaries touch the heap (once). v0.0 uses a uniform DenseGate; specialized representations (const-generic Mat2/Mat4, diagonal-only) arrive with the optimized kernels.

Why separate backends per milestone

Each optimization is its own Backend struct rather than a mutation of one. This keeps the slow, obviously-correct oracle alive next to every fast kernel, so the differential test suite can run an identical circuit through both and diff the result. It also lets the benchmark harness run one circuit through every milestone and plot them head to head — the optimization narrative falls out of the architecture.

How we optimize

This is the heart of the project: a disciplined, evidence-driven path from a naive simulator to a fast one, where every step is measurable and justified by a single principle.

The principle: it's memory-bandwidth-bound

Consider applying one 1-qubit gate to an $N$-qubit state. The kernel must touch all $2^N$ amplitudes, reading and writing each once. Per amplitude it does roughly 2 complex multiply-adds (the $2\times 2$ matrix times a 2-vector). In f64, each amplitude is 16 bytes (real + imaginary).

$$\text{arithmetic intensity}\approx \frac{\sim 8\text{FLOP per pair}}{2\times 16\text{bytes per pair}}\approx 0.13\text{FLOP/byte}.$$

A modern CPU sustains tens of FLOP/byte before it runs out of compute; at 0.13 FLOP/byte the kernel is starved for memory bandwidth, not arithmetic. On the roofline it sits far to the left of the ridge point — its ceiling is the slope (bandwidth), not the flat top (peak FLOPs).

This is confirmed independently across every production simulator we studied (qsim, Qiskit-Aer, QuEST, Yao.jl, cuStateVec, spinoza). It is the single most important fact about the problem, and it dictates everything below.

The governing question for any change: does it reduce bytes moved per gate, or raise arithmetic intensity per byte moved? If neither, it will not move the needle, no matter how clever it looks in isolation.

The clinching evidence

The bit-index generation at the kernel's core (insert-zero-bit, gather, scatter) compiles on x86-64 to the BMI2 PEXT/PDEP instructions, replacing a per-bit loop with a single instruction. qsv ships this as the bmi2 feature, and the measurement is decisive (full numbers in bench/results/SUMMARY-xeon.md):

• In isolation, PEXT gather and PDEP insert are ~4× faster than the scalar loop (3.9 vs 0.9 Gelem/s on a Xeon Gold 6526Y).
• End-to-end on a fused QFT-18, the same bmi2 feature changes throughput by 0.995× — i.e. not at all (3.029 → 3.014 Gelem/s).

It accelerates address computation, not memory traffic, and the kernel is bound by the latter. This is the empirical anchor for treating bandwidth, not arithmetic, as the bottleneck — and for ordering our work accordingly. (The same pattern is well known in the wider community; it was the substance of QuEST's BMI2 work, which this project's author contributed to.)

The priority order

From highest to lowest end-to-end impact:

1. Gate fusion — merge adjacent 1–2 qubit gates into a single 2–5 qubit matrix so that K gates become one pass over the state. Fewer passes = fewer bytes moved. The biggest single win on structured circuits.
2. SoA layout — separate re[]/im[] arrays so SIMD complex multiply needs no de-interleave shuffles.
3. Cache-aware access / blocking — keep the working set hot and minimize passes; non-temporal stores to skip read-for-ownership on write-back.
4. Multithreading — the amplitude pairs are disjoint, so the loop is embarrassingly parallel; the only subtleties are a single-thread size threshold and load balancing.
5. Micro-optimizations — BMI2 index generation, unchecked indexing, ILP/unrolling, prefetch. Real but small; the long tail.

The order matters: fusion and layout change how much memory you move, which dominates; micro-ops only shave the constant on work you are already doing.

The milestone narrative

Rather than ship one optimized backend, qsv evolves through a sequence of separate backends, each a self-contained, benchmarkable diff and one data point on the headline throughput/roofline plot — a naive reference first, then each optimization layered on and measured in turn.

Speedups are regimes, not promises — the real numbers come from the benchmark harness and are reported honestly, including where we lose to mature simulators.

A deliberate sequencing choice

Fusion (v0.8) lands after SIMD and threading on purpose: its value is best measured as a multiplier on top of an already-fast kernel, which is the realistic and most informative framing.

Knowing when to stop

Optimization has diminishing returns, and chasing them is a trap. The stop criterion is built into the methodology: when the dominant kernels sustain ≥ 70–80% of the machine's measured STREAM-triad bandwidth, further micro-optimization is noise — the remaining effort should go to fusion or algorithmic improvements that change the bytes-moved equation, not to shaving the constant. Stating this explicitly is part of the story.

A note on Apple Silicon

The primary development machine is an Apple M3 Pro: 128-bit NEON (only 2× f64 lanes) and a heterogeneous 5 performance + 6 efficiency core layout. Two honest consequences we report rather than hide:

• the SIMD win on NEON is modest (~1.3–1.8×); the 2–4× payoff appears on x86 AVX2/AVX-512, which is why we benchmark on both;
• equal static thread partitioning makes the slow efficiency cores stragglers, so we use work-stealing with many small chunks and report P-only vs all-core scaling separately.

The core kernel

Every production statevector simulator — qsim, Qiskit-Aer, QuEST, Yao.jl, cuStateVec — shares the same indexing trick at its heart. qsv calls it insert_zero_bit. Understanding it is understanding 90% of how a statevector simulator works.

The problem

Applying a 1-qubit gate $G=\left(\begin{array}{cc}{g}_{00}& g_{01}\\ g_{10}& g_{11}\end{array}\right)$ to qubit $q$ couples exactly the amplitude pairs whose basis indices differ only in bit $q$:

$$\begin{array}{rl}{\psi }_{a_0}^{\prime }& =g_{00}{\psi }_{a_0}+g_{01}{\psi }_{a_1}\\ {\psi }_{a_1}^{\prime }& =g_{10}{\psi }_{a_0}+g_{11}{\psi }_{a_1}\end{array}\text{where }{a}_1=a_0\oplus 2^q.$$

There are $2^{N-1}$ such pairs. The naive approach builds a $2^N\times 2^N$ matrix and multiplies — catastrophically wasteful in both memory and compute. The reshape approach (v0.1) avoids the matrix but copies the whole state per gate. The efficient approach visits each pair exactly once, in place.

Enumerating the pairs

We loop a counter i over 0 .. 2^(N-1) and reconstruct the two partner indices from it. The trick is to take i, which has $N-1$ bits, and insert a 0 bit at position q to get $a_0$; flipping that bit gives $a_1$.

#![allow(unused)]
fn main() {
/// Insert a `0` bit at position `bit`, shifting higher bits up by one.
pub fn insert_zero_bit(index: usize, bit: u32) -> usize {
    let left  = (index >> bit) << bit; // bits ≥ bit
    let right = index - left;          // bits < bit
    (left << 1) | right
}

pub fn flip_bit(index: usize, bit: u32) -> usize {
    index ^ (1usize << bit)
}
}

insert_zero_bit splits index at position bit, shifts the high part up by one to open a gap, and ORs the low part back. As i ranges over all $2^{N-1}$ values, $a_0=\text{insert\_zero\_bit}(i,q)$ ranges over exactly the indices with bit q clear — every pair, once.

This indexing trick is the shared heart of every production statevector simulator (it appears under various names — insert-zero-bit, index0, bit-deposit), and on x86-64 it maps directly to the PDEP instruction — see How we optimize for the measured BMI2 path and why it is not the bottleneck.

The in-place kernel

Putting it together (qsv's BitShiftBackend::apply_1q, lightly trimmed):

#![allow(unused)]
fn main() {
let pairs = state.dim() >> 1;
let (re, im) = state.parts_mut();   // SoA: separate real/imag slices
for i in 0..pairs {
    let a0 = insert_zero_bit(i, q);
    let a1 = a0 | (1usize << q);
    let x0 = Cplx::new(re[a0], im[a0]);
    let x1 = Cplx::new(re[a1], im[a1]);
    re[a0] = (g00 * x0 + g01 * x1).re;  im[a0] = (g00 * x0 + g01 * x1).im;
    re[a1] = (g10 * x0 + g11 * x1).re;  im[a1] = (g10 * x0 + g11 * x1).im;
}
}

No $2^N$ matrix, no per-gate allocation, no copy — just one streaming pass over the state, which (per How we optimize) is exactly the bandwidth-bound minimum.

Generalizing

Multi-qubit gates. For an $m$-qubit gate, insert $m$ zero bits at the sorted target positions to anchor each of the $2^{N-m}$ blocks, then enumerate the $2^m$ sub-indices within a block (insert_zero_bits + scatter_bits). qsv's apply_mq gathers a block's $2^m$ amplitudes into a small stack buffer, applies the matrix, and writes back — still in place, still zero-allocation.

Controlled gates. A controlled gate only changes amplitudes where the control bits are 1, so the optimized form iterates just the $2^{N-k}$ active subspace (with $k$ controls), cutting work by $2^k$. (In the current milestone, controlled gates are handled as plain dense gates via apply_mq; the subspace specialization is a later optimization.)

High vs low target qubits. When q is small the pair stride $2^q$ is tiny and both partners fall in the same cache line / SIMD register; when q is large the stride is huge and the access pattern is cache-hostile. qsim and Aer use two code paths — a permute-within- register kernel for low qubits and a blocked streaming kernel for high ones. This is qsv's v0.4 milestone.

Why a separate, different oracle

The optimized kernels all rest on this index arithmetic, so a bug in insert_zero_bit would be invisible to a test that used the same trick. qsv's RefBackend oracle therefore applies gates a structurally different way — gather/scatter per output amplitude — so the differential tests genuinely cross-check the indexing. The helpers themselves are also exhaustively unit-tested (every pair differs in exactly the target bit; the block anchors tile the index space).

Roadmap & milestones

This page tracks status and is updated as each milestone lands. The rationale for the order is in How we optimize.

Status

¹ SIMD measured ~0% on the 1q kernel — it's bandwidth-bound, not arithmetic-bound. See benchmarking. Expected to matter for fused multi-qubit matvecs and x86 AVX-512. ² Promoted to "next" by the v0.8 finding: fusion's win is currently capped at large N by the multi-qubit kernel's scattered access — cache-blocking it is the unlock.

Also planned but out of the v1 critical path: density-matrix / noise simulation, a GPU backend (CUDA/cuTile or Metal) behind the Backend seam, and distributed multi-node support (mirroring QuEST's pairwise rank exchange).

Testing

Robustness rests on differential testing: every optimized backend must reproduce the naive RefBackend oracle amplitude-for-amplitude. Because the oracle is implemented a structurally different way (gather/scatter, not the bit-shift pairing), this genuinely cross-checks the kernels rather than re-running the same logic.

Current suite (grows with each milestone):

• 200 reproducible random circuits (3–8 qubits, depth 20–60) — every optimized backend (Reshape, BitShift, CpuBackend serial and parallel) vs oracle;
• 15 random circuits at 10 qubits, depth 80 — BitShiftBackend vs oracle;
• 6 random circuits at 14 qubits (above the threading threshold) — CpuBackend's rayon path, including the multi-qubit parallel kernel, vs oracle;
• QFT → uniform superposition (n = 1,2,3,5,8) and QFT vs oracle;
• a 3-qubit Toffoli across several qubit orderings, exercising the general apply_mq path;
• unit tests for the index arithmetic (pairs differ in exactly the target bit; block anchors tile the space) and for the RNG/generators.

Every kernel added in a later milestone is wired into this same differential harness, so the optimization can never silently break correctness.

Definition of done for a milestone

1. New backend (or kernel) implemented behind the Backend trait.
2. Differential tests against the oracle pass, including any new path it introduces.
3. cargo clippy --all-targets -- -D warnings and cargo fmt --check clean.
4. A benchmark data point recorded (see Benchmarking & profiling).
5. This roadmap and the relevant design pages updated.

Benchmarking & profiling

Measurement is the whole point: the optimization strategy only means something if each milestone's effect is quantified honestly. The qsv-bench crate holds the criterion benchmarks and a profiling workload.

What we measure

Throughput is reported in amplitude-updates per second — criterion's Throughput::Elements set to $2^n$ per gate. This normalizes across qubit counts and is the natural figure of merit for a bandwidth-bound kernel (it converts directly to effective GB/s: multiply by the bytes touched per amplitude).

Three benchmark groups:

The single_h_gate group runs the slow milestone backends (oracle, reshape) only at small n, so a single run shows the v0.0 → v0.1 → v0.2 progression without taking forever.

Running

cargo bench -p qsv-bench --bench throughput              # everything
cargo bench -p qsv-bench --bench throughput -- bitshift  # filter by name (regex)

# Quick look (short sampling) while iterating:
cargo bench -p qsv-bench --bench throughput -- \
    'single_h_gate/bitshift' --warm-up-time 0.3 --measurement-time 1.0 --sample-size 10

criterion writes HTML reports (with plots) to target/criterion/.

Milestone results so far

A single Hadamard on a mid-range qubit, throughput in amplitude-updates/sec, Apple M3 Pro (11 cores). Higher is better.

¹ below the threading threshold (n < 13), so parallel runs serially — by design.

What each milestone bought, read honestly:

• v0.2 → v0.3/4 (cpu_serial): ~1.9×. Removing per-amplitude bounds checks (via iterator zips that lower to checked-free code — no unsafe needed) and the cache-friendly nested-block walk. Squarely in the predicted 1.3–2× regime.
• v0.3/4 → v0.5 (cpu_parallel): ~2–3.4× more, depending on size.

The kernel just became bandwidth-bound — and we can prove it

Convert the v0.5 throughput to effective memory bandwidth: each update reads and writes a 16-byte complex amplitude = 32 bytes of traffic.

$$\text{n=24:}3.97\text{Gelem/s}\times 32\text{B}\approx 127GB/s\approx 85\%\text{of the M3 Pro’s}\sim 150\text{GB/s peak.}$$

At 256 MB (pure DRAM) the threaded kernel is running into the memory wall: ~85% of peak bandwidth, exactly the regime the thesis predicted. At n=20 (16 MB, partly L2-resident) it reports ~205 GB/s — above DRAM peak — because some traffic is served from cache.

This is the inflection the project was built to surface. Earlier, the single-threaded scalar kernel was a flat ~1 Gelem/s — compute-bound on per-element work, not memory. Removing that work (v0.3) and parallelizing it (v0.5) has now pushed the DRAM-resident case to ~85% of the bandwidth roof. Per the stop criterion (≥ 70–80% of STREAM bandwidth), micro-optimizing this kernel further is noise — the next real win must change the bytes-moved equation, i.e. gate fusion (v0.8). The data wrote the roadmap.

v0.6–v0.8: three findings, reported honestly

v0.6 diagonal fast path — Z/S/T/PHASE/RZ/CZ/RZZ run a single sequential pass (one complex-mul per amplitude, no 2^q stride). Real win for phase-heavy circuits (QFT/QAOA) and high-target-qubit gates; modest at large N where everything is bandwidth-bound.

v0.7 SIMD — a measured null result. wide::f64x4 on the 1-qubit kernel: ~0% at every size (n=12 2.01 vs 2.04, n=20 1.89 vs 1.89, n=24 1.79 vs 1.80 Gelem/s vs the scalar CpuBackend). The 1q gate's arithmetic intensity (~0.13 FLOP/byte) is so low that widening the ALUs changes nothing — the cache-friendly scalar kernel is already bandwidth-bound at L2 and DRAM. This is the roofline taken to its conclusion, not a bug. SIMD is expected to pay where arithmetic intensity is high: the fused multi-qubit matvecs below, and x86 AVX-512 (64-byte loads vs NEON's 16) — flagged for the Intel box (see todo.md).

v0.8 gate fusion — the headline, with a caveat the benchmark surfaced. Fused vs unfused QFT on CpuBackend::parallel():

Fusion cuts passes over memory, so at n=14 it's a clean ~1.8×. At n=18 it washes out — and why is the interesting part: unfused QFT's controlled-phases already use the fast diagonal kernel, whereas fused H+phase blocks run the general apply_mq kernel, whose scattered gather/scatter erodes effective bandwidth at large N enough to cancel the pass reduction. So the fusion win is real but currently gated by the multi-qubit kernel's access pattern — exactly what v0.9 (cache-blocking apply_mq) targets. The benchmark didn't just validate fusion; it located the next bottleneck.

Roofline methodology

1. Measure the machine's empirical peak bandwidth with a STREAM-triad microbench (don't trust the spec sheet).
2. For each kernel, compute arithmetic intensity (FLOP / byte touched) and plot it against achieved performance.
3. As milestones land, watch points move: SIMD/threading pushes the scalar kernel up toward the bandwidth roof; fusion moves points rightward (higher arithmetic intensity) — the visual proof of the central thesis.

Profiling

qsv-bench ships a long-running workload, qsv-profile, for sampling profilers:

cargo build --release -p qsv-bench --bin qsv-profile

# macOS / Linux — samply (Firefox-profiler UI)
samply record ./target/release/qsv-profile 22 40

# macOS — Xcode Instruments
cargo instruments -t "Time Profiler" --release --bin qsv-profile -- 22 40

# Linux — perf with memory-bandwidth counters
perf stat -e cache-misses,mem_load_retired.l3_miss ./target/release/qsv-profile 22 40
likwid-perfctr -g MEM_DP ./target/release/qsv-profile 22 40   # reports achieved GB/s

Arguments are <qubits> <layers> <reps>; the workload runs layers·n random gates reps times so the profiler gets enough samples to attribute time to the kernel.

Cross-simulator comparison (planned)

The fair-comparison harness against Qiskit-Aer, qsim, QuEST, and spinoza lands with the later milestones. The protocol: one circuit definition exported to every tool; time only the statevector evolution (exclude Python import / transpile / JIT); match precision and thread count; same physical box for head-to-heads; report both wall-clock and gate-throughput; and show where we lose, explained by the roofline.

Research notes — overview

Distilled findings from studying production statevector simulators, which justify the design decisions throughout this book. The full reference source is shallow-cloned under _local/ (git-ignored).

The one finding that organizes everything

Statevector simulation is memory-bandwidth-bound, not compute-bound.

A 1-qubit gate streams the entire $2^N$-amplitude array doing only ~2 complex multiplies per 16-byte amplitude → arithmetic intensity ≈ 0.13 FLOP/byte, deep in the bandwidth-bound region of the roofline. Confirmed independently across qsim, Qiskit-Aer, QuEST, Yao.jl, cuStateVec, and spinoza. See How we optimize for the consequences.

Reference repositories studied

cuda-quantum was intentionally not cloned (multi-GB); its cuStateVec approach is captured in GPU, Rust & HPC landscape, and GPU work is deferred behind the Backend seam.

The two pages that follow go into detail on the CPU simulators and the GPU / Rust / HPC landscape.

CPU statevector simulators

Concrete, reimplementable techniques distilled from the production CPU simulators (qsim, Qiskit-Aer, QuEST, Yao.jl). Source lives in the shallow clones under _local/.

The universal core kernel

All of them pair amplitudes via the same insert-zero-bit indexing — explained from scratch in The core kernel. It surfaces under different names (insertZeroBit, index0, bit-deposit/PDEP, controlled iteration over the active subspace), but the arithmetic is identical. It generalizes to m-qubit gates (insert m zero bits) and to controlled gates (iterate only the $2^{N-k}$ active subspace).

Memory layout — SoA vs AoS

• qsim: Structure-of-Arrays (separate real / imaginary blocks, AVX-width aligned) → SIMD complex multiply with no lane shuffles.
• Qiskit-Aer: AoS (interleaved Complex), then works around it with separate real/imag views and de/re-interleave in the AVX2 path.

qsv chooses SoA — see Architecture overview.

SIMD complex multiply

$(a+bi)(c+di)=(ac-bd)+(ad+bc)i$, vectorized with FMA: the gate entry is a broadcast scalar, the amplitudes are SIMD lanes. NEON is 128-bit (2× f64 lanes) so the Apple win is modest; AVX2/AVX-512 give 4×/8× lanes.

High- vs low-qubit dispatch (qsim, Aer)

Low target qubit (small stride) → permute within registers; high target qubit (large stride) → blocked streaming. Two code paths per gate.

Gate fusion — the biggest end-to-end win

Merge adjacent 1–2 qubit gates into a $2^m\times 2^m$ matrix (m ≈ 4–5) so K gates → 1 pass. Cost model under the bandwidth-bound view: each gate ≈ one streaming pass regardless of arity, until the fused matrix is large enough that the $2^m$ matvec overtakes bandwidth. (qsim fuser_mqubit.h; Aer fusion transpiler pass.)

Diagonal-gate fast path (Aer)

Z, S, T, PHASE, RZ, CZ, RZZ are diagonal: a single pass multiplying each amplitude by a phase — no pairing, no second load. QFT and QAOA are phase-heavy, so this is a cheap ~2× there.

Threading

qsim uses a custom ParallelFor with a size threshold (small loops stay single-threaded) and static partitioning; QuEST uses OpenMP collapse(2). The amplitude-pair work is disjoint, so it is embarrassingly parallel.

Profiling-driven lessons qsv adopts

A set of cross-cutting techniques the profiling literature (and our own benchmarks) keep surfacing:

• Type/representation stability is often the single biggest win — mixed or promoted scalar types silently dominate runtime. qsv's analog: keep amplitudes and gate entries the same concrete Real, monomorphized (never dyn).
• Unchecked indexing — a modest gain after representation is fixed (qsv gets it from iterator-shaped loops, no unsafe).
• Per-thread accumulators for expectation values; alias-table $O(1)$ sampling; stack-allocated small gate matrices. qsv implements the last two directly (see crate::sample and the MAX_SUB stack buffer).

QuEST distributed simulation (future multi-node seam)

When the target qubit is higher than the locally-stored qubits, the partners live on different MPI ranks. QuEST pairs ranks by XOR-toggling the relevant bit, exchanges half-buffers with MPI_Sendrecv (deadlock-free, all pairs concurrent), then applies the 2×2 locally. Captured for a future distributed backend behind the same Backend seam.

BMI2 (PEXT/PDEP) — the bandwidth-bound proof

The bit-gather/scatter/insert at the kernel's core map to single PEXT/PDEP instructions on x86-64. qsv's bmi2 feature and index_gen microbenchmark measure this directly: ~4× faster in isolation (PEXT/PDEP vs the scalar loop) but 0.995× end-to-end on a fused QFT-18 — because they accelerate address computation, not memory traffic, and the kernel is bandwidth-bound. This is the empirical anchor for the whole optimization strategy (numbers in bench/results/SUMMARY-xeon.md; the same effect is well documented in the wider community, e.g. QuEST's BMI2 work).

GPU, Rust & HPC landscape

Landscape research for the (deferred) GPU backend and the CPU optimization toolbox.

cuStateVec / cuda-quantum

NVIDIA's cuStateVec applies gates with the same bit-index pairing as the CPU, mapped to the GPU:

• thread → amplitude-pair mapping; the gate matrix is staged into shared memory and reused across a block's pairs;
• coalescing depends on the target-qubit stride — low-order targets coalesce, high-order targets stride badly and are mitigated by qubit reordering (the GPU analog of QuEST's distributed exchange);
• multi-GPU partitions the state by index bits; gates on "global" qubits need GPU↔GPU exchange.

The statevector lives in HBM (the bottleneck); achieved bandwidth is ~60–80% of peak — the same bandwidth-bound story as the CPU. cuda-quantum's nvq++ lowers circuits to custatevecApplyMatrix calls.

cuTile / cuTile-rs — ✅ verified (2026-06)

cuTile is real and shipping: NVlabs/cutile-rs (crate cutile, Rust 1.89+) and NVIDIA/cutile-python, with the paper Fearless Concurrency on the GPU (arXiv 2606.15991). Both clones are under _local/. The full memory-architecture study and the build-vs-buy decision live in the cuTile memory architecture note; the essentials:

• Requirements: GPU sm_80+ (Ada sm_89 ✅), CUDA 13.2+ (13.3 recommended), driver r580+, Linux. Our box is CUDA 12.4 / driver 550 → cuTile cannot run here without a toolkit + driver upgrade (sysadmin-level). cuTile is early-stage ("expect API breakage").
• Memory model: programmer picks tile shapes; the compiler owns coalescing, shared-memory staging, TMA pipelining, and tensor-core lowering. Tiles live in registers; tensors live in HBM; load/store moves between them. There is no native complex type, and strided access is the pattern cuTile's own docs flag as most likely to lose bandwidth — with no manual remedy.
• Decision: for our memory-bound, strided, complex, non-GEMM gate kernel, the performance levers (coalescing, __shared__ staging, __shfl pair exchange, vectorized complex loads) are exactly what cuTile abstracts away. → CudaBackend is built on cudarc + NVRTC (hand-written CUDA C, full memory control, runs on CUDA 12.4 today). cuTile stays as a future validation reference and for the contiguous/low-target-qubit cases.

Rust quantum-sim ecosystem — the gap

No native-Rust statevector simulator combines gate fusion + cache-blocking + SoA + SIMD behind a modern pluggable GPU seam — that is qsv's niche.

amh-code — CPU optimization toolbox

Techniques from Algorithms for Modern Hardware, with rough impact for a bandwidth-bound complex-array kernel:

Net plan implications

• GPU stays behind Backend (associated type State, on-device reductions, single download); the RefBackend second implementation proves the seam today.
• CPU optimization order follows the bandwidth-bound thesis: fusion → SoA+SIMD → cache-blocking → threading → micro-ops, each a benchmarkable milestone.

cuTile memory architecture — investigation

A deep-dive into NVIDIA's cuTile tile-based GPU programming model, written to answer one question for this project: should the GPU backend be written in cuTile (cuTile-rs), or in hand-written CUDA C kernels driven from Rust (cudarc + NVRTC)?

The short answer, derived below: cuTile is real, mature enough to evaluate, and an excellent fit for contiguous/GEMM-shaped work — but for our memory-bound, strided, complex-valued, non-GEMM gate-apply kernel, hand-written cudarc/NVRTC kernels are the better fit, because the workload's performance hinges on exactly the low-level controls (coalescing, shared-memory staging, vectorized complex loads, warp shuffles) that cuTile deliberately abstracts away. cuTile stays on the table as a validation reference and for the contiguous (low-target-qubit) cases.

Sources: the two local clones under _local/cutile-rs and _local/cutile-python, the cuTile-Rust paper Fearless Concurrency on the GPU (arXiv 2606.15991), NVIDIA's cuTile-Python guide, and the cuTile.jl blog.

1. Tiles and the memory hierarchy

cuTile splits state into two objects:

(_local/cutile-rs/cutile-book/guide/tensors-and-tiles.md; mirrored in the cuTile-Python guide.)

A tensor is the addressable, strided object in HBM. A tile is a register-resident, statically shaped fragment that exists only inside the kernel. You load a tile from a tensor, compute on tiles, and store a tile back.

Who decides placement

The programmer chooses tile shapes and access patterns; the compiler owns every level below the tile. From useful-mental-models.md:

Registers — fastest storage; tiles live here during computation. Shared memory — fast on-chip storage shared within a hardware block. L2 cache — hardware-managed, shared across SMs. HBM — large global memory where tensors live.

"In cuTile Rust, you load from tensors in HBM and compute on tiles in registers. The Tile IR compiler and runtime decide how to stage data through shared memory, caches, threads, warps, Tensor Cores, and Tensor Memory Accelerator (TMA) instructions when those mechanisms are useful."

So shared-memory staging, TMA pipelining, and tensor-core (MMA) lowering are compiler-owned. The programmer's knobs are: tile shape, partition shape, dtype, algorithmic access pattern, an opt-in mma call to reach tensor cores, and a few coarse optimization_hints (occupancy, CTA-in-CGA, divisibility). The Triton-migration table in interoperability.md is explicit that "the compiler generates shared memory staging for load_tile operations" and selects TMA automatically.

Load / store / partition API

cuTile-Rust (reference/dsl-api.md):

#![allow(unused)]
fn main() {
tensor.load()                       // -> Tile          (load the output tile)
tensor.store(tile)                  //                  (store a tile to the tensor)
tensor.load_tile(shape, idx)        // -> Tile          (load at a partition index)
load_tile_like(src, dst)            // -> Tile          (load src at dst's tile-block position)

let part_x  = x.partition(const_shape![BM, BK]);   // device-side partition
let tile_x  = part_x.load([pid.0, k_tile]);
}

cuTile-Python (samples/VectorAddition.py):

a_tile   = ct.load(a, index=(bid,), shape=(TILE,))   # HBM -> SMEM/registers (auto-distributed)
sum_tile = a_tile + b_tile
ct.store(c, index=(bid,), tile=sum_tile)             # tile -> HBM

The sample's own comment: "ct.load automatically distributes the load across the threads within the block, bringing the tile into shared memory or registers." You do not write the thread-index math.

2. Execution model vs SIMT

Traditional CUDA C++ is thread-per-element: you compute i = blockIdx.x*blockDim.x + threadIdx.x and write c[i] = a[i] + b[i]. cuTile is tile-per-block: the entry function is written once as straight-line, scalar-looking code over whole tiles, and the compiler performs the parallel decomposition. From the paper:

"The programmer writes sequential code over multi-dimensional tiles, and the compiler maps tile operations to thread blocks, manages shared memory, and performs the parallel decomposition."

"Tile-based programming gives up SIMT-level control (explicit warp primitives, shared memory management) in exchange for the single-threaded semantics that make static safety checking tractable."

The critical caveat for us, from useful-mental-models.md (Coalescing and Strides): "Tile loads are designed to produce coalesced access patterns for regular layouts… Strided access can reduce effective bandwidth because memory requests become scattered." You cannot hand-fix coalescing — you can only restructure the tile/partition shape and re-measure.

3. The Rust ownership / disjoint-partition model

cuTile-Rust extends Rust's aliasing-XOR-mutability across the launch boundary: mutable output tensors are partitioned into disjoint sub-tensors on the host before launch; immutable inputs are shared. The quick-start kernel (_local/cutile-rs/README.md):

#![allow(unused)]
fn main() {
#[cutile::entry()]
fn add<const B: i32>(
    z: &mut Tensor<f32, { [B] }>,   // exclusive mutable output
    x: &Tensor<f32, { [-1] }>,      // shared read-only
    y: &Tensor<f32, { [-1] }>,      // shared read-only
) {
    let tx = load_tile_like(x, z);
    let ty = load_tile_like(y, z);
    z.store(tx + ty);
}
// host:
let z = api::zeros::<f32>(&[1024]).partition([128]);  // disjoint 128-elt chunks
let (_z, _x, _y) = kernel::add(z, x, y).sync()?;      // grid (8,1,1) inferred
}

The invariant (paper): "The mapping from tile programs to sub-tensors is injective: no sub-tensor is assigned to more than one tile program." The store target is the partition view itself, not an index the programmer picks — so the classic index-swap data race is inexpressible. Memory ordering is token-based: mutable references thread a t₀ →load→ t₁ →store→ t₂ token chain establishing happens-before; immutable references emit no tokens and may be freely reordered. Raw *mut T device-pointer entries and unchecked_accesses are the opt-outs.

4. Host execution model

Every host-side call is a lazy DeviceOp (a GPU Future). Three modes over identical kernel code:

Combinators mirror futures (.then, zip!, .map, .shared). CUDA-graph replay swaps inputs and re-launches without recompilation (~0.8 µs/op overhead reported); only non-allocating ops (kernel launches, memcpy) can be graph nodes. cuTile-Python is simpler: ct.launch(stream, grid, kernel, args) plus a CUDA-graph path.

Relevance to qsv: this is genuinely attractive for a long gate sweep — capture the per-gate launches once, replay across the circuit. Our Backend::execute override (batching the whole circuit on one stream) maps cleanly onto either cuTile graphs or cudarc's manual graph API.

5. Requirements (and our box)

Our machine: 2× L40S (Ada sm_89, supported), but CUDA 12.4 / driver 550.127.08. cuTile needs CUDA 13.2+ and driver r580+ — so it cannot run here without a toolkit + driver upgrade (the driver bump is system-level / sysadmin). Note the irony: in early builds Ada is better covered than Hopper (sm_90 only landed in CUDA 13.3). Both projects self-describe as early-stage ("expect bugs, incomplete features, API breakage").

6. Suitability for our gate-apply kernel

Our kernel applies a 1- or 2-qubit gate to the statevector: an elementwise complex multiply combined with a strided amplitude-pair gather/scatter (pairs separated by stride 2^t for target qubit t), bandwidth-bound, goal ~60–80% of HBM peak.

cuTile is not GEMM-only. Memory-bound elementwise is a first-class target: the paper reports 7 TB/s for element-wise ops on B200 (~91% of peak HBM). So 60–80% of peak is a goal cuTile routinely beats for contiguous layouts. It also can express gather/scatter: cuTile-Python ships ct.gather/ct.scatter with index tiles; cuTile-Rust exposes PointerTile, load_ptr_tko / store_ptr_tko, addptr_tile, and atomics. A clean framing is to reshape the statevector so the target qubit splits the flat index into [high, 2, low] and make the size-2 pair axis a tile dimension — then a low/contiguous-target gate becomes a coalesced tile load.

But the real limitations for this workload:

1. No control over coalescing or SMEM staging — the crux. A hand-written CUDA C kernel lets you choose the exact strided index math, use double2/float4 vectorized loads, stage amplitude pairs in __shared__, use __shfl_xor for the in-warp pair exchange, and tune block size to the stride — the levers that take a strided gate kernel from ~40% to 70%+ of peak. cuTile removes those levers by design.
2. Strided access is the one pattern cuTile explicitly flags as a risk. performance.md: "tile loads coalesce well, but algorithmic strides can still reduce effective bandwidth." The bad pattern it names — memory[0], memory[1024], memory[2048], … — is exactly a high-target-qubit gate. The only recourse is to re-tile and re-measure, not hand-optimize.
3. No native complex type. The ElementType set is f16/bf16/f32/f64, ints, FP8/FP4, tf32 — no complex. You carry real/imag yourself (interleaved or a trailing size-2 axis) and write the complex multiply by hand — and the interleaved layout you'd want for coalescing is exactly what you can't force the compiler to honor.
4. Maturity / portability gates: early-stage, Linux-only, CUDA 13.2/13.3 + driver r580.

cuTile vs hand-written cudarc/NVRTC

Verdict

For our memory-bound, strided, complex-valued, non-GEMM gate-apply kernel — targeting 60–80% of HBM peak — hand-written cudarc/NVRTC kernels are the better fit. Not because cuTile is GEMM-only (it delivers ~91% of HBM peak on contiguous elementwise work), but because this workload's performance hinges on exactly the controls cuTile abstracts away, its own docs name strided access as the highest-risk pattern with no manual remedy, and it has no native complex type — and it can't even run on this box's CUDA 12.4 today.

cuTile stays relevant as: (a) a fast-to-write validation reference once a CUDA-13 box is available; (b) a strong option for the contiguous / low-target-qubit cases where a reshape makes the gate axis a clean tile dimension; and (c) interop — cuTile-Rust's borrow_raw / cudarc_interop lets a hand-tuned cudarc kernel handle the hard strided gates while cuTile handles the regular bulk, sharing one CUstream. The structural race-freedom and graph-replay ergonomics are genuinely nice. But to reliably hit the bandwidth roof on the strided cases, we need the low-level memory control that only a hand-written kernel provides — so qsv's CudaBackend is built on cudarc + NVRTC (see the GPU landscape note).