Matmul Ladder I: Give Each Thread A Column Strip

Lesson

This lab runs on the register-blocking harness: the grid is a 1-D flatten of 64x64 block tiles - gridX = (N/64)^2 - and your kernel recovers its tile as bx = blockIdx.x % (N/64), by = blockIdx.x / (N/64). Every test dimension is an exact multiple of 64, so there are no partial tiles to guard; the lesson is register reuse, not edge handling.

Why registers? Shared memory is fast, but a shared-memory read still costs an instruction issue slot and has latency an FMA does not. A one-output-per-thread inner loop does sum += As[row][k] * Bs[k][col]: two shared reads per FMA, so the SM spends more issue slots on loads than on math - and with only 256 threads for 4,096 outputs, the starter must sweep the tile 16 times, restaging the same shared tiles on every sweep. Thread coarsening flips the ratio: give each thread TM = 16 outputs (a column strip of the tile) held in an accumulator array that the compiler keeps in registers. In the k-loop, read one B value from shared into a register, then reuse it against 16 A values: one shared read amortized over 16 FMAs, and the tile is staged exactly once.

This is the correction to the occupancy chapter: registers per thread go UP (the occupancy calculator frowns), yet throughput goes up too, because each resident thread now hides latency with its own instruction-level parallelism. Volkov's classic GTC talk 'Better Performance at Lower Occupancy' made exactly this argument, and Boehm's matmul worklog measures this rung at roughly 3x over the one-output tiled kernel.

Where to go next - 2D register-blocking: the 1-D strip amortizes only B (each A value is still read from shared once per FMA). A 2-D outer-product tile gives each thread a 4x4 micro-tile - copy a 4-element A fragment and a 4-element B fragment into registers, multiply every combination, and each fragment value is reused 4 times, so shared reads per FMA halve again, with float4 (16-byte) staging loads on top. The catch, and why this lab stops at 1-D: on a laptop GPU the 2-D kernel is another ~2x, but on a datacenter GPU (L4-class) the SMs are already saturated by the 1-D kernel, so 2-D buys only ~15-20% - a live reminder that advanced micro-optimizations hit sharply diminishing returns once the machine is full, and that the honest way to know is to measure on the hardware you deploy on.

Distinguish loading a global tile into shared memory from reusing shared values across several outputs per thread. Inspect registers, occupancy and kernel time as a group.

What you will learn

• Use the register-blocking harness contract: 1-D grid of (N/64)^2 fixed 64x64 tiles, block shape yours to choose.
• Apply thread coarsening: one thread computes a TM-row column strip held in a register accumulator array.
• Amortize shared-memory reads: load a B value once per k-step and reuse it across all TM rows.
• Read the Volkov trade in the profiler: more registers per thread, and the kernel runs faster despite lower occupancy.

Starter code

    // Register-blocking harness contract: gridX = (N/64)^2 one-dimensional
    // blocks; this block owns the 64x64 tile of C at (baseRow, baseCol).
    // This starter keeps the one-output-per-thread habit: 256 threads cover
    // the 4,096-output tile in 16 sub-passes - and every sub-pass re-stages
    // the SAME shared tiles it already loaded the sub-pass before.
    const int TILE = 64;
    const int BK = 8;
    __shared__ float As[TILE][BK];
    __shared__ float Bs[BK][TILE];

    int tilesPerRow = N / TILE;
    int baseCol = (blockIdx.x % tilesPerRow) * TILE;
    int baseRow = (blockIdx.x / tilesPerRow) * TILE;
    int tid = threadIdx.x;

    for (int sub = 0; sub < 16; ++sub) {
        int idx = sub * 256 + tid;
        int row = idx / TILE;
        int col = idx % TILE;
        float sum = 0.0f;
        for (int t = 0; t < N; t += BK) {
            for (int i = tid; i < TILE * BK; i += 256) {
                As[i / BK][i % BK] = A[(baseRow + i / BK) * N + (t + i % BK)];
                Bs[i / TILE][i % TILE] = B[(t + i / TILE) * N + (baseCol + i % TILE)];
            }
            __syncthreads();
            for (int k = 0; k < BK; ++k) {
                sum += As[row][k] * Bs[k][col];
            }
            __syncthreads();
        }
        C[(baseRow + row) * N + (baseCol + col)] = sum;
    }