Stop fetching the same numbers

1. The problem: hot rows

Look at who reads what in the coalesced kernel. Row i of C has 4,096 outputs, each computed by its own thread, and every one of those threads reads all of row i of A from global memory. So each float of A is fetched 4,096 times, once per thread that needs it. The same holds for each float of B, read once by every thread in its column of C.

Add it up: A and B are 64 MB each, but the kernel asks global memory for about 550 GB in total, the whole input read four thousand times over. In backend terms, every request re-queries the same hot rows instead of caching them.

A useful way to measure this is arithmetic intensity: floating-point operations done per byte fetched from global memory. Each multiply-add (2 operations) brings in two fresh floats (8 bytes): 0.25 operations per byte. A GPU computes far faster than data can arrive at that rate, so the kernel's speed is capped at intensity × memory bandwidth. Plotted against intensity, that cap is a sloped line that flattens out at the compute peak: the roofline.

Lab: where is the ceiling?

Keep "no tiling" and pick your GPU (Colab's free one is the T4). The red bar is the ceiling that binds; the others are the ones you would hit next.

One honest wrinkle: real untiled kernels run above the "no tiling" line, because the GPU's L2 cache quietly catches some of the reuse. The cache is doing this level's trick for you, partly and by accident. Doing it on purpose is the fix.

2. Where could a copy live?

A GPU has several kinds of memory, from huge and slow to tiny and fast. The latencies below are rough orders of magnitude; they vary by chip.

An SM (streaming multiprocessor) is one of the GPU's core clusters; a block always runs on a single SM. Shared memory is the interesting row: on-chip, fast, and shared by exactly the threads of one block. It is a cache you manage: nothing arrives there unless your code copies it in, and it disappears when the block finishes. You declare it inside the kernel:

__shared__ float As[32][32];   // one copy per block, visible to all its 1,024 threads

3. The plan, one tile at a time

for each step t = 0, 1, …, 127:
    every thread copies ONE float of A's tile and ONE float of B's tile into shared memory
    wait until the whole block has finished copying
    every thread does 32 multiply-adds, reading only shared memory
    wait until the whole block has finished computing
every thread writes its one output

Why one float per thread per tile? A tile has 32 × 32 = 1,024 floats and the block has 1,024 threads: the copy is shared out perfectly. Then each loaded float is used by 32 threads: A's tile element in row r by the 32 threads computing row r of the patch, B's element in column c by the 32 threads computing column c. Watch it:

Lab: walk the tiles

Tap a block of C, then step through t. The block moves along its strip of A (left to right) and down its strip of B (top to bottom), one tile each per step.

Back to the roofline lab above: choose "tile 32x32" and watch the DRAM ceiling move 32 times higher, from 0.25 to 8 operations per byte.

4. The kernel, line by line

#define TILE 32
extern "C" __global__ void matmul_tiled(const float* A, const float* B, float* C, int N) {
    __shared__ float As[TILE][TILE];                                  // (1)
    __shared__ float Bs[TILE][TILE];
    int tx = threadIdx.x, ty = threadIdx.y;
    int row = blockIdx.y * TILE + ty;                               // (2) uniform
    int col = blockIdx.x * TILE + tx;                               // (2) varying
    float acc = 0.0f;
    for (int t = 0; t < N; t += TILE) {                               // (3)
        As[ty][tx] = A[row][t + tx];                         // (4) packed
        Bs[ty][tx] = B[t + ty][col];                         // (4) packed
        __syncthreads();                                              // (5)
        for (int k = 0; k < TILE; ++k)
            acc += As[ty][k] * Bs[k][tx];                     // (6) broadcast · no conflict
        __syncthreads();                                              // (7)
    }
    if (row < N && col < N)
        C[row * N + col] = acc;                                        // (8)
}
// A[r][c] stands for A[r * N + c]; the notebook version also checks bounds on the loads.

1. Two tiles of shared memory, 4 KB each, one pair per block.
2. The output this thread owns, mapped the level-1 way: tx picks the column, so a warp owns a slice of a row of C.
3. The walk along the strips: t is the tile's starting column in A and starting row in B. It is uniform: the whole block is on the same step.
4. The copy. Thread (tx, ty) copies element (ty, tx) of each tile. Read through the lens: A[row][t + tx] is varying only in the last index, packed; B[t + ty][col] too. The loads are coalesced exactly as in level 1.
5. First barrier: nobody reads the tiles before all 1,024 threads have written their element.
6. The work, 32 multiply-adds from shared memory. As[ty][k] is all uniform: one value broadcast to the warp. Bs[k][tx] is varying in the last index. Section 6 shows why both are fast in shared memory too.
7. Second barrier: section 5 is all about it.
8. The write: one output per thread, packed as in level 1.

5. The barrier you already know

__syncthreads() is a barrier for the threads of one block: nobody passes until everybody arrives, the same idea as the barriers of concurrent programming. The first one is obvious. The second is not.

Lab: race two warps

A block shrunk to two warps; each writes half of the tile, and each computes on the whole tile. Step them in any order, switch barriers off, or let W0 run ahead as a fast warp would. With both barriers, no order you choose can produce a wrong step.

Two more rules from the same family. Every thread of the block must reach every barrier: a __syncthreads() inside an if that only some threads take can hang the block. And like any race, the missing-barrier bug shows up only sometimes, depending on timing, which is exactly what makes it dangerous.

6. Your lens, inside shared memory

Shared memory has its own version of sectors. It is split into 32 banks of 4 bytes: float i of a shared array lives in bank i mod 32. Each bank serves one address per step. So a warp's access is fast when its 32 addresses fall in 32 different banks, or are all the same address (a broadcast), and slow when several different addresses land in the same bank: a bank conflict, served one after another.

The level-1 rule carries over almost unchanged: all uniform is a broadcast; varying in the last index lands in 32 different banks; varying in an earlier index, with rows of 32 floats, puts every thread in the same bank.

Lab: count the bank steps

7. The twist

Global-memory traffic is down 32x. Commit to a number before you look.

Only about 1.5x, for 32x less traffic. Two reasons, and both matter for everything that follows:

• The L2 cache was already catching much of the reuse for the coalesced kernel, so the real saving was smaller than the 32x on paper.
• The bottleneck moved. Each multiply-add still reads two numbers, now from shared memory instead of global memory. Shared memory is fast, but it has a bandwidth too. Back in the roofline lab with "tile 32x32", on most GPUs the binding ceiling is now shared memory, not DRAM.

That is exactly the next move: more work per thread, with reused values in registers. It takes the climb from 12.8% to 36.5% and then 68.7% of cuBLAS. That is level 3.

Practice: mark it yourself

Three variations of this kernel. Mark each index uniform or varying and predict the cost: in sectors for global memory, in bank steps for shared memory. In all of them tx is varying and ty, t, k are uniform.

Side quest: tile it yourself

The notebook's level 2 cell has the skeleton with three TODOs: the two tile loads and the two barriers. It compiles as given but produces a wrong answer until you fill them in; the harness checks every result against cuBLAS.

• Open in Colab · download the notebook

What you earned