One matrix multiplication. Seven CUDA kernels. A 70× speedup, documented like a road trip — with the wrong turns left in. By the end you'll know exactly why the naive kernel crawls at 1.3% of what the GPU can do, and how to claw your way to 94%.
GEMM = GEneral Matrix Multiply. The full BLAS ceremony is C = αAB + βC, but the α and β are seasoning. The meal is C = AB: take a matrix A (size M×K), a matrix B (size K×N), and produce C (size M×N), where every cell of C is a dot product — one row of A pointwise-multiplied with one column of B, then summed.
Why does it matter? Because a large language model is, to a first approximation, a machine for doing GEMMs. Attention projections: GEMMs. The feed-forward blocks: two big GEMMs. The output head: a giant GEMM. When you're waiting for a chatbot's next token, you're mostly waiting for matrix multiplications. Make GEMM 2× faster and inference gets roughly 2× cheaper. That's why NVIDIA employs small armies to tune it, and why it's the single best kernel to learn GPU optimization on: simple enough to hold in your head, deep enough to teach you the whole machine.
Here's the property that makes the next seven chapters possible: those 137 billion FLOPs are almost embarrassingly parallel. Every cell of C can be computed independently — no cell needs another cell's answer. 16.7 million independent little jobs (for 4096×4096). That's the shape of workload GPUs were born for.
And here's the property that makes the next seven chapters necessary: the inputs are only A + B + C ≈ 201 MB, but the arithmetic touches those same numbers over and over — each element of A is needed by 4096 different output cells. GEMM is a data-reuse game. Win the reuse game, you win everything. Spoiler: the naive kernel loses the reuse game so badly it re-requests about half a terabyte from memory.
Before writing kernel one, meet the cast. A GPU like the RTX A6000 is a factory with 84 identical workshops called SMs (Streaming Multiprocessors). You don't hire individual workers; you submit a grid of thread blocks, and each block gets assigned to a workshop. Inside a block are up to 1024 threads — the actual workers.
But the threads don't move freely. The hardware bundles them into groups of 32 called warps, and a warp moves in lockstep: 32 threads, one instruction, all together. Think of a warp as a 32-segment caterpillar. It's the single most important creature in this scrapbook, so let's give it a name.
Threads need data, and data lives at very different distances. The numbers below are the whole plot of this scrapbook — every optimization we make is just moving work from a far shelf to a near one:
| Storage | Lives | Latency-ish | Size-ish | Kitchen analogy |
|---|---|---|---|---|
| Registers | per thread | ~1 cycle | 255 per thread | your own two hands |
| Shared memory | per block (on-SM) | ~20–30 cycles | up to 100 KB/SM | the counter you're standing at |
| L2 cache | whole GPU | ~200 cycles | 6 MB | the pantry down the hall |
| Global memory (DRAM) | off-chip | ~400–600 cycles | 48 GB @ 768 GB/s | the warehouse across town |
The A6000 can do ~38.7 TFLOP/s of fp32 math but can only pull ~768 GB/s from DRAM. Divide those and you get the machine's ridge point:
A whole GEMM has arithmetic intensity of roughly 2·4096³ FLOPs / 201 MB ≈ 680 FLOPs/byte — way above 50, so GEMM deserves to be compute-bound. Whether your kernel actually achieves that depends entirely on how much reuse you capture on-chip. That's the game. Let's play badly first.
From here on, every chapter is one kernel. Same problem every time — C = A·B, 4096×4096, fp32 — so the only thing changing is how well we use the machine. We measure GFLOP/s: the 137 billion FLOPs of work divided by wall-clock time. The reference bar is NVIDIA's own cuBLAS at 23,250 GFLOP/s.
A note on honesty, since this will become an article: the numbers on these pages are from Simon Boehm's classic kernel ladder, measured on an RTX A6000 — the canonical modern walk up this mountain. The kernels here follow the same route. A ready-to-run harness ships with this scrapbook so the numbers can be re-measured on any card, and the ratios between rungs are the durable lesson: they look similar on every modern NVIDIA GPU.
The obvious plan, and honestly a beautiful one: C has 16.7 million cells, GPUs have millions of threads, so give every thread one cell. Each thread walks its row of A and its column of B, multiply-accumulates K = 4096 times, writes one float. Done. This is the kernel every CUDA tutorial ends with, and the kernel every fast library begins by deleting.
__global__ void sgemm_naive(int M, int N, int K, const float *A, const float *B, float *C) { // one thread ↔ one output cell of C const uint x = blockIdx.x * blockDim.x + threadIdx.x; // row of C const uint y = blockIdx.y * blockDim.y + threadIdx.y; // col of C if (x < M && y < N) { float acc = 0.0f; for (int k = 0; k < K; ++k) acc += A[x * K + k] * B[k * N + y]; // dot product C[x * N + y] = acc; } } // launch: dim3 block(32, 32); dim3 grid(M/32, N/32);
Each thread does 2K = 8192 FLOPs but requests 2K·4 = 32 KB from global memory. That's the reuse game lost 0–100:
The L2 cache heroically absorbs a lot of that (which is why we're 75× slower than cuBLAS instead of 2700×), but the arithmetic units still spend almost all their time waiting. And it's actually worse than that — the way these threads ask for memory offends Warpy personally. Which brings us to the fastest 6× you will ever earn.
Remember Warpy's warning: memory requests happen per-warp, not per-thread. When 32 threads of a warp load a float each, the hardware looks at the 32 addresses and groups them into as few 32-byte transactions as possible. If the addresses are consecutive, all 32 floats (128 bytes) arrive in a handful of transactions — one bus, everyone rides together. If they're scattered, it's up to 32 separate trips carrying mostly air.
Now look at kernel 1's indexing. Threads that are neighbors in threadIdx.x got consecutive rows — their addresses in A and C are a full 16 KB apart. Same warp, 32 wildly scattered addresses, worst case every single load. The fix is almost insultingly small: remap the indices so that consecutive threads own consecutive columns — neighboring addresses in memory.
// each 32-thread stripe of the block now walks ALONG a row of C, // so warp-neighbors touch neighboring memory. that's the whole trick. const uint x = blockIdx.x * 32 + (threadIdx.x / 32); // row (changes every 32 threads) const uint y = blockIdx.y * 32 + (threadIdx.x % 32); // col (changes every thread) ★ if (x < M && y < N) { float acc = 0.0f; for (int k = 0; k < K; ++k) acc += A[x * K + k] * B[k * N + y]; // B & C now coalesce beautifully C[x * N + y] = acc; }
Nothing about the math changed. Nothing about how much data we touch changed. We only changed the shape of each warp's requests, and throughput jumped from 309 → 1,987 GFLOP/s. This is the recurring theme of GPU work: the hardware rewards you for asking politely.
We're still going to the DRAM warehouse for every ingredient. But notice: the 1024 threads of one block all compute cells in the same 32×32 patch of C — and their dot products read the same 32 rows of A and 32 columns of B. A thousand workers making a thousand separate trips for overlapping shopping lists. The fix: tile the K dimension.
Chop A's rows and B's columns into 32-wide chunks. Each round, the block cooperatively hauls one 32×32 tile of A and one of B from DRAM into shared memory — the fast on-SM counter from chapter 1 — everyone grabbing one element (coalesced, obviously; we learned that lesson). Then __syncthreads(), and everybody computes their partial dot product from the counter at ~20-cycle latency instead of ~500. Slide the tiles along K, repeat, done.
__shared__ float As[32 * 32]; // the counter: 4 KB each __shared__ float Bs[32 * 32]; float acc = 0.0f; for (int tile = 0; tile < K; tile += 32) { // whole block hauls one tile of A and B in together (coalesced) As[row * 32 + col] = A[x * K + (tile + col)]; Bs[row * 32 + col] = B[(tile + row) * N + y]; __syncthreads(); // wait till the counter is stocked for (int k = 0; k < 32; ++k) // shop from the counter, not the warehouse acc += As[row * 32 + k] * Bs[k * 32 + col]; __syncthreads(); // don't restock while people are shopping } C[x * N + y] = acc;
And the payoff for all that engineering is… 1.5×. One point five. We built the beautiful tile choreography and barely moved.
Warpy's diagnosis: too many loads per FLOP. The cure is the same trick as chapter 3, one level down. We tiled DRAM→shared; now we tile shared→registers, the only memory that's actually free. Instead of one output cell, each thread now owns a little column of TM = 8 cells of C, kept in registers the whole time.
Watch what that does to the inner loop. For each k-step, the thread loads one value of B into a register — and reuses it against 8 values of A to update all 8 accumulators:
float acc[8] = {0.0f}; // 8 cells of C live in registers now for (int k = 0; k < BK; ++k) { float bTmp = Bs[k * BN + col]; // ONE shared-mem load… for (int t = 0; t < 8; ++t) acc[t] += As[(row*8 + t) * BK + k] * bTmp; // …amortized over 8 FMAs }
An asymmetry should be bugging you: B values get reused 8× but A values are still loaded fresh every FMA. If reusing along one axis is this good…
…then reuse along both. Each thread now owns an 8×8 patch of C — 64 accumulators living in registers. Per k-step it loads a length-8 sliver of A and a length-8 sliver of B into registers, then computes their outer product: all 64 pairwise products, each feeding its own accumulator. 16 loads buy 128 FLOPs.
This is the moment GEMM's deep structure surfaces: a matmul is just a sum of outer products — for each k, (column k of A) ⊗ (row k of B), all stacked up. Every serious GEMM implementation on earth, cuBLAS and CUTLASS included, is arranged around exactly this shape. You've now derived it from first principles by being annoyed at load counters.
float acc[8][8] = {{0.0f}}; // 64 registers of pure ambition float regA[8], regB[8]; for (int k = 0; k < BK; ++k) { for (int i = 0; i < 8; ++i) regA[i] = As[(tRow*8+i) * BK + k]; // sliver of A for (int j = 0; j < 8; ++j) regB[j] = Bs[k * BN + tCol*8+j]; // sliver of B for (int i = 0; i < 8; ++i) for (int j = 0; j < 8; ++j) acc[i][j] += regA[i] * regB[j]; // 64 FMAs, zero memory ops }
68.7% of cuBLAS. The remaining gap isn't about what we compute anymore — it's about the pipes and the paperwork. Time for two chapters of pure logistics.
Every load instruction has fixed overhead: issue it, track it, retire it. Loading one 4-byte float per instruction means paying that overhead 4× more often than necessary, because the hardware has 128-bit wide load instructions (LDG.128 from global, LDS.128 from shared) that move four floats in one go. In CUDA you unlock them by loading float4:
// four floats, one instruction. compiler emits LDG.128 / LDS.128 float4 tmp = reinterpret_cast<const float4*>(&A[x * K + k])[0]; // while loading, store A's tile TRANSPOSED into shared memory, // so the inner loop can later read slivers of A with LDS.128 too: As[(k+0) * BM + row] = tmp.x; As[(k+1) * BM + row] = tmp.y; As[(k+2) * BM + row] = tmp.z; As[(k+3) * BM + row] = tmp.w;
Two details hiding in that snippet, both classic senior-engineer moves. First, the transpose: storing A's tile column-major in shared memory costs nothing extra during the load phase, but lets the hot inner loop read A-slivers as contiguous 128-bit chunks. Arrange your data at write time so it's cheap at read time — you'll use this trick your whole career. Second, float4 requires 16-byte alignment; feed it a misaligned pointer and it's an instant trip to cudaErrorMisalignedAddress town.
Count our tiling levels: the block owns a big patch of C (tiled from DRAM into shared), and each thread owns an 8×8 patch (tiled from shared into registers). But there's a tier of the hardware we never gave its own tile: Warpy. Kernel 7 fixes the org chart — each warp gets an explicit warp tile, subdivided among its 32 threads.
Why does adding a middle manager make anything faster? Three unglamorous, very real reasons:
1 · Register-cache locality. The warp's fragments of A and B get reused across the warp tile in a tight, predictable pattern, so the compiler schedules loads and FMAs into a smooth pipeline instead of a traffic jam.
2 · Bank-conflict dodging. Shared memory is 32 banks wide; if two threads of a warp hit the same bank, they queue. Warp-aware layouts let you arrange accesses so all 32 segments hit 32 different banks.
3 · It's the shape of the future. A warp computing a small matrix patch together is exactly the contract of tensor-core instructions (mma, wgmma). Warptiling is the fp32 dress rehearsal for how all modern GEMM libraries are built.
One more trick rides along in the fast kernels: double buffering. While the math units chew tile t, the next tile t+1 is already being fetched into a second buffer — compute and delivery overlap, and nobody waits:
Seven kernels, a 70× speedup, and not one line of assembly. Here's the ladder in full — press the button and watch the race we just ran:
| # | Kernel | The one idea | GFLOP/s | ms/matmul | vs cuBLAS |
|---|---|---|---|---|---|
| 1 | Naive | 1 thread = 1 cell | 309 | 444.8 | 1.3% |
| 2 | Coalescing | warps read neighbors | 1,987 | 69.2 | 8.5% |
| 3 | Smem tiling | block shares tiles | 2,980 | 46.1 | 12.8% |
| 4 | 1D tiling | 8 cells per thread | 8,475 | 16.2 | 36.5% |
| 5 | 2D tiling | outer products, 64 cells | 15,972 | 8.6 | 68.7% |
| 6 | float4 | 128-bit loads | 18,237 | 7.5 | 78.4% |
| 7 | Warptiling | a tile per hardware tier | 21,779 | 6.3 | 93.7% |
| — | cuBLAS | several careers of tuning | 23,250 | 5.9 | 100% |
1. GEMM has enormous intrinsic data reuse (~680 FLOPs/byte available); a kernel is fast exactly in proportion to how much of that reuse it captures in fast memory. 2. So you tile at every level the hardware has — DRAM→shared (blocks), shared→registers (threads), and warps in between — and you make every remaining transfer wide, coalesced, and overlapped with compute. 3. Everything else is bookkeeping.
The honest fine print: the real ladder has two rungs we hand-waved — resolving shared-memory bank conflicts and autotuning tile sizes per GPU (that's how you squeeze out 84.8% before warptiling's 93.7%). And above all of this sits a different vehicle entirely: tensor cores — dedicated matmul hardware driven by mma/wgmma instructions, where the A6000 does not do 39 but ~77 TFLOP/s of tf32 (and ~155 of fp16), and libraries like CUTLASS, Triton and DeepGEMM live. The beautiful part: it's all still tiles — the same hierarchy you just learned, with the thread-tile FMA swapped for a warp-wide matrix instruction. You now speak the language; the tensor-core chapter is a vocabulary lesson.
This scrapbook ships with a benchmark harness — all seven kernels as runnable CUDA files with timing, cuBLAS verification, and a Makefile — so every number here can be re-measured on whatever GPU you get your hands on. The reference numbers and route follow Simon Boehm's worklog; for the visual-explanation genre this page lives in, see bbycroft's LLM visualization and Wattenberger's essays; for what lies above kernel 7, start with CUTLASS and the CUDA C++ Programming Guide.