Learn: Optimizing GEMM

What you're learning

How to take the naive matmul from 04-first-cuda-kernels — which runs at ~1% of cuBLAS — and climb the classic optimization ladder to ~90%+, understanding the mechanism and the roofline consequence of every rung. By the end you should be able to write each kernel from memory, predict its arithmetic intensity, explain why each jump happens, and say precisely why the last ~10% belongs to vendor libraries. This is the canonical kernel-optimization case study, and it mirrors Simon Boehm’s superb worklog “How to Optimize a CUDA Matmul Kernel for cuBLAS-like Performance” — read it alongside this.

This lesson assumes 04-first-cuda-kernels (the naive SGEMM you’re optimizing), and leans hard on Tier 1: 01-gpu-architecture-and-simt (warps, occupancy), 02-gpu-memory-hierarchy (coalescing, SMEM, bank conflicts), and 03-performance-modeling-roofline (arithmetic intensity, the ridge point). We now put all of it to work on one kernel.

All performance numbers below are Boehm’s, measured on an RTX A6000 (Ampere, GA102), FP32, 4092×4092 square matmul, with cuBLAS = 100% ≈ 23.2 TFLOP/s.¹ The ladder is hardware-independent; only the constants move.

---

1. Learning map

graph TD
    A["Naive SGEMM<br/>(lesson 04)<br/>I = 1/4 FLOP/byte"] --> B["1. Global memory<br/>coalescing"]
    B --> C["2. Shared-memory<br/>tiling (blocking)"]
    C --> D["3. 1D register<br/>blocking (TM)"]
    D --> E["4. 2D register<br/>blocking (TM×TN)<br/>outer product"]
    E --> F["5. Vectorized float4 +<br/>transposed SMEM +<br/>bank-conflict-free"]
    F --> G["6. Warptiling"]
    G --> H["cuBLAS / CUTLASS<br/>double buffering,<br/>tensor cores, ptx"]

    R["Roofline (lesson 03)<br/>ridge = peak / BW"] -.classify.-> A
    R -.-> C
    R -.-> E

    style A fill:#a44,color:#fff
    style C fill:#44a,color:#fff
    style E fill:#4a4,color:#fff
    style H fill:#666,color:#fff

Each rung attacks a specific bottleneck exposed by the previous one. The whole story is one idea applied recursively down the memory hierarchy: do more math per byte you pull from the slow tier — HBM first (coalescing, SMEM tiling), then SMEM itself (register blocking), then the register file / instruction stream (vectorization, warptiling).

---

2. Why this matters

GEMM is the kernel. In a transformer the QKV projection, attention output projection, and the two MLP matmuls are all GEMMs, and they dominate FLOPs — for a compute-bound training step they are ~$>90\%$ of the useful work. A frontier training run’s cost (the $6ND$ estimate from 03-performance-modeling-roofline) is, at bottom, “how efficiently does your GEMM run?” A kernel at 40% of peak literally doubles your GPU-hours bill versus one at 80%.

You will almost never ship a hand-written GEMM — cuBLAS² and CUTLASS³ exist and are excellent. So why learn to build one? Three reasons:

1. It is the complete curriculum in one kernel. Coalescing, SMEM tiling, register blocking, bank conflicts, vectorization, occupancy — every technique in Tier 1 shows up, with an immediately measurable payoff. Nothing teaches the memory hierarchy like watching a matmul go $70\times$ faster.
2. The moment cuBLAS doesn’t fit, you’re on your own. Fused epilogues, unusual shapes (tall-skinny decode GEMVs), novel dtypes (FP8/MX with custom scaling), fused attention — the frontier is full of matmuls the vendor library doesn’t cover. CUTLASS is this exact ladder, templated.
3. It is the mental model for everything else. FlashAttention (06-memory-wall-and-flashattention) is “tiled GEMM + online softmax.” Once you see GEMM optimization, you see every IO-aware kernel.

The one sentence

Optimizing GEMM is the art of increasing data reuse: each byte loaded from a slow tier must feed as many FMAs as possible before you throw it away. Every rung on the ladder increases reuse at one level of the hierarchy.

The ladder at a glance

#	Kernel	GFLOP/s	% of cuBLAS	Bottleneck it fixes	New bottleneck
0	Naive	309	1.3%	—	uncoalesced HBM
1	GMEM coalescing	1987	8.5%	wasted HBM bandwidth	HBM traffic (no reuse)
2	SMEM tiling	2980	12.8%	HBM traffic (reuse in SMEM)	SMEM bandwidth
3	1D register blocking	8475	36.5%	SMEM bandwidth (reuse in regs)	SMEM bandwidth (still)
4	2D register blocking	15972	68.7%	SMEM bandwidth (outer product)	SMEM/GMEM load efficiency
5	Vectorized + transposed SMEM	18237	78.4%	load width, bank conflicts	warp-level data movement
6	Warptiling	21779	93.7%	register/warp scheduling	the last ~6%
—	cuBLAS	23250	100%	double-buffering, ptx, assembly	—

Watch the second column: coalescing is a $6\times$ jump, SMEM tiling $1.5\times$, but register blocking is the big one ($2.8\times$ then $1.9\times$) that takes you from a toy to a serious kernel.¹ Keep that in mind — it’s the least obvious and most important rung.

---

3. Rung 0 — Naive (recap): memory-bound and uncoalesced

Recall the naive kernel from 04-first-cuda-kernels. One thread per output element, dot product straight out of global memory:

