Chapter 11 · Section 11.5 · 60 min

Tiled matrix multiplication

Matrix multiplication is the operation a GPU is built for. It is also the easiest one to write badly. In a well-known worklog, Simon Boehm (2022) wrote the obvious CUDA matmul, one thread per output element, and measured it on an A6000 at 1.3 percent of the speed of NVIDIA's cuBLAS library on the same $4092\times4092$ float32 problem. Fixing the memory access pattern took it to 8.5 percent, caching tiles in shared memory to 12.8, having each thread compute a short column of outputs to 36.5, a small square tile of outputs to 69, and a series of further refinements to 94.

The arithmetic was the same in every version, $2n^3$ FLOPs. What changed was how many times each number was fetched from memory. Each output element is a reduction of the kind §11.4 studied, the dot product of a row of the first matrix with a column of the second, so a $4092\times4092$ matmul is 16.7 million reductions running at once. Unlike the sums of §11.4, they overlap almost completely in their inputs, since every row of the first matrix feeds a whole row of outputs and every column of the second a whole column. Everything in this section is about making those reductions share what they load, and its measure of success is the number of loads, which we derive at every step. We start with the naive kernel and its traffic, derive exactly how much tiling with shared memory saves, push the same idea down to registers, prove that no schedule can do fundamentally better, and end with the tensor-core instructions that real kernels are built from. A NumPy simulation counts every read along the way, and the exercises ask you to write it yourself, a common interview problem.

Throughout, $\mathbf{C} = \mathbf{A}\mathbf{B}$ with $\mathbf{A} \in \R^{M\times K}$ , $\mathbf{B} \in \R^{K\times N}$ , and $\mathbf{C} \in \R^{M\times N}$ , all row-major, with $s$ bytes per element. (These are the standard names in GEMM libraries, so in this section $\mathbf{A}$ is a matrix operand, not attention weights, and $K$ is the inner dimension, not a number of classes.)

The naive kernel and its traffic

The obvious kernel gives each thread one element $C_{ij}$ and lets it compute a dot product of length $K$ .

// One thread per element of C. A is M x K, B is K x N, C is M x N, all row-major float32.
// Launch: dim3 block(32, 8); dim3 grid((N + 31) / 32, (M + 7) / 8);
__global__ void matmul_naive(const float* A, const float* B, float* C, int M, int N, int K) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    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];     // two global loads per multiply-add
        C[row * N + col] = acc;
    }
}

The mapping already follows §11.3. x is the column, so the 32 threads of a warp read 32 consecutive elements of a row of $\mathbf{B}$ (coalesced) and the same element of $\mathbf{A}$ (a broadcast). Map x to rows instead and the loads of $\mathbf{A}$ become strided by $K$ , which is most of the difference between Boehm's first two kernels.

Now count. Each thread loads $K$ elements of $\mathbf{A}$ and $K$ of $\mathbf{B}$ , so the kernel issues

Q_{\text{naive}} = 2MNK \text{ element loads}

against the $MK + KN$ elements that actually exist. Every element of $\mathbf{A}$ is loaded $N$ times, once by each thread in its row of $\mathbf{C}$ , and every element of $\mathbf{B}$ is loaded $M$ times. The arithmetic intensity, if every load went to HBM, would be $2MNK/(2MNK\,s) = 1/s$ , a quarter of a FLOP per byte in float32, the bottom of the roofline (§11.2). For $n = 4096$ in float32 that is about 550 GB of loads for 134 MB of distinct input data.

The caches soften this in practice, since a warp's neighbors often want the same row of $\mathbf{A}$ and L1 and L2 catch some of the repeats. But caches are small and shared by thousands of threads, and the kernel gives them no help. The fix is to make the reuse explicit.

Tiling with shared memory

Picture the output divided into $b\times b$ tiles, one per thread block, and one thread per output element. Tile $(I, J)$ of $\mathbf{C}$ needs rows $I$ of $\mathbf{A}$ (a $b \times K$ strip) and columns $J$ of $\mathbf{B}$ (a $K\times b$ strip). Cut both strips into $b\times b$ tiles along $K$ . Then

