CUDA C++ · inline PTX · AVX-512 intrinsics · NVIDIA RTX 4000 Ada · SM-level performance engineering
A tiled SGEMM kernel taken from 27.12 ms to 7.19 ms — a 3.77× speedup on a 3072×3072×3072 matrix multiply — by moving the working set from DRAM into shared memory and then into registers. Alongside it: GPU cache-latency and coalescing microbenchmarks written in inline PTX, a three-kernel parallel prefix sum, and a run-length encoder built on top of it.
Mandelbrot set, rendered by the CUDA kernel — 311 ms for the scalar GPU version down to 14.1 ms once the kernel used the full 32-wide SIMT width.
A matrix multiply is four lines of code and is almost never limited by arithmetic. On an RTX 4000 Ada, a 3072³ SGEMM is 58 GFLOP against 113 MB of data — compute-bound in theory, which means every millisecond above the roofline is a memory-system or scheduling failure rather than a math failure.
So the work is entirely in data movement: which values live in registers versus shared memory versus L1, whether a warp's 32 loads collapse into one memory transaction or scatter into 32, and how many times each byte fetched from DRAM gets reused before it is evicted. Each section here isolates one of those and measures it.
| # | Component | Techniques | Result |
|---|---|---|---|
| 1 | Mandelbrot: SIMD and CUDA | AVX-512 intrinsics with mask registers; CUDA SIMT kernel | 29.7 → 2.99 ms (CPU), 311 → 14.1 ms (GPU) |
| 2 | Memory-hierarchy microbenchmarks | Inline PTX cache operators, pointer chasing, coalescing | 84 / 5 cycle latencies; 1.7× from coalescing alone |
| 3 | Tiled matmul | 32×32 shared-memory tiles, 4×4 register microtiles | 27.12 → 7.19 ms, 8.06 TFLOP/s |
| 4 | Parallel scan + RLE | Three-kernel reduce-then-scan; scan-based compression | 9.70 GB/s scan; 14.95 ms compression — parallel across blocks only, see below |
| Kernel | Time | TFLOP/s | Key change |
|---|---|---|---|
| Shared-memory tiling | 27.12 ms | 2.14 | Each value loaded from DRAM reused 32× out of shared memory |
| + register microtiling | 7.19 ms | 8.06 | 4×4 accumulator per thread held entirely in registers |
Each thread block computes a 128×128 tile of C; within it each thread owns a 4×4 microtile, giving 1024 threads per block. The inner loop loads 4 values of a and 4 of b from shared memory into registers, then issues all 16 FMAs from registers before touching shared memory again. SASS confirms this: LDS.U.128 vector loads followed by uninterrupted runs of FFMA.
Instruction-level parallelism helped the CPU and actively hurt the GPU. Processing two pixels per thread instead of one is textbook ILP, and on the CPU it delivered 1.6×. The identical change on the GPU made it 2.3× slower. Two causes compounded: the two pixels diverge — one escapes the Mandelbrot iteration before the other, leaving the thread idle on the finished pixel — and the extra live state raised register pressure enough to cut occupancy, undermining latency hiding. The penalty got worse with more threads (2.66× slower in the multithreaded configuration), which is the signature of divergence amplifying rather than a fixed overhead. The fastest GPU configuration used no ILP at all.
One microbenchmark returned a physically impossible result, and that was the finding. Measured latencies came out at 84 cycles for global memory, 5 for L2 — and 7 for L1. L1 cannot be slower than L2; it is closer to the SM. The explanation is that the .ca/.cg cache operators the benchmark relies on are hints the compiler may ignore, so the "L1" measurement was most likely never served from L1 at all. Worth more than a clean number: the benchmark was measuring something other than what it claimed, and reporting the figure at face value would have been wrong.
Coalescing is worth 1.7× for free. Two loops moving byte-for-byte identical data, differing only in indexing arithmetic: 0.2627 ms versus 0.1497 ms. When consecutive threads read consecutive addresses, a warp's 32 loads fall in one cache line and collapse into a single transaction; when each thread walks its own contiguous run, they scatter into 32. No algorithmic change, no extra memory — just which lane touches which address.
The scan is correct but leaves most of its parallelism on the table, and the measured throughput says so. The decomposition is right — reduce per block, scan the block totals, then re-walk each tile with its offset — but every kernel begins if (threadIdx.x != 0) return, so one thread per block does a sequential loop over all 256 elements while the other 255 exit immediately. The result is parallel across blocks and serial within them, which is exactly what 9.70 GB/s reflects: roughly an order of magnitude under what this GPU's memory system can sustain. The fix is a standard block-level scan — each thread scans ITEMS_PER_THREAD elements in registers, then a warp-shuffle scan combines the per-thread totals. Reported as measured rather than quietly omitted, because the gap between the design and the implementation is the interesting part.
Register reuse mattered nearly 4× more than getting into shared memory did. Shared-memory tiling is the well-known first move and delivered 2.14 TFLOP/s. But holding a 4×4 accumulator in registers on top of that was worth another 3.77×. Shared memory is fast relative to DRAM and slow relative to registers, and once the DRAM traffic is handled, shared-memory bandwidth becomes the next wall. The reuse hierarchy has more than two levels, and stopping at shared memory leaves most of the performance unclaimed.
The same SGEMM continued through two further stages that are not included in this repository — see the note on scope below. Reporting the results since they are the outcome of the same line of work:
| Stage | Time | TFLOP/s | Technique |
|---|---|---|---|
| Occupancy tuning + split-K | 4.98 ms | 11.65 | 256 threads/block and an 8×8 microtile to fit 2 blocks per SM; splitting the K dimension across blocks for small problems |
| TF32 tensor cores | 3.03 ms | 19.12 | Warp-cooperative mma.sync.aligned.m16n8k8 issued through inline PTX |
Two findings from those stages worth stating:
- The single largest speedup in the whole project — 28.87× — came from changing the decomposition, not optimizing the fast path. A 128×128×32768 matmul produces exactly one thread block, leaving 47 of 48 SMs idle no matter how well-tuned the inner loop is. Splitting K across blocks and reducing afterward turned a single-SM job into a whole-GPU one. No amount of register or shared-memory tuning could have found that.
- Tensor cores cost roughly 750× in accuracy, which the speed number hides. RRMSE went from 1.03e-06 with FMA to 7.78e-04 with TF32, a direct consequence of TF32's 10-bit mantissa versus FP32's 23. The speedup is a deliberate precision trade, not a free win.
Languages: CUDA C++, C++17
GPU: inline PTX (ld.global.{ca,cg,cv}, mma.sync), shared memory, register blocking, warp-level primitives, SASS/PTX inspection, NVIDIA RTX 4000 Ada (compute capability 8.9)
CPU: AVX-512 intrinsics with __mmask16 mask registers, pthreads
Method: roofline modeling, Little's Law, occupancy analysis
Written from scratch rather than pulled from a library: the shared-memory and register tiling schedules, the three-kernel parallel scan, the RLE compressor, the AVX-512 Mandelbrot kernel with mask-register control flow, and the PTX cache-operator latency probes.
01-mandelbrot-simd-and-cuda/ AVX-512 and CUDA SIMT kernels
02-memory-hierarchy-microbenchmarks/ cache latency + coalescing, inline PTX
03-tiled-matmul/ shared-memory and register tiling
04-parallel-scan-and-rle/ prefix sum and GPU compression
figures/
All measurements and analysis are in this README; each directory holds the kernels for that section.
This repository contains only kernels I wrote. They were built against an existing benchmark harness — timing loops, correctness checkers, data loaders and reference implementations — which is not reproduced here. The sources are therefore excerpts rather than buildable programs, and are here to be read.
The occupancy/split-K and tensor-core kernels are described above but not published.
See RUNNING.md for the hardware and toolchain the measurements depend on.