Tiling and Blocking for Reuse

the short version

A naive matrix multiply reads each input element once per output element that needs it, which for an N-by-N problem means reading the same values N times from device memory. The arithmetic is unavoidable; the re-reading is not.

Tiling loads a small block of each input into on-chip scratchpad once, has every thread in the block reuse it many times, then moves on. The FLOP count is identical and the memory traffic drops by roughly the tile dimension. That is the entire difference between a naive kernel and a fast one.

why naive is slow

Consider the obvious kernel: one thread per output element, looping over the shared dimension.


__global__ void matmulNaive(const float *A, const float *B, float *C, int N)
{
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    if (row < N && col < N) {
        float acc = 0.0f;
        for (int k = 0; k < N; ++k)
            acc += A[row * N + k] * B[k * N + col];
        C[row * N + col] = acc;
    }
}

Each thread performs N multiply-adds and loads 2N values from device memory. That is one FLOP per four bytes before you even count the write, which puts the kernel firmly on the bandwidth-limited part of the roofline despite being, in principle, the most compute-heavy operation there is.

The waste is that the thread next door reads almost exactly the same row of A, and the thread in the next block reads it again later. Every value of A is pulled out of HBM N times.

the tiled version

The fix is cooperative: threads in a block load a tile of A and a tile of B into shared memory, synchronize, then every thread computes partial results from those tiles before moving to the next pair.


#define TILE 32

__global__ void matmulTiled(const float *A, const float *B, float *C, int N)
{
    __shared__ float As[TILE][TILE];
    __shared__ float Bs[TILE][TILE];

    int row = blockIdx.y * TILE + threadIdx.y;
    int col = blockIdx.x * TILE + threadIdx.x;
    float acc = 0.0f;

    for (int t = 0; t < N / TILE; ++t) {
        As[threadIdx.y][threadIdx.x] = A[row * N + t * TILE + threadIdx.x];
        Bs[threadIdx.y][threadIdx.x] = B[(t * TILE + threadIdx.y) * N + col];
        __syncthreads();

        for (int k = 0; k < TILE; ++k)
            acc += As[threadIdx.y][k] * Bs[k][threadIdx.x];
        __syncthreads();
    }
    C[row * N + col] = acc;
}

Each tile is loaded from device memory once and used TILE times by the block. Device traffic falls by a factor of TILE, arithmetic intensity rises by the same factor, and the kernel moves up the sloped roof toward the compute ceiling. The two barriers are both required: one so every thread sees the full tile before computing, one so nobody overwrites the tile while others are still reading it.

choosing the tile size

Bigger tiles mean more reuse and less traffic, which argues for the largest tile you can manage. Three things push back.

Constraint	Effect of a larger tile
Shared memory per unit	Two tiles must fit; larger tiles mean fewer resident blocks
Registers per thread	More accumulators per thread means fewer resident groups
Tail effects	Problem dimensions that are not a multiple of the tile waste work at the edges

So tile size trades reuse against occupancy, and the optimum depends on the architecture, the precision and the problem shape. There is no formula that beats measuring, which is exactly why the tuned libraries and template frameworks treat tile size as their central autotuned parameter.

the next level: register blocking

The version above still performs one multiply-add per shared memory read, and shared memory, while fast, is not free. The next step is to have each thread compute a small patch of output rather than a single element, holding several accumulators in registers.

A thread computing a 4-by-4 patch loads 4 values from each tile and performs 16 multiply-adds with them, so the ratio of arithmetic to scratchpad traffic improves fourfold. This is how hand-tuned kernels get close to the vendor libraries, and it is the point at which register pressure becomes the binding constraint rather than shared memory.

reference

CUDA C++ Best Practices Guide
NVIDIA CUTLASS
AMD Composable Kernel