The CUDA Handbook ▸ Chapter 17. Matrix Multiplication

17.3 Cache Blocking

The naive kernel threw out its reuse; drawing on rich tradition, a CPU developer’s first instinct would be to block the loops so the working set stays in cache.¹ A GPU has real L1 and L2 caches too, so the same strategy applies. Give each thread a small TM×TN microtile of C, keep those accumulators in registers, and nest the loops so the strip of A and B touched by a thread stays hot in cache – with no shared memory at all.

template<int TM, int TN>__global__ void sgemm_cache( int M, int N, int K, const float *A, const float *B, float *C ){    for ( int row0 = (blockIdx.y*blockDim.y + threadIdx.y)*TM; row0 < M; row0 += gridDim.y*blockDim.y*TM ) {        for ( int col0 = (blockIdx.x*blockDim.x + threadIdx.x)*TN; col0 < N; col0 += gridDim.x*blockDim.x*TN ) {            float acc[TM][TN] = {};            for ( int k = 0; k < K; k++ ) {                float a[TM], b[TN];                for ( int i = 0; i < TM; i++ ) a[i] = A[(row0+i)*K + k];                for ( int j = 0; j < TN; j++ ) b[j] = B[k*N + (col0+j)];                for ( int i = 0; i < TM; i++ ) {                    for ( int j = 0; j < TN; j++ ) acc[i][j] += a[i]*b[j];                }            }            // ... store acc[TM][TN] back to C ...        }    }}

The block’s footprint – blockDim multiplied by the microtile MT – specifies how much of A and B stays resident, so the best value of MT is an empirical question; on this GPU, sweeping the two tile parameters independently, the best microtile size is 8×8. The cache-blocked kernel reaches 3762 GFLOP/s, 4.2× the naive kernel, with only register accumulation and the hardware caches capturing the reuse the threads of a block share. The lesson to carry forward is that what matters is reuse per thread: a thread that computes 64 outputs from a strip of A and B amortizes each global load across 64 multiply-adds.

Although blocking for the GPU caches is a real and nearly free first win, as we shall see in the next section, shared memory offers even more performance.

Figure 17-2. A thread block computes a BM×BN tile of C, marching the K dimension in BK-deep slabs of A and B.