// C = alpha * (A @ B) + beta * C ;  A: MxK, B: KxN, C: MxN  (all row-major)
__global__ void sgemm_naive(int M, int N, int K, float alpha,
                            const float *A, const float *B,
                            float beta, float *C) {
  const uint row = blockIdx.x * blockDim.x + threadIdx.x;  // maps to M
  const uint col = blockIdx.y * blockDim.y + threadIdx.y;  // maps to N
  if (row < M && col < N) {
    float acc = 0.0f;
    for (int k = 0; k < K; ++k)
      acc += A[row * K + k] * B[k * N + col];   // K reads of A, K reads of B
    C[row * N + col] = alpha * acc + beta * C[row * N + col];
  }
}
// launch: dim3 block(32,32); dim3 grid(ceil(M/32), ceil(N/32));

Arithmetic intensity. Every output element does a length-$K$ dot product, reading $K$ elements of $A$ and $K$ of $B$ from HBM — nothing is cached across output elements. Total HBM read traffic is $2MNK$ elements for $2MNK$ FLOPs:

$$I_{\text{naive}}=\frac{2MNK\text{ FLOP}}{2MNK\cdot 4\text{ B}}=\frac{1}{4}\frac{\text{FLOP}}{\text{byte}}$$

That is $0.25$ FLOP/byte, wildly below the A6000’s FP32 ridge point of $\approx 50$ FLOP/byte (38.7 TFLOP/s FP32 $÷$ 768 GB/s).⁴ Deeply memory-bound — the SM’s FP32 ALUs sit idle almost the entire time waiting on HBM.

And it’s worse than that, because the access pattern is uncoalesced. With block(32,32), threadIdx.x is the fast index within a warp, so a warp holds 32 consecutive row values (same col). Look at the two loads:

• A[row*K + k]: 32 threads → 32 different rows, stride $K$ → 32 separate cache lines per load. ~3% bandwidth efficiency.
• B[k*N + col]: col is the same for all 32 threads → all read the same address → broadcast (fine).
• C[row*N + col] write: 32 different rows, stride $N$ → uncoalesced.

So the dominant traffic (A and C) is uncoalesced. You’re not even getting the full memory-bound performance — you get $\sim 1/32$ of it. That’s why the naive kernel lands at a miserable 1.3%.¹

Two different sins, don't conflate them

“Low arithmetic intensity” (memory-bound) and “uncoalesced access” are independent problems. The naive kernel has both. Coalescing (rung 1) fixes only the second — it recovers the bandwidth you were wasting but does not change $I$; the kernel is still memory-bound at $I=1/4$. Raising $I$ needs reuse, which is rungs 2+.

---

4. Rung 1 — Global memory coalescing (1.3% → 8.5%)

The fix is pure thread-index remapping — no algorithm change, no extra memory. We want consecutive threads in a warp to hit consecutive columns so their loads/stores coalesce into single cache lines.

Use a 1D block of BLOCKSIZE*BLOCKSIZE threads and map:

template <const int BLOCKSIZE>
__global__ void sgemm_coalesced(int M, int N, int K, float alpha,
                                const float *A, const float *B,
                                float beta, float *C) {
  // consecutive threadIdx.x -> consecutive column (fast), same row within a warp
  const uint row = blockIdx.x * BLOCKSIZE + (threadIdx.x / BLOCKSIZE);
  const uint col = blockIdx.y * BLOCKSIZE + (threadIdx.x % BLOCKSIZE);
  if (row < M && col < N) {
    float acc = 0.0f;
    for (int k = 0; k < K; ++k)
      acc += A[row * K + k] * B[k * N + col];
    C[row * N + col] = alpha * acc + beta * C[row * N + col];
  }
}
// launch: dim3 block(BLOCKSIZE*BLOCKSIZE); grid(ceil(M/BS), ceil(N/BS)); BS=32

Now within a warp of 32 threads (threadIdx.x = 0..31), col = ...% 32 runs $\mathrm{0..31}$ and row is constant. Re-examine:

• B[k*N + col]: consecutive col → 32 contiguous floats = 1 cache line. Coalesced.
• C[row*N + col] write: contiguous → coalesced.
• A[row*K + k]: same row for the whole warp → broadcast (one address). Fine.

Every HBM transaction is now fully used. Same FLOPs, same $I=1/4$, but effective bandwidth jumps ~$6\times$ → 8.5% of cuBLAS.¹ On the roofline you haven’t moved right (intensity unchanged); you’ve moved up onto the memory roof you were sitting far below.

Coalescing is the cheapest big win in GPU programming

A thread-index permutation, zero extra memory, ~$6\times$. This is why “make threadIdx.x map to the contiguous memory dimension” is the first thing to check in any kernel. In Nsight Compute, look at global load/store efficiency — the naive kernel shows ~3–12%, the coalesced one ~100%.

Quiz: In the naive kernel, A[row*K+k] was uncoalesced but B[k*N+col] was a harmless broadcast. After the coalescing remap, which of A/B/C is a broadcast and which are coalesced — and why did just swapping the row/col mapping fix the dominant traffic?

Answer

After the remap, a warp fixes row and varies col across its 32 lanes. So:

• B (B[k*N+col]) is now coalesced — consecutive col → contiguous addresses.
• C write (C[row*N+col]) is coalesced — same reason.
• A (A[row*K+k]) is now the broadcast — row and k are identical across the warp, so all 32 lanes read one address (served in a single transaction).

It worked because the naive kernel had it exactly backwards: the two large, strided accesses (A-read and C-write) were the uncoalesced ones and the reused access (B) was already a broadcast. Swapping which thread index maps to the contiguous dimension turned the two costly strided patterns into a coalesced pattern + a broadcast. The total bytes fetched from HBM barely changed conceptually, but the wasted bytes (lines fetched but mostly unused) collapsed — recovering the ~$32\times$ efficiency you were throwing away, hence the ~$6\times$ measured speedup.

---

5. Rung 2 — Shared-memory tiling / blocking (8.5% → 12.8%)

Coalescing fixed how we read HBM; it did nothing about how much. We still stream $2MNK$ elements. The problem is no reuse: element A[row][k] is read $N$ times (once per output column in that row), each time from HBM. The whole point of shared memory (from 02-gpu-memory-hierarchy) is to load a chunk from HBM once, then let every thread in the block reuse it from the ~20 TB/s SMEM scratchpad.