\mathbf{C}_{IJ} = \sum_{p=0}^{K/b - 1} \mathbf{A}_{Ip}\,\mathbf{B}_{pJ},

a sum of products of small tiles. The block computes it in $K/b$ phases. In each phase, its $b^2$ threads cooperatively load one tile of $\mathbf{A}$ and one tile of $\mathbf{B}$ into shared memory, one element each, wait at a barrier, then each thread does $b$ multiply-adds reading only shared memory, and waits again before the next phase overwrites the tiles.

In plain words, a block brings two small tiles on chip, squeezes $b$ uses out of every element, and only then fetches the next two.

Click a tile of $\mathbf{C}$ in the widget and step through its phases to see which tiles of $\mathbf{A}$ and $\mathbf{B}$ each phase loads, and watch the global traffic counter against the naive kernel's.

#define TILE 32
// Launch: dim3 block(TILE, TILE); dim3 grid((N + TILE - 1) / TILE, (M + TILE - 1) / TILE);
__global__ void matmul_tiled(const float* A, const float* B, float* C, int M, int N, int K) {
    __shared__ float As[TILE][TILE];
    __shared__ float Bs[TILE][TILE];
    int tx = threadIdx.x, ty = threadIdx.y;
    int row = blockIdx.y * TILE + ty, col = blockIdx.x * TILE + tx;
    float acc = 0.0f;
    for (int p = 0; p < (K + TILE - 1) / TILE; ++p) {
        int ka = p * TILE + tx, kb = p * TILE + ty;
        As[ty][tx] = (row < M && ka < K) ? A[row * K + ka] : 0.0f;   // coalesced, zero-padded edges
        Bs[ty][tx] = (kb < K && col < N) ? B[kb * N + col] : 0.0f;
        __syncthreads();                                              // tiles are complete
        for (int k = 0; k < TILE; ++k)
            acc += As[ty][k] * Bs[k][tx];                             // reads only shared memory
        __syncthreads();                                              // done before overwriting
    }
    if (row < M && col < N) C[row * N + col] = acc;
}

Both barriers are needed. The first guarantees the tiles are complete before anyone reads them, and the second that everyone has finished reading before anyone overwrites them. In the inner loop, As[ty][k] is the same word for every thread of a warp (a broadcast) and Bs[k][tx] is a row read across 32 banks, so there are no bank conflicts.

How much traffic tiling saves

Count the global loads. Each phase loads $2b^2$ elements, there are $K/b$ phases per block, and $MN/b^2$ blocks, so

Q_{\text{tiled}} = \frac{MN}{b^2}\cdot\frac{K}{b}\cdot 2b^2 = \frac{2MNK}{b}.

In plain words, tiling divides global memory traffic by the tile size $b$ . Each element brought on chip is used by $b$ threads instead of one.

The same count works for rectangular tiles of $b_M \times b_N$ outputs with an inner tile length $b_K$ . Each phase loads $b_M b_K$ elements of $\mathbf{A}$ and $b_K b_N$ of $\mathbf{B}$ , there are $K/b_K$ phases and $MN/(b_Mb_N)$ blocks, and $b_K$ cancels,

Q_{\text{tiled}} = MNK\left(\frac{1}{b_N} + \frac{1}{b_M}\right), \qquad I_{\text{block}} = \frac{2MNK}{s\,Q_{\text{tiled}}} = \frac{2\,b_Mb_N}{s\,(b_M + b_N)} .

The first term is the traffic for $\mathbf{A}$ , read once per column of blocks, and the second for $\mathbf{B}$ , read once per row of blocks. Only the output tile's shape matters, and square tiles are best for a given area. The shared memory each block needs is $s\,b_K(b_M + b_N)$ , which is why $b_K$ is kept small (32 or 64) while $b_M$ and $b_N$ grow to 128 or 256.

