Making a CPU Multiply Fast

The idea in one sentence

The same two billion multiply-adds take 2.6 seconds or 2.7 milliseconds on one chip, and every step between the two is a decision about the order memory is walked and how much of the chip is allowed to work.

One line, six rungs

Lecture 23 ended by promising this line:

with $A$, $B$ and $C$ all $1024 \times 1024$ and all single precision. Written as three nested loops it is four lines of C. Compiled the obvious way it takes two and a half seconds. Six changes later it takes under three milliseconds, and not one of the six changes which numbers get multiplied together or in what precision.

Every number in this chapter was produced on the machine the chapter was written on. That machine is an Apple M3 Max: twelve performance cores and four efficiency cores, 128 KiB of L1 data cache per performance core, 16 MiB of L2 shared by each group of six of them, 128-byte cache lines, 48 GiB of memory. The compiler is Apple clang 21.0.0. Nothing here is borrowed except one table near the end, which is labelled.

Count the work first

Before making anything faster it helps to know exactly how much work there is, because that number never changes and everything else is measured against it.

Every entry of $C$ is a dot product of a row of $A$ with a column of $B$. Both have 1024 numbers in them, so one entry costs 1024 multiplications and 1024 additions. There are $1024^2 = 1048576$ entries. Multiply:

Call it 2.15 billion operations, or 2.15 GFLOP. That is the arithmetic, all of it, for every rung below.

Now count the memory, at both extremes. Suppose the machine read each of $A$, $B$ and the old $C$ once and wrote $C$ once. That is $4 \times 1024^2 \times 4$ bytes, or 16 MiB, and the ratio of work to traffic would be

operations for every byte moved. A machine asked to do 128 things with each byte it fetches is being asked a very easy question.

The obvious loop asks a much harder one. Each of the 1048576 entries reads a whole row of $A$ and a whole column of $B$, which is 8192 bytes, and does 2048 operations on them. That is a quarter of an operation per byte, and 8.59 GB of requests for a problem that holds 12 MiB of data. The distance between 0.25 and 128 is what the rest of this chapter closes.

What the chip can do

The other end of the measurement is the ceiling. Two things in the hardware set it, and both are worth one paragraph.

The first is the fused multiply-add. The inner statement of a matrix multiply is always the same shape, a product added to a running sum, so the hardware provides it as a single instruction that does both and rounds once. The second is the vector unit. One register holds four single-precision numbers side by side, and one instruction applies the same operation to all four.

The measured ceiling is better than a datasheet, so here it is measured. Sixteen independent accumulators, four floats wide, nothing but arithmetic, no memory in the loop at all:

/* 16 independent accumulators, four floats each, nothing but arithmetic */
float32x4_t a[16], x = vdupq_n_f32(1.000001f), y = vdupq_n_f32(0.999999f);
for (int i = 0; i < 16; i++) a[i] = vdupq_n_f32((float)i);

for (long t = 0; t < 400000000L; t++)
    for (int i = 0; i < 16; i++)
        a[i] = vfmaq_f32(a[i], x, y);         /* one fmla.4s each */

threads=1     0.453 s     113.1 GFLOP/s  (14.1 G vector-FMA/s per thread)
threads=4     0.453 s     452.2 GFLOP/s  (14.1 G vector-FMA/s per thread)
threads=12    0.477 s    1288.3 GFLOP/s  (13.4 G vector-FMA/s per thread)
threads=16    0.650 s    1260.7 GFLOP/s  (9.8 G vector-FMA/s per thread)

One performance core issues 14.1 billion vector multiply-adds per second, which is 113 GFLOP/s. Twelve of them reach 1288 GFLOP/s, and the scaling out to twelve is nearly perfect. Sixteen threads is slower than twelve, which is the first appearance of a fact that returns at the last rung: four of this chip’s sixteen cores are small ones, and asking them to keep up with the big ones slows everyone down.

So 2.15 GFLOP at 1288 GFLOP/s would take 1.7 milliseconds. That is the target the loops below are walking towards, and none of them will reach it.

The loops as anybody writes them

Here is the multiply, written the way the definition reads. Rows of $C$ outside, columns next, the dot product inside. Simon Boehm calls this order RCI, for row, column, inner, and the name is useful enough to borrow.