The tiling scheme

Each block computes one BM×BN output tile of $C$. It marches across the $K$ dimension in steps of BK. At each step it cooperatively loads a BM×BK tile of $A$ and a BK×BN tile of $B$ into SMEM, syncs, then every thread accumulates its partial dot product from SMEM. Loop, accumulate, write once.

Loop over K/BK tiles: load As,Bs into SMEM → __syncthreads() → accumulate As·Bs into C → __syncthreads()

Figure 1 — Shared-memory tiling. Each block owns one BM×BN output tile of C (corner). It sweeps the shared K-dimension in steps of BK: at each step it stages a BM×BK tile of A and a BK×BN tile of B from HBM into shared memory (As, Bs), then every thread reads those tiles many times to accumulate its dot-product partials. The highlighted A-tile × B-tile is one term of the block’s accumulation.

The kernel

template <const int BLOCKSIZE>   // BM = BN = BK = BLOCKSIZE (square, e.g. 32)
__global__ void sgemm_smem(int M, int N, int K, float alpha,
                           const float *A, const float *B,
                           float beta, float *C) {
  const uint cRow = blockIdx.x;   // block's output-tile row
  const uint cCol = blockIdx.y;   // block's output-tile col
 
  __shared__ float As[BLOCKSIZE * BLOCKSIZE];
  __shared__ float Bs[BLOCKSIZE * BLOCKSIZE];
 
  const uint threadRow = threadIdx.x / BLOCKSIZE;
  const uint threadCol = threadIdx.x % BLOCKSIZE;   // coalesced: fast index = col
 
  // advance the base pointers to this block's tiles
  A += cRow * BLOCKSIZE * K;                     // row-block of A
  B += cCol * BLOCKSIZE;                         // col-block of B
  C += cRow * BLOCKSIZE * N + cCol * BLOCKSIZE;  // the output tile
 
  float acc = 0.0f;
  for (int bkIdx = 0; bkIdx < K; bkIdx += BLOCKSIZE) {
    // cooperative, coalesced load: each thread stages one element of A and B
    As[threadRow * BLOCKSIZE + threadCol] = A[threadRow * K + threadCol];
    Bs[threadRow * BLOCKSIZE + threadCol] = B[threadRow * N + threadCol];
    __syncthreads();                 // (1) all tiles present before anyone reads
 
    A += BLOCKSIZE;                  // slide right along K
    B += BLOCKSIZE * N;             // slide down along K
 
    // dot-product this tile from SMEM
    for (int dotIdx = 0; dotIdx < BLOCKSIZE; ++dotIdx)
      acc += As[threadRow * BLOCKSIZE + dotIdx] *
             Bs[dotIdx * BLOCKSIZE + threadCol];
    __syncthreads();                 // (2) don't overwrite SMEM until all done reading
  }
  C[threadRow * N + threadCol] = alpha * acc + beta * C[threadRow * N + threadCol];
}

The two __syncthreads() are not optional and not interchangeable

Barrier (1) guarantees every thread has finished writing its As/Bs element before any thread starts reading the tile — without it, a fast warp reads garbage (a read-before-write race). Barrier (2) guarantees every thread has finished reading this tile before any thread overwrites As/Bs with the next tile — without it you get a write-before-read race that corrupts slow warps’ results. Both are races on __shared__ memory across warps in the block. A classic bug is dropping barrier (2) “because it seemed to still pass” — it passes until a scheduler timing change silently corrupts results. Also: __syncthreads() must be reached by all threads in the block; putting it inside a divergent if (row < M) branch deadlocks.

How much HBM traffic did we save?

This is the quantitative payoff. With a square tile of side $S=$ BLOCKSIZE:

• Number of output tiles $=\frac{M}{S}\cdot \frac{N}{S}$.
• Each tile loops $\frac{K}{S}$ times, loading $S\times S$ of $A$ and $S\times S$ of $B$ each iteration.
• HBM loads per tile $=\frac{K}{S}(S^2+S^2)=2KS$.

Total HBM load traffic (elements):