Put numbers on it, because the result is surprising. With $128\times128$ output tiles in bfloat16, $I_{\text{block}} = 2\cdot128\cdot128/(2\cdot256) = 64$ FLOPs per byte, and even $128\times256$ tiles give only 85, both well below the H100's ridge point of 295. If every tile load went to HBM, a large GEMM would be memory-bound. It is not, because the blocks running at the same time share panels of $\mathbf{A}$ and $\mathbf{B}$ , and most of their tile loads hit the 50 MB L2 cache, so the true HBM traffic sits much closer to the minimum $s(MK + KN + MN)$ . That is why the order in which blocks are launched matters (the grouped ordering of §11.6), and why Hopper added thread block clusters, in which the TMA can multicast one tile into the shared memory of several blocks at once, cutting L2 traffic as well.

When the dimensions are not multiples of the tile, the edge tiles are partial. A careful kernel does not load out-of-range elements (our CUDA kernel writes zeros instead), so the exact count is that each of the $\lceil N/b\rceil$ columns of blocks reads all of $\mathbf{A}$ , and each of the $\lceil M/b\rceil$ rows of blocks reads all of $\mathbf{B}$ ,

Q_{\text{tiled}} = MK\left\lceil\frac{N}{b}\right\rceil + KN\left\lceil\frac{M}{b}\right\rceil .

Interview

Global traffic of naive and tiled matmul

Implement matmul_traffic(M, N, K, b) returning a dictionary of global element reads (not bytes) with keys "naive" (one thread per output, every operand from global memory), "tiled" (square $b\times b$ output tiles that read only in-bounds elements, as derived above), and "minimum" (each input element read once). Use exact integer arithmetic.

