Notes From Simon Boehm’s CUDA Matmul Blog

These are my notes from reading Simon Boehm’s How to Optimize a CUDA Matmul Kernel for cuBLAS-like Performance. I am going through it one kernel at a time and trying to understand not only what changes, but why each change makes the GPU faster.

We are going to optimize SGEMM: single-precision general matrix multiplication.

$$C=\alpha AB+\beta C$$

For $A\in {\mathbb{R}}^{M\times K}$, $B\in {\mathbb{R}}^{K\times N}$, and $C\in {\mathbb{R}}^{M\times N}$, the literal arithmetic in the displayed kernel is approximately:

$$2MNK+3MN\text{ FLOPs}$$

The $3MN$ comes from the final $\alpha \cdot \text{sum}+\beta \cdot C$ operation. In normal GEMM benchmarks we ignore this comparatively small part and count the workload as $2MNK$ FLOPs.

Now let us calculate the minimum amount of data that must move. Assume square matrices with width $W$. We must read $W^2$ FP32 values from A, read $W^2$ from B, read the old $W^2$ values of C, and write the new $W^2$ values of C. That is four matrix-sized transfers, and every FP32 value occupies 4 bytes:

$$4W^2\text{ values}\times 4\text{ bytes}=16W^2\text{ bytes}$$

For $W=4092$:

$$16\times 4092^2=267,911,424\text{ bytes}\approx 267.91\text{ MB}$$

This is the minimum data volume, not memory bandwidth. Data volume tells us how many bytes must move. Bandwidth tells us how many bytes the GPU can move per second.

That amount of mathematical work is fixed. If two kernels compute the same SGEMM, we cannot magically remove the required FLOPs. What we can optimize is how quickly the GPU performs them.

$$\text{achieved FLOP/s}=\frac{2MNK}{\text{execution time}}$$

So the goal of every optimization in these notes is simple: perform the same $2MNK$ useful FLOPs in less time.

Is The Kernel Memory-Bound Or Compute-Bound?

For the $4092\times 4092$ square SGEMM, the ideal workload is:

• approximately $137.04$ GFLOPs of useful arithmetic;
• at least $267.91$ MB of global-memory traffic, assuming perfect reuse.

Assume the GPU provides:

• FP32 compute throughput: $30$ TFLOP/s;
• global-memory bandwidth: $768$ GB/s.

If arithmetic were the only limitation, the fastest possible compute time would be:

$$t_{\text{compute}}\ge \frac{137.04\text{ GFLOP}}{30,000\text{ GFLOP/s}}\approx 4.57\text{ ms}$$

If memory movement were the only limitation, the fastest possible transfer time would be:

$$t_{\text{memory}}\ge \frac{0.26791\text{ GB}}{768\text{ GB/s}}\approx 0.349\text{ ms}$$

The kernel cannot finish before both the arithmetic and the required memory movement finish. In the ideal case they can overlap, so the lower bound is controlled by the slower one:

$$t_{\text{kernel}}\ge \max (4.57,0.349)\text{ ms}=4.57\text{ ms}$$

The compute work takes about $13.1\times$ as long as the minimum memory movement. Therefore an ideally reused implementation of this SGEMM is compute-bound.

The reuse assumption is everything. The $267.91$ MB calculation assumes A and B are fetched from global memory once and then reused on-chip. A naive kernel repeatedly requests the same values for different outputs, so its memory traffic can be vastly larger. The mathematical operation may be ideally compute-bound while a bad implementation of it is still memory-bound.

What About Tensor Cores?

The $30$ TFLOP/s roof above refers to ordinary FP32 arithmetic. Tensor cores are separate matrix-multiply hardware and can provide much higher throughput when we use their supported matrix instructions, data types, and layouts.

The kernels in this article perform scalar FP32 fused multiply-add instructions on the regular FP32 pipelines. They do not issue tensor-core MMA instructions. Therefore using a tensor-core peak in this calculation would give us the wrong roof for the code we are actually optimizing. Warp tiling later in the article organizes work in a tensor-core-like hierarchy, but it still does not turn these kernels into tensor-core kernels.

Naive Implementation

The naive implementation assigns one thread to each output element of C. A thread that owns $C[row,col]$ must calculate a length-$K$ dot product, so it reads the complete required row of A and the complete required column of B. It also reads the previous value of C for the $\beta C$ term.

Therefore each output thread logically reads:

$$K+K+1=2K+1\text{ FP32 values}$$

There are $M\times N$ output elements, and therefore $M\times N$ output threads. The total logical read demand is:

$$MN(2K+1)\times 4\text{ bytes}$$

For the article’s square case, $M=N=K=W=4092$:

$$4092^2(2\times 4092+1)\times 4\approx 548.21\text{ GB of logical reads}$$

Including the final write of every C element gives approximately $548.28$ GB of logical memory operations. The ideal perfectly reused implementation needed only $267.91$ MB, so the naive demand is roughly $2,046\times$ the ideal minimum.

This is a logical request count, not measured DRAM traffic. Caches, warp-wide reuse, broadcasts, coalescing, and transaction overfetch change how many bytes physically cross global memory.

The FLOP count has not changed. The naive kernel and every optimized kernel still perform the same $2MNK$ useful FLOPs. We are not removing matrix-multiplication work. We are removing redundant data movement so the GPU can finish those same FLOPs in less time and therefore achieve a higher FLOP/s rate.