$$\underset{\#\text{tiles}}{\underbrace{\frac{MN}{S^2}}}\cdot \underset{\text{per tile}}{\underbrace{2KS}}=\frac{2MNK}{S}$$

Compare to naive’s $2MNK$. Tiling cuts HBM traffic by a factor of $S$.⁵ For $S=32$ that’s a 32× reduction. Equivalently, each element loaded into SMEM is reused $S$ times before eviction. Arithmetic intensity rises:

$$I_{\text{tiled}}=\frac{2MNK}{\frac{2MNK}{S}\cdot 4}=\frac{S}{4}\Rightarrow I_{S=32}=8\frac{\text{FLOP}}{\text{byte}}$$

On the roofline you’ve now moved right by $32\times$ — from $I=0.25$ to $I=8$.

So why only 12.8%, not more? ¹

Because $I=8$ is still below the FP32 ridge ($\approx 50$): the kernel is now bound by SMEM bandwidth, not HBM. Look at the inner loop: for every FMA (2 FLOPs) each thread issues two SMEM loads (one As, one Bs). That’s a SMEM arithmetic intensity of ~1 FLOP per SMEM load — the shared-memory pipe is now the bottleneck, exactly as HBM was before. The fix is the same idea one level down: reuse SMEM values in registers. That’s rung 3, and it’s where the real speed lives.

Quiz: You run the SMEM-tiled SGEMM with BLOCKSIZE = 64 on a 4096×4096×4096 FP32 problem. (a) By what factor does HBM load traffic drop vs. naive? (b) What is the arithmetic intensity, and is it above the A6000 FP32 ridge (~50)? (c) How many bytes of SMEM does the block need, and is that a problem?

Answer

(a) Traffic drops by the tile side $S=64$ → 64× less HBM traffic ($2MNK\to 2MNK/64$). Each loaded element is reused 64 times.

(b) $I_{\text{tiled}}=S/4=64/4=16$ FLOP/byte. Still below the ~50 ridge → the HBM roofline says memory-bound at the HBM level, but in practice the binding constraint has already shifted to SMEM bandwidth (2 SMEM loads per FMA). Either way you are not yet compute-bound; you need register blocking.

(c) SMEM = As[64*64] + Bs[64*64] = $2\cdot 64^2\cdot 4\text{ B}=32,768$ B = 32 KB per block. That fits (A6000 allows up to 100 KB SMEM/block with opt-in), but it hurts occupancy: fewer blocks co-resident per SM. Bigger tiles reuse more but cost more SMEM and registers — the central tension of GEMM tuning. (This is also why the ideal traffic — reading each operand ~once, $I=2n/3b$, hundreds of FLOP/byte — is only approached with the help of the L2 cache capturing reuse across blocks, not by SMEM alone.)

---

6. Rung 3 & 4 — Register blocking / thread coarsening (12.8% → 68.7%)

This is the pivotal rung. The insight: the register file is the tier below SMEM — ~256 KB/SM, effectively unlimited bandwidth to the ALUs, ~1-cycle latency. If SMEM bandwidth is the bottleneck, stage a small tile of As/Bs into registers and reuse those across many FMAs. Instead of one thread computing one output, each thread now computes a small micro-tile of outputs, amortizing every SMEM load over many FLOPs. This is thread coarsening: fewer, fatter threads.

1D register blocking (TM): 12.8% → 36.5%

First step: each thread computes a column of TM output elements (e.g. TM=8). Tile params: BM=64, BN=64, BK=8, TM=8. The trick is in the inner loop — load one value of Bs into a register, reuse it across all TM rows:

template <const int BM, const int BN, const int BK, const int TM>
__global__ void sgemm_1d(int M, int N, int K, float alpha,
                         const float *A, const float *B, float beta, float *C) {
  const uint cRow = blockIdx.y, cCol = blockIdx.x;
 
  __shared__ float As[BM * BK];
  __shared__ float Bs[BK * BN];
 
  // each thread owns a TM-tall column of the output tile
  const uint threadRow = threadIdx.x / BN;
  const uint threadCol = threadIdx.x % BN;
 
  A += cRow * BM * K;
  B += cCol * BN;
  C += cRow * BM * N + cCol * BN;
 
  // separate indices for the (coalesced) cooperative loads
  const uint innerRowA = threadIdx.x / BK, innerColA = threadIdx.x % BK;
  const uint innerRowB = threadIdx.x / BN, innerColB = threadIdx.x % BN;
 
  float threadResults[TM] = {0.0f};   // <-- TM accumulators live in registers
 
  for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
    As[innerRowA * BK + innerColA] = A[innerRowA * K + innerColA];
    Bs[innerRowB * BN + innerColB] = B[innerRowB * N + innerColB];
    __syncthreads();
 
    A += BK;  B += BK * N;
 
    for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
      float tmpB = Bs[dotIdx * BN + threadCol];         // 1 SMEM load...
      for (uint resIdx = 0; resIdx < TM; ++resIdx)      // ...reused TM times
        threadResults[resIdx] +=
            As[(threadRow * TM + resIdx) * BK + dotIdx] * tmpB;
    }
    __syncthreads();
  }
 
  for (uint resIdx = 0; resIdx < TM; ++resIdx)
    C[(threadRow * TM + resIdx) * N + threadCol] =
        alpha * threadResults[resIdx] +
        beta  * C[(threadRow * TM + resIdx) * N + threadCol];
}

The inner loop now does TM FMAs (2*TM FLOPs) per iteration while loading 1 + TM SMEM values (tmpB once, plus TM from As). SMEM arithmetic intensity goes from ~1 to $\frac{2TM}{1+TM}\approx 1.8$ for TM=8. Small change in the ratio, but it nearly triples throughput → 36.5%, because we were badly SMEM-bound.¹

2D register blocking (TM×TN): the outer product — 36.5% → 68.7%

Now do it in both dimensions. Each thread computes a TM×TN micro-tile (e.g. 8×8 = 64 outputs). The inner kernel becomes an outer product: load a length-TM column of As and a length-TN row of Bs into registers, then take their outer product into the TM×TN accumulator. This is the formulation that maps directly onto tensor-core MMAs and is what CUTLASS generalizes.

Figure 2 — The register-blocking micro-tile (outer product). Each thread loads a length-TM slice of column dotIdx from As into registers (regM) and a length-TN slice of row dotIdx from Bs (regN), then does all TM×TN FMAs of their outer product into a register accumulator. TM+TN SMEM loads feed TM·TN FMAs — with TM=TN=8, 16 loads feed 64 FMAs, an 8× jump in SMEM reuse over rung 2.

