The SGEMM Ladder

Canonical source: Simon Boehm, How to Optimize a CUDA Matmul Kernel for cuBLAS-like Performance — https://siboehm.com/articles/22/CUDA-MMM Repo: https://github.com/siboehm/SGEMM_CUDA

This blog post is the single most valuable thing you will read in Phase 2. It is a 10-step ladder from a naive matmul to ~93–95% of cuBLAS on an A6000 in fp32. Every step introduces exactly one optimization, isolates its effect with a measurement, and inspects the Nsight numbers. You are going to walk this ladder yourself.

Why fp32 and not BF16 tensor cores? Because the ladder isolates memory-hierarchy and tiling skills without hiding them behind tensor-core wall-clock. Once you can get 90% of cuBLAS on fp32, you know your tiling is right and adding BF16 tensor cores is one more layer, not a rewrite.

Ladder overview — what you will build

Step

Kernel

Key idea

Expected % of cuBLAS (fp32, A6000-class)

Ncu metric to verify

1

Naïve

1 thread per output element, all reads from GMEM

~1–2%

dram__throughput.avg.pct_of_peak_sustained_elapsed low, sm__throughput low

2

Global memory coalescing

Reorder threads so warp reads contiguous 128 B

~8–10%

l1tex__t_sectors_pipe_lsu_mem_global_op_ld.sum drops; L1 hit rate up

3

Shared memory blocking

Tile A, B into SMEM (32×32)

~13–15%

l1tex__data_pipe_lsu_wavefronts_mem_shared.sum high; GMEM traffic reduced ~32×

4

1D block tiling

Each thread computes multiple output elements

~30–35%

Register file usage up; arithmetic intensity up

5

2D block tiling

Each thread computes an 8×8 tile in registers

~55–65%

Even higher register usage, ~4 KB/thread; sm__inst_executed_pipe_fp32 up

6

Vectorized SMEM loads

Use float4 for SMEM → RF loads

~65–75%

Fewer LDS instructions; SMEM bandwidth utilization higher

7

Auto-tuning

Sweep BLOCK, BK, TM, TN, WM, WN

~80–85%

Occupancy vs tile size sweet spot

8

Warp tiling

Explicit warp-level tile alongside block/thread tile

~88–92%

Warp partitioning — 4 warps × explicit shape

9

Double buffering

Overlap next-tile GMEM load with current-tile compute

~93–95%

Fewer stall cycles (stall_long_scoreboard down)

10

(Advanced) Tensor cores

wmma / mma.sync in fp16 or TF32

> cuBLAS_fp32 (different regime)

sm__inst_executed_pipe_tensor_op.sum non-zero

The % numbers are approximate and hardware-dependent — Boehm ran on A6000. On RTX 4090 or L4 you’ll see slightly different ratios, on H100 fp32 tensor performance is very different, but the shape of the curve is identical.

Step-by-step commentary

Step 1 — Naïve kernel

__global__ void sgemm_naive(int M, int N, int K, float alpha,
                            const float* A, const float* B, float beta, float* C) {
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    if (x < M && y < N) {
        float acc = 0.0f;
        for (int i = 0; i < K; ++i) acc += A[x*K + i] * B[i*N + y];
        C[x*N + y] = alpha*acc + beta*C[x*N + y];
    }
}

Why it’s slow: each output does 2K flops but 2K global loads. Intensity ≈ 1 flop/byte-ish. Far below ridge. HBM-bound.

Ncu check: dram__throughput.avg.pct_of_peak_sustained_elapsed will be low (< 20%) and sm__throughput low — you’re neither memory nor compute bound, you’re stalled on latency.

Step 2 — GMEM coalescing

Swap the mapping: x = blockIdx.x*blockDim.x + (threadIdx.x % 32) and let the 32 threads in a warp span the contiguous dimension of C. Now a warp issues one 128 B load, not 32 scattered loads.

Speedup: ~7×. This is often the single biggest “free” win in kernel work — and the reason you spend time on Ch. 6 of PMPP.

Step 3 — Shared memory blocking