void mm_naive(const float *A, const float *B, float *C, int n) {
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            for (int k = 0; k < n; k++)
                C[i * n + j] += A[i * n + k] * B[k * n + j];
}

Compiled with clang -O0, which is what you get when nobody remembers to turn the optimiser on, it runs in 2.575 seconds. That is 0.83 GFLOP/s, seven tenths of one percent of what a single core was doing in the listing above.

Compiled with clang -O3 the same source runs in 1.494 seconds. The optimiser is worth 1.7 times and then it stops, which should be surprising. It had the whole function in front of it and it did not find the factor of a thousand sitting there. Optimisers rewrite instructions. The problem here is not the instructions.

The rung that does almost nothing

The obvious next move is to stop touching memory in the inner loop. The statement C[i * n + j] += ... reads and writes the same address 1024 times in a row, and the compiler is not allowed to keep that address in a register, because for all it knows C overlaps A or B. So say it explicitly:

void mm_acc(const float *A, const float *B, float *C, int n) {
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++) {
            float acc = 0.0f;
            for (int k = 0; k < n; k++)
                acc += A[i * n + k] * B[k * n + j];
            C[i * n + j] += acc;
        }
}

This is a real improvement to the code and it buys 1.02 times. Measured: 1.464 seconds against 1.494. Within the run-to-run spread it is nothing at all.

The rung is here because failing is informative. Removing 1024 loads and 1024 stores per entry changed nothing, so those loads and stores were not what the loop was waiting for. Something else is, and it has to be $B$.

The swap

In the loop above, k moves fastest. Watch what each array does as k advances by one. A[i * n + k] moves forward four bytes, so it walks along a row. C[i * n + j] does not move. B[k * n + j] moves forward $n$ floats, 4096 bytes, so it walks down a column.

That is Figure 5 of Lecture 23, arriving in real code. A cache line here is 128 bytes, which holds 32 floats. Walking a row, 32 consecutive values ride in on one line and 31 of them were free. Walking a column with a stride of 4096 bytes, every value needs a line of its own, and the other 31 numbers on each line are fetched and thrown away.

The fix is to make j the fast index instead of k. Pull the one value of $A$ out into a scalar, because it no longer changes in the inner loop, and let $B$ and $C$ both be walked along their rows. The loops now read row, inner, column, which is RIC.

void mm_ric(const float *A, const float *B, float *C, int n) {
    for (int i = 0; i < n; i++)
        for (int k = 0; k < n; k++) {
            float a = A[i * n + k];
            for (int j = 0; j < n; j++)
                C[i * n + j] += a * B[k * n + j];
        }
}

Every multiplication that happened before still happens, and every sum still gets the same 1024 terms in the same order. Only the order of the visits changed. The time goes from 1.464 seconds to 75.9 milliseconds, a factor of 19.3, and it is by a long way the largest jump on this ladder.

Two things happened at once there, and separating them is worth doing. Compile the same two functions again with the vectoriser switched off, -fno-vectorize -fno-slp-vectorize, and time them:

The RCI loop does not care, because it never vectorised. The RIC loop cares a great deal. So the reordering on its own is worth 4.7 times, from 1.464 down to 0.3119, and the vector unit is worth another 4.1 on top of that, down to 0.0759. Multiply the two and you get the 19.3.

The vector unit did not appear because a flag was added. It appeared because the reordering made the inner loop vectorisable at all: 32 consecutive values of $B$, 32 consecutive values of $C$, and one value of $A$ that stays put. The generated assembly says so plainly. The RIC inner loop is built out of fmla.4s v5, v1, v0[0], a four-lane multiply-add whose second operand is one lane broadcast across the vector, which is the bottom half of Figure 2 exactly. The RCI inner loop contains no vector multiply-add anywhere, only the scalar fmadd s0, s1, s2, s0.

That is the deeper reason loop order matters. A bad order costs you memory traffic, and then it quietly costs you the arithmetic units too, because the instructions that would use them need neighbouring data to work on.

A tile of the answer

At 75.9 milliseconds the kernel runs at 28.3 GFLOP/s, a quarter of what one core proved it could do. The inner loop now does one multiply-add per element of $B$ and pays a load and a store on $C$ for each one, so most of its instructions still move numbers rather than multiply them.

The cure is to make one trip through memory do more arithmetic. Give the innermost work a rectangle of $C$ instead of a single row of it.