def matmul_traffic(M, N, K, b):
    return {"naive": M * N * K, "tiled": 2 * M * N * K // b, "minimum": M * K + K * N}   # placeholder

Hidden tests, hints, and AI help, from $1.99 a month.

A tiled matmul, simulated

The simulation follows the kernel's loop structure exactly. For each output tile and each phase it copies the in-bounds parts of one tile of $\mathbf{A}$ and one of $\mathbf{B}$ into "shared memory", counts those reads, and accumulates their product, counting the shared reads of the inner loop (each output element reads one row of the $\mathbf{A}$ tile and one column of the $\mathbf{B}$ tile).

Python · cell 1

The global reads fall by exactly the tile size (with $b = 32$ , $80/32$ is not an integer, so the edge tiles make the factor a little less than 32), while the shared reads stay at $2MNK$ for every $b$ . Tiling moved the traffic from HBM to shared memory. It did not remove it.

Interview

A tiled matmul that reports its traffic

Write tiled_matmul(A, B, b) that computes $\mathbf{A}\mathbf{B}$ with explicit loops over $b\times b$ output tiles and over phases along $K$ , as a thread block would, and returns (C, stats) with stats["global_reads"], the number of elements copied from A and B into tile buffers (in-bounds elements only), and stats["shared_reads"], the number of element reads from the tile buffers by the inner products (two per multiply-add). Inside a phase you may multiply the two small tile buffers with @. It must work for any shapes, including ones that are not multiples of $b$ .

import numpy as np

def tiled_matmul(A, B, b):
    C = A @ B                                       # placeholder: right answer, no tiling, no counts
    return C, {"global_reads": 0, "shared_reads": 0}

Hidden tests, hints, and AI help, from $1.99 a month.

Register tiling

Tiling made HBM traffic manageable, but look at the new bottleneck. The inner loop does one multiply-add per two shared memory reads. On an H100, shared memory delivers 128 bytes per cycle per SM (§11.1), 32 float32 values, enough for 16 multiply-adds per cycle. The SM's 128 float32 lanes can do 128. So the shared-tiled kernel runs at most at about one eighth of the float32 peak, which is why Boehm's shared memory version was still far below cuBLAS.

The fix is the same idea one level down. Give each thread a small $t_M\times t_N$ tile of outputs, held in registers. At each step $k$ of the inner loop, the thread loads $t_M$ values from a column of the $\mathbf{A}$ tile and $t_N$ from a row of the $\mathbf{B}$ tile into registers, and performs all $t_Mt_N$ multiply-adds of their outer product.

// Inner loop of a register-tiled kernel. Each thread owns an 8 x 8 accumulator tile.
float acc[8][8] = {0.0f};
float a_reg[8], b_reg[8];
for (int k = 0; k < BK; ++k) {
    for (int i = 0; i < 8; ++i) a_reg[i] = As[k][thread_row * 8 + i];   // As stored transposed (BK x BM)
    for (int j = 0; j < 8; ++j) b_reg[j] = Bs[k][thread_col * 8 + j];
    for (int i = 0; i < 8; ++i)
        for (int j = 0; j < 8; ++j)
            acc[i][j] += a_reg[i] * b_reg[j];                              // 64 FMAs from 16 loads
}

Shared memory reads per multiply-add fall from 2 to $(t_M + t_N)/(t_Mt_N)$ , which for $8\times8$ is $16/64 = 1/4$ , an eightfold reduction that exactly closes the gap above. The price is registers. The accumulators alone take $t_Mt_N = 64$ per thread, plus 16 for the operands and some for addresses, so such kernels use well over 100 registers per thread and run at low occupancy, by design (the pitfall in §11.1). This is why the occupancy calculator is a guide and not a goal.

In plain words, reuse happens at every level. A block reuses each element it loads from HBM $b$ times, and a thread reuses each value it loads from shared memory $t_N$ times (for $\mathbf{A}$ ) or $t_M$ times (for $\mathbf{B}$ ).

level	data held	typical size	reuse per element brought in
block tile in shared memory	$b_M\times b_K$ of $\mathbf{A}$ , $b_K\times b_N$ of $\mathbf{B}$	$128\times32$ and $32\times128$	$b_N$ for $\mathbf{A}$ , $b_M$ for $\mathbf{B}$
warp tile	a $64\times64$ piece of the block's output
thread tile in registers	$t_M\times t_N$ accumulators	$8\times8$	$t_N$ for $\mathbf{A}$ , $t_M$ for $\mathbf{B}$
tensor core instruction	a $16\times8\times16$ product (below)

Python · cell 2

The table hides an exact count worth deriving once. With $b\times b$ block tiles and $t\times t$ register tiles, a block has $(b/t)^2$ threads, and at each of the $K$ steps of the inner loop each thread loads $t$ values of $\mathbf{A}$ and $t$ of $\mathbf{B}$ from shared memory, so the whole matmul reads $K \cdot (M/t)(N/t) \cdot 2t = 2MNK/t$ elements from shared memory, while the global reads stay at $2MNK/b$ . In plain words, the block tile divides HBM traffic by $b$ , and the register tile divides shared-memory traffic by $t$ , independently.

Interview

A two-level tiled matmul

Extend Exercise 11.5b to two levels. Write two_level_matmul(A, B, b, t) that computes $\mathbf{A}\mathbf{B}$ with $b\times b$ block tiles of the output, phases of width $b$ along $K$ , and, inside each block tile, $t\times t$ register tiles owned by one thread (with $t$ dividing $b$ ). Each thread, at each step $k$ of a phase, loads its slice of column $k$ of the $\mathbf{A}$ tile and its slice of row $k$ of the $\mathbf{B}$ tile into registers and adds their outer product to its accumulators. Return (C, stats) with stats["global_reads"] (elements copied from A and B into the block's tile buffers, in-bounds only, as in 11.5b) and stats["shared_reads"] (elements loaded from the tile buffers into registers). At ragged edges a thread's register tile is clipped to the matrix, and it loads only the rows and columns that exist.

import numpy as np

def two_level_matmul(A, B, b, t):
    M, K = A.shape
    N = B.shape[1]
    g = M * K * -(-N // b) + K * N * -(-M // b)
    return A @ B, {"global_reads": g, "shared_reads": 2 * M * N * K}   # placeholder: no register reuse

Hidden tests, hints, and AI help, from $1.99 a month.

No schedule can do much better

Is tiling just a good trick, or is it optimal? Hong and Kung (1981) answered this for the classical $O(n^3)$ algorithm. With a fast memory that holds $S$ elements, any order of computing the $MNK$ products must move $\Omega(MNK/\sqrt{S})$ elements between fast and slow memory. This is a theorem.

Tiling meets the bound. Three $b\times b$ tiles fit in fast memory when $3b^2 \le S$ , so $b \approx \sqrt{S/3}$ , and the traffic $2MNK/b \approx 2\sqrt{3}\,MNK/\sqrt{S}$ is within a constant factor of the lower bound. The square-root dependence is why bigger on-chip memories matter so much. Doubling shared memory cuts the minimum traffic of every matmul by only $\sqrt 2$ , but that is still a 30 percent cut. FlashAttention's IO analysis has the same $\sqrt{S}$ structure (§12.1).

∂ The lower bound, through the Loomis-Whitney inequality

We follow the argument of Irony, Toledo & Tiskin (2004), which makes Hong and Kung's bound explicit. The computation consists of the $MNK$ multiply-adds $C_{ij} \mathrel{+}= A_{ik}B_{kj}$ , one for each triple $(i, j, k)$ . Split the execution into consecutive segments that each contain exactly $S$ transfers between fast and slow memory (the last segment may contain fewer).

How many products can one segment do? During a segment, the elements of $\mathbf{A}$ it can use are those in fast memory when the segment starts (at most $S$ ) or loaded during it (at most $S$ ), so at most $2S$ . The same holds for $\mathbf{B}$ . Each partial sum of $\mathbf{C}$ that the segment touches must survive it, either by staying in fast memory at the end (at most $S$ of them) or by being written out during the segment (at most $S$ ), so at most $2S$ of them are involved.

Counting triples by their shadows. Let $V$ be the set of triples $(i, j, k)$ the segment performs. Its projections onto the $(i,k)$ , $(k,j)$ , and $(i,j)$ planes are the elements of $\mathbf{A}$ , $\mathbf{B}$ , and $\mathbf{C}$ it used. The Loomis-Whitney inequality says a finite set of lattice points in three dimensions is at most the square root of the product of the sizes of its three projections,

|V| \le \sqrt{|V_{ik}|\,|V_{kj}|\,|V_{ij}|} \le \sqrt{(2S)^3} = (2S)^{3/2}.

Adding up. All $MNK$ triples must be performed, so there are at least $MNK/(2S)^{3/2}$ segments, and all but the last contain $S$ transfers. The total traffic is therefore

Q \;\ge\; S\left(\frac{MNK}{(2S)^{3/2}} - 1\right) = \frac{MNK}{2\sqrt{2}\,\sqrt{S}} - S.

For large matrices the $-S$ is negligible, and $Q = \Omega(MNK/\sqrt{S})$ . In plain words, fast memory of size $S$ lets each transferred element take part in at most about $\sqrt S$ products, so the products force the traffic. Note that the bound is for the classical algorithm. Strassen-like algorithms do fewer multiplications and have a different bound.

Tensor cores and MMA shapes

On a modern GPU the innermost level is not a thread doing FMAs. It is a warp issuing matrix multiply-accumulate (MMA) instructions to the tensor cores of §11.1. The instruction names its shape as $m\times n\times k$ , a product of an $m\times k$ tile by a $k\times n$ tile added to an $m\times n$ accumulator.

Ampere and later, warp-level mma.sync.aligned.m16n8k16 for bfloat16 or float16 inputs with float32 accumulation. One warp computes a $16\times8$ output from a $16\times16$ tile of $\mathbf{A}$ and a $16\times8$ tile of $\mathbf{B}$ . The operands are spread across the warp's registers in a fixed layout, so that each of the 32 threads holds 8 elements of the $\mathbf{A}$ tile, 4 of the $\mathbf{B}$ tile, and 4 of the accumulator.
Hopper, warpgroup-level wgmma.mma_async. Four warps (a warpgroup of 128 threads) issue one instruction of shape $64\times N\times16$ for 16-bit inputs, with $N$ any multiple of 8 up to 256, and $k = 32$ for FP8. The $\mathbf{B}$ operand is read directly from shared memory, $\mathbf{A}$ from registers or shared memory, and the instruction is asynchronous, so the warps can issue loads for the next tile while the tensor cores work.

So a fast GEMM is a hierarchy of tiles. The grid tiles $\mathbf{C}$ into block tiles, a block's warps split its tile into warp tiles, and each warp tile is a grid of MMA instructions. The simplest way to use tensor cores from CUDA C++ is the WMMA API, which hides the register layout behind opaque fragments. The kernel below gives each warp one $16\times16$ output tile. It is correct but naive (operands come straight from global memory), and real kernels add every trick of this section around the same instructions.

#include <mma.h>
#include <cuda_fp16.h>
using namespace nvcuda;

// C (M x N, float) = A (M x K, half) * B (K x N, half), all row-major, M, N, K multiples of 16.
// One warp per 16 x 16 tile of C. Launch with blockDim.x a multiple of 32.
__global__ void wmma_gemm(const half* A, const half* B, float* C, int M, int N, int K) {
    int warp = (blockIdx.x * blockDim.x + threadIdx.x) / 32;
    int tiles_n = N / 16;
    int tm = warp / tiles_n, tn = warp % tiles_n;
    if (tm >= M / 16) return;                                   // the whole warp leaves together

    wmma::fragment<wmma::matrix_a, 16, 16, 16, half, wmma::row_major> a;
    wmma::fragment<wmma::matrix_b, 16, 16, 16, half, wmma::row_major> b;
    wmma::fragment<wmma::accumulator, 16, 16, 16, float> acc;
    wmma::fill_fragment(acc, 0.0f);
    for (int k = 0; k < K; k += 16) {
        wmma::load_matrix_sync(a, A + (tm * 16) * K + k, K);    // leading dimension K
        wmma::load_matrix_sync(b, B + k * N + tn * 16, N);      // leading dimension N
        wmma::mma_sync(acc, a, b, acc);                         // acc += a @ b on the tensor cores
    }
    wmma::store_matrix_sync(C + (tm * 16) * N + tn * 16, acc, N, wmma::mem_row_major);
}

The cell decomposes a matmul into $16\times8\times16$ MMA instructions exactly as the hardware would, padding the edges with zeros, and checks the result. It also shows why the alignment advice of §11.1 is not a superstition. Every dimension is rounded up to the instruction shape, so a dimension of 17 costs as much as 32.

Python · cell 3

Warm-up

Count tensor core instructions

Implement mma_count(M, N, K, m=16, n=8, k=16) returning (instructions, utilization), the number of $m\times n\times k$ MMA instructions needed when each dimension is padded up to a multiple of the instruction shape, and the fraction of the instructions' multiply-adds that are real, $MNK$ divided by instructions * m * n * k.

def mma_count(M, N, K, m=16, n=8, k=16):
    instructions = (M * N * K) // (m * n * k)          # placeholder: ignores padding
    return instructions, 1.0

Hidden tests, hints, and AI help, from $1.99 a month.

The rest of the road to cuBLAS

Tiling at three levels is the core of every fast GEMM. The remaining distance is a list of refinements, each fixing one specific inefficiency, and knowing which fixes what is what interviewers probe.

technique	what it fixes
coalesced global loads, `x` on columns	wasted sectors (§11.3)
shared memory tiles, $b_M\times b_N$	HBM traffic, divided by the tile size
register tiles, $t_M\times t_N$ per thread	shared memory bandwidth, divided by the register tile size
vectorized 16-byte loads (`float4`)	instruction count of loads
swizzled or padded shared layouts	bank conflicts on tile reads
double buffering, multi-stage pipelines (`cp.async`, TMA)	the load latency of each phase, by overlapping the next tile's load with this tile's math
warp tiling and tensor core MMA	the FP32 lanes' 15 times lower throughput
split-K, stream-K, persistent kernels	wave quantization and too few tiles for skinny shapes
grouped (swizzled) block order	L2 hit rate between blocks that share panels (§11.6)
epilogue fusion (bias, activation, cast)	an extra round trip of $\mathbf{C}$ through HBM (§11.7)

Split-K deserves a sentence, because it is §11.4 again. When $M$ and $N$ give too few output tiles to fill the GPU but $K$ is long, several blocks each sum one slice of $K$ for the same output tile, and their partial tiles are combined by a second pass or by atomic adds, with the same trade between an extra launch and nondeterministic rounding as any cross-block reduction.

Libraries such as cuBLAS and CUTLASS combine all of these, with tile shapes chosen per GPU and per problem size, and newer layers such as ThunderKittens, CuTe, and Triton make the same structure easier to write. You will almost never need to write a production GEMM yourself. You will constantly need to reason about one, to predict whether a shape is efficient, to fuse an epilogue, or to write an attention kernel that has exactly this structure inside.

⚠️ Benchmarking a matmul honestly

Report FLOP/s as $2MNK/t$ for the shape you care about, after warm-up and with synchronization, and compare against the dense peak for the format, not the sparse headline number. Small or skinny shapes are often memory-bound or wave-quantized and will never approach peak, and that is a property of the shape, not of the kernel. Also check the result against a float64 reference. Each output is a sum of $K$ products accumulated phase after phase, so its rounding error grows with $K$ as in the sequential bound of §11.4, and an accumulator narrower than float32 loses accuracy silently. This is not hypothetical. DeepSeek-V3's FP8 matmuls copy the tensor cores' partial sums into float32 registers every 128 elements of $K$ , because Hopper's tensor cores accumulate FP8 products with only about 14 bits of precision (DeepSeek-AI, 2024).

In PyTorch

PyTorch's @ calls cuBLAS on a GPU and a tuned CPU library otherwise. The block checks that the blocked algorithm is the same computation as the library call, then times a tiled loop in PyTorch against one library call, a reminder that tiling only pays when the tile loop runs on the hardware, not in Python.

🔥 In PyTorchfor a machine with PyTorch

🎯 Interview questions on this section

How much global memory traffic does a naive matmul generate, and how much does tiling save? Derive it.

With one thread per output, each thread reads a row of $\mathbf{A}$ and a column of $\mathbf{B}$ , $2K$ elements, so the kernel reads $2MNK$ elements, against $MK + KN$ that exist. With $b\times b$ output tiles, each block loads $2b^2$ elements per phase for $K/b$ phases, and there are $MN/b^2$ blocks, so the total is $2MNK/b$ , a reduction by the tile size. For rectangular tiles it is $MNK(1/b_M + 1/b_N)$ , and the intensity per block is $2b_Mb_N/(s(b_M + b_N))$ , which is why output tiles are large and the inner tile $b_K$ only needs to be large enough to amortize the barriers.

Why do you need two __syncthreads() per phase in a shared-memory tiled matmul?

The first, after loading the tiles, ensures every element of both tiles is in shared memory before any thread reads them, since each thread reads elements loaded by others. The second, after the inner product, ensures every thread has finished reading the current tiles before any thread starts overwriting them with the next phase's data. Remove the first and threads read stale or missing data. Remove the second and fast warps corrupt tiles that slow warps are still using. Double-buffered kernels use two sets of buffers so loads of the next tile can overlap the current math.

Your shared-memory tiled kernel is still far from peak. What is the bottleneck, and what is the fix?

The inner loop reads two operands from shared memory per multiply-add. Shared memory delivers about 128 bytes per cycle per SM, 16 float32 multiply-adds' worth, while the SM can do 128 per cycle, so the kernel is bound by shared memory bandwidth at about an eighth of peak. Register tiling fixes it. Each thread computes a $t_M\times t_N$ block of outputs from $t_M + t_N$ loaded values per step, so shared reads per multiply-add fall to $(t_M+t_N)/(t_Mt_N)$ , a quarter for $8\times8$ . The cost is many registers per thread and therefore lower occupancy, which is fine because each thread has abundant independent work.

State the I/O lower bound for matrix multiplication and explain why tiling achieves it.

Hong and Kung proved that the classical algorithm on a machine with fast memory of $S$ words needs $\Omega(MNK/\sqrt{S})$ transfers. A clean proof splits execution into segments of $S$ transfers, notes each segment can touch at most $2S$ elements of each of $\mathbf{A}$ , $\mathbf{B}$ , and $\mathbf{C}$ , and uses Loomis-Whitney to bound the products per segment by $(2S)^{3/2}$ . Tiling with $b \approx \sqrt{S/3}$ moves $2MNK/b = O(MNK/\sqrt S)$ elements, matching the bound up to a constant.

What do tensor core instructions look like, and why do dimensions need to be multiples of 8 or 16?

A warp (Ampere) or a warpgroup of four warps (Hopper) issues a matrix multiply-accumulate on fixed tile shapes, such as mma.sync.m16n8k16 (a $16\times16$ by $16\times8$ product into a $16\times8$ float32 accumulator) or Hopper's asynchronous wgmma of shape $64\times N\times16$ . Every operand is consumed in these units, so dimensions are padded up to them, and a dimension of 17 costs as much as 32. Operands are also loaded in 16-byte vectors, which need row lengths that are multiples of 8 elements in 16-bit formats.

Name the main optimizations between a naive CUDA matmul and cuBLAS, and what each fixes.

Coalesced loads (wasted sectors), shared memory tiling (HBM traffic, divided by the tile size), register tiling (shared memory bandwidth), vectorized loads (load instruction count), padded or swizzled layouts (bank conflicts), multi-stage pipelining with cp.async or TMA (exposed load latency), warp tiling onto tensor core MMAs (the throughput of the FP32 lanes), split-K or stream-K (wave quantization for awkward shapes), grouped block ordering (L2 reuse), and epilogue fusion (an extra pass over $\mathbf{C}$ ). Boehm's worklog shows the first several taking a kernel from about 1 percent to over 90 percent of cuBLAS.

🧭 What to remember

The naive matmul loads $2MNK$ elements for $MK + KN$ distinct ones, an intensity of $1/s$ , and reaches about 1 percent of cuBLAS.
Shared memory tiling divides global traffic by the tile size, $2MNK/b$ , or $MNK(1/b_M + 1/b_N)$ for rectangular tiles, with two barriers per phase.
Register tiling applies the same reuse one level down, cutting shared reads per multiply-add from 2 to $(t_M + t_N)/(t_Mt_N)$ , at the price of many registers and low occupancy.
No schedule of the classical algorithm moves fewer than $\Omega(MNK/\sqrt S)$ elements with fast memory $S$ (Hong and Kung), and tiling achieves this bound.
Tensor cores execute fixed MMA shapes ( $16\times8\times16$ per warp on Ampere, $64\times N\times16$ per warpgroup on Hopper), so a fast GEMM is a hierarchy of block, warp, and instruction tiles, and shapes should be multiples of those tiles.
Next, Triton lets you write this whole tiled structure in a few lines of Python, with the compiler handling threads, shared memory, and tensor cores.

Boehm (2022). How to Optimize a CUDA Matmul Kernel for cuBLAS-like Performance: a Worklog. Ten kernels from naive to near cuBLAS, each measured and explained, the best single reading on this topic.
Kirk, Hwu & El Hajj (2022). Programming Massively Parallel Processors, 4th edition. Chapter 5 develops the shared memory tiled matmul and its traffic analysis.
Hong & Kung (1981). I/O Complexity: The Red-Blue Pebble Game. The original lower bound on data movement for matrix multiplication and other computations.
Irony, Toledo & Tiskin (2004). Communication lower bounds for distributed-memory matrix multiplication. The segment and Loomis-Whitney argument we used for the lower bound.
NVIDIA. Matrix Multiplication Background User's Guide. Tile and wave quantization, tensor core alignment requirements, and how shapes map onto the hardware.
NVIDIA. CUTLASS. The open-source template library that implements the full tiling hierarchy, with documentation of each level.
NVIDIA. Parallel Thread Execution ISA, warp-level matrix instructions. The exact MMA shapes, operand layouts, and the Hopper wgmma instructions.

Zero to Frontier is made for larger screens

Tiled matrix multiplication

The naive kernel and its traffic

Tiling with shared memory

How much traffic tiling saves

Global traffic of naive and tiled matmul

A tiled matmul, simulated

A tiled matmul that reports its traffic

Register tiling

A two-level tiled matmul

No schedule can do much better

Tensor cores and MMA shapes

Count tensor core instructions

The rest of the road to cuBLAS

In PyTorch

Read more

Watch next