gpu-kernels-explained
← /learn · 02

The roofline

Arithmetic intensity decides whether a kernel waits on memory or on the ALUs; the ridge point is where they meet.

Loading the animation…

Concept

Chapter 1 found each kernel's slowest level. The roofline (Williams, Waterman and Patterson, 2009) turns the most common case, HBM against the ALUs, into one picture. Put a kernel's arithmetic intensity II, the flops it does per byte it moves from HBM, on the horizontal axis. Then the best it can do is either limited by the bandwidth (the sloped line: each byte delivered allows II flops) or by the peak (the flat roof). The two meet at the ridge point: 12.5 flop/byte for the A100's FP32 units, 20 flop/byte for the H100's.

The animation steps a 4096 × 4096 × 4096 FP32 matrix multiply through tile sizes from 1 to 128. Every block of b×bb \times b outputs loads its bb rows of A and bb columns of B once and reuses each value bb times, so the intensity grows with the tile: about b/4b/4 flops per byte. With 1×1 tiles the kernel manages 389 GFLOP/s; with 32×32 tiles, 12.4 TFLOP/s, still on the slope; at 64×64 it passes the A100's ridge and meets the FP32 roof. On the H100, whose ridge sits further right, it takes 128×128.

The dashed roof is the BF16 tensor cores, about 16 (A100) or 15 (H100) times higher than FP32. Its ridge is at 201 flop/byte on the A100: tensor-core kernels need far more reuse before they stop waiting on memory. That is why the GEMM chapter's last step, and FlashAttention, are about keeping tiles on chip.

Slide the intensity yourself: below the ridge, every doubling of II doubles the attainable speed; above it, nothing changes.

Maths

For peak PP and bandwidth WW, the attainable rate at intensity II is

P(I)=min⁡(P, WI),I∗=PW.P(I) = \min(P,\ W I), \qquad I^{*} = \frac{P}{W}.

Below the ridge I∗I^{*} the kernel is memory-bound and P(I)=WIP(I) = WI; above it, compute-bound and P(I)=PP(I) = P. For the A100, I∗=19.5×1012/1.555×1012≈I^{*} = 19.5 \times 10^{12} / 1.555 \times 10^{12} \approx 12.5.

For the square GEMM (M=N=K=nM = N = K = n) with b×bb \times b tiles and no L2 reuse,

F=2n3,BHBM=4(2n3b+n2),I=2n34(2n3/b+n2)=b4⋅11+b/(2n)≈b4.F = 2n^3, \qquad B_{\text{HBM}} = 4\left(\frac{2n^3}{b} + n^2\right), \qquad I = \frac{2n^3}{4\left(2n^3/b + n^2\right)} = \frac{b}{4}\cdot\frac{1}{1 + b/(2n)} \approx \frac{b}{4}.

The n2n^2 term is writing C once. Setting I=I∗I = I^{*} gives the smallest tile that can reach the roof: b≈4I∗b \approx 4 I^{*}, about 50 on the A100 and 80 on the H100, so 64 and 128 in powers of two.

Code

The sweep, cut from src/lib/gpu/model.ts:

  for (const b of SWEEP_TILES) {
    const t = b > 16 ? idiv(b, 16) : 1;
    const g = gemmTraffic(p, size, size, size, b, b, t, t);
    const ai = g.flops / g.bytes.hbm;
    const pt = rooflinePoint(p, ai, "fp32");

Where real kernels sit

KernelIntensity (flop per HBM byte)On the A100's FP32 roofline
Vector add, FP320.0833far down the slope
GEMM 4096³, 16×16 tiles3.99on the slope (and shared memory binds first, chapter 1)
GEMM 4096³, 128×128 tiles31.5on the roof

Batch-1 decoding in an LLM is a matrix-vector product: each weight is read once and used for one multiply-add, about 1 flop per byte at BF16, which is why decode is memory-bound (LLM Inference Explained, chapter 3).