Roofline: what actually limits your kernel

Before you optimize a kernel you should be able to state, in one sentence, what is stopping it. There are only three candidate answers, and a single quantity decides between the first two. Getting this right early saves you from the classic failure mode: spending a week on tensor-core utilization for a kernel that was never going to be limited by arithmetic.

Arithmetic intensity is the number of floating-point operations a kernel performs per byte it moves between the GPU and its memory:

The machine has a matching number, the ridge point: peak arithmetic throughput divided by peak memory bandwidth. It is the intensity at which the two ceilings cross.

Below the ridge you are memory-bound: the arithmetic units are idle waiting for data, so you win by moving fewer bytes or by using more of the available bandwidth. Above it you are compute-bound: memory keeps up, so you win by doing the arithmetic more efficiently (better instructions, better tiles, lower precision).

Two properties of this plot are worth stating explicitly, because they are the whole reason the model is useful:

• Horizontal moves are optimizations that change what you touch — tiling, fusion, caching, quantizing weights. They move a kernel right, toward the ridge.
• Vertical moves are optimizations that change how well you use the ceiling you are already under — coalescing, bank-conflict removal, better instruction mix. They move a kernel up toward the roof it is already beneath, and they cannot move it above it.

A kernel at 4 FLOP/byte running at 90% of the bandwidth ceiling is done. There is nothing more to win without restructuring the algorithm.

The plot has a pedigree: it is the Roofline model of Williams, Waterman and Patterson (CACM 2009), drawn originally for multicore CPUs and the GPUs of that era. It has since been promoted into Nsight Compute, which will draw this exact chart for any kernel you point it at.

Horace He's taxonomy, from his post Making Deep Learning Go Brrrr From First Principles, is the practical version of the model: every slow kernel is slow for one of three reasons, and each has a different fix.

The third regime — overhead-bound — is invisible to the roofline model because it is not about the GPU at all. A kernel launch costs a few microseconds of CPU and driver work. A transformer layer with 20 small unfused operations, each taking 3 µs of GPU time and 5 µs of launch overhead, spends most of its wall clock in Python and the driver. This is the normal state of an unoptimized PyTorch model at small batch size, and no kernel tuning fixes it — CUDA graphs, torch.compile, and fusion do.

Here is the whole table for a transformer, worked out. Assume BF16 (2 bytes) unless stated, and count only HBM traffic — the intensity of a well-tiled kernel is set by what crosses the DRAM boundary, not by what crosses the SMEM boundary.

Read the last three rows together, because they are the central fact of LLM serving. Long-sequence or well-batched prefill is typically a compute problem; low-batch decode is typically a bandwidth problem, and they are the same mathematics with different shapes. A single-sequence decode step reads the entire weight matrix to do two FLOPs per element. Its arithmetic intensity is 1. The ridge point is 295. You are using 0.34% of the machine's arithmetic, and that is the optimal outcome for that shape.

For this weight-stationary shape, the main reuse lever is batching — sharing that weight read across more queries. Quantizing the weights is the other direct way to reduce the compulsory bytes.

This is why inference servers use continuous batching, and why FP8 weights improve both sides of the model: half the weight bytes help the memory-bound side, while twice the compute ceiling helps after the knee. Both effects are about 2× in this idealized model, not a compounded 4×; because FP8 doubles both arithmetic intensity and the ridge, the crossover batch stays roughly unchanged. An H200 — the same peak compute as H100 with 43% more bandwidth — is dramatically better at low-batch decode and can also help bandwidth- or capacity-limited training phases, although it does not raise the ceiling for compute-bound training kernels.

A naive matmul computes each output element independently: it reads a full row of A and a full column of B for every one of the MN outputs, so its A/B operand reads move 2MNK elements. Its intensity is 1 FLOP per operand element read — hopeless. The output write is omitted for this first reuse derivation and added back when counting total DRAM traffic.

Now give each thread block a T×T output tile. The block reads TK elements of A and KT of B, and there are (M/T)(N/T) blocks:

In intensity terms, for a BM × BN tile of depth BK and element size s:

Note what is missing from that expression: BK cancels out of this SM-inbound A/B operand model. It still affects shared-memory footprint, pipeline depth, edge handling, and how output traffic is amortized in a complete kernel. The harmonic-mean form is maximised for a fixed accumulator area when BM = BN, so square-ish tiles get more reuse from that fixed area. Real CUTLASS shapes also answer to tensor-instruction geometry, registers, SMEM, occupancy and problem shape; intensity is one constraint, not the entire explanation.

