Systems · GPU performance — 2026
Thirty-four CUDA kernels written from scratch and taken past NVIDIA's own library — one optimization at a time, with every stage measured so the value of each technique is visible rather than asserted.
Shipped — public, CI green across four GPU architectures
Image slot
warpsmith — SGEMM throughput across nine optimization stages
1600×1000 · the stage bar chart from docs/charts/sgemm-progression.svg
Role
Solo build — kernels, measurement harness, tooling, documentation
Year
2026
Stack
- CUDA C++
- cuBLAS
- CUB
- WMMA / Tensor Cores
- CMake
- Python
- GitHub Actions
The problem
A GPU is almost never slow because it runs out of arithmetic. It is slow because the arithmetic units are waiting on memory, and the material that explains this tends either to stay abstract or to hand over a finished kernel with no account of how it got there. I wanted the account: the same operation implemented repeatedly, each version differing from the last by exactly one idea, so that the gap between two adjacent versions measures that idea instead of describing it.
Approach
Rebuilt single-precision matrix multiply in nine stages — uncoalesced naive, then global-memory coalescing, shared-memory tiling, 1D and 2D register tiling with an outer-product inner loop, float4 vectorization over a transposed tile, warp-level tiling, and finally tensor cores through WMMA in TF32.
Wrote a fused attention kernel that never materializes the S×S score matrix, carrying a running maximum and denominator through the accumulator with the online-softmax correction — the trick that makes FlashAttention possible — plus a causal variant that skips whole key tiles above the diagonal before issuing a single dot product.
Covered the bandwidth-bound half of the problem as well: six reduction formulations from global atomics to vectorized warp shuffles, transpose kernels that isolate coalescing from shared-memory bank conflicts, online softmax, and RMSNorm with the residual fusion that production code always does.
Built the measurement harness first — CUDA events on the stream, untimed warm-up launches, trimmed-median statistics with p95 and a coefficient of variation on every row, and correctness validated against cuBLAS, CUB or a double-precision host reference before anything is timed.
Made the documentation a build artifact: the report, the charts and the README's own tables are generated from the results file, and CI regenerates them and fails if they disagree with what was measured.
Results
- of cuBLAS FP32 throughput at 4096³
114%
of cuBLAS FP32 throughput at 4096³
- naive to fastest, same problem
48×
naive to fastest, same problem
- fused causal attention
2.83×
fused causal attention
materialized, and 64 MiB of scratch removed
- of theoretical bandwidth on reduction
88%
of theoretical bandwidth on reduction
- correctness failures across 79 measurements
0
correctness failures across 79 measurements
Measured on a thermally constrained 4 GB laptop GPU, so absolute throughput sits well below a datacentre card and the ratios are the portable result. The tensor-core kernel is the weakest number in the repository — roughly a fifth of the TF32 ceiling — and the README says so, along with the specific pipelining it lacks.
Image slot
Roofline plot of every kernel against the memory and compute ceilings
1600×1000
Image slot
Throughput and register pressure side by side across stages
1600×1000
What it came to
The finding I did not expect was how much of the total came from a single index swap: making a warp read 128 contiguous bytes instead of 32 scattered ones was worth 5.6× on its own, more than any other stage contributed. The second was that occupancy moved the wrong way — the fastest kernels hold so many accumulators in registers that fewer warps fit on an SM, and they win anyway, because what the hardware needs is independent work in flight rather than resident threads.