template <const int BM, const int BN, const int BK, const int TM, const int TN>
__global__ void sgemm_2d(int M, int N, int K, float alpha,
                         const float *A, const float *B, float beta, float *C) {
  const uint cRow = blockIdx.y, cCol = blockIdx.x;
  const uint totalResults = BM * BN;
  const uint numThreads = totalResults / (TM * TN);   // e.g. 128*128/64 = 256
 
  __shared__ float As[BM * BK];
  __shared__ float Bs[BK * BN];
 
  // this thread's position in the grid of TM×TN micro-tiles
  const uint threadCol = threadIdx.x % (BN / TN);
  const uint threadRow = threadIdx.x / (BN / TN);
 
  A += cRow * BM * K;
  B += cCol * BN;
  C += cRow * BM * N + cCol * BN;
 
  // cooperative-load indices (each thread loads several elems; strided across the tile)
  const uint innerRowA = threadIdx.x / BK, innerColA = threadIdx.x % BK;
  const uint strideA   = numThreads / BK;
  const uint innerRowB = threadIdx.x / BN, innerColB = threadIdx.x % BN;
  const uint strideB   = numThreads / BN;
 
  float threadResults[TM * TN] = {0.0f};   // 64 accumulators in registers
  float regM[TM] = {0.0f};
  float regN[TN] = {0.0f};
 
  for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
    // load As and Bs (multiple elements per thread, all coalesced)
    for (uint off = 0; off < BM; off += strideA)
      As[(innerRowA + off) * BK + innerColA] = A[(innerRowA + off) * K + innerColA];
    for (uint off = 0; off < BK; off += strideB)
      Bs[(innerRowB + off) * BN + innerColB] = B[(innerRowB + off) * N + innerColB];
    __syncthreads();
 
    A += BK;  B += BK * N;
 
    // the outer-product inner loop
    for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
      for (uint i = 0; i < TM; ++i) regM[i] = As[(threadRow * TM + i) * BK + dotIdx];
      for (uint i = 0; i < TN; ++i) regN[i] = Bs[dotIdx * BN + threadCol * TN + i];
      for (uint m = 0; m < TM; ++m)
        for (uint n = 0; n < TN; ++n)
          threadResults[m * TN + n] += regM[m] * regN[n];   // TM×TN FMAs, 0 SMEM loads
    }
    __syncthreads();
  }
 
  // write the micro-tile back
  for (uint m = 0; m < TM; ++m)
    for (uint n = 0; n < TN; ++n) {
      uint r = threadRow * TM + m, c = threadCol * TN + n;
      C[r * N + c] = alpha * threadResults[m * TN + n] + beta * C[r * N + c];
    }
}
// launch: block(BM*BN/(TM*TN)); grid(ceil(N/BN), ceil(M/BM));
//         BM=BN=128, BK=8, TM=TN=8  ->  256 threads/block

Why this is the key jump. Count the inner loop’s SMEM traffic per dotIdx: TM + TN SMEM loads feed TM·TN FMAs. The SMEM arithmetic intensity is

$$\frac{\text{FLOPs}}{\text{SMEM loads}}=\frac{2TM\cdot TN}{TM+TN}\stackrel{TM=TN=8}{\to }\frac{2\cdot 64}{16}=8\frac{\text{FLOP}}{\text{load}}$$

versus ~1 for the plain SMEM kernel — an 8× improvement in SMEM reuse, achieved by moving the innermost reuse into the register file. Now the ALUs, not SMEM, are close to the bottleneck. This single idea takes you from 12.8% to 68.7% — more than the entire rest of the ladder combined.¹ The outer-product structure is also exactly the shape of a tensor-core MMA (a small matrix-multiply-accumulate held in registers), which is why this formulation is the bridge to CUTLASS and tensor cores.³

The recursion made explicit

Coalescing + SMEM tiling raised FLOPs per HBM byte. Register blocking raises FLOPs per SMEM byte. It’s the same move — cache-and-reuse — applied at the next tier down. Each tier is ~10–100× faster than the one above it, so pushing reuse down the hierarchy is always the lever. GEMM is the purest demonstration of this principle.

Register blocking trades occupancy for reuse — profile the sweet spot

TM×TN = 64 accumulators + regM/regN + addressing pushes register usage up (often 100+ regs/thread). Per 02-gpu-memory-hierarchy, that caps resident threads and can lower occupancy. Usually worth it — a coarsened thread hides latency with instruction-level parallelism (many independent FMAs) instead of thread-level, so it tolerates low occupancy. But push TM,TN too far and you spill to local memory (HBM!), which is catastrophic. Check ptxas -v for spills and sweep BM,BN,BK,TM,TN — don’t guess. (See rung 5’s autotuning.)

Quiz: The plain SMEM kernel and the 2D register-blocked kernel move the same amount of data from HBM (same tiling). Yet register blocking is ~5× faster. If HBM traffic is unchanged, where does the speedup come from — and why couldn't you get it just by making the SMEM tile bigger?

Answer

The speedup comes from relieving shared-memory bandwidth, not HBM bandwidth. In the plain SMEM kernel every FMA needs 2 SMEM loads, so the kernel is bottlenecked on the SMEM pipe (~1 FLOP/load). Register blocking stages a TM-column and TN-row into registers once per dotIdx and reuses them across TM·TN FMAs, giving $\frac{2TMTN}{TM+TN}=8$ FLOP per SMEM load — an 8× cut in SMEM traffic per FLOP. The ALUs (and instruction issue) become the limiter, which is much closer to peak.

A bigger SMEM tile would not fix this: it reduces HBM traffic (more reuse of HBM loads), but the inner loop still issues 2 SMEM loads per FMA regardless of tile size. The SMEM-bandwidth bottleneck is a property of the inner loop’s reuse pattern, not the tile dimensions. Only moving reuse into registers (a faster tier) removes it. This is why the roofline lens must be applied per tier: you can be “HBM compute-bound” yet “SMEM memory-bound” at the same time.

---

7. Rung 5 — Vectorized loads, transposed SMEM, bank conflicts (68.7% → 78.4%) + warptiling (→ 93.7%)

At ~69% the remaining losses are in how data moves between HBM/SMEM/registers, not in reuse. Three refinements:

(a) Vectorized float4 loads. Instead of 4 separate 32-bit loads, issue one 128-bit LDG.E.128 / LDS.128. Wider transactions mean fewer instructions and better bus utilization. Reinterpret the pointer:

// load 4 contiguous floats from HBM in one instruction
float4 tmp = reinterpret_cast<const float4 *>(&A[(innerRowA)*K + innerColA*4])[0];
// ...and store them into SMEM (see transpose below)

Requirement: addresses must be 16-byte aligned and contiguous. float4 also lets you halve instruction count in the SMEM read of the inner loop.

