optimized-GEMM

Optimized CUDA matmul kernel for cuBLAS-like performance.

Every kernel computes the same thing — C = alpha * (A @ B) + beta * C for row-major A (M x K), B (K x N), C (M x N) — and each one removes exactly one bottleneck that the previous one exposed. Read in order, the sequence is the argument for why a production SGEMM looks the way it does.

Developed and profiled on an RTX 5080 (sm_120).

How to study this

Read them in the order below. For each kernel:

1. Read the previous kernel's closing bottleneck first. Each file's header comment states what the last one was stalling on. That claim is the reason this kernel exists.
2. Diff against the previous file — git diff --no-index prev.cu next.cu. The steps are deliberately small; almost every kernel changes one idea.
3. Work out the index arithmetic by hand before trusting the comments. Pick a concrete threadIdx.x, a concrete tile, and trace which element of A, B, C it touches. This is where the actual learning is — the performance ideas are short, the indexing is what takes time.
4. Profile it. ncu --set full ./kernel and look at the stall reasons, not just the runtime. The whole ladder is only motivated if you can see the bottleneck move.

Build any single kernel with:

nvcc -O3 -arch=native -o kernel <file>.cu

sgemm_autotuned.cu additionally needs -lcublas.

Every file shares the same harness — M = N = K = 4096, alpha = 1, beta = 0, device-only allocation, and the signature (M, N, K, alpha, A, B, beta, C) — so a diff between two kernels shows only the idea that changed, never boilerplate. The individual files launch once and exit without timing or checking anything; they are meant to be read and profiled. Use bench_all.cu below to actually measure them.

---

Benchmarking

One command compiles every kernel into a single binary and prints a table against cuBLAS:

nvcc -O3 -arch=native -lcublas -o bench bench_all.cu && ./bench

bench_all.cu #includes each kernel file with its main macro'd to a unique name, so only the __global__ template survives and each kernel is launched with the exact template arguments and grid/block shape from its own main(). Nothing is duplicated — editing a kernel and rebuilding measures the edit.

On an RTX 5080 at M = N = K = 4096:

  kernel                          time (ms)      GFLOP/s  % cuBLAS  max rel err
  ------------------------------------------------------------------------------
  cuBLAS (reference)                  3.727      36878.3    100.0%            -
  1. naive                          305.838        449.4      1.2%     4.13e-06
  2. coalesced (2D block)            43.966       3126.0      8.5%     4.13e-06
  3. coalesced (1D block)            43.780       3139.3      8.5%     4.13e-06
  4. shared-memory tiling            30.744       4470.5     12.1%     4.13e-06
  5. register tiling 1D               9.811      14008.3     38.0%     4.13e-06
  6. register tiling 2D               6.522      21073.9     57.1%     4.13e-06
  7. vectorized (float4)              5.232      26268.8     71.2%     4.13e-06
  8. warp tiling                      4.607      29834.9     80.9%     4.13e-06
  9. padded smem (no conflicts)       4.673      29413.1     79.8%     4.13e-06
  10. double buffered                 4.300      31964.4     86.7%     4.13e-06
  11. final (tuned + beta0)           3.924      35023.0     95.0%     4.13e-06

Two files are deliberately not in the table. sgemm_autotuned.cu is a tuning harness whose kernel is a copy of kernel 8's (see below). sgemm_reg_pipelined.cu is a documented negative result — it adds SMEM→register double buffering on top of kernel 10 and is measurably slower. Both build and run on their own.

Options:

./bench 1024                 # square: M = N = K = 1024
./bench 8192 4096 2048       # M N K
./bench --reps 10            # timing repetitions (default 5)
./bench 2>/dev/null | ...    # progress goes to stderr, table to stdout

Reading the numbers

• max rel err is against cuBLAS, not against an exact result. Both accumulate FP32 in a different order, so ~1e-6 is agreement. ~1e-1 or worse means a real bug.
• The reference is default-math cublasSgemm — FP32 CUDA cores, not tensor cores. On this card that is ~39 TFLOP/s against a 56.3 TFLOP/s FP32 peak. For scale, CUBLAS_TF32_TENSOR_OP_MATH reaches ~60 TFLOP/s and FP16-in/FP32-accumulate via cublasGemmEx reaches ~113 TFLOP/s. None of the kernels here use tensor cores.
• % cuBLAS is a ratio, and the denominator moves with size. At 1024 neither implementation fills the GPU and the score flatters you (~94%); at 8192 the working set exceeds the 64 MB L2 and it drops (~78%). Quote absolute GFLOP/s first and never quote a percentage without the size attached. Run-to-run clock variation alone moves it ±2%.
• Sizes must be multiples of the tile dimensions. The tiled kernels have no bounds checking, so ./bench 4092 runs but returns garbage (max rel err ~1.2). This is why the repo standardizes on 4096 — it is a constraint, not a free choice.

sgemm_autotuned.cu is deliberately excluded from the table: its kernel is a copy of sgemm_warp_tiling.cu's with the same template signature, so including both is a redefinition error. It is a tuning harness, not a rung on the ladder — build and run it on its own.

---

The ladder

1. naive-kernel.cu — the baseline

One thread per output element, each walking the full K dimension out of global memory. Establishes the shape of the problem: 2*M*N*K FLOPs against M*N*K loads of A plus M*N*K loads of B. The arithmetic intensity is ~1 FLOP per byte, so this is entirely memory bound and nowhere near peak.