For memory-bound work, the optimization is always the same: touch DRAM fewer times. A chain of k elementwise operations, each launched as its own kernel, reads and writes the whole tensor k times. Fused into one kernel, it reads and writes once.

Concretely, for a Llama-style block at hidden size 8192, batch 8192 tokens, BF16 — one activation tensor is 8192×8192×2 = 134 MB:

The same arithmetic as an audit you can run — four materialized PyTorch operations on one side, one Triton kernel on the other:

# Four high-level eager operations. Profile your PyTorch/CUDA version:
# a native reduction may use more than one internal kernel.
h = torch.nn.functional.rms_norm(x, (x.shape[-1],), weight=None, eps=eps)
s = torch.nn.functional.silu(h)
t = s * g
y = t + x

# DRAM audit, bf16, one activation tensor = 8192×8192×2 B = 134 MB:
#   rmsnorm   read x, write h          268 MB
#   silu      read h, write s          268 MB
#   mul       read s and g, write t    402 MB
#   add       read t and x, write y    402 MB
#   10 passes = 1342 MB ⇒ floor 400 µs at 3.35 TB/s

@triton.jit
def rms_silu_gate_add(X, G, Y, stride, H, eps, BLOCK: tl.constexpr):
    row  = tl.program_id(0)
    cols = tl.arange(0, BLOCK)
    mask = cols < H
    x = tl.load(X + row*stride + cols, mask=mask, other=0.).to(tl.float32)
    g = tl.load(G + row*stride + cols, mask=mask, other=0.).to(tl.float32)
    inv = tl.math.rsqrt(tl.sum(x*x, axis=0) / H + eps)
    h  = x * inv
    tl.store(Y + row*stride + cols, (h * tl.sigmoid(h) * g + x).to(tl.bfloat16), mask=mask)

# 3 passes: read x and g, write y = 402 MB ⇒ floor 120 µs.

The same logic drives GEMM epilogue fusion: the accumulator is already in registers at the end of the mainloop, so applying a bias, an activation, a residual add or a quantization step there costs nothing but a few instructions. Writing it out and reading it back in another kernel costs two full passes over the output. Part 5 shows where the epilogue sits in a real kernel.

And it is the reason FlashAttention exists at all: standard attention writes the N×N score matrix to HBM and reads it back, then softmax writes and reads an equally large probability matrix — four one-way N×N transfers, or two round trips, about 17 GB for seq 8192 and 32 heads, none of it useful output. Part 7.

The ridge point is the single most useful summary of a GPU generation, and it has a surprising history: it roughly doubled from A100 to H100, meaning Hopper made it substantially harder to be compute-bound. Blackwell then held it steady for BF16 while pushing it far out for FP4.

Two consequences worth carrying:

1. Lower precision raises both ceilings. FP8 weights halve the bytes, so bandwidth-bound decode gets about 2×; FP8 tensor cores also double peak compute, so the compute-bound side gets about 2×. Arithmetic intensity and ridge both double, however, so in this idealized weight-stationary model the batch-size crossover stays near the same point rather than doubling.
2. The gap between generations is mostly not in the ridge. B200 gives you 2.3× the BF16 FLOPs and 2.3× the bandwidth of an H100. If your kernel was memory-bound on H100 it is still memory-bound on B200 — just 2.3× faster. Restructuring the algorithm is the only thing that changes which ceiling you are under.

The procedure, before writing any code:

1. Count the FLOPs. For a matmul, 2MNK. For attention, 4N²d per head forward, ×2.5 for backward.
2. Count the compulsory bytes — the inputs and outputs that must cross the DRAM boundary at least once. Ignore intermediates you intend to keep on chip.
3. Divide, and compare to the ridge. Below it, your target is bytes / bandwidth. Above it, your target is FLOPs / peak.
4. Now measure. If you are at 85% of the relevant ceiling, stop. If you are at 20%, the gap is a Parts 1–2 problem: coalescing, bank conflicts, occupancy, or an empty pipeline.

That last step is the discipline the model buys you. "This kernel takes 400 µs" is not actionable. "This fused kernel must move about 402 MB, so it cannot beat roughly 120 µs at peak H100 bandwidth, and it takes 400" is a bug report.

Next: the kernels that live at the bottom of the roofline — elementwise, reductions, softmax, RMSNorm — and how to actually reach the bandwidth ceiling rather than merely aiming at it.