Here is the tile. It holds a 4 by 16 patch of $C$ in local variables for the entire sweep of $k$, so those 64 running sums live in registers and touch memory exactly twice, once at the start and once at the end. Sixteen floats is four vector registers wide, and four rows means four accumulator vectors per column group, so 16 vector registers in total out of the 32 the instruction set provides.

static void tile(const float *A, const float *B, float *C, int n,
                 int i0, int j0, int k0, int k1) {
    float acc[MR][NR] = {{0.0f}};
    for (int k = k0; k < k1; k++) {
        float a[MR], b[NR];
        for (int ii = 0; ii < MR; ii++) a[ii] = A[(i0 + ii) * n + k];
        for (int jj = 0; jj < NR; jj++) b[jj] = B[k * n + j0 + jj];
        for (int ii = 0; ii < MR; ii++)
            for (int jj = 0; jj < NR; jj++)
                acc[ii][jj] += a[ii] * b[jj];
    }
    for (int ii = 0; ii < MR; ii++)
        for (int jj = 0; jj < NR; jj++)
            C[(i0 + ii) * n + j0 + jj] += acc[ii][jj];
}

static void rows(const float *A, const float *B, float *C, int n, int r0, int r1) {
    for (int k0 = 0; k0 < n; k0 += KC)
        for (int i = r0; i < r1; i += MR)
            for (int j = 0; j < n; j += NR)
                tile(A, B, C, n, i, j, k0, k0 + KC);
}

With MR 4, NR 16 and KC 256 that runs in 24.5 milliseconds, 87.7 GFLOP/s, a factor of 3.1 over the reordered loop. It is 78 percent of the 113 GFLOP/s one core managed with no memory in the loop at all, which is close enough to stop.

Where the working set has to fit

The outer loop above cuts $k$ into panels of KC, which is the classic cache blocking move, and it deserves an honest accounting, because at this size it does almost nothing.

Read the figure as a list of things that must fit. The 4 by 16 tile is 256 bytes and lives in registers. The slices the tile reads at each step, 4 by 256 of $A$ and 256 by 16 of $B$, come to 20 KiB against 128 KiB of L1. The panel of $B$ that the whole i loop sweeps through is $256 \times 1024$ floats, 1 MiB, against 16 MiB of L2.

And there is the problem. All three matrices together are 12 MiB, and L2 is 16 MiB. At $N = 1024$ the entire problem already lives in cache, so cutting $k$ into panels has nothing to save. Sweeping KC over a sixteen-fold range confirms it:

A KC equal to $N$ means no blocking at all, so the last filled cell in each row is the unblocked run. Along the top row nothing happens: sixteen-fold changes in the panel depth move the time by less than the spread between repeats of a single setting.

The bottom row is where blocking earns its reputation. At $N = 4096$ the three matrices are 192 MiB, which fits in nothing, and the unblocked run takes 10.74 seconds against 2.15 for KC = 512. Five times, from the same instructions, for choosing how much of $B$ to keep hot.

So the rung is real and this size is the wrong place to see it. That is a common shape in performance work: an optimisation is a bet about which store the data is spilling out of, and if it is not spilling yet the bet pays nothing. The code above keeps the panel loop because it costs nothing here and it is what stops the same kernel falling over at 4096.

Threads

One core is now doing 78 percent of what one core can do. The remaining factor has to come from the other eleven.

Matrix multiply is the friendliest parallel problem there is. Split $C$ into horizontal bands. Each band needs the matching rows of $A$ and all of $B$, both read-only, and it writes only its own rows of $C$. No two threads ever touch the same output byte, so there is nothing to lock and nothing to synchronise except the join at the end.

static void *worker(void *p) {
    job_t *j = p;
    int nb = j->n / BAND, t;
    while ((t = atomic_fetch_add(&next_band, 1)) < nb)
        rows(j->A, j->B, j->C, j->n, t * BAND, (t + 1) * BAND);
    return NULL;
}

The counter is the interesting part. The obvious thing is to hand each thread an equal slice up front, and on a chip whose cores are all the same that would be right. This chip has twelve fast cores and four slow ones, so an equal split finishes when the slowest slice finishes. Cutting $C$ into 32 bands of 32 rows and letting each thread take the next free one whenever it is idle means the small cores simply take fewer bands.

