Heads up: posts on this site are drafted by Claude and fact-checked by Codex. Both can still get things wrong — read with care and verify anything load-bearing before relying on it.
why → how

Why GPU kernels are still hand-tuned

A modern GPU can do tens of teraflops of matrix math. A naive, correct implementation of the same math leaves most of that on the floor. Here's why moving the bytes — not doing the FLOPs — is the actual job.

AI & ML intermediate Apr 29, 2026 · updated Aug 25, 2026 · 13 min read
On this page

The picture version

Six pictures for a reader who has never written a line of GPU code. The prose below fills in the seams the pictures skip.

1 · The problem

The chip is not broken. The maths is not wrong. You get 5% of the box number.

what the box says 100% hundreds of trillions of multiplications a second what you measure ~5% the textbook matrix multiply, written correctly Nothing is wrong with the answer. Everything is wrong with where the numbers were when they were needed.
Three nested loops, an operation older than computers, run on a modern data-centre chip. Correct, and leaving almost all the machine idle — that same matmul is the running example, and fixing it is the whole reason a profession exists around writing these by hand.

2 · The ratio that explains it

You must do about 295 operations on every byte you fetch. The naive loop does half of one.

989 trillion operations/s ÷ 3.35 TB/s ≈ 295 operations per byte 295 below here you are waiting for bytes above here the silicon is actually busy the naive matmul, at about half an operation per byte operations done per byte fetched → speed
Divide what the chip can compute by what its memory can deliver and you get a break-even point. Below it you are bandwidth-limited and the arithmetic units idle; above it you are actually using the machine. The naive loop sits two orders of magnitude below the line.

3 · Fix 1

Use each byte many times before you let go of it.

before: fetch, use once, drop main memory one multiply two fetches per multiplication, for every multiplication in the job half an operation per byte after: fetch a block, use it exhaustively main memory once the fast on-chip shelf every pair inside the block gets multiplied before anything is fetched again now it scales with the block size — over the line but the best block size depends on the matrix shape, the number format, and the chip — there is no single right answer
Loading a block once and doing every multiplication inside it before fetching again makes the operations-per-byte ratio scale with the block. Which converts “write a fast matmul” into “write a family of them and pick well” — and that conversion is most of a kernel engineer’s job.

4 · Fix 2

Then someone hands you an operation that fights back.

why blocking works on a matmul the partial results just add up blocks can be done in any order, independently why it fails on attention to finish one entry you need a total over the entire row so no block can be finished on its own The fix: rewrite the step to carry running totals as blocks arrive. the enormous intermediate table then never has to be written to main memory at all — same answer, far fewer bytes moved
Blocking is only easy when the partial results simply add. Attention’s normalisation step needs a whole row before any entry can be finished — and rewriting it to carry running totals is what let the same trick apply, which is the single most famous kernel in the field.

5 · And then the ground moves

Each new chip generation adds arithmetic faster than it adds bandwidth.

one generation ago the next the one you just bought the bar you must clear keeps moving right arithmetic has grown faster than bandwidth, generation after generation your kernel, written to clear the old bar — now below the new one, unchanged operations per byte → Which is why this is a permanent job, not a one-time cost.
Compute throughput has grown faster than memory bandwidth generation over generation, so the bar for “busy” keeps rising. A kernel that cleared it last year can fall below it without a line changing — and each generation also adds new hardware whose whole purpose is to hide a delay the last generation’s kernels were carefully built around.

6 · Keep this card

The whole thing on one index card.

a good kernel = the correct maths + bytes in the right place at the right moment + and “right” depends on shape, format, chip which is why you ship hundreds of kernels and a rule for choosing if good placement were a fixed property of an algorithm, you would solve it once and ship a library forever
Picture to keep: a kitchen where the chef chops faster than anyone can carry ingredients from the walk-in freezer. Hand-tuning is not sharpening the knife — it is stacking the right ingredients on the counter so the chef never stops. Where it breaks: a real chef can walk to the freezer while something simmers, and a naive kernel mostly can’t.

Why it exists

You replaced a five-year-old laptop with one whose spec sheet is better in every line — more cores, faster clock, more memory bandwidth — and the project you build every day takes almost as long as it did before. You check that it’s really using the new hardware. It is. Somewhere between the number on the box and the work getting done, most of the improvement went missing.

On a GPU that gap is measurable, and it’s enormous. Here’s the running example for this post: the textbook matrix multiply — three nested loops, the operation that’s been in linear algebra books since the 1800s — written correctly and run on a modern data-center GPU advertising hundreds of teraflops. You measure. You get a small fraction of the advertised number. Maybe 5%, maybe 10%. The GPU is not broken. The math is not wrong. The kernel is wrong — in a sense that has nothing to do with correctness and everything to do with where the bytes were when you needed them. We’ll fix that same three-loop matmul step by step, and the fixes turn out to be the entire reason a specialist industry exists around writing kernels by hand.