(b) Transposed As in SMEM. In the 2D kernel, the inner loop reads a column of As (As[(threadRow*TM+i)*BK + dotIdx] for varying i) — a strided SMEM access. Store As transposed (As[BK][BM] layout) when loading, so the inner loop reads a contiguous row and can use float4:

float4 t = reinterpret_cast<const float4 *>(&A[innerRowA*K + innerColA*4])[0];
// transposed store: column of A-tile becomes a row in As
As[(innerColA*4 + 0) * BM + innerRowA] = t.x;
As[(innerColA*4 + 1) * BM + innerRowA] = t.y;
As[(innerColA*4 + 2) * BM + innerRowA] = t.z;
As[(innerColA*4 + 3) * BM + innerRowA] = t.w;

The GMEM read stays coalesced (float4 over contiguous A), and the transposed layout makes the inner-loop SMEM read contiguous and vectorizable.

(c) Bank-conflict-free SMEM. Recall from 02-gpu-memory-hierarchy: SMEM has 32 banks; if a warp’s 32 lanes hit the same bank with different addresses, accesses serialize. In the outer-product inner loop, the way the 32 threads of a warp index As/Bs determines conflicts. The float4 + transposed layout, plus laying out the warp’s threads so their regM/regN reads span distinct banks (and padding where a stride shares a factor with 32), removes the serialization. Boehm’s vectorized kernel resolves the remaining conflicts and lands at 78.4%.¹

reinterpret_cast<float4> will fault on misaligned pointers

LDG.128/LDS.128 require 16-byte alignment. If your tile offset makes &A[...] not a multiple of 16 bytes (e.g. odd BK, or a leading-dimension not divisible by 4), the load is illegal — silent corruption or an unaligned-address fault. Vectorization constrains your tile sizes to multiples of 4; that constraint is one reason autotuners search a discrete grid of shapes.

Warptiling (→ 93.7%). The final rung adds a middle level between the block-tile and the thread-tile: an explicit warp-tile. The output tile is partitioned block → warp → thread. Why bother? Because the warp is the real unit of scheduling and of SMEM/register access. Giving each warp a contiguous sub-tile (a) keeps a warp’s operands in a fixed register/SMEM region for better locality and instruction scheduling, (b) enables register reuse within the warp across its threads’ micro-tiles, and (c) matches how tensor-core MMAs are issued per-warp. Warptiling is the structure CUTLASS uses (block-tile → warp-tile → MMA/thread-tile).³ It buys the last ~15 points to 93.7%.¹

Quiz: Storing As transposed in shared memory adds work at load time (a scatter of 4 elements). Why is it a net win — what does it buy in the inner loop, and how does it relate to both float4 and bank conflicts?

Answer

The inner loop reads a column of the A-tile (regM[i] = As[... ][dotIdx] for i=0..TM-1). In the natural As[BM][BK] layout those TM elements are BK apart — a strided SMEM read: not vectorizable, and (since the stride can share a factor with 32) prone to bank conflicts. Transposing to As[BK][BM] makes those same TM elements contiguous, so the inner loop can issue a single float4/LDS.128 and the 32 lanes of a warp hit 32 distinct banks — conflict-free.

The one-time transpose at load costs little because the load itself stays coalesced (you read contiguous A from HBM with float4), you just store it scattered into SMEM. You pay a cheap scatter once per tile to make the hot inner loop — executed BK times per tile, over all tiles — both vectorized and conflict-free. Net: fewer SMEM instructions, no serialization, ~10% faster. It’s the same “spend a little at the boundary to speed up the hot loop” logic as padding a transpose tile [32][33].

---

8. Rung 6 — Where this lands vs cuBLAS / CUTLASS, and the last ~10%

At warptiling you’re at ~94% of cuBLAS with a ~600-line kernel.¹ The remaining ~6–10% is real engineering that vendor libraries have and you (probably) don’t:

• Double buffering / software pipelining. Overlap the next tile’s HBM→SMEM load with the current tile’s compute, so the ALUs never stall at the __syncthreads(). On Ampere+ this uses cp.async (asynchronous global→shared copies) and a multi-stage SMEM ring buffer. This is the single biggest remaining lever and the reason cuBLAS hides memory latency almost perfectly.
• Tensor cores. For FP32 this ladder uses CUDA cores; the huge wins come from TF32/BF16/FP8 tensor-core MMAs (mma.sync / wgmma on Hopper), 8–16× the FP32 throughput (03-performance-modeling-roofline). Your outer-product micro-tile is already the right shape to feed them, but issuing MMAs correctly (fragment layouts, ldmatrix) is intricate.
• Hand-tuned PTX/SASS. cuBLAS ships per-architecture assembly with optimal instruction scheduling, register allocation, and dual-issue that nvcc won’t match from C++.
• Shape-specialized kernels + heuristics. cuBLAS dispatches among hundreds of pre-tuned kernels by (M,N,K, dtype, layout).² CUTLASS is the open-source expression of this: templated block/warp/thread tiling + pipelining + MMA atoms — literally this ladder, parameterized and productionized.³

Know when to stop

The economics are stark (Boehm’s own note): reaching ~80% took two weekends; the next 14% took four more.¹ For almost all work, call cuBLAS/CUTLASS and move on — you will not beat them on standard shapes, and your time is better spent elsewhere (fusion, attention, parallelism). Build the ladder once to understand it; reach for it in anger only when your shape/dtype/epilogue is one the library doesn’t serve well.

Premature optimization is the real trap

The biggest mistake isn’t a slow kernel — it’s optimizing the wrong one, or optimizing before measuring. Always nsys to find the kernel that dominates runtime, then ncu to learn why it’s slow (HBM-bound? SMEM-bound? occupancy? spills?), then pick the matching rung. Guessing which rung you need without a profiler is how you spend a weekend on a 2%-of-runtime kernel, or add register blocking to a kernel that was actually uncoalesced. The roofline + profiler tell you which bottleneck you’re on; the ladder tells you the fix for that bottleneck.