Both were measured in the same program, three passes each:

Twelve threads through the counter give 2.71 milliseconds. Against the single-threaded rung of the previous section that is 9.0 times, and against the one-thread entry in the table above, which pays for the atomic counter it does not need, it is 9.4. The same twelve threads on fixed slices give 3.55, so the counter is worth 1.3 times for three lines of code. Adding the four efficiency cores moves 2.71 to 2.69, which is nothing, and moves the fixed-slice version from 3.55 to 3.43, which is also nothing. Those four cores exist for background work at low power, and this is not that.

Nine times on twelve cores rather than twelve is the usual story. The bands are not all equally fast, the last threads finish alone, and the six cores in a cluster share one L2 while every one of them streams the same $B$ through it.

The whole staircase

Here is the ladder in one place. Rung 1 and rung 2 are the same source file compiled twice; the rest differ only in which function is called.

clang -O0 -o ladder0 ladder.c -framework Accelerate
clang -O3 -o ladder3 ladder.c -framework Accelerate

./ladder0 naive  7        # 1. naive, row-column-inner, no optimiser
./ladder3 naive 15        # 2. the same source, -O3
./ladder3 acc   15        # 3. accumulate in a register
./ladder3 ric   25        # 4. swap the two inner loops
./ladder3 block 25        # 5. hold a 4x16 tile of C in registers
./ladder3 par   41 12     # 6. hand row bands to 12 threads
./ladder3 blas  41        # -- Apple's Accelerate, for scale

naive  reps=7    median   2.575446 s       0.83 GFLOP/s   spread 1.7%   max|err| 0.00e+00 (|C|max 13.47)
naive  reps=15   median   1.494318 s       1.44 GFLOP/s   spread 20.0%   max|err| 0.00e+00 (|C|max 13.47)
acc    reps=15   median   1.499144 s       1.43 GFLOP/s   spread 31.1%   max|err| 1.34e-05 (|C|max 13.47)
ric    reps=25   median   0.075923 s      28.29 GFLOP/s   spread 2.3%   max|err| 0.00e+00 (|C|max 13.47)
block  reps=25   median   0.023807 s      90.20 GFLOP/s   spread 22.7%   max|err| 1.72e-05 (|C|max 13.47)
par    reps=41   median   0.002711 s     792.14 GFLOP/s   spread 6.6%   max|err| 1.72e-05 (|C|max 13.47)
blas   reps=41   median   0.000948 s    2265.28 GFLOP/s   spread 23.2%   max|err| 0.00e+00 (|C|max 13.47)

That is one pass. Five of them, with the median taken across passes, give the table the figure below is drawn from.

Four things in that table are worth saying out loud.

The max|err| column says the fast versions are not sloppier than the slow ones. Every rung lands within $2 \times 10^{-5}$ of Apple’s sgemm on entries of size up to 13.5, which is float32 noise. Two of them, the naive loop and the reordered loop, match sgemm bit for bit, and all three differ from a float64 computation of the same product by the same $1.8 \times 10^{-5}$. Reordering a floating point sum does change the last bits. It changed nothing here that anybody would care about.

The largest jump is the loop swap, which changes no arithmetic at all. The second largest is threads, which changes no arithmetic either. The rung that looks most like optimising, hoisting a value into a register, bought two percent.

Twelve threads reach 791 GFLOP/s against the 1288 the cores proved they could do, so the finished kernel runs at 61 percent of its own arithmetic ceiling. Almost 40 percent of the machine is still spent waiting for memory, after all six rungs.

And the same ladder, on very different hardware, has the same shape. Simon Boehm built these six rungs on an Intel i7-6700, a four-core Haswell at 3.4 GHz, at the same size of 1024, and published the times in his CPU matrix multiplication worklog: 4481 ms naive, 1621 with compiler flags, 1512 with the register accumulate, 89 after the loop reorder to RIC, 70 after L1 tiling, 16 with multithreading, against 8 ms for NumPy on Intel’s MKL. His flags were -O3 -march=native -ffast-math where the flags rung here is a plain -O3, and his sizes were compile-time constants where these are not. Line the two ladders up anyway: the register accumulate is worth 1.07 times on his machine and 1.02 on this one, and the loop reorder is worth 17 times on his and 19.3 on this one. Different decade, different vendor, different instruction set, and the same two rungs carry the chapter.