The deep reason is the same reason memory hierarchy exists at all: compute got faster much faster than main memory got faster, and the gap is huge. The job of a “good” kernel is mostly not to do the FLOPs faster — you can’t, the silicon picks the rate. The job is to keep the FLOP units fed, by carefully arranging which bytes live in which level of the memory hierarchy at which moment. That arranging is what hand-tuning is.

Why it matters now

Every LLM you use runs on top of a stack of these hand-tuned kernels. The canonical example is FlashAttention (Dao, Fu, Ermon, Rudra, Ré, May 2022, arXiv:2205.14135). It computes the exact same attention output as the textbook formulation — not an approximation — and reports a 15% end-to-end speedup on BERT-large at sequence length 512, a 3x speedup on GPT-2 at sequence length 1K, and 2.4x on long-range arena at 1K–4K. The mathematics didn’t change. Only where the intermediate values lived changed.

That pattern repeats everywhere. cuBLAS ships many GEMM kernels and uses runtime heuristics (via cuBLASLt) to pick one based on the matrix shapes, GPU model, and data type. Compiler stacks like Triton exist so people other than NVIDIA’s own kernel team can write kernels that land in the same ballpark. When a new GPU generation ships (Ampere → Hopper → Blackwell), the top-end kernels often need retuning and sometimes substantial rewrites, because the hierarchy changed shape.

This is also why “just port it to the GPU” rarely gets you the headline number. The headline number assumes someone hand-tuned the kernel.

The short answer

good kernel = correct math + data placement that keeps the FLOP units fed

Picture to keep: a kitchen where the chef chops faster than anyone can carry ingredients from the walk-in freezer. Hand-tuning is not sharpening the knife — it’s stacking the right ingredients on the counter so the chef never stops. (Where the picture breaks: a real chef can walk to the freezer while something simmers. A naive GPU kernel mostly can’t overlap the two, which is why the fetching time shows up as idle time.)

A GPU’s compute units run far faster than its main memory can deliver bytes. A naive kernel reads each operand from main memory once per use, which starves the compute. A hand-tuned kernel rearranges the work so each byte loaded from main memory gets used many times before being evicted — by tiling, fusing adjacent operations, and carefully choosing what lives in registers vs. on-chip cache vs. main memory. Same FLOPs, very different wall-clock.

How it works

The whole story is captured by one ratio.

The bandwidth wall, in one number

A modern GPU has a memory hierarchy roughly like this, fastest and smallest at the top (SRAM here means the on-chip scratchpad each streaming multiprocessor owns):

registers           : tiny, per-thread, ~instant
shared memory/SRAM  : ~tens to hundreds of KB per SM, very fast (on-chip)
L2 cache            : ~tens of MB, fast
HBM (main memory)   : tens of GB, slow (relatively)

HBM is fast in absolute terms — an H100 SXM has roughly 3.35 TB/s of HBM bandwidth, per public NVIDIA specs. The problem is that the compute side is even faster. The same spec page puts H100 SXM’s dense BF16 tensor-core throughput at 989 TFLOPS (the headline sparse figure is double that, but it requires weights pruned to the 2:4 structured-sparsity pattern, which dense LLM matmuls are not). Divide dense compute by bandwidth and you get the arithmetic intensity break-even point: about 295 FLOPs per byte loaded. Below that ridge, you’re memory-bound and the tensor cores are idle waiting for bytes. Above it, you’re compute-bound and the silicon is actually busy.

That number — “you must do hundreds of FLOPs per byte you load” — is the whole game. The roofline model just draws this as a graph: there’s a sloped ceiling on the left (bandwidth-limited) and a flat ceiling on the right (compute-limited), and your kernel sits somewhere underneath. Naive kernels live on the sloped part. Good kernels claw their way to the flat part.

Attempt 1: the three-loop matmul, and exactly how far it misses

Multiply two N×N matrices the textbook way. Total FLOPs: 2N³. If every multiply-add reloads both inputs from HBM, the ratio is constant — and miserably small. About 0.25 FLOP/byte for fp32, 0.5 FLOP/byte for fp16 or bf16. Two orders of magnitude below the 295-FLOP/byte ridge. Memory-bound isn’t even the right word; you’d be HBM-bound to the point of caricature. That’s where the 5%-of-peak number in the opening comes from — not from slow arithmetic, but from arithmetic that spends its life waiting.

The fix: reuse each byte before you drop it. Load a B×B block of A and a B×B block of B into shared memory once. Do B³ multiply-adds against those blocks. Now arithmetic intensity scales with B, the tile size. Make B big enough — limited by how much shared memory you have — and you cross the roofline ridge. cuBLAS does this. So does every reasonable matmul library. siboehm’s CUDA-MMM worklog walks through twelve progressively more careful versions and lands at ~94% of cuBLAS for square matrices — “within 95% on a good day,” in his words (siboehm.com/articles/22/CUDA-MMM).

