Lab 01 — GPU Roofline, Occupancy & Tiled-GEMM Performance Model

Phase: 18 — C++/CUDA & GPU Performance Engineering Difficulty: ⭐⭐⭐☆☆ (the arithmetic is small; the kernel intuition is ⭐⭐⭐⭐⭐) Time: 3–4 hours

You will not write a single line of CUDA in this lab — and that is the entire point. A senior GPU engineer can tell you, before touching a profiler, whether a kernel is compute- or memory-bound, what its occupancy will be and which resource caps it, how badly a strided access pattern bleeds bandwidth, how much tiling will cut HBM traffic, and what fusion buys. All of that is arithmetic over a mental model of the hardware. This lab builds that model: the roofline (deepening Phase 00 to the kernel level), the occupancy calculation that names the binding resource, memory coalescing (the #1 perf lever), tiling for data reuse (the canonical shared-memory GEMM optimization), operator fusion (the FlashAttention idea), and warp divergence.

What you build

• roofline_attainable / ridge_point / is_compute_bound — the roofline at the kernel level: min(peak, bandwidth·intensity), the ridge = peak/bandwidth, and the classifier that proves decode (I≈1) is far left of an A100's ridge (~156). Same min() as Phase 00, now aimed at one kernel.
• gemm_arithmetic_intensity — 2·M·N·K FLOPs over the bytes a naive GEMM moves, and the reason tiling raises it.
• occupancy — model the SM's four limits (register file, shared memory, warp slots, block slots), return the active-warp fraction and the limiting resource string. This is the soul of the lab: a register-heavy kernel is register-limited, a shared-mem-heavy one is smem-limited.
• coalescing_efficiency — bytes requested ÷ bytes transferred for a strided warp access; stride 1 ≈ 1.0, a large stride is proportionally worse. The single biggest memory lever.
• tiled_gemm_hbm_traffic — HBM bytes for a tile×tile blocked matmul; a larger tile reuses more, so it moves less data. The classic CUDA optimization, quantified.
• fused_vs_unfused_traffic — an element-wise chain: 2·n_ops array passes unfused vs 2 fused; the ratio is the n_ops× bandwidth-bound speedup ceiling (the FlashAttention idea).
• warp_divergence_factor — the execution-time multiplier when warp lanes branch differently (both paths run, masked) — 1 for uniform, >1 for divergent.

Key concepts

Files

Run

pip install -r requirements.txt
pytest test_lab.py -v                       # against your lab.py
LAB_MODULE=solution pytest test_lab.py -v   # against the reference
python solution.py                          # the worked example

Success criteria

• All tests pass against your implementation (and LAB_MODULE=solution).
• You can state the A100 ridge (~156 FLOP/byte) and explain why decode (I≈1) is memory-bound and a big GEMM (I≈500) is compute-bound — at the kernel level, not just the model level.
• You can explain why test_occupancy_register_limited and test_occupancy_shared_memory_limited return different limiter strings, and what you would change in each kernel to raise occupancy.
• You can explain why coalescing_efficiency(1) ≈ 1.0 and coalescing_efficiency(32) ≈ 1/32, and why coalescing is the first thing you check on a slow memory-bound kernel.
• You can explain why a larger tile in tiled_gemm_hbm_traffic moves less data, and tie it to the shared-memory GEMM in every CUDA tutorial.
• Given any "this kernel is slow" prompt you can, on paper, name the bound (roofline), the occupancy limiter, the coalescing efficiency, and the fix (tile/fuse/cut registers).

How this maps to the real stack

Limits of the miniature (be honest in the interview): the occupancy model uses round A100-ish SM limits and ignores register-allocation granularity (registers are rounded to a quantum), the warp/block-allocation quanta, and the static-vs-dynamic shared-memory split; the coalescing model counts aligned segments but ignores L2 caching, partial-sector loads, and the read/write asymmetry; the GEMM traffic model assumes a perfectly blocked schedule with full reuse and ignores the C-accumulation reads; high occupancy is necessary-ish but not sufficient (a low-occupancy, high-ILP kernel can beat a high-occupancy one — the occupancy-vs-ILP tradeoff). Real numbers come from a profiler on real silicon; this lab makes you predict them.

Extensions (build these on real hardware)

These are the bridge from "I can model a kernel" to "I wrote one." Each is a real nvcc/PyTorch build; do them on a machine with a GPU.

A) A coalesced vs strided copy kernel (the coalescing soul test, in CUDA)

