a CUDA matmul — where
“utilization” finally means something

00 Start here — why a matrix multiply?

In Walk 01 each node of the graph was a scalar +=. In Walk 02 the heavy nodes — attention, the MLP — were matrix multiplies. That's not a coincidence: on a real model, the overwhelming majority of the GPU's time is spent multiplying matrices. Get matmul right and you've got the machine right.

So here we drop one level below Python and PyTorch, to the actual code the GPU runs: akernel. The task is SGEMM — Single-precision GEneral Matrix Multiply:C = α·(A·B) + β·C. It's the "hello world" of GPU programming and also the single most important operation in machine learning.

01 The mental model: threads, blocks, and the memory cliff

A GPU runs one small function — the kernel — across thousands of threads at once. You don't loop over the output; you launch a thread per output element and they all run in parallel. Threads are grouped into blocks; blocks tile the whole output.

The one fact that decides everything below is the memory hierarchy — and how steep the cliff is between its levels:

Every kernel below computes the exact same answer. The only thing that changes is how many times each number is dragged across that slow global-memory cliff. That is the whole game.

02 Kernel 1 — the naive version. Correct, and honest.

Here it is, verbatim. One thread computes one output element by walking a full row of A against a full column of B:

__global__ void sgemm_naive(int M, int N, int K, float alpha,
                            const float *A, const float *B,
                            float beta, float *C) {
  // this thread owns exactly one output element C[x][y]
  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];   // one row of A · one col of B
    }
    // C = α*(A@B) + β*C
    C[x * N + y] = alpha * tmp + beta * C[x * N + y];
  }
}

And this is how you launch it — the line that fires one thread per output cell:

dim3 gridDim(CEIL_DIV(M, 32), CEIL_DIV(N, 32));  // enough blocks to cover C
dim3 blockDim(32, 32);                           // 1024 threads per block
sgemm_naive<<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
//          ^^^^^^^^^^^^^^^^^^^^^^ launch: fire one thread per output element

In plain words:

• x, y — this thread figures out which output element it owns from its block and thread index.
• the for loop — the dot product: march across row x of A and column y of B, accumulating into tmp.
• the last line — write the result. Nothing wrong here. It compiles, it's correct, it'll pass any test.

Launch it on a 4096×4096 matmul on an A6000 and it delivers 309 GFLOP/s. The GPU's real capacity is ~23,000. So this correct kernel runs at 1.3% of what the chip can do — whilenvidia-smi shows it pinned at 100%.

Before we ask why, feel the operation in your hands. Every cell of C is one independent dot product — hover a cell to see its row × column, click to accumulate, and drag the dimensions to watch the FLOP count explode. The kernel above is exactly this: one thread per cell of C.

Want the full linear-algebra build-up behind this widget? It lives inStudy · Module 1 — Linear algebra & matmul.

03 Why so slow? Follow one number.

The GPU here isn't idle — it's starving. Every iteration of that inner loop reaches all the way out to global memory for a fresh A and B value. The arithmetic units can multiply far faster than memory can feed them, so they spend almost all their time waiting.

Worse: adjacent threads read A and B in a scattered pattern, so the hardware can't bundle their reads together. Each thread pays for its own slow trip. The chip looks fully occupied because threads are scheduled and running — they're just running the instruction that means "wait for memory."

04 Kernel 2 — coalescing: let neighbors read neighbors (→ 8.5%)

The code barely changes. All that changes is which thread computes which element, so that threads standing next to each other read memory addresses next to each other. When they do, the hardware fuses 32 separate reads into one transaction — "coalescing."

const int cRow = blockIdx.x * BLOCKSIZE + (threadIdx.x / BLOCKSIZE);
const int cCol = blockIdx.y * BLOCKSIZE + (threadIdx.x % BLOCKSIZE);

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

Same math, same memory hierarchy, one reindexing. Result: 309 → 1986 GFLOP/s, a6.4× speedup from nothing but access order. We went from 1.3% → 8.5% of peak without touching the arithmetic. The lesson: on a GPU, how you touch memory matters more than how much you compute.

05 Kernel 3 — shared memory: stop re-reading the same data (→ 12.8%)

The naive kernel reads every element of A and B from slow global memory many times — once for every output element that needs it. The fix: each block first copies a tile of A and B into fast on-chip shared memory, then does all its math against that. Read slow memory once, reuse it many times.

__shared__ float As[BLOCKSIZE * BLOCKSIZE];   // a tile of A, in fast on-chip memory
__shared__ float Bs[BLOCKSIZE * BLOCKSIZE];   // a tile of B

float tmp = 0.0;
for (int bkIdx = 0; bkIdx < K; bkIdx += BLOCKSIZE) {
  // every thread cooperatively loads ONE element of each tile from slow global memory
  As[threadRow * BLOCKSIZE + threadCol] = A[threadRow * K + threadCol];
  Bs[threadRow * BLOCKSIZE + threadCol] = B[threadRow * N + threadCol];
  __syncthreads();                       // wait until the whole tile is loaded

  A += BLOCKSIZE;  B += BLOCKSIZE * N;    // slide to the next tile

  // now do the dot-product against FAST shared memory, not global memory
  for (int dotIdx = 0; dotIdx < BLOCKSIZE; ++dotIdx) {
    tmp += As[threadRow * BLOCKSIZE + dotIdx] *
           Bs[dotIdx * BLOCKSIZE + threadCol];
  }
  __syncthreads();                       // don't overwrite the tile until everyone's done
}

The new machinery — and it's the pattern behind nearly every fast GPU kernel:

• __shared__ — a scratchpad on the chip itself, shared by every thread in the block.
• cooperative load — each thread grabs one element of the tile; together they fill it.
• __syncthreads() — a barrier: "nobody starts the math until the whole tile is loaded," and later "nobody overwrites the tile until everyone's finished with it."
• the inner loop now multiplies against As / Bs — fast memory.

Result: 2980 GFLOP/s, 12.8% of peak. Still far from the ceiling — but now the bottleneck hasmoved, and that's the real skill: each kernel fixes the current limiter and exposes the next one.

06 The whole ladder — one answer, 72× of speed on the table

Boehm keeps going: blocktiling (each thread computes many outputs), vectorized loads, resolving bank conflicts, warptiling, double-buffering. Every step targets the memory system, not the math. Here's the full climb on an A6000, 4096×4096 — every row computes the identical matrix:

Top to bottom: a 70× speedup on the same operation, same GPU, same correct answer. Every one of these kernels would light up nvidia-smi at ~100%. Only one number tells them apart, and it's not the one on the dashboard.

07 What actually separated them — the taste to take away

You don't need to memorize twelve kernels. You need the one mental model they all share: