Chapter 1 gave you a formula, time ≈ max(FLOPs / peak, bytes / bandwidth) + overhead. This chapter turns it into a picture you can draw for any kernel in thirty seconds: the roofline. One axis is how much arithmetic a kernel does per byte it moves, the other is how fast it runs, and the hardware draws a roof over both that no kernel can pass. Where a kernel sits under that roof tells you what limits it and which optimizations can possibly help. You will learn to count FLOPs and bytes honestly, see why tiling and batching both work by moving a point to the right, read the several roofs a real GPU has, recognize the regime the chart does not show (overhead), and diagnose a kernel that is running “at the roof” while wasting most of its time. Lab 3 traces your own card’s roofline, one point at a time.
A kernel reads and writes some number of bytes and performs some number of floating-point operations. Divide the two and you get its arithmetic intensity I, in FLOPs per byte. The hardware offers two rates: a peak FLOP rate P and a memory bandwidth β. If the kernel moved no bytes it could run at P. If it did no math it would be limited to moving bytes at β, which at intensity I means doing I × β FLOPs per second. It gets whichever is smaller:
attainable FLOP/s = min( P, I × β )
Plot that on log-log axes, intensity across and FLOP/s up, and you get the shape that gave the model its name: a slope of bandwidth rising from the left, meeting a flat roof of peak compute. Williams, Waterman and Patterson introduced it in 2009 for multicore CPUs. It became the standard way to reason about GPU kernels because it compresses the whole of chapter 1 into one glance.
Three things about the log-log axes make the chart useful and not just pretty:
10 × 672 GB/s = 6.7 TFLOP/s, about 22% of 30.9. It is 4.6× to the left of the ridge, so it is memory-bound: no amount of arithmetic optimization gets it past 6.7. If it runs at 6 TFLOP/s it is already excellent, even though “19% of peak” sounds bad. This is the single most useful reframing the roofline gives you: percent of peak is the wrong question for memory-bound work; percent of the roof above your point is the right one.
The same formula works as a machine: a feed that can only move so many bytes a second, and a press that can only strike so many times. Turn the knob through the kernels of this chapter and watch which of the two waits.
The chart is only as good as the two numbers you feed it, and both are easy to get wrong. The conventions:
exp runs on the special-function units at a fraction of the FMA rate. Counting it as one FLOP is the usual convention for placing a point, but remember that a kernel heavy in exp has a lower private roof.There are two different byte counts, and the difference between them is the most important idea in this chapter.
dram__bytes in Nsight Compute).Achieved intensity can never be higher than algorithmic intensity, only equal or lower. So a memory-bound kernel has exactly two ways to get faster: move its bytes closer to the bandwidth roof (coalescing, bytes in flight: chapter 3), or close the gap between achieved and algorithmic intensity (reuse, fusion). If both are already done, the only remaining move is changing the math itself, for example storing the data in a smaller type.
| Operation (per element unless noted) | FLOPs | Bytes | Algorithmic intensity |
|---|---|---|---|
| vecadd, FP32 | 1 | 12 | 0.083 |
SAXPY y = a·x + y, FP32 | 2 | 12 | 0.17 |
| ReLU, FP16 | ~1 (a compare) | 4 | ~0.25 |
| softmax over a row, fused, FP32 | ~5 | 8 | ~0.6 |
| layernorm, fused, FP32 | ~8 | 8 | ~1 |
| matrix × vector (GEMV) M×K, FP16 | 2MK | ≈ 2MK | ≈ 1 |
| square GEMM N×N×N, FP32 | 2N³ | 12N² | N / 6 |
| square GEMM, FP16 | 2N³ | 6N² | N / 3 |
Look at the shape of that table. Elementwise and normalization kernels, which are most of the kernels in a transformer by count, sit between 0.1 and 1, two orders of magnitude left of the FP32 ridge. Only matrix multiplication has an intensity that grows with the problem size, which is why it is the one operation that can use tensor cores fully, and why the rest of the network is about moving bytes efficiently.
It saves 128 MB of DRAM traffic, about 0.2 ms at the practical bandwidth, and one kernel launch. The FLOPs are irrelevant: at 10 FLOPs per 4 bytes GELU has an intensity of 2.5, deep in memory-bound territory, so the unfused kernel spends essentially all its time moving bytes. Fused into the matmul’s epilogue, the same 10 FLOPs happen on values that are already in registers, for free. This is why every serious GEMM library offers “epilogue fusion”, and why compilers like torch.compile (Part V) spend most of their effort on fusion.
A 4096×4096×4096 FP32 matrix multiply does 2 × 4096³ ≈ 137 GFLOPs and must touch 3 matrices of 64 MB: an algorithmic intensity of 4096 / 6 ≈ 683, fifteen times past the ridge. On paper, it is compute-bound. Whether it runs compute-bound depends entirely on the achieved intensity, and that depends on how often each byte fetched from DRAM gets used before it is thrown away.
That is the roofline reading of Part II’s whole matmul progression, and of Simon Boehm’s famous blog post: the kernels that close the gap to cuBLAS are, one after another, ways to raise achieved intensity. Real kernels don’t reach 683 (tiles are finite, and the caches have their own limits), but getting past the ridge is enough, because past the ridge the remaining work is about the math units, not the memory.
Figure 3 plots each kernel at its attainable performance: the best it could possibly do at that intensity. Real kernels land somewhere below their point, and the vertical gap is the second dimension of optimization (coalescing, occupancy, bank conflicts, instruction mix). Horizontal moves change what is possible; vertical moves get closer to it. Chapter 5’s profiler tells you which one you need.
Chapter 1 argued that LLM decode for one user is memory-bound and that batching fixes it. The roofline shows exactly how. A linear layer applies a weight matrix W (K×N) to the activations of every sequence in the batch.
So batch size is a dial that slides the point to the right along the slope, and throughput rises in proportion, for free, until the point reaches the ridge of the FP16 tensor roof at 92. That is chapter 1’s “compute catches up at batch ≈ 92”, now as a picture. Prefill, where the model processes a whole prompt of thousands of tokens at once, behaves like a huge batch: it is a big GEMM and sits far right, compute-bound.
No. The weights of a linear layer are shared by every sequence, so a bigger batch reuses them more. The KV cache is private to each sequence: 64 sequences read 64 different caches, and each byte is still used once. Batching raises the total FLOPs and the total bytes equally, and the point stays at intensity ≈ 1. This is why decode attention remains memory-bound at any batch size, why long contexts hurt decode so much, and why so much serving research (paged KV caches, KV quantization, multi-query and grouped-query attention, which share keys and values across heads) is about shrinking or sharing those bytes. Chapters 13 and 15 pick this up.
Put the pieces together below. Every preset uses the algorithmic intensity worked out in this chapter; try the batch slider on the decode layer, then compare attention decode and prefill.
The “≈88% DRAM” column uses a typical achievable fraction of the spec bandwidth. Lab 0 measured your real one, and it replaces the 88% in every calculation from now on.
The basic chart has one slope and one roof. A real GPU has several of each, and a kernel is limited by the lowest one that applies to it.
The FP32 roof assumes every lane does an FMA every cycle. Plenty of kernels can’t, and each reason draws a lower flat line:
Likewise, 672 GB/s is a spec, not a promise. The practical DRAM ceiling, the one lab 0 measured, sits some 10–20% lower. Uncoalesced access lowers it further: a warp that uses 4 bytes of every 32-byte sector (chapter 3) turns 672 GB/s into 84 GB/s of useful data. Too few bytes in flight lowers it again (chapter 3, section 5).
Bytes don’t only cross the DRAM interface. They also flow from L2 to the SMs and from L1 or shared memory into registers, and each of those links has its own, much higher, bandwidth. The hierarchical roofline draws one slope per level and places the kernel once per level, each time using the bytes that cross that level. The same kernel has a different intensity at each level: every byte from DRAM passes through L2, but many L2 and L1 accesses never touch DRAM, so the intensities decrease as you go up the hierarchy.
The L1/shared slope follows from chapter 3: 32 banks × 4 bytes per cycle per SM, × 48 SMs × 2.5 GHz ≈ 15 TB/s. Nsight Compute can draw the hierarchical version for you (its “Hierarchical Roofline” sections, per precision), using hardware counters for the bytes at each level. That’s chapter 5.
The roofline has no time axis. It assumes the kernel runs long enough that fixed costs don’t matter, and for small kernels that assumption is false. Launching a kernel costs a few microseconds whatever it does, so below a certain size every kernel takes the same time, and its FLOP/s and GB/s collapse towards zero. The kernel hasn’t moved on the chart; it has left the chart.
You recognize the overhead regime by two symptoms: time that doesn’t scale with size (halve the work and the time barely changes), and gaps between kernels on a timeline. The fixes live outside the kernel: fewer, larger kernels (fusion), capturing sequences of launches in a CUDA graph so the CPU submits them all at once (lab 3 measures the difference), or batching more work into each launch. For LLM decode with a small model, this is often the bottleneck, as chapter 1’s GPT-2 estimate suggested.
One more effect hides between the roofline’s assumptions. A grid runs in waves: if 6 blocks fit on each of your 48 SMs, the GPU runs 288 blocks at a time. A grid of 300 blocks needs two waves, and the second one keeps only 12 blocks busy while the other 276 slots sit idle. The kernel does 4% more work than one wave and takes nearly twice as long. Grids that are large compared with one wave hide this; grids of a few waves, common in inference with small batches, show it. When a kernel runs at a strange fraction of its roof, count its waves.
A worked example, with made-up but realistic numbers. You are handed a model that implements layer normalization by hand, out of elementary tensor operations, on a 16,384 × 4,096 FP32 activation (67 million elements, 268 MB):
mean = x.mean(-1, keepdim=True) # kernel 1: read x
xc = x - mean # kernel 2: read x, write xc
var = (xc * xc).mean(-1, keepdim=True) # kernels 3-4: read xc, write xc², read xc²
y = xc * torch.rsqrt(var + 1e-5) # kernel 5: read xc, write y
out = y * gamma + beta # kernel 6: read y, write outThe profiler reports 4.5 ms in total, spread evenly over six kernels, and each kernel achieves about 590 GB/s of DRAM throughput: 88% of 672. Every kernel is sitting on its roof.
No, and the roofline shows why once you compare the two intensities. Count the achieved bytes: ten passes of 268 MB (read x, read x + write xc, read xc + write xc², read xc², read xc + write y, read y + write out), about 2.7 GB. The algorithmic bytes are one read of x and one write of out: 0.54 GB. Achieved intensity is five times lower than algorithmic. Each kernel is efficient at what it does, but the program does five times too much of it.
One read and one write, 0.54 GB, at the same 590 GB/s: about 0.9 ms, five times faster than the six kernels. It is still memory-bound, and that’s fine: it now sits on the roof at the algorithmic intensity, which is the best any implementation can do. The fused kernel handles one row per block: 4,096 floats is 16 KB, which fits in registers and shared memory, so the mean, the variance and the normalization all happen on data that is already on chip. This is exactly what PyTorch’s built-in layer_norm kernel does, and what torch.compile would generate from the code above. You will write it yourself in chapter 9.
1. Work out the algorithmic FLOPs and bytes, and place the ideal point. 2. Measure time and achieved bytes, and place the real point. 3. If the real point is left of the ideal one, you are moving extra bytes: look for reuse or fusion. 4. If it is below its roof, you are using the hardware badly: look at access patterns, occupancy, instruction mix (chapter 5). 5. If it is on the roof at the ideal intensity, stop: the only improvements left change the algorithm or the data type.
Counting an FMA as one FLOP halves your numbers. Using spec bandwidth makes a kernel at 85% of the practical ceiling look like it has 15% to gain. Measuring with small, cache-resident data (lab 0’s L2 rows) puts points above the DRAM slope, which is impossible and means you measured L2. Comparing across precisions with the wrong roof: an FP16 kernel on the CUDA cores is judged against 30.9, not 61.7.
Two programs, three experiments. As always under WSL2, work inside the Linux filesystem and close games and video before benchmarking.
cd ~/unibo-lessons/mlsys/labs/lab03-roofline
make # ptxas: dial<1> ... dial<2048>, the two mma_peak kernels
make test # every dial kernel vs a CPU reference, the tensor-core tile, the launch counters
make bench # A intensity dial · B tensor-core peak · C launch overhead| Exp. | What it does | Prediction to write down |
|---|---|---|
| A | Loads one float, does F FMAs, stores one float: intensity = F/4, from 0.25 to 512 (sections 1–2) | the F where GB/s stops being flat; the GFLOP/s plateau |
| B | Tensor-core loop with operands in registers, FP32 vs FP16 accumulators (section 5) | the two TFLOPS numbers, and their ratio |
| C | Empty kernels: back-to-back, launch + synchronize, and from a CUDA graph (section 6) | three times in microseconds |
Paste the PASTE: lines back (two from make test, three from make bench). Experiment A gives you your card’s measured ridge point; it replaces the 46 in this chapter.
The ridge is where F/4 equals peak / practical bandwidth: roughly 30,900 / 590 ≈ 52 FLOPs per byte, so F ≈ 210. Below that, from F = 1 up to about 128, every run moves the same 537 MB and should take the same time, with GB/s flat near your lab 0 ceiling and GFLOP/s doubling with each step. From F = 256 on, time should roughly double with each step and GFLOP/s flatten. At 2048 expect something near the FP32 peak. Don’t be surprised by more than 100% of 30.9: GeForce cards often boost above their rated clock when cool, and the rated boost clock is what the spec uses.
About 2×. NVIDIA’s spec sheet gives 61.7 TFLOPS for FP16 inputs with FP32 accumulation and 123.5 with FP16 accumulation: on GeForce cards, FP32 accumulation runs at half rate (datacenter parts don’t have this limit). The same halving applies to FP8. Most training and inference uses FP32 accumulation for accuracy, so 61.7 is the tensor roof that matters in practice, which is why this course’s charts use it. WMMA is a convenience API, so reaching 85–95% of either peak would already be good.
Three details keep experiment A honest. The array is 2 × 256 MB, far bigger than the 48 MB L2, so the bytes really come from DRAM. The FMA count is a template parameter, so the loop is fully unrolled into straight-line FMAs: the build proves it (the 2048-FMA kernel contains exactly 2,048 FFMA instructions in its SASS). And the FMAs are split over four independent chains per thread, so the compute plateau measures the FMA units, not the latency of one dependent chain (section 5’s ILP ceiling).
Make a single-chain version of the dial (CHAINS = 1) and rerun experiment A. Where does the compute plateau land now? You have just measured the “not enough independent work” ceiling from section 5, and the ratio to the four-chain plateau tells you roughly how many cycles an FMA takes to produce its result.
Attainable FLOP/s = min(P, I × β): P is the peak compute rate of the units the kernel uses, β the memory bandwidth, I the kernel’s arithmetic intensity in FLOPs per byte moved. The ridge point P/β separates memory-bound (left) from compute-bound (right).
Algorithmic intensity uses the bytes the computation must move (inputs once, outputs once); achieved intensity uses the bytes the kernel actually moves across DRAM, including re-reads, wasted sector bytes and round trips of intermediates. Achieved can only be equal or lower. The gap between them is what reuse and fusion close.
It raises achieved intensity: each element fetched from DRAM is used T times from shared memory (and more with register tiling), so the point moves right along the bandwidth slope, raising the attainable FLOP/s, until it passes the ridge and the compute roof becomes the limit.
Weights are shared by all sequences in a batch, so a batch of B uses each weight byte B times: intensity ≈ B. Each sequence’s KV cache is private, so batching multiplies FLOPs and bytes equally and attention intensity stays ≈ 1.
You can’t tell without its intensity. At intensity 10 its roof is 6.7 TFLOP/s (22% of peak), so 25% of peak would be impossible over DRAM (it must be hitting cache). At intensity 100 its roof is the full peak, and 25% means real inefficiency. Always compare with the roof above the point.
That the kernel moved fewer DRAM bytes than you counted: its data was served from L2 or L1 (for example, a benchmark small enough to be cache-resident). Either recount the bytes that really crossed DRAM, or place it on the hierarchical roofline against the level that served it.
No FMA (half the rate: 15.5 TFLOP/s), divergence (only the active fraction of lanes does work), and insufficient independent instructions (too few warps or too little ILP to cover FMA latency). Using the CUDA cores for matrix math instead of the tensor cores is a fourth, relative to the tensor roofs.
The overhead-bound regime: kernels so short that launch cost dominates. Symptoms: time that doesn’t scale with problem size, and gaps between kernels on a timeline. Fixes: fusion, CUDA graphs, larger launches.
It runs in two waves; the second has only 12 blocks and leaves almost the whole GPU idle, so the kernel takes nearly twice as long as a 288-block grid for 4% more work. Wave quantization matters when a grid is only a few waves long.