Load a 32×32 tile of A and B into SMEM once per block. All threads in the block reuse those tiles across the inner-K dimension.

Consequence: GMEM traffic drops by ~32× (tile size). Intensity goes up. You are now within roofline touch.

Step 4 — 1D block tiling

One thread computes 8 output elements (a 1×8 strip). Reduces load pressure, increases register reuse.

Step 5 — 2D block tiling

Each thread computes an 8×8 tile in registers. This is the big jump. Arithmetic intensity per thread now ≈ 64 flops per pair of loads. You’re crossing into compute-bound territory.

Watch: register pressure. If you use too many registers, occupancy tanks. Check sm__warps_active.avg.pct_of_peak_sustained_active and register-per-thread count.

Step 6 — Vectorized SMEM loads

Use float4 (reinterpret_cast<float4*>) for SMEM → RF loads. Cuts LDS instruction count by 4.

Step 7 — Auto-tuning

Expose block size (BM, BN, BK), thread tile (TM, TN), and sweep. Boehm shipped a small script; write your own. Autotune once per (M, N, K, dtype, GPU) tuple.

Step 8 — Warp tiling

Explicit warp-tile in the hierarchy: block-tile (SMEM) → warp-tile (register + tensor core prep) → thread-tile (registers). This is exactly the CUTLASS layered abstraction — you are hand-rolling what CUTLASS gives you for free.

Step 9 — Double buffering

While computing on tile k, load tile k+1 into a second SMEM buffer (__pipeline_memcpy_async on Ampere+, TMA on Hopper+). Overlap latency with compute.

Ncu check: stall_long_scoreboard (memory-stall) drops.

Step 10 — Tensor cores (bonus)

Switch A, B to fp16 (or TF32), issue wmma fragments. This changes regime entirely — you’re no longer fp32; you’re using the tensor cores for real. On H100+ this means going through wgmma.

Modern companion resources

  • Lei Mao’s blog — https://leimao.github.io — has parallel writeups on matmul optimization with slightly different framing. Read after Boehm.

  • Salykova’s article on beating cuBLAS on consumer HW (RTX 4090) — more recent than Boehm, worth the second pass.

  • Curated meta-list: https://www.abhik.ai/articles/best-resources-cuda-matmul-optimization — aggregates the good writeups.

  • CUTLASS gemm walkthrough — once you’ve done Boehm’s 10 steps by hand, read the equivalent CUTLASS kernel to see how the same ideas are expressed at industrial scale. Start with include/cutlass/gemm/collective/* for Hopper.

What “done” looks like

You can:

  • Draw the memory hierarchy of your kernel from block tile → warp tile → thread tile, labeling which lives in SMEM vs registers.

  • Explain, kernel-by-kernel, which Nsight metric moved and why.

  • Hit ≥ 70–80% of cuBLAS fp32 on your target GPU. If you can hit 90%, you are in genuinely competent territory.

  • Reproduce Boehm’s kernel 8 without looking at the blog.

If you cannot, do not proceed to Triton yet — fix your understanding first.

Common failure modes

  • Optimizing without profiling. If a step doesn’t show up in ncu, you probably didn’t change what you thought you changed.

  • Skipping steps. Each step teaches an isolated lesson. Skipping step 4 (1D tiling) makes step 5 (2D tiling) confusing.

  • Chasing the last 5%. Getting from 90% → 95% of cuBLAS involves subtle ordering + PTX-level tricks and is not a good ROI unless you’re targeting kernel-writing as a job title. 80% is fine for phase-2 exit.

  • Overfitting to one shape. Your ladder should measure at M=N=K in {512, 1024, 2048, 4096} and small non-square shapes. Fast at 4096 but 5× slower at 512 is a real problem for real inference (batch=1 decode).

Path to Blackwell / Hopper GEMM

After you finish Boehm’s fp32 ladder, the next question is BF16/FP8 GEMM on real modern hardware. Do not roll your own Hopper GEMM from scratch — read CUTLASS’s Hopper kernels and Colfax Research’s technical notes:

References