What to take away: the ratio above. Everything that follows is an attempt to raise it.

2. coalesced-kernel-2d.cu — coalescing

Identical body to kernel 1 with row and col swapped relative to threadIdx.x. Now consecutive threads in a warp read consecutive columns of B and write consecutive elements of C, so a 32-thread access collapses into a small number of 128-byte transactions instead of 32 separate ones.

What to take away: which index varies with threadIdx.x determines whether a warp's loads coalesce. Compare this file to kernel 1 line by line — the change is that small, and the speedup is large.

3. coalesced-kernel-1d.cu — flattening the block

Same access pattern as kernel 2, but the block is launched as a flat 32*32 array and the 2D position is recovered inside the kernel with / and %. Functionally a no-op; it exists because every kernel from here on maps a 1D thread index onto a tile by hand, and it's worth seeing that mapping in isolation before it's tangled up with tiling.

What to take away: threadIdx.x / BLOCKSIZE is the row, % BLOCKSIZE is the column — and the % has to be the fast-moving one, or you undo kernel 2.

4. smem-tiling.cu — shared memory

The first real change in arithmetic intensity. Each block loads a 32x32 tile of A and of B into shared memory, syncs, and every thread does its 32 MACs out of SMEM. Each global element is now read once per block instead of once per thread.

What to take away: the __syncthreads() pair. One after the load (nobody computes on a half-filled tile), one after the compute (nobody overwrites a tile others are still reading). Understand why removing either is a race.

5. register-tiling-1d.cu — one thread, TM outputs

The block tile becomes rectangular (BM x BN, walked over K in chunks of BK), and each thread now owns TM=8 rows of the output instead of one. A single value of B loaded into a register feeds 8 FMAs. Thread count drops to (BM*BN)/TM.

What to take away: the inner loop hoists tmpB out and reuses it TM times. That one hoist is the whole kernel — SMEM traffic per FMA drops by nearly half. Also note blockIdx.x now walks N so neighbouring blocks share rows of A and hit in L2.

6. register-tiling-2d.cu — the outer product

Extends kernel 5 to TM x TN outputs per thread: load TM values of A and TN values of B into registers, then do TM*TN FMAs. 8+8 loads for 64 FMAs. The load loops become strided because there are now fewer threads (256) than tile elements.

What to take away: the load-to-FMA ratio as a function of TM/TN, and why you can't just keep raising them — registers per thread cap occupancy. This kernel is the structural core; 7–9 are refinements on top of it.

7. sgemm_vectorized_smem_gemm.cu — 128-bit access

Same tiling, every memory access widened to float4. Same bytes moved, ~4x fewer memory instructions, which is what clears the MIO-pipe congestion showing up as mio_throttle in Nsight. Requires storing As transposed so a thread's TM values are contiguous.

What to take away: the transpose happens on the store side, in registers, because the global read has to stay coalesced along K. Trace that scatter carefully — it's the subtlest indexing in the repo, and it's also the thing that makes kernel 8 necessary.

8. sgemm_warp_tiling.cu — the warp tile

A third tiling level between block and thread. Each warp owns one contiguous WM x WN rectangle, so a warp's SMEM reads stay in a narrow address range and one regM/regN fetch feeds WMITER*WNITER sub-tiles. Tile shape came from the sweep in sgemm_autotuned.cu, not from reasoning.

What to take away: bigger block tiles are not automatically better. 128x64 moves more bytes than 128x128 and is still faster, because L2 absorbs the reuse and what actually binds is blocks resident per SM.

9. sgemm_padded_smem.cu — bank conflicts

Pads the As row stride so the transposed stores stop colliding on 8 of 32 banks. Nsight reported 4.2-way conflicts on every shared store in kernel 8.

What to take away: it makes no measurable difference, and kernel 11 later removes it for a real gain. Fixing what the profiler flags loudest is not the same as fixing what costs time — the conflicts were real and cheap.

10. sgemm_double_buffered.cu — software pipelining

Two SMEM stages. While the FMAs for tile i run, tile i+1 is already loading — B straight to shared via cp.async, A through registers because it has to be transposed. One __syncthreads() per K-tile instead of two.

What to take away: this lowers occupancy (5 blocks/SM → 3) and is still 6 points faster. That is the moment the repo stops being occupancy-bound and starts being latency-bound, and it invalidates the tuning of kernels 8 and 9 without changing a line of them.

11. sgemm_final.cu — the tuned final kernel · 95% of cuBLAS

No new technique at all. A re-tune of kernel 10 against the pipeline it now has (thread tile 4x4 → 16x4, BN 64 → 128, padding removed), a beta == 0 fast path that skips a pointless 67 MB read of C, and two compiler-facing declarations — if constexpr instead of a runtime branch, and the __launch_bounds__ min-blocks argument. The last two are worth ~10 points between them.

What to take away: read the header comment; it is the longest in the repo and the most useful. At 227 registers per thread the register allocator is the optimizer that matters, and most of the work is arranging for it to make a better choice. Also note that half of this kernel consists of undoing correct decisions from kernels 8 and 9 — nothing about those files changed, the constraint they were tuned against did.

Side branch: sgemm_reg_pipelined.cu — a negative result

Adds the second CUTLASS pipelining level (SMEM → registers) on top of kernel 10. It is slower, and the header explains why: the dotIdx loop is fully unrolled, so ptxas already schedules those loads ahead of the FMAs. Kept because knowing which plausible optimization doesn't work is worth as much as the ones that do.