Roofline summary of the whole climb

graph LR
    N["Naive<br/>I≈0.25, uncoalesced<br/>≪ HBM roof"] --> C1["Coalesced<br/>I≈0.25, ON HBM roof<br/>(recovered BW)"]
    C1 --> C2["SMEM tiled<br/>I≈S/4, moved RIGHT<br/>now SMEM-bound"]
    C2 --> C3["Register blocked<br/>SMEM reuse ↑8×<br/>near compute roof"]
    C3 --> C4["Vectorized+warptiled<br/>≈ compute roof"]
    style N fill:#a44,color:#fff
    style C2 fill:#44a,color:#fff
    style C4 fill:#4a4,color:#fff

The through-line: coalescing lifts you onto the memory roof; SMEM tiling and register blocking walk you right toward the ridge (raising intensity at each tier); vectorization/warptiling/double-buffering squeeze you up onto the compute roof.

---

9. Practice problems

Practice problem 1: Pick tile sizes from the SMEM and register budget

Problem: You’re writing a 2D register-blocked FP32 SGEMM (rung 4) for an SM with 100 KB usable SMEM/block and 65,536 registers/SM (max 255 regs/thread, max 2048 threads/SM). You want BM=BN=128, BK=8, TM=TN=8. (a) How many threads per block? (b) How much SMEM per block, and does it fit? (c) Roughly how many registers per thread just for accumulators + operands, and what does that imply for occupancy? (d) What single change would you try first to raise occupancy, and what’s the risk?

Worked solution:
(a) Threads $=\frac{BM\cdot BN}{TM\cdot TN}=\frac{128\cdot 128}{8\cdot 8}=\frac{16384}{64}=\text{256}$ threads/block.

(b) SMEM $=(BM\cdot BK+BK\cdot BN)\cdot 4=(128\cdot 8+8\cdot 128)\cdot 4=2048\cdot 4=\text{8192 B}=8\text{ KB}$. Fits easily in 100 KB (so you could raise BK or use double-buffering with room to spare).

(c) Accumulators alone: threadResults[TM*TN] = 64 registers. Plus regM[8] + regN[8] = 16, plus pointers/indices/loop vars (~10–20). So ~90–110 regs/thread realistically. At ~100 regs: resident threads $\le 65536/100\approx 655$ → about $655/2048\approx$ 32% occupancy from registers, and only $\lfloor 655/256\rfloor =2$ blocks/SM. Low occupancy — but 2D register blocking hides latency via ILP (64 independent FMAs), so it’s usually fine.

(d) First lever: shrink TM,TN (e.g. to 8×4 → 32 accumulators) to cut registers and raise occupancy — or raise occupancy via __launch_bounds__. Risk: too few accumulators reduces SMEM reuse (back toward the SMEM-bound regime), and forcing occupancy up can trigger register spills to local memory (HBM), which is far worse than low occupancy. The only way to know is to sweep and profile (ptxas -v for spills, ncu for achieved occupancy and SM/mem throughput) — this is exactly what an autotuner does.

Practice problem 2: Prove the traffic reduction and place each rung on the roofline

Problem: For a square FP32 SGEMM ($M=N=K=n$) on a machine with FP32 ridge $I^{\ast }=50$ FLOP/byte: (a) derive HBM traffic and $I$ for the naive kernel; (b) for SMEM tiling with side $S$; (c) for the ideal kernel that reads each operand exactly once (the asymptote cuBLAS approaches via L2). For $n=4092$, evaluate (b) at $S=32,128$ and (c), and state which are HBM-compute-bound.

Worked solution:
(a) Naive. Each of $n^2$ outputs reads $n$ from A and $n$ from B → $2n^3$ element-reads. FLOPs $=2n^3$. $I=\frac{2n^3}{2n^3\cdot 4}=\frac{1}{4}=0.25$. $\ll 50$ → memory-bound (and uncoalesced on top).

(b) SMEM tiled, side $S$. Traffic $=\frac{2n^3}{S}$ elements (derived in §5). $I=\frac{2n^3}{(2n^3/S)\cdot 4}=\frac{S}{4}$.

• $S=32$: $I=8$ FLOP/byte $\ll 50$ → still HBM-memory-bound (and SMEM-bound in practice).
• $S=128$: $I=32$ FLOP/byte, closer but still below 50 → not yet HBM-compute-bound by the SMEM-only model.

(c) Ideal (each operand read once). Traffic $=(MK+KN+MN)=3n^2$ elements. $I=\frac{2n^3}{3n^2\cdot 4}=\frac{n}{6}$. For $n=4092$: $I\approx 682$ FLOP/byte $\gg 50$ → firmly compute-bound.

The lesson: SMEM tiling alone (even with big tiles) doesn’t push $I$ above the ridge for FP32 — yet real GEMM is compute-bound. The gap is closed by the L2 cache capturing reuse across blocks (an A block-row is reused by all blocks in that row of C; a B block-col likewise), so effective HBM traffic falls toward the ideal $3n^2$. This is why the actual kernel reaches the compute roof even though the per-block SMEM model predicts otherwise — and why measuring dram__bytes.sum in ncu (real HBM traffic) matters more than a hand model.

---

10. Practice & resources

Hands-on — build the whole ladder

