# Matrix multiplication · Kernel optimization A naive matrix-multiply kernel assigns one output element to one thread: each thread walks the corresponding row of A and column of B and accumulates a dot product straight out of global memory. On an RTX A6000 multiplying two 4092×4092 float32 matrices, that kernel manages about 309 GFLOP/s, roughly 1.3% of what cuBLAS reaches on the same GPU and the same problem. The gap is not arithmetic: neighboring threads in a warp end up reading rows of A that sit nowhere near each other in memory, so almost none of those loads can be combined into a single wide transaction. The worklog’s first kernel is the whole algorithm in fifteen lines, launched with one thread per entry of C: > Figure. sgemm_naive, from How to Optimize a CUDA Matmul Kernel Reassigning which thread owns which output element so that threads in a warp read consecutive addresses, letting the hardware coalesce those reads, already lifts throughput to about 1986.5 GFLOP/s with no other change. The next step is **tiling** (Copying a block-sized chunk of data from slow memory into fast on-chip memory once, then reusing that on-chip copy across many arithmetic operations before moving on to the next chunk.): the kernel stages a block-sized tile of A and a tile of B into shared memory once, and every thread in the block reads back out of that on-chip copy instead of returning to global memory for each partial product, which pushes throughput to about 2980.3 GFLOP/s. The heart of that shared-memory kernel is the loop every later version elaborates: stage a tile of A and a tile of B, synchronize, accumulate, synchronize, advance to the next tile along the reduction dimension: > Figure. Shared-memory tiling, from How to Optimize a CUDA Matmul Kernel The two barriers are the cooperation cost the transpose exercise introduced, now guarding both directions: the first keeps any thread from computing against a tile another thread has not finished staging, and the second keeps a fast thread from overwriting the tile with the next chunk while a slower one is still reading the current chunk. At a block size of 32 the two tiles occupy 8KB of the 48KB of shared memory a block can address on this GPU. Even with tiling, each thread is still computing exactly one output element, so most instructions in the inner loop are shared-memory loads rather than the fused multiply-adds actually doing the work; a profiler shows warps repeatedly stalling in the “Stall MIO Throttle” state, waiting on the memory pipe rather than on arithmetic. **Register blocking** (Giving each thread several output elements to accumulate in its own registers instead of one, so a single value read from shared memory gets reused across many fused multiply-adds rather than just one.) fixes that ratio directly: giving each thread a small 1D tile of outputs held in registers moves throughput to 8474.7 GFLOP/s, and a 2D tile of outputs per thread reaches 15971.7 GFLOP/s, 68.7% of cuBLAS. Vectorized memory instructions, autotuned tile sizes, and tiling at the warp level close most of what remains, reaching 21779.3 GFLOP/s, 93.7% of cuBLAS, without a tensor core in sight. 309.0GFLOP/s Naive kernel (1.3% of cuBLAS) 21779.3GFLOP/s Warptiled kernel (93.7% of cuBLAS) 23249.6GFLOP/s cuBLAS (same GPU, same problem) None of these kernels is limited by **occupancy** (The fraction of the warps an SM could run at once that are actually resident, capped by whichever per-block resource, threads, shared memory, or registers, runs out first.) in the way a first guess might suggest: the register-blocked kernel above fits only one block per SM and still reaches 66% occupancy, and pushing that number higher would not by itself close the remaining gap to cuBLAS. What separates a kernel like this from a vendor library is mostly the same tiling idea applied recursively, at the block, warp, and instruction level. Libraries such as CUTLASS describe those nested tiles as compositions of [layouts](https://learn-kernels.com/chapters/reading#matrix-multiplication), a shape paired with a stride, and talk about a matrix being “K-major” when it is stride-1 along the reduction dimension, rather than relying on the row-major and column-major vocabulary borrowed from BLAS. The case against occupancy-first reasoning is older than any GPU in this chapter. In 2008, Volkov and Demmel benchmarked dense linear algebra across four NVIDIA GPUs and built an SGEMM that sustained 58 to 60% of each chip’s peak where NVIDIA’s own CUBLAS 1.1 sustained 36 to 44%, and the design contradicted the official guidance of the era point by point: instead of many threads, shared memory as the primary storage, and long vectors, their kernel kept each 64×16 output block of C entirely in registers, staged only B’s block through shared memory, and ran short 64-element vector threads. The paper’s method is what survived: measure, then reason. Varying the thread count showed the code reaching 32, 49, 58, and 59% of peak at one to four threads per core, and on the GTX280 those four threads correspond to 25% occupancy, from which the authors conclude that one should not over-optimize for occupancy, though extremely low occupancy can also hurt. Cycle counting on the disassembled binaries located the real bound in instruction throughput: CUBLAS ran twice as many warps yet was 1.6× slower, because keeping both input blocks in shared memory forced an extra register move for every two multiply-adds, diluting the multiply-add share of its inner loop to 56% of instructions against 82% in theirs. > Figure. Where a matmul tile lives. One thread block computes one tile of C. The band of A rows and B columns it needs is staged tile by tile into shared memory, and each thread accumulates a small register tile of outputs. Shaded regions mark the data for the highlighted C tile.Illustrative numbers Source Kernel-by-kernel numbers, the optimization sequence, and both code listings from [How to Optimize a CUDA Matmul Kernel for cuBLAS-like Performance](https://siboehm.com/articles/22/CUDA-MMM). The register-blocked SGEMM, its occupancy measurements, and the CUBLAS instruction-mix analysis from Volkov and Demmel’s [Benchmarking GPUs to Tune Dense Linear Algebra](https://mc.stanford.edu/cgi-bin/images/6/65/SC08_Volkov_GPU.pdf). The layout and tiler vocabulary is from NVIDIA’s [CuTe GEMM tutorial](https://docs.nvidia.com/cutlass/latest/media/docs/cpp/cute/0x_gemm_tutorial.html). More on tile programming and layout algebra in [the reading list](https://learn-kernels.com/chapters/reading#matrix-multiplication).