Why that fix doesn’t finish the job: the right tile size depends on the matrix shape, the GPU model, the data type, and how much shared memory the surrounding kernels are using. There is no value of B that is correct everywhere, so there is no single best kernel — which is why cuBLAS is hundreds of kernels plus a heuristic to choose among them, rather than one good one. The first fix converted “write a fast matmul” into “write a family of matmuls and pick well,” and that conversion is most of what a kernel engineer’s job actually is.

Attempt 2: the same trick on something that resists it

Tiling a matmul is easy because a matmul is associative across blocks — the partial results just add. The interesting question is what happens when the operation you want to tile isn’t. Attention is the famous case, and it’s worth walking through because the answer generalizes.

Standard attention does softmax(QK^T / √d) V. Implemented naively: compute the N×N matrix QK^T, write it to HBM, read it back to apply softmax, write it again, read it again to multiply by V. The N×N intermediate doesn’t fit anywhere except HBM, and you traverse it multiple times.

FlashAttention’s insight is that you don’t need the whole N×N matrix at any one moment. You can tile Q, K, V into blocks, hold a block in SRAM, and compute a chunk of the output incrementally — provided you can do softmax incrementally too. That last bit is the online softmax trick (which predates FlashAttention; the contribution was wiring it into a tiled GPU kernel). Result: the N×N matrix never gets written to HBM at all. The kernel still moves bytes — there are just far fewer of them to move, which is the only thing that was ever in the way.

Same output. Same numerical answer (modulo non-associative floating-point reordering, which is a real but small caveat). Several times faster.

Attempt 3: ship it — and then the hardware changes underneath you

Each GPU generation moves the cliff. Hopper (H100) added the TMA, which extends Ampere’s async-copy model with a hardware unit that can move whole tiles in the background, freeing registers and simplifying the producer/consumer pipelines kernels rely on. Ampere kernels could overlap data movement with compute too — Hopper just makes it cheaper, easier, and faster to express. FlashAttention-3 (Shah, Bikshandi, Zhang, Thakkar, Ramani, Dao, July 2024) is largely a rewrite of FlashAttention-2 to exploit TMA and warp-specialization, and reports lifting H100 attention utilization from FlashAttention-2’s roughly 35% to around 75% in FP16 (PyTorch blog; the later NeurIPS 2024 version of the paper reports a higher BF16 figure, so check which version a number came from before quoting it).

The kernel didn’t get smarter about the math. It got smarter about the machine — and this is the part that makes hand-tuning a permanent job rather than a one-time cost. Every generation moves the ridge, changes the capacities, and adds a unit whose whole purpose is to hide a latency the previous generation’s kernels were carefully structured around.

The seam: this is genuinely hard, and getting harder

A few honest caveats.

You started with good kernel = correct math + data placement that keeps the FLOP units fed. After three attempts on that three-loop matmul, what did the post add? — + and the right placement is a function of the shape, the dtype, and the chip you're on. That last clause is the whole answer to the title. If good placement were a fixed property of an algorithm, you’d solve it once and ship a library forever. It isn’t, so you ship hundreds of kernels, a heuristic, and a team that rewrites them when the next generation lands.

Check yourself

Before you go — someone ports a workload to a new GPU with 2× the tensor-core throughput and the same memory bandwidth as the old one, and reports that their kernel got no faster at all. Is that a bug?

Answer

Almost certainly not — it’s the roofline doing exactly what it says. Doubling compute while holding bandwidth fixed doubles the break-even arithmetic intensity: a kernel that used to sit right at the ridge is now on the memory-bound side of it, and memory-bound kernels don’t care how fast the tensor cores are. The speedup was available only to code above the new ridge. This is the mechanism behind the post’s claim that each generation makes ignoring the memory hierarchy more costly, not less.

And a design one: you have a chain of small elementwise operations — add a bias, apply an activation, scale the result — each one a separate kernel. Nobody’s arithmetic is expensive. Where does the time go, and what’s the fix?

Answer

Every one of those kernels reads the whole tensor from HBM and writes the whole tensor back, to do roughly one FLOP per element. Read a value, write a value, do one operation: that’s a fraction of a FLOP per byte, hundreds of times under the ridge no matter how you count the exact accounting — you are paying pure bandwidth for almost no work, three times over. The fix is fusion: do all three in one kernel while the values are still in registers, so the tensor crosses the HBM boundary once instead of six times. This is the same move as FlashAttention, applied to a boring case; it’s also why fused elementwise chains are among the first things a compiler like Triton or torch.compile goes after.

Going deeper