1. Implement the ladder and benchmark each rung vs cuBLAS. Write all six kernels (naive → coalesced → SMEM → 1D → 2D → vectorized) plus a cublasSgemm baseline. Time each with CUDA events on a 4096³ FP32 problem, report GFLOP/s and % of cuBLAS, and reproduce a table like §2. You should see roughly $1\%\to 8\%\to 13\%\to 37\%\to 69\%\to 78\%$. Boehm’s repo (github.com/siboehm/SGEMM_CUDA) is the reference implementation — build it first, then rewrite from scratch without looking.
2. Profile every rung with Nsight Compute. ncu --set full ./sgemm <kernel_id>. For each rung read: global load/store efficiency (jumps to ~100% at rung 1), SM vs memory throughput % (memory-bound until register blocking, then SM climbs), achieved occupancy, shared-memory bank conflicts (l1tex__data_bank_conflicts_pipe_lsu_mem_shared, → 0 after rung 5), and spill loads/stores (must stay 0). Confirm each rung fixes the bottleneck the previous one exposed — this is the whole pedagogical point.
3. The SGEMM tuning exercise. Take the 2D kernel and sweep (BM,BN,BK,TM,TN) over a valid grid (multiples of 4, threads = BM*BN/(TM*TN) ≤ 1024, SMEM ≤ limit). Plot GFLOP/s vs each parameter; find the sweet spot; watch it move if you change $n$ or the GPU. This is a hand-rolled autotuner — the same idea Triton/CUTLASS automate.
4. Compare against CUTLASS. Run a CUTLASS SGEMM at the same shape; note it matches cuBLAS and inspect how its template parameters (ThreadblockShape, WarpShape, InstructionShape) map onto the block/warp/thread tiling you just built by hand.

Real resources

• The worklog (must-read, this lesson mirrors it): Simon Boehm, “How to Optimize a CUDA Matmul Kernel for cuBLAS-like Performance: a Worklog.” siboehm.com/articles/22/CUDA-MMM. Code: github.com/siboehm/SGEMM_CUDA.
• PMPP — Programming Massively Parallel Processors, 4th ed. (Hwu, Kirk, El Hajj). Ch. 3–5 build tiled matmul from scratch; the canonical textbook derivation of the reuse/traffic math.
• CUTLASS docs & the CUTLASS GEMM “efficient GEMM” guide. github.com/NVIDIA/cutlass — the productionized, templated version of this exact ladder (block→warp→thread tiling, pipelining, MMA atoms). Read media/docs/efficient_gemm.md.
• GPU MODE lectures (formerly CUDA MODE). github.com/gpu-mode/lectures — sessions on GEMM, CUTLASS, and profiling; pairs directly with this lesson.
• CUDA C++ Programming Guide — §on shared memory, __syncthreads, float4, and (for double buffering) cp.async/pipelines. docs.nvidia.com/cuda/cuda-c-programming-guide
• NVIDIA “CUTLASS: Fast Linear Algebra in CUDA C++” dev blog — the design rationale for the tiling hierarchy and tensor-core MMAs.
• Roofline & profiling background: 03-performance-modeling-roofline and Nsight Compute docs — for reading the per-kernel roofline chart you’ll generate in exercise 2.

---

11. What’s next

Two threads, in order:

1. 06-memory-wall-and-flashattention (next) — the payoff. FlashAttention is “tiled GEMM + online softmax”: it applies exactly the SMEM-tiling + register-blocking machinery from this lesson to the attention computation, fusing softmax so the $s\times s$ score matrix never touches HBM. Everything you just learned about tiles, __syncthreads, and reuse transfers directly.
2. Back to 04-first-cuda-kernels — reread the naive matmul now that you can name every reason it’s slow and every rung that fixes it. Then build the ladder (exercise 1) against the kernel you started from.

And the frame that made it all predictable: 03-performance-modeling-roofline — every rung was a move on the roofline (onto the memory roof, then right toward the ridge, then up to the compute roof). The roofline is why; this ladder is how.

---

Reference topic: gpu-systems-for-llms | Related concepts: gemm, memory-coalescing, shared-memory-bank-conflicts, occupancy | Filed: 2026-09-02

---

References

Footnotes

1. Simon Boehm, “How to Optimize a CUDA Matmul Kernel for cuBLAS-like Performance: a Worklog.” siboehm.com/articles/22/CUDA-MMM (2022). siboehm.com/articles/22/CUDA-MMM. Reference implementation: github.com/siboehm/SGEMM_CUDA. All GFLOP/s and %-of-cuBLAS figures, the per-rung speedups, the RTX A6000 / FP32 / 4092³ / cuBLAS≈23.2 TFLOP/s measurement setup, and the “~80% in two weekends, next 14% in four more” economics note are taken from this worklog. ↩ ↩² ↩³ ↩⁴ ↩⁵ ↩⁶ ↩⁷ ↩⁸ ↩⁹ ↩¹⁰

2. NVIDIA cuBLAS Library documentation (cublasSgemm and dispatch of shape/dtype/layout-specialized kernels). docs.nvidia.com/cuda/cublas. ↩ ↩²

3. NVIDIA CUTLASS — “Efficient GEMM in CUDA” (media/docs/efficient_gemm.md) and the templated block→warp→thread tiling + pipelining + MMA-atom hierarchy. github.com/NVIDIA/cutlass. See also NVIDIA dev blog, “CUTLASS: Fast Linear Algebra in CUDA C++.” ↩ ↩² ↩³ ↩⁴

4. NVIDIA RTX A6000 spec: peak FP32 ≈ 38.7 TFLOP/s, memory bandwidth 768 GB/s (GDDR6), Ampere GA102. NVIDIA, “NVIDIA RTX A6000 Datasheet.” nvidia.com/en-us/design-visualization/rtx-a6000/. Ampere GA10x per-block shared-memory opt-in up to 100 KB: NVIDIA, CUDA C++ Programming Guide (Compute Capability 8.x). docs.nvidia.com/cuda/cuda-c-programming-guide. ↩

5. W. Hwu, D. Kirk, I. El Hajj, Programming Massively Parallel Processors: A Hands-on Approach, 4th ed. (Morgan Kaufmann, 2022), ch. 3–5 — canonical derivation of tiled matmul and the data-reuse/HBM-traffic math. ↩