What a library does instead

The last row of the table is the one nobody wrote here. Apple’s Accelerate does the same multiply in 0.947 milliseconds, 2268 GFLOP/s, another 2.9 times beyond six rungs of work.

That number is larger than the whole chip’s vector units can produce. The measured NEON ceiling across all sixteen cores was 1288 GFLOP/s, and Accelerate is running at 1.8 times it, so it is not doing this arithmetic on the vector units at all. Apple silicon carries matrix hardware that C source cannot reach, and a call into sgemm is how you reach it.

NumPy sits on that same library, so the same thing is visible from Python:

import numpy as np
from timeit import repeat

n = 1024
rng = np.random.default_rng(0)
A = rng.standard_normal((n, n), dtype=np.float32)
B = rng.standard_normal((n, n), dtype=np.float32)

t = min(repeat(lambda: A @ B, number=50, repeat=11)) / 50
print("%.6f s   %.0f GFLOP/s" % (t, 2 * n**3 / t / 1e9))

# 0.000741 s   2899 GFLOP/s

Which makes the other Python experiment worth running. Transliterate the same three loops into Python, shrink the problem to $256 \times 256$ so it finishes, and compare against @:

import numpy as np, time
from timeit import repeat

n = 256
rng = np.random.default_rng(0)
A = rng.standard_normal((n, n))
B = rng.standard_normal((n, n))
Al, Bl = A.tolist(), B.tolist()

def triple(A, B, n):                     # the C loops, transliterated
    C = [[0.0] * n for _ in range(n)]
    for i in range(n):
        for j in range(n):
            s = 0.0
            for k in range(n):
                s += A[i][k] * B[k][j]
            C[i][j] = s
    return C

t0 = time.perf_counter()
C = triple(Al, Bl, n)
t_py = time.perf_counter() - t0
t_np = min(repeat(lambda: A @ B, number=20, repeat=9)) / 20

print("triple loop  %8.4f s  %8.4f GFLOP/s" % (t_py, 2 * n**3 / t_py / 1e9))
print("A @ B        %8.6f s  %8.1f GFLOP/s" % (t_np, 2 * n**3 / t_np / 1e9))
print("factor       %8.0f" % (t_py / t_np))
print("same answer  %s" % np.allclose(np.array(C), A @ B))

# triple loop    0.4492 s    0.0747 GFLOP/s
# A @ B        0.000089 s     375.9 GFLOP/s
# factor           5032
# same answer  True

Five thousand times, for the identical answer. Read that carefully, because it is not five thousand times of interpreter overhead. The same three loops in C at -O3, timed at this same size, take 14.2 milliseconds, which is only 32 times faster than the Python. The remaining factor of about 160 is the six rungs of this chapter plus the matrix hardware, all of it hiding inside one @.

That is the honest reason every machine learning framework hands its matrix multiplies to a library. A forward pass through a linear layer is this call and nothing else, and the distance between writing the loop and calling the library is not a matter of taste. It is the factor of a thousand you have just watched being assembled, one decision at a time.

Where this is going

Look back at what the six rungs actually were. One was a compiler flag. One did nothing. The other four are the same instruction: move the data less, and reuse it more once it has arrived. Reordering the loops reused a cache line 32 times instead of once. The tile reused a loaded value across 64 multiplications instead of one. Threads reused a chip that was sitting idle. Cache blocking, when the problem is big enough to need it, reuses a panel of $B$ across a whole band of $C$.

That instruction does not depend on the hardware, which is why the next three chapters can repeat it on a machine that looks nothing like this one. A GPU trades twelve clever cores for thousands of simple ones, the arithmetic ceiling goes up by more than an order of magnitude, and every rung of the ladder has to be climbed again from the bottom. Lecture 25 draws the model that organises those threads. Lecture 26 rebuilds this exact ladder inside it, ten rungs deep, and the largest jump there is the same one it was here: a swap that changes which index runs fastest.

If you want to write these rather than read about them, the Library’s practice grounds are all on the GPU side, which is where the next two chapters go: GPU Puzzles and Triton Puzzles for the language, LeetGPU for a browser with real hardware behind it, and the GPU MODE reference kernels for problems people compete on. For the CPU, the honest exercise is the one this chapter is. Take the four lines, and go and get your own factor of a thousand.