// copy.cu — compile: nvcc -O3 copy.cu -o copy
__global__ void copy_coalesced(const float* __restrict__ in,
                               float* __restrict__ out, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;   // adjacent threads → adjacent addrs
    if (i < n) out[i] = in[i];                        // stride 1 → fully coalesced
}
__global__ void copy_strided(const float* __restrict__ in,
                             float* __restrict__ out, int n, int stride) {
    int i = (blockIdx.x * blockDim.x + threadIdx.x) * stride;  // scattered
    if (i < n) out[i] = in[i];                                 // bleeds bandwidth
}

Time both with CUDA events, compute achieved GB/s, and confirm the strided version's bandwidth drops by roughly the factor coalescing_efficiency(stride) predicts. Then run ncu and read the "sectors per request" — that is your efficiency number.

B) A tiled (shared-memory) matmul (the tiling soul test, in CUDA)

// tiled_gemm.cu — load TILE×TILE blocks of A and B into __shared__, reuse TILE times
#define TILE 32
__global__ void gemm_tiled(const float* A, const float* B, float* C, int M, int N, int K) {
    __shared__ float As[TILE][TILE], Bs[TILE][TILE];
    int row = blockIdx.y * TILE + threadIdx.y;
    int col = blockIdx.x * TILE + threadIdx.x;
    float acc = 0.f;
    for (int t = 0; t < (K + TILE - 1) / TILE; ++t) {
        As[threadIdx.y][threadIdx.x] = A[row * K + t * TILE + threadIdx.x];
        Bs[threadIdx.y][threadIdx.x] = B[(t * TILE + threadIdx.y) * N + col];
        __syncthreads();
        for (int k = 0; k < TILE; ++k) acc += As[threadIdx.y][k] * Bs[k][threadIdx.x];
        __syncthreads();
    }
    C[row * N + col] = acc;
}

Benchmark naive vs tiled vs cuBLAS; confirm the tiled version moves roughly the HBM bytes tiled_gemm_hbm_traffic(M, N, K, TILE) predicts, and that bigger TILE (within register/smem limits) helps until occupancy collapses — that is the occupancy-vs-reuse tradeoff this lab modeled.

C) Build it as a PyTorch C++/CUDA extension (the role of C++)

# setup.py — turn your .cu into a torch op you can call from Python
from torch.utils.cpp_extension import CUDAExtension, BuildExtension
from setuptools import setup
setup(name="mygemm",
      ext_modules=[CUDAExtension("mygemm", ["mygemm.cpp", "tiled_gemm.cu"])],
      cmdclass={"build_ext": BuildExtension})

Then import mygemm; mygemm.gemm(a, b) from Python and benchmark it against torch.matmul. This is the real "custom op" path. Compare the effort to writing the same kernel in Triton (@triton.jit, autotuned tiles) — the middle path that gives you most of the speedup for a fraction of the C++.

D) Profile and read the roofline

Run ncu --set roofline ./your_kernel and find where prefill, decode, and your GEMM land on the chart. Confirm decode sits far left (memory-bound) and the tiled GEMM sits near the ridge. Then run nsys profile to see kernel-launch overhead and gaps — the thing CUDA Graphs and stream overlap exist to kill (and why launch overhead bites decode hardest).

Interview / resume

• Talking points: "Is this kernel compute- or memory-bound — prove it from the roofline." "Your occupancy is 25% — which resource is binding and how do you raise it?" "Walk me through why AoS→SoA (coalescing) is a 10× win on a memory-bound kernel." "Why does tiling speed up GEMM, and what limits the tile size?" "What does FlashAttention actually fuse, and why does it help?"
• Resume bullet: Built a from-scratch GPU performance model — kernel-level roofline, SM occupancy with binding-resource analysis, memory-coalescing efficiency, tiled-GEMM HBM traffic, and operator-fusion traffic ratios — to predict and triage CUDA kernel performance without hardware, mirroring Nsight Compute's roofline and occupancy analyses.

Next: Phase 19 — Capstone: Production Multimodal Agentic Serving Platform.