Why Matmul & the Roofline
GEMM is the canonical "show me you can optimize a GPU kernel" exercise. First, the numbers you must know cold: FLOPs, arithmetic intensity, and what "% of peak" actually means.
Can you compute GEMM FLOPs, decide memory- vs compute-bound at a given size, and explain why we measure against cuBLAS? This whole course is "walk me through optimizing a matmul" — the highest-yield GPU systems question.
The Numbers
Each output is a K-length dot product = K multiplies + K adds → the factor of 2. For a 4092³ FP32 matmul: 2·4092³ ≈ 137 GFLOP; minimum traffic (read A,B,C + write C) ≈ 268 MB.
Roofline check on an A6000 (~30 TFLOP/s FP32, ~768 GB/s): compute takes 137e9/30e12 ≈ 4.5 ms; minimum memory 268e6/768e9 ≈ 0.34 ms. Compute is ~10× longer → a well-written GEMM at this size should be compute-bound. The whole course is a march to make that true.
The Optimization Ladder (Preview)
Every step raises arithmetic intensity, marching rightward toward the compute roof. Real measured numbers (Boehm, A6000, FP32, 4092³):
The Estimation Round
1. 4096³ FP32 GEMM on a 20 TFLOP/s, 1.5 TB/s GPU. Bound (well-tiled)?
2. Why measure % of cuBLAS instead of % of silicon peak?
If the roofline model is new, do GPU Performance Engineering first — this course is that model made concrete, one measured optimization at a time.
• Don't forget the factor of 2 in 2·M·N·K (multiply and add).
• "Peak" is precision-specific — FP32 and tensor-core FP16 differ >10×.
• Small/skinny matrices are latency/launch-bound, not roofline-bound (Module 8 revisits the decode-phase GEMV regime).
The Naive Kernel — 1.3% of Peak
One thread per output element. It does the exact same FLOPs as cuBLAS yet runs ~70× slower — because it loses on data movement, not math. This is the "why is this slow?" starting point.
That identical FLOPs can be 70× slower — the loss is operand re-reads from global memory plus uncoalesced access. This reframes the whole problem as data movement.
One Thread Per Output
int x = blockIdx.x*blockDim.x + threadIdx.x;
int y = blockIdx.y*blockDim.y + threadIdx.y;
if (x < M && y < N) {
float tmp = 0;
for (int i = 0; i < K; ++i)
tmp += A[x*K+i] * B[i*N+y]; // re-read from GMEM
C[x*N+y] = alpha*tmp + beta*C[x*N+y];
}Conceptual Round
1. The naive kernel does the same 137 GFLOP as cuBLAS. Why 70× slower?
2. Does adding __restrict__ / const fix the naive kernel?
• Mapping threadIdx.x → row (as above, A[x*K+i]) is the bad layout for coalescing — fixed next module.
• blockDim(32,32) = 1024 threads is the max per block; fine here but no headroom.
• The insight to carry forward: GEMM is a data-reuse problem.
Global Memory Coalescing — 1.3% → 8.5%
The single biggest "free" speedup: 6.4× with no change to FLOPs or reuse. Just remap which thread computes which output so a warp reads consecutive addresses.
That a 6× speedup can come purely from access pattern — same bytes, same math. And which index should map to threadIdx.x for a row-major output (the contiguous one).
Remap for Consecutive Addresses
Threads run in warps of 32. The naive layout made adjacent threads in a warp stride across non-consecutive addresses. Remap so consecutive threadIdx.x produce consecutive output columns:
// consecutive threadIdx.x -> consecutive col
int x = blockIdx.x*BS + (threadIdx.x / BS);
int y = blockIdx.y*BS + (threadIdx.x % BS);Spot the Bug: Which Line Kills Coalescing?
A candidate "optimized" the indexing but throughput stayed at ~15 GB/s. Click the line that fragments memory access.
int bx = blockIdx.x, by = blockIdx.y;int row = by*BS + (threadIdx.x % BS), col = bx*BS + (threadIdx.x / BS);C[row*N + col] = tmp;For row-major C[row*N + col], map threadIdx.x to col (the contiguous, innermost axis), so a warp writes consecutive addresses.
Conceptual Round
1. Same FLOPs and bytes, one 6× faster. Why?
2. Which axis maps to threadIdx.x for row-major C?
• Coalescing needs alignment, not just contiguity.
• A single swapped / vs % silently destroys it (the bug above).
• This bought 8.5% — we're still memory-bound because there's no reuse yet. Tiling is next.
Shared-Memory Tiling — 8.5% → 12.8%
The core GEMM interview idea: stage tiles of A and B in shared memory and reuse them. Explaining the tiling arithmetic — and why the gain here is modest — is a top signal.
The tiling mechanism (load tile → __syncthreads → reuse), the arithmetic-intensity improvement it buys, and the subtle point: why SMEM tiling alone only gives +50% (the bottleneck shifts to SMEM bandwidth).
Stage in Shared Memory, Reuse
for (int bk = 0; bk < K; bk += BS) {
As[tr*BS+tc] = A[tr*K + tc];
Bs[tr*BS+tc] = B[tr*N + tc];
__syncthreads(); // loads land
A += BS; B += BS*N;
for (int d = 0; d < BS; ++d)
tmp += As[tr*BS+d] * Bs[d*BS+tc];
__syncthreads(); // before overwrite
}Why Only +50%? Watch the Bottleneck Move
Interviewers love this: tiling cut GMEM traffic 32× but only gave +50%. Step through why.
Conceptual Round
1. What breaks if you delete the second __syncthreads()?
2. SMEM tiling only gave +50%. Why so little?
• SMEM per block caps occupancy — but here occupancy isn't the bottleneck (the MIO/SMEM stall is).
• Tile dim must divide K cleanly or you need boundary guards.
• Both __syncthreads are load-bearing — one guards reads-after-writes, one guards overwrites.
Register Blocking — 12.8% → 68.7%
The crux of the whole course. Each thread computes a TM×TN micro-tile of outputs held in registers, so each SMEM load feeds many FMAs. This is where the kernel becomes genuinely compute-bound.
Why a TM×TN square beats a TM×1 column (reuse in both dimensions), and what limits tile size (register pressure → spilling). This is the highest-signal GEMM optimization.
The Outer Product: TM+TN Loads → TM×TN FMAs
for (int d = 0; d < BK; ++d) {
for (i: 0..TM) regM[i] = As[...]; // TM loads
for (i: 0..TN) regN[i] = Bs[...]; // TN loads
for (m: 0..TM)
for (n: 0..TN)
acc[m][n] += regM[m]*regN[n]; // TM*TN FMAs
}The Crux Round
1. Why does a TM×TN square beat a TM×1 column for the same thread count?
2. What limits how big TM×TN can go?
• Register spilling is silent — check -Xptxas -v for spill stores/loads.
• TM/TN must be compile-time so the compiler unrolls and keeps accumulators in registers; dynamic bounds defeat it.
• Occupancy drops as regs/thread rises, but ILP compensates (this is the Volta lesson from GPU Perf Engineering Module 5).
Advanced Memory Optimizations — → 78.4%
Vectorized float4 loads, transposing the A tile, and bank-conflict avoidance. Mostly engineering practice, but the float4 alignment insight and bank-conflict fixes are real interview material.
Why casting to float4* speeds things up (it's the alignment promise, not just fewer instructions), and how to detect + fix a bank conflict. Framed as "how would you squeeze the last 10%?"
Vectorized Loads via float4
// one 128-bit transaction = 4 floats
float4 t = reinterpret_cast<float4*>(
&A[row*K + col*4])[0];
// transposing store into As
As[(col*4+0)*BM + row] = t.x;
As[(col*4+1)*BM + row] = t.y; // ...z, .wfloat* is 16-byte aligned, so it won't emit 128-bit loads.float4* cast asserts alignment → unlocks LDG.E.128, quartering memory instructions.LDS.128). Net: 15972 → 18237 GFLOP/s.Match the Optimization to Its Mechanism
reinterpret_cast<float4*>pad SMEM leading dimtranspose A tile in SMEMdouble bufferingConceptual Round
1. Why does casting to float4* help beyond fewer instructions?
2. How do you detect and fix a 2-way bank conflict?
• float4 needs 16-byte-aligned base pointers and tile dims divisible by 4 — a hard constraint on legal BM/BN/BK.
• Transposing A adds SMEM store complexity and can introduce bank conflicts if not padded.
• Double buffering doubles SMEM usage, pressuring occupancy.
Tensor Cores — Why They Dominate Peak
Dedicated matrix-multiply units on a separate, ~17× higher FLOP ceiling. Low-precision inputs, FP32 accumulate. Rising fast in interviews as FP8/BF16 take over.
Why FP16 inputs but FP32 accumulator, and that a WMMA call is executed by the whole warp cooperatively (not per-thread). The peak-FLOP gap over CUDA cores is the "why bother" answer.
A Separate, Much Higher Ceiling
Tensor Cores compute a small matrix-multiply-accumulate in a few instructions. H100 peaks (why libraries target them):
using namespace nvcuda;
wmma::fragment<matrix_a, 16,16,16, half, row_major> a;
wmma::fragment<matrix_b, 16,16,16, half, col_major> b;
wmma::fragment<accumulator, 16,16,16, float> acc;
wmma::load_matrix_sync(a, A+..., lda);
wmma::mma_sync(acc, a, b, acc); // D = A*B + CConceptual Round
1. Why FP16 inputs but FP32 accumulator?
2. WMMA is written per-thread — what actually executes it?
• The WMMA fragment layout is opaque — don't index it manually; use load/store/mma.
• Every active lane must participate — no divergence.
• BF16 keeps FP32's exponent range (safer than FP16 for training); TF32 truncates the mantissa — wrong for full-precision HPC.
• Below WMMA, mma.sync (PTX) + ldmatrix and Hopper's wgmma give finer control — what libraries use.
Reaching cuBLAS & the Skinny-GEMM Regime
The last ~15% comes from autotuning + warptiling — exactly what CUTLASS codifies. Plus the regime that dominates LLM inference: skinny/decode-phase GEMV, which is memory-bound.
Why there's no universal tile size (autotune per GPU), how CUTLASS uses software pipelining to hide latency despite low occupancy, and — the LLM link — why decode-phase matmul is memory-bound.
Warptiling + Autotuning → 93.7%
BM,BN,BK,TM,TN. A6000 optimum 128/128/16/8/8; A100's was 64/64/16/4/4. Tiles are not portable.
A hierarchical tiled loop nest mirroring hardware: threadblock tile (GMEM→SMEM), warp tile (SMEM→registers via mma.sync), thread/instruction tile, plus software pipelining (double-buffered SMEM + register fragments), threadblock swizzle for L2, split-K, and — on Hopper — warp specialization (producer TMA warps + consumer wgmma warps).
The Regime That Runs LLM Inference
Everything so far assumed large square matrices (compute-bound). But LLM decode multiplies a weight matrix by a single token vector each step — a tall-skinny GEMV. That flips the roofline:
AI = N/8, grows with N → compute-bound → tensor cores shine. This is training / prefill.
Batch=1 → low AI → memory-bound: you re-read all the weights to produce one token. This is why decode is bandwidth-bound and why quantization (fewer weight bytes) helps.
This ties directly to vLLM & Distributed Kernels, where "why is decode memory-bound?" is a core question. Same roofline, different regime.
The Staff-Level Round
1. Why isn't there one universal tile config?
2. Accumulators eat most registers, killing occupancy. How does CUTLASS still hide latency?
Naive (1.3%) → coalescing (8.5%) → SMEM tiling (12.8%) → register blocking (68.7%) → vectorized (78.4%) → tensor cores → warptiling + autotune (93.7%). The through-line, and the answer to "walk me through optimizing a matmul": every step raises arithmetic intensity to march up the roofline — quantify the reuse at each stage.
For standard dense GEMM, call cuBLAS/CUTLASS — matching them by hand took the reference author ~6 weekends. Write your own only for fusion (epilogue/activation), unusual shapes/dtypes, or learning. And know: swizzling for L2 only helps if L2 hit rate is low (a documented negative result), and the last 14% costs far more effort than the first 80%.