Given an M×K matrix A and a K×N matrix B, the product C = AB is the M×N matrix whose entries are
\[C_{ij} = \sum_{k=0}^{K-1} A_{ik} B_{kj}.\]
The examples keep M, N, and K distinct – production GEMMs are rarely square, whatever the textbooks assume – and store A, B, and C in row-major order.
Computing C takes 2MNK floating point operations – MNK multiply-adds – while the three matrices together occupy only MK + KN + MN floats. The ratio of arithmetic to data therefore grows with the matrix dimensions. That is what makes GEMM special: unlike reduction or scan, whose one-pass-over-the-data nature limits their performance to memory bandwidth (Chapters 12 and 13), matrix multiplication has enough arithmetic per byte that a well-written kernel can be made compute-bound.
The catch is that the arithmetic intensity is only available. Realizing it means moving each input element from DRAM into registers once and reusing it as much as possible. With the notable exception of the introduction of Tensor Cores, the journey below amounts to a sequence of increasingly effective reuse of the input data as the output is computed.