The level in one paragraph. You write the most obvious GPU matrix multiply there is, one thread per output, and it runs at about 1% of what the hardware can do. The whole level turns on one distinction: inside a warp, 32 threads that run in lockstep, which values are the same for all of them and which differ. Once you mark that in the code, the reason for the 1% is visible, and so is the fix: move the varying index to a different place, swap two letters, and get about 6x. Then a new problem shows up, and that is level 2.
Throughout this level: uniform = the same for all 32 threads of a warp (fixed) · varying = different for each thread (mobile).
C = A × B. Each cell of C is one row of A "dotted" with one column of B: multiply them element by element and add up.
C[i][j] = A[i][0]·B[0][j] + A[i][1]·B[1][j] + … + A[i][N-1]·B[N-1][j]
The loop variable k walks along row i of A and down column j of B at the same time. Tap cells of C:
Two facts to keep: neighbours in the same row of C share their row of A; neighbours in the same column of C share their column of B. Our matrices are 4096×4096 floats: 16.7 million outputs, 4,096 multiply-adds each, 137 billion operations with nothing depending on anything else. The best library for it, cuBLAS, does it in about 6 milliseconds on an RTX A6000.
The obvious GPU plan: 16.7 million workers (threads), one per cell of C. Every worker runs the same code; the only difference between them is their ID, a pair (x, y) with both values from 0 to 4095. The one design decision is the mapping: which cell each worker computes. The natural first choice is "x picks the row":
// every worker runs this; its ID is (x, y) row = x; // naive mapping: x picks the row of C col = y; acc = 0; for (k = 0; k < N; ++k) acc += A[row][k] * B[k][col]; // the dot product from section 1 C[row][col] = acc;
It is correct. Written for a real GPU, it is also a textbook case of how to waste one. Commit to a number before reading on.
The GPU does not run workers one at a time. It runs them in warps: 32 workers with neighbouring x and the same y, such as (0, y) to (31, y), executing the same line at the same moment. At tick 0 they all do their k = 0 step, at tick 1 their k = 1 step, and so on. When the line says A[row][k], the hardware issues one request for all 32 of them.
Only x differs inside a warp; y, k and N are shared. Push that through the code:
row = x; // varying: 32 different rows col = y; // uniform: one column for the whole warp for (k = 0; k < N; ++k) // uniform: everyone is at the same step acc += A[row][k] * B[k][col];
So at every tick the warp asks for A[0..31][k], 32 different rows at the same column, and for B[k][col], one single float.
Here is warp y = 0 for its first two ticks. Read each row of the table as one request:
| tick | worker 0 | worker 1 | worker 2 | … | the warp asks for |
|---|---|---|---|---|---|
k = 0 | A[0][0], B[0][0] | A[1][0], B[0][0] | A[2][0], B[0][0] | 32 rows of A · one float of B | |
k = 1 | A[0][1], B[1][0] | A[1][1], B[1][0] | A[2][1], B[1][0] | the same, one step on |
Down a column is one worker over time: the ordinary loop. Across a row is what memory sees. Whether that request is cheap depends on how memory is laid out.
Memory is one long line. A 4096×4096 matrix is stored row after row: A[0][0] … A[0][4095], then A[1][0] … So A[r][c] sits at position r × 4096 + c. Each float takes 4 bytes: the next c is 4 bytes away, the next r is 4,096 × 4 bytes = 16 KB away.
Memory is delivered in chunks. When a thread needs one float, the memory system sends the whole 32-byte chunk around it, a sector: 32 ÷ 4 = 8 neighbouring floats. You pay for all 8 even if you use one, like a database reading a whole page from disk to hand you one row. A warp's request costs as many sectors as its 32 addresses touch. (Two different 32s: 32 threads per warp, 32 bytes per sector. A coincidence.)
Put the two facts together and the cost of any request follows from where its varying part sits:
| the address | its 32 addresses are | cost |
|---|---|---|
all uniform: B[k][col] | the same float | broadcast: 1 sector feeds everyone |
varying in the last index: B[k][col] | 32 neighbours, 128 bytes | packed: 4 sectors, all used |
varying in an earlier index: A[row][k] | 16 KB apart | scattered: 32 sectors, 1 float of 8 used in each |
The naive kernel's A[row][k] is the worst row of the table: 1,024 bytes moved to deliver 128, at every tick.
The same thing computed for every thread: each cell is a thread (t0–t31) with the sector its float lands in (s0, s1 …); the bottom row is what memory actually fetches.
x go? (think, then open)Into the last index of something. Let x pick the column and y the row: the warp now owns a slice of a row of C instead of a slice of a column.
row = y; col = x; acc += A[row][k] * B[k][col]; // broadcast · packed
A becomes all uniform, a broadcast; B becomes varying in the last index, packed. Per tick the warp goes from 33 sectors (32 + 1) to 5 (1 + 4): a 6.6x cut in memory traffic. Before you look at the measurement, guess.
The memory throughput Boehm measured went from 15 GB/s to 110 GB/s. The rule, neighbouring threads touch neighbouring addresses, is called coalescing, and it comes back in every level of this campaign.
"But the swapped kernel still jumps 16 KB through B at every tick." True, and so did the naive one: over time, both kernels do exactly the same thing. The swap only changes which workers stand next to each other.
| over time (one worker, tick after tick) | across the warp (one tick) | |
|---|---|---|
| naive | A: +4 bytes per tick · B: +16 KB per tick | A scattered (32) · B broadcast (1) |
| swapped | A: +4 bytes per tick · B: +16 KB per tick | A broadcast (1) · B packed (4) |
If your instinct says "each thread should walk memory in order", that is the CPU rule: one thread streaming sequentially is what CPU caches love, and the naive kernel does exactly that along its row of A. On a GPU the rule is transposed: what must be in order is the 32 threads side by side at one moment.
B were stored column by column, so B[k][c] sits at c × 4096 + k. Is the swapped kernel's read of B still packed? (think, then open)No. Now the first index is the contiguous one, and the varying col is the 16 KB-stride index: scattered again. You would want x to drive k, but k is uniform, the loop variable, not something a worker owns. The right mapping depends on the layout; this is why libraries such as cuBLAS take a "transpose" flag, and one reason real kernels stage data in shared memory (level 2).
A sector is 8 floats, so why not have each thread consume all 8, reading A[row][k] to A[row][k+7] in one go? It works: the warp still touches 32 sectors, but every byte of each is used, and the loop needs 8 times fewer requests. This is vectorised loading, a later rung of this ladder (level 5, where Boehm goes from 68.7% to 78.4%). The widest load one thread can do is float4, 16 bytes, so a thread takes two loads per sector. The swap fills a sector with 8 threads reading one float each, which the hardware merges for free, so it comes first; vectorised loads fill it from one thread, and stack on top later.
16.7 million workers are too many for the hardware to treat as one group, so a real kernel cuts them into blocks of 32 × 32 = 1,024 threads, arranged in a grid of 128 × 128 blocks. Each thread gets its block's position in the grid, blockIdx.x and blockIdx.y (0 to 127), and its own position inside the block, threadIdx.x and threadIdx.y (0 to 31). A worker ID is rebuilt the way you rebuild an offset from a page: page × page size + position on the page.
A warp always lives inside one block, as 32 threads with neighbouring threadIdx.x and the same threadIdx.y. So blockIdx, blockDim and threadIdx.y are uniform, and only threadIdx.x is varying: the same picture as before.
// naive: threadIdx.x picks the row extern "C" __global__ void matmul_naive(const float* A, const float* B, float* C, int N) { int row = blockIdx.x * blockDim.x + threadIdx.x; // varying int col = blockIdx.y * blockDim.y + threadIdx.y; // uniform if (row < N && col < N) { // edge guard: idle threads past the matrix float acc = 0.0f; for (int k = 0; k < N; ++k) acc += A[row * N + k] * B[k * N + col]; // A[row][k], B[k][col], stored flat C[row * N + col] = acc; } } // swapped: threadIdx.x picks the column int col = blockIdx.x * blockDim.x + threadIdx.x; // varying int row = blockIdx.y * blockDim.y + threadIdx.y; // uniform
Two details. The matrices arrive as flat arrays, so A[row * N + k] is exactly A[row][k] in the "one long line" layout of section 4. And which block axis you pair with rows or columns does not matter for a warp: all 32 threads share their block, so only the threadIdx choice decides uniform or varying.
Choose any block and thread. The arithmetic shows which terms are uniform and varying; the patch shows the thread's output in red and the rest of its warp in blue. Switch the mapping and the warp's strip turns from vertical to horizontal.
The swapped kernel sits at 8.5% of cuBLAS. Nothing is wasted per request any more, so what is left? Every output still reads a full row of A and a full column of B from global memory, and each element of A is read by all 4,096 workers that compute its row of C. You are fetching the same hot rows four thousand times. Any backend developer knows the next move. That is level 2.
Three new kernels. For each access, mark every index uniform or varying and predict what one warp request costs. This is the whole skill of level 1, applied where the page has not done it for you.
The notebook runs the naive kernel, gives you the same kernel to fix, times both against cuBLAS on your GPU and prints a score line for the ladder below. Colab's free T4 is enough.
pip install cupy-cuda12x)GPU programmers call the distinction uniform vs varying (per-lane) values, and recent NVIDIA GPUs even keep uniform values in separate registers, stored once per warp. Nsight Compute, NVIDIA's kernel profiler, reports what the lab counts: in its memory workload analysis the naive kernel shows many more sectors per request than the coalesced one. The full climb this campaign follows, from 1.3% to 93.7% of cuBLAS, is Simon Boehm's How to Optimize a CUDA Matmul Kernel for cuBLAS-like Performance (RTX A6000, 4092×4092 FP32).