Matrix multiplication

A naive matrix-multiply kernel assigns one output element to one thread: each thread walks the corresponding row of A and column of B and accumulates a dot product straight out of global memory. On an RTX A6000 multiplying two 4092×4092 float32 matrices, that kernel manages about 309 GFLOP/s, roughly 1.3% of what cuBLAS reaches on the same GPU and the same problem. The gap is not arithmetic: neighboring threads in a warp end up reading rows of A that sit nowhere near each other in memory, so almost none of those loads can be combined into a single wide transaction.

The worklog’s first kernel is the whole algorithm in fifteen lines, launched with one thread per entry of C:

__global__ void sgemm_naive(int M, int N, int K, float alpha, const float *A,
                            const float *B, float beta, float *C) {
  // compute position in C that this thread is responsible for
  const uint x = blockIdx.x * blockDim.x + threadIdx.x;
  const uint y = blockIdx.y * blockDim.y + threadIdx.y;

  if (x < M && y < N) {
    float tmp = 0.0;
    for (int i = 0; i < K; ++i) {
      tmp += A[x * K + i] * B[i * N + y];
    }
    C[x * N + y] = alpha * tmp + beta * C[x * N + y];
  }
}

Reassigning which thread owns which output element so that threads in a warp read consecutive addresses, letting the hardware coalesce those reads, already lifts throughput to about 1986.5 GFLOP/s with no other change. The next step is tiling: the kernel stages a block-sized tile of A and a tile of B into shared memory once, and every thread in the block reads back out of that on-chip copy instead of returning to global memory for each partial product, which pushes throughput to about 2980.3 GFLOP/s.

The heart of that shared-memory kernel is the loop every later version elaborates: stage a tile of A and a tile of B, synchronize, accumulate, synchronize, advance to the next tile along the reduction dimension:

// advance pointers to the starting positions
A += cRow * BLOCKSIZE * K;
B += cCol * BLOCKSIZE;
C += cRow * BLOCKSIZE * N + cCol * BLOCKSIZE;

float tmp = 0.0;
for (int bkIdx = 0; bkIdx < K; bkIdx += BLOCKSIZE) {
  As[threadRow * BLOCKSIZE + threadCol] = A[threadRow * K + threadCol];
  Bs[threadRow * BLOCKSIZE + threadCol] = B[threadRow * N + threadCol];

  // block threads in this block until cache is fully populated
  __syncthreads();

  // advance pointers onto next chunk
  A += BLOCKSIZE;
  B += BLOCKSIZE * N;

  // execute the dotproduct on the currently cached block
  for (int dotIdx = 0; dotIdx < BLOCKSIZE; ++dotIdx) {
    tmp += As[threadRow * BLOCKSIZE + dotIdx] *
            Bs[dotIdx * BLOCKSIZE + threadCol];
  }
  // need to sync again at the end, to avoid faster threads
  // fetching the next block into the cache before slower threads are done
  __syncthreads();
}
C[threadRow * N + threadCol] =
    alpha * tmp + beta * C[threadRow * N + threadCol];

The two barriers are the cooperation cost the transpose exercise introduced, now guarding both directions: the first keeps any thread from computing against a tile another thread has not finished staging, and the second keeps a fast thread from overwriting the tile with the next chunk while a slower one is still reading the current chunk. At a block size of 32 the two tiles occupy 8KB of the 48KB of shared memory a block can address on this GPU.

Even with tiling, each thread is still computing exactly one output element, so most instructions in the inner loop are shared-memory loads rather than the fused multiply-adds actually doing the work; a profiler shows warps repeatedly stalling in the “Stall MIO Throttle” state, waiting on the memory pipe rather than on arithmetic. Register blocking fixes that ratio directly: giving each thread a small 1D tile of outputs held in registers moves throughput to 8474.7 GFLOP/s, and a 2D tile of outputs per thread reaches 15971.7 GFLOP/s, 68.7% of cuBLAS. Vectorized memory instructions, autotuned tile sizes, and tiling at the warp level close most of what remains, reaching 21779.3 GFLOP/s, 93.7% of cuBLAS, without a tensor core in sight.

None of these kernels is limited by occupancy in the way a first guess might suggest: the register-blocked kernel above fits only one block per SM and still reaches 66% occupancy, and pushing that number higher would not by itself close the remaining gap to cuBLAS. What separates a kernel like this from a vendor library is mostly the same tiling idea applied recursively, at the block, warp, and instruction level. Libraries such as CUTLASS describe those nested tiles as compositions of layouts, a shape paired with a stride, and talk about a matrix being “K-major” when it is stride-1 along the reduction dimension, rather than relying on the row-major and column-major vocabulary borrowed from BLAS.

The case against occupancy-first reasoning is older than any GPU in this chapter. In 2008, Volkov and Demmel benchmarked dense linear algebra across four NVIDIA GPUs and built an SGEMM that sustained 58 to 60% of each chip’s peak where NVIDIA’s own CUBLAS 1.1 sustained 36 to 44%, and the design contradicted the official guidance of the era point by point: instead of many threads, shared memory as the primary storage, and long vectors, their kernel kept each 64×16 output block of C entirely in registers, staged only B’s block through shared memory, and ran short 64-element vector threads.

The paper’s method is what survived: measure, then reason. Varying the thread count showed the code reaching 32, 49, 58, and 59% of peak at one to four threads per core, and on the GTX280 those four threads correspond to 25% occupancy, from which the authors conclude that one should not over-optimize for occupancy, though extremely low occupancy can also hurt. Cycle counting on the disassembled binaries located the real bound in instruction throughput: CUBLAS ran twice as many warps yet was 1.6× slower, because keeping both input blocks in shared memory forced an extra register move for every two multiply-adds, diluting the multiply-add share of its inner loop to 56% of instructions against 82% in theirs.