from a C loop to a GPU kernel
The hook I wanted was clean:
512x512 matmul
one C loop
460x slower than BLASThe machine refused to give me that exact number.
(The 460x came from cross-build arithmetic: the unoptimized naive time divided by an optimized build’s Accelerate sample, 385.4 ms over ~0.84 ms. Two numbers that never lived in the same binary.)
If I compile like a beginner, without optimization, it is worse:
Accelerate SGEMM 0.152 ms 1767.00 GFLOP/s checksum -88.858457
naive ijk 385.396 ms 0.70 GFLOP/s checksum -88.858457
max abs diff vs BLAS 0, slowdown 2536.9xIf I compile like someone who knows enough to type -O3 -ffast-math, it is less dramatic:
Accelerate SGEMM 0.787 ms 341.09 GFLOP/s checksum -88.858404
naive ijk 140.085 ms 1.92 GFLOP/s checksum -88.858162
max abs diff vs BLAS 2.19345e-05, slowdown 178.0xAnd if I use a separate optimized build without -ffast-math, I get another honest number:
Accelerate SGEMM 0.808 ms 332.14 GFLOP/s checksum -88.858457
naive ijk 107.875 ms 2.49 GFLOP/s checksum -88.858457
max abs diff vs BLAS 0, slowdown 133.5xSo this is not the 460x post. It is the better version: the one where the number moves when the toolchain moves, and the lesson survives.
A tiny C loop can be hundreds or thousands of times slower than a tuned matrix multiply. A first Metal GPU kernel can still be much slower than the library. Tiling can make code worse. The compiler can emit NEON and FMA (arm64’s vector instructions and its fused multiply-add, both dissected below) for a loop you wrote in ordinary C, but still fail to make the whole algorithm good. The GPU can have enormous arithmetic throughput and still spend its time waiting for bytes.
Assumed background: beginner C or C++, beginner Python-level array intuition, and enough comfort with nested loops to read this:
for (int i = 0; i < N; ++i) {
for (int j = 0; j < N; ++j) {
float sum = 0.0f;
for (int k = 0; k < N; ++k) {
sum += A[i * N + k] * B[k * N + j];
}
C[i * N + j] = sum;
}
}The triple loop is matrix multiplication in row-major memory, and row-major means rows are contiguous: A[row][col] lives at A[row * N + col].
For N = 512, each matrix has:
512 * 512 floats = 262144 floats
262144 * 4 bytes = 1048576 bytes = 1 MiBThe computation does:
512 * 512 * 512 multiplies
512 * 512 * 512 adds
= 268435456 floating-point operationsThat number is large enough to be interesting and small enough that we can run bad versions without making tea.
Local machine:
Apple M4
Apple clang version 15.0.0 (clang-1500.3.9.4)
Target: arm64-apple-darwin24.6.0Accelerate SGEMM is the local BLAS baseline. BLAS means Basic Linear Algebra Subprograms, and SGEMM means single-precision general matrix multiply. BLAS is an interface and a tradition: you call cblas_sgemm, and the vendor implementation tries very hard to make matrix multiply fast on that machine. I’m not claiming Accelerate is “just C with good loops.” On Apple Silicon it’s a vendor library, and the exact path can involve implementation details Apple doesn’t expose. That is the point of using it as the adult in the room.
The benchmark is a trap too
Before touching the loops, look at the measurement itself. The probe prints a time, a GFLOP/s figure, a checksum, a max-abs-diff against the BLAS result, and a slowdown ratio. The checksum is a cheap smoke alarm, not a proof of correctness, it samples the output so a totally broken kernel (the hot inner function doing the actual math; later in this post, literally the GPU function) can’t look fast by writing zeros or skipping work.
The real correctness check is:
static double max_abs_diff(const std::vector<float>& a,
const std::vector<float>& b) {
double m = 0.0;
for (size_t i = 0; i < a.size(); ++i) {
m = std::max(m, double(std::fabs(a[i] - b[i])));
}
return m;
}Every candidate gets compared with the BLAS result. Even this isn’t perfect. Floating-point matrix multiply can produce small differences when the accumulation order changes, and -ffast-math gives the compiler permission to reassociate some floating-point operations, which can change low bits. A maximum absolute difference around 1e-5 for this fp32 probe (fp32, ordinary 32-bit floats) isn’t surprising. A difference of 10 would be.
The diff that actually made me stop isn’t the 2.19345e-05. It’s the 0. The ikj and tiled loops below match Accelerate bit-for-bit, max abs diff exactly zero across 262,144 floats, and so does the naive loop in both builds without -ffast-math. A hand-written triple loop agreeing with a vendor BLAS to the last bit is stranger than it looks, because any difference in summation order or rounding anywhere would surface. My unverified hypothesis: for this size, Accelerate accumulates each output element in plain increasing-k order with fused multiply-adds, and that happens to match what clang emits for my loops under its default fp-contract. The candidates that match are exactly the ones where each output still receives its contributions in strict k order, vectorizing across j doesn’t disturb that; vectorizing the k accumulation does, which fits the guess. I haven’t chased this to the bottom.
The benchmark takes the best time out of repeated runs, best of 8 for the Accelerate baseline, best of 3 for the loop variants, and a single run for the naive loop, which is too slow to repeat politely:
double best = 1e100;
for (int r = 0; r < reps; ++r) {
auto t0 = Clock::now();
f();
auto t1 = Clock::now();
best = std::min(best, elapsed);
}Best-of timing is common for small kernels because it tries to capture what the code can do when the OS mostly leaves it alone. Median timing answers a different question, and worst timing mostly measures interruptions.
The BLAS number itself refuses to sit still, and this matters more than any other caveat here. Across the builds above, Accelerate’s best-of-8 came out at 0.152 ms, 0.787 ms, 0.808 ms, and, in the Metal probe further down, 0.184 ms. I re-ran the same fast-math binary ten times in a row while editing this post: nine runs landed at 0.148 ms and one at 0.677 ms, and an eleventh full run caught 0.176 ms. That isn’t noise around a mean. It’s two modes about 4-5x apart, and a whole process lands in one mode or the other, even best-of-8 inside the slow run never saw a fast sample.
So every “Nx slower than BLAS” ratio in this post depends on which Accelerate mode that run happened to catch. The main fast-math table below caught a slow sample, 0.787 ms; against the fast mode the same naive loop reads about 950x, not 178x, and today’s 0.176 ms re-run printed it as 687.7x. My best guess for the bimodality, and it is a guess, is warm-up or dispatch state around the machine’s matrix hardware: some processes get the fast path from the first call, and some never do. I haven’t traced it.
I did try to pin it down. I wrote a probe that does best-of-20 SGEMM inside one process and prints every sample, then launched it as 30 fresh processes back to back:
best-of-20 across 30 processes, sorted and counted:
0.151(4) 0.153(2) 0.155 0.156(3) 0.157(2) 0.159(2)
0.161(4) 0.163 0.166 0.167 0.168 0.169
0.178 0.179 0.188(2) 0.189 0.191 0.206Every single one landed in the fast mode, between 0.151 and 0.206 ms, and inside one process the twenty samples were tight too, 0.151 to 0.198, no split. So today I could not reproduce the slow mode at all. That doesn’t unsay the earlier 0.787 and 0.808 samples, they were real, printed by the same code on the same machine. It just means the slow mode isn’t a stable property of the binary or the input; it’s a property of some machine state, boot, thermal, scheduler, matrix-unit residency, that was present when I first drafted this and absent when I went back to chase it. The mystery got deeper, not shallower: whatever causes it, I can’t summon it on demand, which is the worst kind of bug to own.
That’s why I don’t build the post around one sacred slowdown ratio. The fact that survives the noise isn’t 460x, or 178x. It’s that simple source-level changes move performance by orders of magnitude, and the tuned library stays far away, whichever mode it’s in.
The first loop
Here is the full naive function:
static void matmul_ijk(const float* A, const float* B, float* C) {
for (int i = 0; i < N; ++i) {
for (int j = 0; j < N; ++j) {
float sum = 0.0f;
for (int k = 0; k < N; ++k) {
sum += A[i * N + k] * B[k * N + j];
}
C[i * N + j] = sum;
}
}
}The math is correct. For every output element C[i][j], we walk along row i of A and column j of B. For a 3x3 example:
C[1][2] = A[1][0] * B[0][2]
+ A[1][1] * B[1][2]
+ A[1][2] * B[2][2]The CPU doesn’t see “row” and “column.” It sees addresses. For A[i * N + k], when k increases by 1, the address moves by one float:
A[i*N + 0]
A[i*N + 1]
A[i*N + 2]Good. Adjacent memory.
For B[k * N + j], when k increases by 1, the address jumps by an entire row:
B[0*N + j]
B[1*N + j]
B[2*N + j]For N = 512, that jump is 512 floats * 4 bytes = 2048 bytes. The loop is doing perfect textbook matrix multiplication and terrible memory traversal at the same time; the math expression and the execution are different objects.
The 2048-byte problem
Main memory is large. Registers are tiny. Between them lives a cache hierarchy. The simplified picture:
registers
current values inside the core
L1 cache
very small, very fast, private to a core
L2 cache
bigger, slower, often private or shared by a small cluster
last-level / system cache
bigger again, shared more widely
DRAM / unified memory
the main-memory chips: huge compared with caches, slow compared with coresA cache is a bet: if you asked for this address, you’ll probably ask for nearby addresses soon, or the same address again. So the machine doesn’t usually fetch one float. It fetches a cache line, a small fixed-size chunk such as 64 or 128 bytes; exact line sizes vary by machine. The important part is that one cache miss buys you a neighborhood, not one value.
I stopped hand-waving and asked this machine what its numbers actually are:
hw.cachelinesize 128 bytes
hw.pagesize 16384 bytes
perflevel0 (P-core) l1dcachesize 131072 = 128 KiB
perflevel0 l2cachesize 16777216 = 16 MiB, shared by 4 P-cores
perflevel1 (E-core) l1dcachesize 65536 = 64 KiBSo on the M4, a cache line is 128 bytes, which is 32 floats. One B[k*N + j] load pulls in B[k*N + j .. j+31], 32 contiguous columns of one row of B, and the naive loop’s next k iteration jumps to a different row and never touches 31 of them before eviction pressure builds. A whole 512-float row of B is 2048 bytes, exactly 16 cache lines. And the number that quietly reorganizes the rest of this post: all three matrices together are 3 MiB, and the P-cores’ shared L2 is 16 MiB. The entire problem fits in L2 with room to spare. That single fact turns out to explain why my tiling regressed and why, later, threading the loop helped when I bet it wouldn’t.
Our A access uses the neighborhood:
A[i*N + 0] A[i*N + 1] A[i*N + 2] A[i*N + 3] ...Our B access mostly throws the neighborhood away:
B[0*N + j]
jump 2048 bytes
B[1*N + j]
jump 2048 bytes
B[2*N + j]The hardware prefetcher (logic that watches your access stream and fetches the next lines before you ask) may notice some regular strides, and the caches may rescue some data because the whole 512x512 matrix isn’t huge. But the access pattern is still awkward for the memory system. When people say “cache-friendly,” they mean spatial locality, after you read one address, you soon read its neighbors, and temporal locality, reading the same address again soon. The naive loop has decent locality for A, bad spatial locality for B, and mediocre reuse overall: it computes just one output value at a time, so it throws away opportunities to reuse the same loaded A and B values across many C entries.
The CPU is fast. This loop makes it beg for bytes.
Two different words hide in that begging, and they get mixed together constantly. Latency is how long one request takes to come back; bandwidth is how much data per second can be delivered once the system is flowing. A water pipe analogy is overused because it works: latency is how long until the first water reaches you, bandwidth is how much water per second the pipe can carry. The naive ijk loop has trouble with both. The stride through B creates loads that are far apart, and far-apart loads are harder for caches and prefetchers to turn into a smooth stream. The core can have multiple cache misses in flight, but it can’t have infinity of them. If the arithmetic depends on the loaded value, the instruction stream eventually waits.
The ikj loop turns the inner work into streams: for each fixed i and k, it reads B[k*N + j] and reads/writes C[i*N + j] for j = 0..511. Streaming is no guarantee of speed, but hardware is good at it. Adjacent addresses make cache lines useful, predictable access helps prefetch, vector loads become easier, and the store buffers (the queue of pending writes between the core and the cache) behave better.
There’s another detail hiding here: stores matter. The naive ijk loop accumulates into a local sum and stores one C value at the end. The ikj loop updates C[i*N + j] many times:
C[i * N + j] += a * B[k * N + j];so it reads and writes C repeatedly. It’s still faster because those reads and writes are contiguous and repeatedly touch the same row of C while it’s hot in cache; the bad B traversal in ijk costs more than the extra C traffic in this particular probe. That’s what this matrix size, this compiler, and this machine did, not a universal law.
Two lines swapped, 17x
Matrix multiplication has three loops, and you can legally reorder them if you preserve the arithmetic result. The naive order is i, j, k: for each row, for each column, sum across k. Try i, k, j:
static void matmul_ikj(const float* A, const float* B, float* C) {
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];
}
}
}
}Now the innermost loop moves across j. For each fixed i and k:
B[k*N + 0] B[k*N + 1] B[k*N + 2] ...
C[i*N + 0] C[i*N + 1] C[i*N + 2] ...Both are contiguous. The same number of multiplies and adds, in a different order, with very different machine behavior. Measured in the main -O3 -ffast-math probe (Accelerate sample here: 0.787 ms, the slow mode, catch the fast mode and every ratio in this table grows about 5x):
Accelerate SGEMM 0.787 ms 341.09 GFLOP/s checksum -88.858404
naive ijk 140.085 ms 1.92 GFLOP/s checksum -88.858162
max abs diff vs BLAS 2.19345e-05, slowdown 178.0x
loop ikj 8.167 ms 32.87 GFLOP/s checksum -88.858404
max abs diff vs BLAS 0, slowdown 10.4xThat is the kind of result that feels like cheating the first time you see it.
We did not use threads. We did not use intrinsics. We did not use the GPU. We only moved the loops around, and the slowdown fell from 178x to 10.4x.
Low-level performance has a strange moral character this way. Sometimes the machine rewards deep hardware knowledge. Sometimes it rewards noticing that the innermost loop is walking down a column of row-major memory.
All six orders, and the one worse than naive
Two orders felt like enough to make the point, but three loops have six legal orderings, and I’d been quietly assuming ijk was the bad one. It isn’t even close to the bad one. I wrote all six into one probe and ran them at -O3 -ffast-math:
ijk 128.226 ms 2.09 GFLOP/s
ikj 8.342 ms 32.18 GFLOP/s
jik 129.992 ms 2.07 GFLOP/s
jki 441.067 ms 0.61 GFLOP/s
kij 10.187 ms 26.35 GFLOP/s
kji 458.897 ms 0.58 GFLOP/sThere’s a clean rule hiding in that mess, and it’s the innermost loop variable, nothing else. Sort the six by their inner index:
inner k (ijk, jik): ~128 ms B walked down a column
inner j (ikj, kij): ~9 ms B and C both contiguous
inner i (jki, kji): ~450 ms A and C both walked down columnsThe two fast orders (ikj, kij) are exactly the two whose innermost loop runs along j, where B[k*N+j] and C[i*N+j] are both contiguous. The two middling-awful orders (ijk, jik) put k innermost, so B strides by a full row each step, which is the original sin from the naive loop. And the two catastrophic orders (jki, kji) put i innermost, so both A[i*N+k] and C[i*N+j] stride down columns at once. That doubles the bad traversal, and it costs almost exactly what you’d fear: jki at 441 ms is 3.4x slower than the naive ijk I built the whole post around.
That reframes the opening. The naive ijk loop I keep calling slow isn’t the worst thing a beginner can type; it’s a middle-of-the-pack accident. If your first instinct had been “sum down a column, it’s a column of the answer,” you’d have written jki and been three times slower still, and you’d have blamed the CPU instead of the stride. The 512x512 matmul has a factor of about 55x hiding in nothing but which of three characters you put on the innermost for.
Transposing B
Another way to fix the column walk is to store a transposed copy of B. If BT[j][k] = B[k][j], then the dot product for C[i][j] becomes dot(A row i, BT row j), and both rows are contiguous. Code:
static void transpose(const float* B, float* BT) {
for (int i = 0; i < N; ++i) {
for (int j = 0; j < N; ++j) {
BT[j * N + i] = B[i * N + j];
}
}
}
static void matmul_transposed_b(const float* A, const float* BT, float* C) {
for (int i = 0; i < N; ++i) {
for (int j = 0; j < N; ++j) {
float sum = 0.0f;
const float* a = A + i * N;
const float* b = BT + j * N;
for (int k = 0; k < N; ++k) {
sum += a[k] * b[k];
}
C[i * N + j] = sum;
}
}
}Main probe:
B transposed 6.780 ms 39.59 GFLOP/s checksum -88.858396
max abs diff vs BLAS 1.90735e-05, slowdown 8.6xStill not BLAS. Much less embarrassing. There’s a subtle trade, though: transposing costs time and memory. If you multiply by the same B once, the transpose may or may not pay. If you multiply by the same B many times, it probably does.
Performance is often amortization wearing a fake mustache.
Tiling, badly
At this point people learn the word “tiling.”
A tile is a small block of the matrix. Instead of trying to multiply the whole matrices in one sweep, you multiply blocks that fit better in cache.
The innocent idea:
do C block (0..31, 0..31)
using A block (0..31, 0..31)
and B block (0..31, 0..31)
then the next blockI wrote this:
static void matmul_tiled_ikj(const float* A, const float* B, float* C) {
constexpr int BS = 32;
for (int ii = 0; ii < N; ii += BS) {
for (int kk = 0; kk < N; kk += BS) {
for (int jj = 0; jj < N; jj += BS) {
for (int i = ii; i < ii + BS; ++i) {
for (int k = kk; k < kk + BS; ++k) {
float a = A[i * N + k];
for (int j = jj; j < jj + BS; ++j) {
C[i * N + j] += a * B[k * N + j];
}
}
}
}
}
}
}And it got worse:
tiled 32 26.885 ms 9.98 GFLOP/s checksum -88.858404
max abs diff vs BLAS 0, slowdown 34.2xThat regression is worth keeping.
The ikj loop was already using contiguous B and C in its innermost loop. My tiled version added loop overhead and chopped the contiguous stream into smaller runs. It did not pack panels. It did not use a register-blocked microkernel. It did not unroll enough to keep the vector units (the hardware lanes that chew several floats per instruction, next section) fed. It did not choose tile sizes based on actual cache and register pressure.
My first instinct was that 32 was just the wrong tile size, so I swept it. Same ikj body, block size from 8 up to 256, everything else identical:
untiled ikj 8.342 ms
tiled 8 67.497 ms
tiled 16 15.818 ms
tiled 32 14.178 ms
tiled 64 11.447 ms
tiled 128 9.082 ms
tiled 256 8.555 msThere is no good tile size. The curve is monotonic: every block smaller than the whole matrix is slower than not tiling at all, and it gets worse the smaller the tile, with tiled 8 an eightfold catastrophe. tiled 256 finally crawls back to within a hair of untiled, which makes sense, a 256-block is half the matrix, so it barely chops the j stream. The tiling never helps because it was never solving a problem I had. Cache blocking exists to keep a working set resident when the whole thing doesn’t fit in cache. But all three of my matrices are 3 MiB and this machine’s L2 is 16 MiB. The naive ikj stream already lives in L2. Blocking it into smaller pieces can only add loop overhead and shorten the contiguous inner run that the vectorizer and the prefetcher both love. I was packing a suitcase to cross a room.
That’s the honest version of “tiling made it worse”: not a subtle cache-associativity story, just a size mismatch. Tiling would start earning its keep at, say, N = 4096, where one matrix alone is 64 MiB and nothing stays resident. At N = 512 on this machine, the textbook optimization is a pure tax.
Real GEMM tiling is not:
put ii, jj, kk loops around the old loopIt is more like:
pack A panel into a friendly layout
pack B panel into a friendly layout
compute a small C block in registers
unroll the inner loop
use SIMD/FMA
write C back once
repeat at L1, L2, and maybe last-level-cache scalesWhat clang actually emitted
SIMD means single instruction, multiple data. Instead of adding one float to one float:
a0 + b0a SIMD instruction can operate on several lanes:
[a0 a1 a2 a3] + [b0 b1 b2 b3]On arm64 Apple machines, the common vector instruction set is NEON. A 128-bit NEON register can hold:
4 float32 values
2 float64 values
16 uint8 valuesYou can write NEON intrinsics (C functions that map one-to-one onto the vector instructions) yourself, but the first thing to check is whether the compiler already did the obvious thing. I compiled this plain C++ dot product:
extern "C" float dot_scalar(const float* a, const float* b, size_t n) {
float s = 0.0f;
for (size_t i = 0; i < n; ++i) {
s += a[i] * b[i];
}
return s;
}With:
clang++ -std=c++20 -O3 -ffast-math -S simd_probe.cpp -o simd_probe.sThe assembly did not stay scalar:
ldp q4, q5, [x9, #-32]
ldp q6, q7, [x9], #64
ldp q16, q17, [x10, #-32]
ldp q18, q19, [x10], #64
fmla.4s v0, v16, v4
fmla.4s v1, v17, v5
fmla.4s v2, v18, v6
fmla.4s v3, v19, v7
subs x11, x11, #16
b.ne LBB0_5Reading arm64 here: ldp loads a pair of registers at once; q4 and v4 are two names for the same 128-bit vector register (q when moved whole, v when used lane-wise, .4s meaning four single-precision lanes); x registers hold 64-bit integers such as addresses and counters; subs subtracts and sets flags; b.ne branches back if the count hasn’t hit zero.
fmla.4s is a vector fused multiply-add on four single-precision lanes. It does:
v0 = v0 + v16 * v4for four float lanes at once.
There are four accumulators: v0, v1, v2, v3. That is not cosmetic. A CPU pipeline can often start a new instruction before the previous one has fully finished, but a single accumulator creates a dependency chain:
s = s + a0*b0
s = s + a1*b1 waits on previous s
s = s + a2*b2 waits againMultiple accumulators let the compiler keep independent work in flight:
s0 accumulates one stream
s1 accumulates another
s2 accumulates another
s3 accumulates anotherThen it adds them together at the end.
So by the time it runs, the original C loop isn’t “one instruction at a time.” A compiler may turn it into a small schedule of vector loads, vector FMAs, scalar cleanup code, and horizontal reductions (adding the lanes of one vector together into a single value).
That word fused deserves its ten seconds. FMA means fused multiply-add: a * b + c as one instruction, rounded once at the end instead of after the multiply and again after the add, faster, and very slightly different numerically. In dot products and matrix multiply the hot operation is sum += a * b, which is exactly what FMA wants. The scalar cleanup path in the same assembly used:
fmadd s0, s2, s1, s0The vector path used:
fmla.4s v0, v16, v4So even before we write a GPU kernel, the CPU compiler is already trying to map the loop onto the hardware’s multiply-add machinery.
That raises an annoying question. If clang can do this, why is the simple matrix multiply still far behind BLAS? Because SIMD is only one layer. The compiler can vectorize a local pattern; it doesn’t automatically become a GEMM implementation with packing, register blocking, cache blocking, architecture-specific kernels, prefetch choices, and years of tuning.
The fast-math bargain
The assembly above came from -O3 -ffast-math. -O3 asks for aggressive optimization; -ffast-math asks for a different floating-point contract. Without fast-math, the compiler has to respect more of the source program’s floating-point behavior: NaNs, signed zero, infinities, rounding, operation order, floating-point exceptions. With fast-math, it can treat floating-point arithmetic more like real-number algebra. That unlocks transformations such as reassociation, (a + b) + c becoming a + (b + c), which is boring for integers and can change the result for floating point.
In a performance post, fast-math is useful because it lets us see the machine’s vector path more clearly. In numerical software, it’s a decision: it can be exactly right, and it can break a carefully written stability trick. Serious libraries often ship multiple kernels or code paths for exactly this reason. They don’t merely ask what’s fastest; they ask what contract the caller requested, what accuracy is acceptable, and what hardware path exists.
Where the vectorizer switches on
I’d been compiling everything at -O3 -ffast-math and quietly crediting fast-math for the speed. So I compiled the same ikj loop at each optimization level to see where the win actually lives:
ikj at -O0 127.728 ms 2.10 GFLOP/s
ikj at -O1 35.781 ms 7.50 GFLOP/s
ikj at -O2 8.322 ms 32.26 GFLOP/s
ikj at -O3 8.310 ms 32.30 GFLOP/s
ikj at -O3 -ffast-math 8.319 ms 32.27 GFLOP/sThe step that matters is -O1 to -O2: 35.8 ms collapses to 8.3 ms, a 4.3x jump, which is suspiciously close to the four-float NEON width. That’s the auto-vectorizer turning on. -O3 adds nothing here, and -ffast-math adds nothing either, because the ikj inner loop is C[j] += a * B[j], a plain axpy with no reduction to reassociate; the compiler is allowed to vectorize it without any fast-math permission. Fast-math was a red herring for this loop the whole time. I’d been paying for a floating-point contract I wasn’t using.
The other thing that ladder kills is the beginner’s mental model where -O0 is “the real code” and -O3 is “the same code but the CPU tries harder.” At -O0 the ikj loop is 15x slower than at -O2. That isn’t the CPU trying harder; it’s an entirely different instruction stream, scalar loads and scalar FMAs one at a time, no vectors at all. The C you wrote is a suggestion. The binary is a negotiation, and the optimization level is one of the terms.
Chasing the zero
Earlier I flagged the strange bit-exact 0 diffs and admitted I hadn’t chased them. That optimization ladder handed me a way in, so I went back and chased.
The guess was that the loops matching Accelerate bit-for-bit are exactly the ones that still add up each output element in strict increasing-k order, and that reassociating the k-sum, which is what a vectorizer does when it splits a reduction across lanes, would move the low bits. The transposed-B loop is the cleanest test, because its inner loop is a genuine dot product over k, exactly the reduction a compiler wants to vectorize. So I watched its diff against Accelerate across the same optimization ladder:
B transposed at -O0 diff 0
B transposed at -O1 diff 0
B transposed at -O2 diff 8.58307e-06
B transposed at -O3 diff 8.58307e-06There it is. Bit-exact at -O0 and -O1, and then at exactly the level where the vectorizer switches on, -O2, the diff jumps off zero. Nothing in my source changed. The only thing that changed is that clang started accumulating that dot product in four NEON lanes instead of one running sum, and four partial sums added at the end is a different order of floating-point additions than one sum built left to right. The ikj loops stay bit-exact at every level because their vectorization is across j, which never touches the order in which a single C[i][j] accumulates its k terms.
To make sure it was the reassociation and not some fast-math side effect, I hand-wrote a transposed-B variant with four explicit accumulators, the tree the vectorizer builds, and compiled it with fast-math off so the compiler couldn’t reorder anything on its own:
transp split4 at -O3 (no fast-math) diff 2.00272e-05Nonzero, and off by a slightly different amount than the compiler’s own four-lane version, 2.0e-05 versus 8.6e-06. Three separate knobs, the optimization level, the fast-math flag, and a hand-written accumulator tree, all move the diff off zero the moment the k-sum stops being one strict running total, and they move it by different amounts because they build different trees.
That’s as far as the evidence takes me, and it’s worth being precise about where it stops. I’ve shown that my strict-k-order loops match Accelerate and every reassociated order doesn’t, which strongly suggests Accelerate is also accumulating each output in strict-k order with a single fused multiply-add chain. But “strongly suggests” is not “I read the kernel.” The fact that my split4 tree lands 2.0e-05 away rather than 0 even says Accelerate is not doing my four-wide tree; whatever order it uses, it happens to line up with the plain sequential one for this size. I can see the shape of the answer. I still can’t see inside Accelerate.
The paranoia in the prologue
There’s one more assembly detail worth seeing. Consider SAXPY (single-precision a*x plus y, the classic BLAS level-1 loop):
extern "C" void saxpy(float* y, const float* x, float a, size_t n) {
for (size_t i = 0; i < n; ++i) {
y[i] = a * x[i] + y[i];
}
}The assembly contains this before the vector loop:
lsl x8, x2, #2
add x9, x1, x8
cmp x9, x0
b.ls LBB1_8
add x8, x0, x8
cmp x8, x1
b.ls LBB1_8(lsl is a left shift, the count times 4 to get bytes, and b.ls branches if lower-or-same.)
That is a runtime overlap check.
The compiler is asking:
do x and y point to overlapping memory?If they overlap in a dangerous way, vectorizing might change the meaning. C and C++ pointers are powerful enough to make optimization harder.
You can sometimes help the compiler with restrict in C, or compiler-specific equivalents in C++, or by designing APIs that make aliasing impossible. But that’s a promise, if you lie, the program can become wrong. “The compiler will optimize it” is only half a sentence. Optimize under which promises?
The roofline, with this machine’s numbers in it
The roofline model is a simple performance picture:
attainable FLOP/s = min(peak FLOP/s, memory_bandwidth * arithmetic_intensity)Arithmetic intensity means floating-point operations per byte moved. If a loop does very little math per byte loaded, it’s memory-bound: faster arithmetic units won’t help much, because they’re already waiting for data. If a loop does lots of math per byte, it can become compute-bound: memory delivers enough, and the arithmetic units become the limit. The trick is that arithmetic intensity is a property of the implementation, not only of the formula.
Matrix multiply has high potential intensity. In the ideal fantasy, for 512x512:
read A once 1 MiB
read B once 1 MiB
write C once 1 MiB
do 268M FLOPsThat would be almost:
268M FLOPs / 3 MiB = about 85 FLOPs per byteReal life isn’t that clean. You don’t magically read each matrix once; you load pieces into caches, evict them, reload them, write partial results, and fight cache associativity (how many places a given address is allowed to land in the cache, no relation to the floating-point associativity coming up later), prefetch behavior, vector width, register count, and loop overhead.
Now plug in this machine. The M4’s reported spec is about 120 GB/s of memory bandwidth. The ikj inner loop does 2 FLOPs per iteration and streams roughly 8 to 12 bytes for them, 4 bytes of B plus a read and a write of C, minus whatever the caches absorb. Call it an arithmetic intensity around 0.2 FLOPs per byte, which against 120 GB/s puts the roof at roughly 20-30 GFLOP/s. The measured ikj number is 32.87 GFLOP/s, sitting on that roof, maybe slightly through it because the C row stays hot in cache and stops costing DRAM traffic. The caches blur the edges of this arithmetic, but the conclusion doesn’t move: the loop lives on the memory roof, nowhere near a compute roof.
Accelerate lives somewhere else entirely. Its fast-mode samples in these tables run between about 1460 and 1820 GFLOP/s, far more than one core’s NEON pipes can plausibly deliver. That throughput points at the matrix unit: the hardware Apple long shipped as the undocumented AMX coprocessor and now exposes architecturally as SME (the Scalable Matrix Extension) on the M4. The naive loop and the library aren’t even fighting over the same silicon.
Roofline thinking asks a blunt question: am I short on arithmetic, or short on delivered bytes? For matrix multiply, the best implementations work hard to make the answer “arithmetic.” The bad implementations accidentally make the answer “bytes.”
What BLAS is buying
A good SGEMM implementation doesn’t merely change the order of three loops. It usually has layers: outer blocking that picks big panels to fit a higher cache level; packing, which copies A and B panels into contiguous layouts the microkernel likes; a microkernel that computes a small C tile, often entirely in registers; aggressive SIMD/FMA; enough unrolling to expose independent work and hide latency; prefetching so data arrives before the core stalls; and clean handling of matrix edges that aren’t multiples of the tile sizes.
Packing is the part that feels suspicious at first. Why copy data before computing? Because a copy can be cheaper than a bad access pattern repeated millions of times. Suppose a microkernel wants a little rectangle of B in a layout where vector loads are easy: if the original matrix layout makes those loads awkward, the implementation copies a panel into a packed buffer once, then reuses it many times. Packing is the same amortization idea as transposing B, but more local and more tuned.
The microkernel then tries to keep a small block of C in registers:
C00 C01 C02 C03
C10 C11 C12 C13
C20 C21 C22 C23
C30 C31 C32 C33For each k, it loads a few values from A, a few from B, performs many FMAs, and only writes C back after accumulating enough work. Registers are the fastest storage you get; if partial sums can live there, you avoid repeatedly loading and storing C.
The simple ikj loop improved locality. BLAS changes the shape of the work so the machine can chew continuously.
The CPU plateau
Here is the main CPU table again, same slow-mode 0.787 ms Accelerate sample:
Accelerate SGEMM 0.787 ms 341.09 GFLOP/s checksum -88.858404
naive ijk 140.085 ms 1.92 GFLOP/s checksum -88.858162
max abs diff vs BLAS 2.19345e-05, slowdown 178.0x
loop ikj 8.167 ms 32.87 GFLOP/s checksum -88.858404
max abs diff vs BLAS 0, slowdown 10.4x
tiled 32 26.885 ms 9.98 GFLOP/s checksum -88.858404
max abs diff vs BLAS 0, slowdown 34.2x
B transposed 6.780 ms 39.59 GFLOP/s checksum -88.858396
max abs diff vs BLAS 1.90735e-05, slowdown 8.6xThe sequence is not monotonic.
bad loop order:
awful
better loop order:
huge win
naive tiling:
regression
transposed B:
another win
BLAS:
still far awayThat is a useful emotional correction. Performance work is not climbing stairs. It is probing a machine that has many cliffs and shelves.
And on this machine, the shelf I’m standing on looks close to the top for a loop this shape on one core: ikj sits pinned to the memory roof, the transpose bought about 20% more, and the thing running 9x to 46x faster than me, depending on which Accelerate mode woke up, is the AMX, matrix hardware inside the CPU that I can’t program directly.
The lever I’d been ignoring: more cores
Every number so far ran on one core, and this machine has ten, four performance and six efficiency. Before reaching for the GPU I did the obvious thing and split the ikj loop’s rows across std::threads. And I made a prediction I was fairly sure of: it would barely help. I’d just spent a whole section calling ikj memory-bound, pinned to the roofline. If it’s already waiting on memory, handing the same memory system to more cores shouldn’t buy much. Four cores splitting one starved pipe stay starved.
I was wrong, and being wrong told me something I’d gotten confused about:
T= 1 8.529 ms 31.47 GFLOP/s
T= 2 4.770 ms 56.27 GFLOP/s
T= 4 3.380 ms 79.43 GFLOP/s
T= 8 2.525 ms 106.32 GFLOP/s
T=10 2.120 ms 126.59 GFLOP/sIt scaled. Not perfectly, four cores gave 2.5x rather than 4x, and all ten gave about 4x rather than ten, but it moved, from 31 to 127 GFLOP/s. My prediction was built on a bad reading of my own roofline. The “memory roof” I’d drawn at ~120 GB/s is the DRAM roof, the road to the main memory chips. But I’d already established that all 3 MiB of this problem sit inside the 16 MiB L2. The ikj loop isn’t waiting on DRAM; it’s mostly being fed by cache, and every core brings its own private L1 while the four P-cores share that big L2. More cores means more aggregate cache bandwidth, so the work spreads instead of contending on one DRAM straw. The same fact that made tiling pointless, the whole thing fits in L2, is what let threading scale. I’d used the fact once and then forgotten it one section later.
The sublinear part is honest too. Four P-cores giving 2.5x, not 4x, is what you’d expect when they share one 16 MiB L2 and start competing for its bandwidth. And the jump from T=4 to T=10 is small because those extra six are efficiency cores, slower and on a different cache, so they add throughput but not in proportion to their count.
And here’s the part that keeps the GPU section honest: even all ten cores flat out reach 127 GFLOP/s, and Accelerate’s fast mode was doing 1600. Ten cores of hand-written NEON is still 12x behind the library. Parallelism was never the missing 40x. The missing 40x is the matrix unit, and I can’t get there by adding cores or moving loops. Which leaves the last lever that isn’t a library: the other processor in the package.
The other processor
A CPU function usually runs because one thread calls it:
matmul_ijk(A, B, C);A GPU kernel runs because the CPU asks the GPU to launch many instances of a function. The shape is:
CPU:
create buffers
compile or load GPU function
encode command
dispatch many threads
wait or continue
GPU:
run the same kernel function for many thread positionsIn Metal, a tiny vector-add kernel looks like this:
#include <metal_stdlib>
using namespace metal;
kernel void add(device const float* a [[buffer(0)]],
device const float* b [[buffer(1)]],
device float* c [[buffer(2)]],
uint gid [[thread_position_in_grid]]) {
c[gid] = a[gid] + b[gid];
}device means the pointer refers to device-visible memory, [[buffer(0)]] means the first buffer argument the CPU binds, and thread_position_in_grid is the global thread index, dispatch a million threads and each instance sees a different gid.
For matrix multiply, use a 2D grid:
kernel void matmul_naive(device const float* A [[buffer(0)]],
device const float* B [[buffer(1)]],
device float* C [[buffer(2)]],
constant uint& N [[buffer(3)]],
uint2 gid [[thread_position_in_grid]]) {
uint col = gid.x;
uint row = gid.y;
if (row >= N || col >= N) return;
float sum = 0.0f;
for (uint k = 0; k < N; ++k) {
sum += A[row * N + k] * B[k * N + col];
}
C[row * N + col] = sum;
}The kernel is the same idea as the C loop, but each GPU thread computes one output element. The CPU no longer loops over i and j. The grid does:
thread (0,0) computes C[0][0]
thread (1,0) computes C[0][1]
thread (2,0) computes C[0][2]
...
thread (511,511) computes C[511][511]So we launch 512 * 512 = 262144 GPU threads, and each thread loops over k = 0..511.
The host-side Metal ceremony
The kernel is only half the program. The CPU side has to create a device:
id<MTLDevice> device = MTLCreateSystemDefaultDevice();then compile or load the kernel:
id<MTLLibrary> library = [device newLibraryWithSource:ns_source
options:nil
error:&error];
id<MTLFunction> fn = [library newFunctionWithName:@"matmul_naive"];
id<MTLComputePipelineState> pso =
[device newComputePipelineStateWithFunction:fn error:&error];The pipeline state is the compiled GPU function plus the state Metal needs to dispatch it efficiently; you generally don’t want to rebuild it inside a hot loop. Then create buffers:
id<MTLBuffer> a = [device newBufferWithBytes:A.data()
length:A.size() * sizeof(float)
options:MTLResourceStorageModeShared];and encode a command:
id<MTLCommandBuffer> command = [queue commandBuffer];
id<MTLComputeCommandEncoder> enc = [command computeCommandEncoder];
[enc setComputePipelineState:pso];
[enc setBuffer:a offset:0 atIndex:0];
[enc setBuffer:b offset:0 atIndex:1];
[enc setBuffer:c offset:0 atIndex:2];
[enc setBuffer:nbuf offset:0 atIndex:3];
MTLSize grid = MTLSizeMake(N, N, 1);
MTLSize threads = MTLSizeMake(16, 16, 1);
[enc dispatchThreads:grid threadsPerThreadgroup:threads];
[enc endEncoding];
[command commit];
[command waitUntilCompleted];grid is the total logical problem shape, 512 by 512 output elements. threadsPerThreadgroup chunks that grid into 16 by 16 pieces, so Metal creates (512 / 16) * (512 / 16) = 1024 threadgroups of 16 * 16 = 256 threads each.
In this probe I call waitUntilCompleted because I’m measuring one command. In a real app or inference runtime, waiting after every kernel is often a performance bug, you want queues full enough that the GPU doesn’t starve.
Unified memory
On many CUDA systems, the CPU and GPU have physically separate memory pools: CPU RAM on one side, GPU VRAM or HBM (high-bandwidth memory, DRAM stacked next to the GPU die) on the other, and PCIe or NVLink (the interconnects between host and device) in between. A beginner CUDA program carries real ceremony: allocate host memory, allocate device memory, copy host to device, launch the kernel, copy device to host.
Apple Silicon has unified memory. The CPU and GPU share the same physical memory system. In my Metal probe I used:
MTLResourceOptions opts = MTLResourceStorageModeShared;
id<MTLBuffer> a = [device newBufferWithBytes:A.data()
length:A.size() * sizeof(float)
options:opts];That makes the experiment feel almost indecently easy. The CPU fills a vector, the GPU reads the buffer, and the CPU reads the result after the command completes.
It doesn’t mean memory is free. You still pay for command submission, synchronization, cache effects, bandwidth limits, and occupancy limits (occupancy, how much work the GPU keeps resident and ready to run; more on it below).
The first Metal result
There was a toolchain stumble first. This Mac has only the Command Line Tools, not full Xcode, and that installation doesn’t ship the metal shader compiler:
xcrun -find metal
xcrun: error: unable to find utility "metal", not a developer tool or in PATHSo there’s no offline compile to a .metallib on this machine. Instead the probe hands Metal the shader source at runtime, the newLibraryWithSource: call in the host code above, and the framework compiles it on the spot. That’s the honest reason the probe is one self-contained binary.
The run produced this (Accelerate sample: 0.184 ms, the fast mode; the checksums differ from the CPU tables because the Metal probe is a separate binary with its own fill):
Accelerate SGEMM 0.184 ms 1461.20 GFLOP/s checksum -15.717581
Metal naive 5.079 ms 52.85 GFLOP/s checksum -15.717581
max abs diff vs BLAS 0
Metal tiled 16 3.312 ms 81.05 GFLOP/s checksum -15.717581
max abs diff vs BLAS 0This is the right kind of humiliating.
The naive Metal kernel is much faster than the worst C loop. It is much slower than the tuned library. The tiled Metal kernel is better than naive Metal and still nowhere near the library.
The GPU did not save the bad algorithm. It gave the bad algorithm more parallel lanes to be bad on.
(And look at the diffs: both GPU kernels match Accelerate bit-for-bit too. Each thread accumulates its output in strict increasing-k order, the same order as my ikj loop, one more data point for the zero-diff pattern from the CPU section, and still not an explanation.)
For scale: the 10-core M4 GPU is commonly estimated at around 4 TFLOPS of fp32, Apple doesn’t publish the number, so treat it as an estimate, which puts the tiled kernel’s 81 GFLOP/s at roughly 2% of the machine.
The distinction that matters: mine is a learning kernel, one output per thread with simple threadgroup tiling, while a vendor GPU matrix library brings tuned, architecture-specific kernels and possibly specialized matrix hardware. The result doesn’t say GPU matmul is slower than CPU BLAS; it says a first GPU kernel is not a BLAS implementation. If I wanted the production Metal path, I’d compare against Metal Performance Shaders or another tuned matrix library, not against my 40-line kernel.
(If you learned this vocabulary from CUDA, as most of the internet did, the translation is mechanical: a CUDA block is a Metal threadgroup, __shared__ memory is the threadgroup address space, __syncthreads() is threadgroup_barrier, a warp is a SIMD-group, and global memory is the device address space. The names differ; the pressure, adjacent threads reading adjacent addresses, reuse on chip, enough work in flight, is identical. I’ll stay in Metal words.)
What I had to unlearn about the GPU
Sitting with that 52.85 GFLOP/s, my first instinct was wrong in an instructive way. I’d launched 262,144 threads; my CPU has ten cores; where’s my enormous ratio? The mistake is thinking a GPU is a CPU with more threads. A CPU core is good at one complicated thread: it predicts branches (guessing which way an if goes before the answer arrives), executes out of order (running whichever instructions are ready instead of waiting in program order), and keeps large caches close for low-latency scalar work. A GPU is good at many simpler lanes: lots of arithmetic units, lots of resident threads (threads parked on the GPU with their state loaded, ready to run), high memory bandwidth, throughput over latency.
When a CPU load misses cache, the core tries hard to hide it with out-of-order work. When a GPU load waits, the GPU often hides it by switching to other ready groups of threads. This is occupancy: how much work is resident and ready to run. Occupancy isn’t the same as performance, you can have high occupancy and do useless work, or lower occupancy and excellent reuse, but if occupancy is too low, the GPU can’t hide memory latency.
GPU performance is a balancing act between registers per thread, threadgroup memory per threadgroup, threads per threadgroup, resident threadgroups, memory coalescing (adjacent lanes combining their loads into one transaction, next section), and instruction mix. If each thread uses too many registers, fewer threads fit. If each threadgroup uses too much threadgroup memory, fewer threadgroups fit. If your memory accesses are scattered, the memory system does more transactions than necessary. If your branches diverge, some lanes sit idle while others run.
GPU threads also aren’t fully independent little CPU threads. Hardware executes lanes in groups, a SIMD-group in Metal, a warp in CUDA, and lanes in a group usually execute the same instruction together. If half the lanes take one branch of an if and half take the other, the hardware may run both paths with lanes masked off. That’s divergence. My kernels only branch on the edge check (if (row >= N || col >= N) return;), which every lane in a full tile answers the same way, so divergence isn’t what’s costing me here, but it becomes real in the reduction and attention shapes later, where row lengths and masks make lanes disagree. The ideal GPU loop looks boring: neighboring lanes, same instruction, neighboring addresses, many times.
Coalescing
Memory coalescing means adjacent lanes combine their memory requests. If a group of neighboring GPU threads reads:
B[k*N + col + 0]
B[k*N + col + 1]
B[k*N + col + 2]
B[k*N + col + 3]the hardware can fetch a contiguous chunk. If they read:
B[0]
B[2048]
B[4096]
B[6144]that’s a pile of separate transactions.
Our naive Metal kernel isn’t completely hopeless. For a fixed row and k, neighboring columns read neighboring B values, B[k * N + col] as col changes across adjacent threads. But those same threads all need related A values, many threads in a tile need overlapping pieces of A and B, and the naive kernel leaves that reuse mostly to caches.
Tiling, again, on a machine where it works
Metal calls shared scratch memory threadgroup memory. Inside a threadgroup, threads can cooperate:
load a tile of A into threadgroup memory
load a tile of B into threadgroup memory
barrier
compute using the loaded tile
barrier
move to the next k tileThe tiled Metal kernel I ran used 16x16 tiles:
constant uint TILE = 16;
kernel void matmul_tiled(device const float* A [[buffer(0)]],
device const float* B [[buffer(1)]],
device float* C [[buffer(2)]],
constant uint& N [[buffer(3)]],
uint2 tid [[thread_position_in_threadgroup]],
uint2 tg [[threadgroup_position_in_grid]]) {
threadgroup float As[TILE][TILE];
threadgroup float Bs[TILE][TILE];
uint row = tg.y * TILE + tid.y;
uint col = tg.x * TILE + tid.x;
float sum = 0.0f;
for (uint kk = 0; kk < N; kk += TILE) {
uint a_col = kk + tid.x;
uint b_row = kk + tid.y;
As[tid.y][tid.x] =
(row < N && a_col < N) ? A[row * N + a_col] : 0.0f;
Bs[tid.y][tid.x] =
(b_row < N && col < N) ? B[b_row * N + col] : 0.0f;
threadgroup_barrier(mem_flags::mem_threadgroup);
for (uint k = 0; k < TILE; ++k) {
sum += As[tid.y][k] * Bs[k][tid.x];
}
threadgroup_barrier(mem_flags::mem_threadgroup);
}
if (row < N && col < N) {
C[row * N + col] = sum;
}
}The barrier matters: without it, one thread could start using As before another thread has finished writing its part. The barrier says every thread in this threadgroup must arrive here before any thread continues.
The arithmetic for a single 16x16 output tile is the whole point. Naively, for one 16-wide k chunk:
256 output elements
each output reads 16 A values and 16 B values
256 * 32 = 8192 float readsWith threadgroup tiling:
load one 16x16 tile of A = 256 floats
load one 16x16 tile of B = 256 floats
total = 512 float readsTiling cut global reads 16x for that tile-shaped chunk of work. The measured speedup was not 16x:
Metal naive 5.079 ms 52.85 GFLOP/s
Metal tiled 16 3.312 ms 81.05 GFLOP/sbecause nothing is ever only one effect. The tiled kernel adds barriers, uses threadgroup memory, keeps the same one-output-per-thread structure, and still doesn’t use a serious register-blocked GEMM strategy. But it moved in the right direction.
The CPU tiled loop got worse. The GPU tiled kernel got better, from the same idea.
The bubble between tiles
The toy tiled kernel has a simple rhythm:
load tile
barrier
compute tile
barrier
load next tile
barrier
compute next tileThat creates bubbles: during the load phase the arithmetic units aren’t doing the main multiply-add work, and during the compute phase the next tile isn’t being loaded. More serious kernels overlap the phases, while computing tile 0, prefetch tile 1 into a second staging buffer, then swap roles, which usually goes by the name double buffering. The idea isn’t unique to GPUs; CPUs also prefetch and pipeline work. GPUs just make the staging explicit, because threadgroup memory is visible in the programming model.
There’s a cost: two buffers use more threadgroup memory, and more staging can mean fewer resident threadgroups.
The Metal tile sweep I said I hadn’t done
In the first pass I confessed, near the end, that I hadn’t tested multiple Metal tile sizes. That bothered me, so I went back and swept them: naive, then threadgroup tiles of 8, 16, and 32, all in one probe with ten warm reps each. The first thing the sweep showed had nothing to do with tiles:
first (cold) run Metal naive 3.995 ms
warm runs Metal naive 1.841 msThe original probe ran five reps and reported the naive kernel at 5.079 ms. This one, warmer, with ten reps, sees the same kernel at 1.841 ms, less than half. Nothing about the kernel changed. The GPU did. On the first dispatch the GPU is at a low clock and its caches are cold, and it takes a few kernels to ramp. This is the GPU cousin of the Accelerate bimodality: the number you print depends on what state the hardware was in when you asked, and a benchmark that doesn’t warm up is measuring the ramp, not the kernel. So the numbers below are the warm ones, and they are not comparable to the 5.079 ms in the original table, which was cold. Both are real. They answer slightly different questions.
With that caveat, the tile sweep:
Metal naive 1.841 ms 145.79 GFLOP/s
Metal tiled 8 1.199 ms 223.90 GFLOP/s
Metal tiled 16 1.237 ms 217.06 GFLOP/s
Metal tiled 32 1.135 ms 236.43 GFLOP/sOn the GPU the tile sweep is nearly flat, which is the opposite of the CPU. Every tile size beats naive, and 8, 16, and 32 all land within about 10% of each other. The CPU punished small tiles because they chopped a contiguous stream; the GPU barely cares, because the win here isn’t stream length, it’s that the threadgroup loads each value from global memory once and then reuses it out of on-chip memory. Any tile that’s big enough to get that reuse captures most of the benefit. The 32x32 tile uses 1024 threads per threadgroup, which is exactly this device’s maxThreadsPerThreadgroup, so that’s the ceiling for this shape, and it happens to be the best of the three by a whisker.
One thread, four outputs
Flat is a hint you’re doing the wrong kind of tiling. The one-output-per-thread structure is the real ceiling, not the tile size, because every thread still does one multiply-add per loaded pair and immediately needs the next pair. The next move, the one the CPU version couldn’t do because it had no on-chip scratch to exploit, is to make each thread own more than one output, so a single loaded value feeds several multiply-adds before it’s thrown away.
So I wrote one: a 32x32 output tile staged into threadgroup memory, but only 16x16 = 256 threads, and each thread computes a 2x2 block of C. In the inner loop each thread loads two A values and two B values and does four FMAs into four register accumulators:
float a0 = As[ty*2+0][k], a1 = As[ty*2+1][k];
float b0 = Bs[k][tx*2+0], b1 = Bs[k][tx*2+1];
acc[0][0] += a0*b0; acc[0][1] += a0*b1;
acc[1][0] += a1*b0; acc[1][1] += a1*b1;Four multiply-adds for four loaded values, versus one multiply-add per two loaded values in the one-output kernel. The reuse doubled, and it doubled in registers, which are the fastest storage on the chip:
Metal tiled 32 1.135 ms 236.43 GFLOP/s
Metal tiled16 rb2 0.661 ms 405.82 GFLOP/sThat is the biggest single GPU jump in the whole post, and it came from the one idea the toy kernels kept refusing: do more arithmetic per byte you bring on chip. At 406 GFLOP/s the register-blocked kernel is roughly 10% of the machine’s estimated ~4 TFLOPS, up from the 2% the toy tiled kernel was stuck at. The output still matches Accelerate bit-for-bit, maxdiff 0, because every thread still sums its four outputs in strict increasing-k order, the same pattern as everywhere else.
It’s also still not a serious GEMM. The block is 2x2; real kernels go to 4x4, 8x8, or wider, and they unroll the k loop, prefetch the next tile while computing the current one, and often hand the innermost fragment to simdgroup matrix instructions. But the shape of the climb is now visible in numbers I ran: naive 146, tiled 236, register-blocked 406, and Accelerate in this same probe still up near 1300 doing something I haven’t reached. Each step was the same lesson pointed one level deeper.
Why this still is not GEMM
The simplest Metal tiled kernel gives each thread one output element. That’s easy to understand, and it’s not how you chase peak GEMM. The 2x2 register-blocked kernel above took exactly one step past it, and a more serious GPU GEMM kernel keeps going: each thread computes many output values, or a small group of lanes cooperates on a fragment, to increase reuse in registers. Push my 2x2 to 4x4 and a thread might accumulate:
C[row + 0][col + 0]
C[row + 0][col + 1]
C[row + 1][col + 0]
C[row + 1][col + 1]Now one loaded A value feeds multiple B values and vice versa, and the partial sums stay in registers.
At higher performance tiers, kernels use specialized matrix instructions where available, Tensor Cores in CUDA land; SIMD-group and matrix-style operations in Metal, depending on GPU family and OS/toolchain support. I’m not walking through simdgroup matrix operations here: I haven’t run them, and they deserve their own post.
The conceptual stack stays:
global memory:
large, expensive
threadgroup memory:
smaller, on-chip, shared by cooperating threads
registers:
private, fastest, limited
matrix/SIMD execution units:
want carefully shaped fragmentsThe kernel writer’s job is to move data down that stack, reuse it hard, and move results back up as rarely as possible.
Two on-chip budgets
Registers are wonderful until you use too many. If each thread computes more output elements, it needs more accumulators, one C value needs one, sixteen need sixteen, plus address variables, loop counters, loaded fragments, and temporaries. The compiler has a finite register budget, and if the kernel wants too much, it spills: values that should have lived in registers get stored and loaded from slower memory. So “compute more per thread” helps until it hurts.
Threadgroup memory is the same trap one level up. A 16x16 tile of float is 16 * 16 * 4 = 1024 bytes, so the A and B tiles together cost 2 KiB, modest. At 32x32 the pair costs 8 KiB; at 64x64, 32 KiB. Larger tiles increase reuse, and they also consume more on-chip memory per threadgroup, which can reduce how many threadgroups are resident, demand more registers, and introduce bank conflicts (multiple threads hitting the same slice of on-chip memory at once, serializing the access) depending on the hardware.
Kernel tuning feels empirical for exactly this reason. You have a model, increase reuse, improve coalescing, raise arithmetic intensity, hide latency, and then a tile size crosses a register threshold and performance falls off a ledge.
The other kernel shape
Everything the matmul taught me, tile, stage, reuse, respect the memory stack, covers half the kernels that matter in ML. The other half compress instead of combining: a reduction takes many values and produces fewer, a sum, a max, an argmax (the index of the largest value). A reduction was hiding inside my matmul all along: each output element is a small dot-product reduction over k. The difference is that the matmul had 262,144 independent reductions to hand out, one per thread, while “sum this array” has one result and a million inputs, so the parallelism has to come from somewhere else. The beginner CPU version:
float sum = 0.0f;
for (int i = 0; i < n; ++i) {
sum += x[i];
}One dependency chain. One accumulator.
The GPU version can’t have every thread do sum += x[i] on one shared variable; that would serialize or race. Instead, use a tree:
thread 0 sums x[0], x[256], x[512], ...
thread 1 sums x[1], x[257], x[513], ...
...
store partial sums in threadgroup memory
barrier
halve active threads
add pairs
barrier
halve again
...
one result per threadgroupThen launch another pass, or use atomics (operations the hardware guarantees happen indivisibly, so concurrent updates can’t corrupt each other) carefully to combine the threadgroup results.
The tree changes the order of floating-point addition:
(((a + b) + c) + d)is not necessarily equal to:
(a + b) + (c + d)The GPU isn’t wrong here. Floating-point addition isn’t associative, so a faster reduction may produce a slightly different low-order result.
A stripped-down Metal-shaped reduction looks like this:
constant uint GROUP = 256;
kernel void reduce_sum(device const float* x [[buffer(0)]],
device float* partials [[buffer(1)]],
constant uint& N [[buffer(2)]],
uint tid [[thread_position_in_threadgroup]],
uint tg [[threadgroup_position_in_grid]]) {
threadgroup float scratch[GROUP];
uint base = tg * GROUP * 2;
uint i = base + tid;
float s = 0.0f;
if (i < N) s += x[i];
if (i + GROUP < N) s += x[i + GROUP];
scratch[tid] = s;
threadgroup_barrier(mem_flags::mem_threadgroup);
for (uint stride = GROUP / 2; stride > 0; stride >>= 1) {
if (tid < stride) {
scratch[tid] += scratch[tid + stride];
}
threadgroup_barrier(mem_flags::mem_threadgroup);
}
if (tid == 0) {
partials[tg] = scratch[0];
}
}It’s not a production reduction, the fixed scratch[256] assumes a maximum group size, it doesn’t use SIMD-group operations, and it still needs another pass to reduce partials, but it shows the shape: every thread loads a couple of elements and writes one partial, then the active threads halve each round with a barrier after every shared-memory update.
Once you can see that shape, softmax stops looking like a mysterious neural-network function. It’s a couple of reductions with exponentials in the middle.
Softmax
Softmax turns a row of scores into probabilities:
softmax(x_i) = exp(x_i) / sum_j exp(x_j)The numerically stable version subtracts the maximum:
m = max(x)
softmax(x_i) = exp(x_i - m) / sum_j exp(x_j - m)Subtracting m doesn’t change the mathematical result; it prevents overflow when scores are large.
A straightforward GPU softmax row kernel has three phases, reduce the row to its max, compute exp(x_i - max) and reduce to a sum, then write the normalized values, which means the same row may be read multiple times. For a short row, launch overhead and underfilled hardware can dominate. For a long row, reductions and memory traffic dominate. For a batch of rows, you must decide how many rows per threadgroup, how many threads per row, whether to vectorize loads, and where to keep partial results.
Softmax looks like a formula. As a kernel, it’s a max reduction, an exponential, a sum reduction, and a normalization write, and every one of those has a memory shape and a parallelism shape.
Attention
If you’ve never touched a transformer, here’s the minimum. A transformer processes a sequence of L positions (tokens), and for each position it computes three learned vectors of length D: a query Q, a key K, and a value V. Attention scores every pair of positions by dotting queries against keys, pushes each row of scores through softmax, and uses the resulting weights to mix the values together.
With that, transformer attention is often written as:
S = Q K^T
P = softmax(S)
O = P VThe formula is wonderfully clear and dangerously literal. Suppose sequence length is L and head dimension (the length of each query/key/value vector) is D. Shapes:
Q: L x D
K: L x D
V: L x D
S: L x L
P: L x L
O: L x DThe scary part is S and P: they’re L x L. For long context, materializing the whole score matrix can become the bottleneck. You do a lot of math, yes, but you also write a large intermediate matrix to memory, read it back for softmax, write probabilities, read probabilities, then multiply by V.
FlashAttention’s core idea is not “approximate attention.” It is exact attention reorganized around IO: tile Q, K, and V; compute blocks of scores; maintain the softmax normalization online; never write the full S and P to global memory. It’s the 512x512 story with a nastier middle, instead of computing a huge intermediate, writing it, reading it back, transforming it, and writing another one, you keep partial state on chip, stream through the K/V blocks, and write the final O once. Attention kernels are where GEMM, softmax, reduction, numerical stability, and fusion stop being separate chapters.
The online softmax trick
Here is the problem inside tiled attention.
Softmax wants the maximum of the whole row:
m = max(all scores in row)
sum = sum(exp(score - m))But a tiled kernel only sees part of the row at a time.
Suppose we have already processed some blocks. We keep:
m_old:
maximum score seen so far
s_old:
sum of exp(score - m_old) for scores seen so far
o_old:
partially accumulated output, scaled to m_oldNow a new block arrives. It has:
m_block:
maximum score in this block
s_block:
sum of exp(score - m_block) in this block
o_block:
block contribution to the output, scaled to m_blockThe new maximum is:
m_new = max(m_old, m_block)The old sum was scaled relative to m_old. The block sum was scaled relative to m_block. To combine them, rescale both to m_new:
s_new =
s_old * exp(m_old - m_new)
+ s_block * exp(m_block - m_new)The output accumulator is rescaled the same way:
o_new =
o_old * exp(m_old - m_new)
+ o_block * exp(m_block - m_new)At the end:
O = o_final / s_finalThat is the heart of the IO-aware attention trick. You do not need the whole score row in memory at once. You need enough state to combine tiles as if you had seen the whole row.
Fusion
Kernel fusion means combining operations that would otherwise be separate kernels.
Without fusion:
kernel 1:
C = A @ B
write C to memory
kernel 2:
D = C + bias
read C, write D
kernel 3:
E = gelu(D)
read D, write EWith fusion:
kernel:
compute A @ B tile
add bias while C is still hot
apply gelu
write final E onceHere bias is a learned per-output offset added after the matrix multiply, and gelu is a smooth activation function applied elementwise, the details don’t matter for the memory story. Fusion saves memory traffic and launch overhead, and it can improve locality because intermediate values stay in registers or threadgroup memory.
But fusion can go too far. If the fused kernel uses too many registers, occupancy drops. If it becomes branchy, lanes diverge. If it combines operations with incompatible tiling needs, one part gets fast and another part suffers. Good fusion isn’t “put everything in one kernel”; it’s “don’t write and reread an intermediate unless you have a reason.” The reason might be debuggability, reuse by multiple consumers, register pressure, different optimal tile shapes, a library boundary, or an autograd boundary.
The other thing fusion buys is fewer launches. The GPU doesn’t run a kernel the moment the C++ reaches the call site, there’s command encoding, queueing, scheduling. For a big GEMM that overhead is a rounding error; for a tiny elementwise op it can dominate, and a model is full of tiny elementwise ops. If every add, activation, and mask launches its own kernel and round-trips full tensors through memory, the GPU spends its day on bookkeeping around trivial math. Framework compilers fuse what they can see; hand-written kernels pick up what they can’t, dynamic shapes, custom ops, layout constraints.
LayerNorm, same story
Layer normalization is usually written as:
mean = average(x)
variance = average((x - mean)^2)
y = gamma * (x - mean) / sqrt(variance + eps) + betawhere gamma and beta are learned scale-and-shift parameters. As separate kernels, that becomes a reduction to the mean, a reduction to the variance, and an elementwise normalize-and-affine pass, multiple reads of the same row and multiple launches. A fused-ish row kernel keeps the row’s work together:
load row chunk
reduce to mean
barrier
reduce to variance
barrier
write normalized valuesFor small hidden sizes, one threadgroup might handle one row. For larger hidden sizes, multiple threadgroups may cooperate and require an extra reduction pass. If the row fits, threadgroup memory can hold values between phases; if it doesn’t, the kernel may reload from global memory.
Variance also has a numerics problem: computed naively as
E[x^2] - E[x]^2it can suffer cancellation (subtracting two nearly equal numbers, which wipes out the leading digits and leaves mostly rounding error, the failure mode dissected in the numerics post). So the implementation choice isn’t just one pass or two, it’s memory traffic, reduction count, accumulation precision, and stability all at once. The same kernel can be fast and numerically sketchy, or stable and slower, or stable and fast if someone did the hard work.
The CPU is still part of the GPU program
The host code matters. A Metal program is host code and kernel code: buffers, pipeline state, queues, and synchronization on one side; thread indexing, address spaces, barriers, and arithmetic on the other. If the CPU waits after every tiny kernel, it kills overlap:
[command commit];
[command waitUntilCompleted];Waiting is fine for a benchmark probe and often bad in a real pipeline. You want the CPU to enqueue work while the GPU runs previous work, with data dependencies expressed through command buffer ordering, events, or resource usage, not through unnecessary full waits.
On a unified-memory machine it’s easy to forget this, because the buffers are so convenient. The CPU and GPU can see the same memory, but they aren’t the same executor. You still have to reason about when writes become visible and when work is complete.
What I did not measure
I did not profile with Xcode GPU counters.
I did not measure memory bandwidth directly.
I did not use SIMD-group matrix operations.
I did not compare against Metal Performance Shaders matrix multiplication.
I did not write a production GEMM.
Since the first draft I did knock two items off this list, and I’m leaving that visible rather than pretending they were always here: I swept the Metal tile sizes (8, 16, 32, all within 10% of each other) and I made each thread compute a 2x2 block, which was the biggest single GPU jump in the post. So the list is shorter than it was, which is the only honest direction for a list like this to move.
The rest of the omissions are not small. They are the difference between “I understand the first layers” and “I have a tuned kernel.” This post is the walk from a C loop to a GPU kernel, not from a C loop to beating Accelerate.
The useful thing is that the experiment has teeth now: loop order gave 17x on the CPU and there’s an order 3.4x worse than the naive one I started from, bad tiling regressed across every tile size because the whole problem already fits in L2, clang emitted NEON and FMA and only at -O2, threading scaled to 4x when I bet it wouldn’t, naive Metal wasn’t enough, threadgroup tiling helped, register blocking helped more, and BLAS stayed far ahead in whichever mode it woke up in. I even chased the two open threads part of the way down: the bit-exact diffs track the vectorizer’s reassociation of the k-sum exactly, off zero at -O2 and beyond, and the Accelerate slow mode simply refused to reappear across 30 fresh processes. Neither is closed. I can see the shape of both answers and I still can’t see inside Accelerate, and I can’t summon its slow mode on command. Those are the two threads I’m shipping open.
The open kernel
One admission first: most of the time, the correct move is to not write the kernel. Check whether a vendor library already does the thing, Accelerate and MPS here, cuBLAS and cuDNN (NVIDIA’s tuned linear-algebra and neural-net kernel libraries) elsewhere, count bytes and FLOPs to learn which wall you’re against, and measure against the library rather than a strawman. The honest reasons to write one anyway: the operation doesn’t exist in a library, the library boundary forces extra memory traffic, you can fuse neighbors, you need a weird layout, or you’re learning. “I suspect I can beat the vendor’s GEMM on my first try” is also a reason, but mostly for character development.
I’m in the learning column, and the next experiment isn’t mysterious. I already took the working Metal probe two steps past toy, the tile sweep and the 2x2 register block, so the remaining list is the part that still scares me a little:
done: try TILE = 8, 16, 32
done: make each thread compute a 2x2 C block
next: push the block to 4x4 and 8x8, watch for register spill
unroll the k loop
measure register pressure
measure threadgroup occupancy
try simdgroup reductions
compare against MPSMatrixMultiplication
profile with Xcode countersThe 4x4 step is where I expect this to stop being free. A 2x2 block needs four accumulators; 4x4 needs sixteen, plus the loaded fragments and addresses, and somewhere between there and 8x8 the compiler runs out of registers and starts spilling to memory, and the throughput curve that’s been going up, 146, 236, 406, turns back down. I haven’t found that ledge yet. It’s the next thing I’d measure.
Then stop treating matmul as the endpoint:
write a row softmax
fuse max/sum/normalize where possible
write a tiled attention forward pass
keep the online softmax state in registers
compare materialized attention vs tiled attentionSources and artifacts
Local scratch directory:
/private/tmp/kernel-notesLocal tools and versions:
Apple M4
Apple clang version 15.0.0 (clang-1500.3.9.4)
Target: arm64-apple-darwin24.6.0Local artifacts used:
matmul_probe.cpp: naiveijk, loop-interchangedikj, naive 32x32 tiled CPU matmul, transposed-B matmul, Acceleratecblas_sgemmbaselinematmul_probe_ext.cpp: all six loop orders, the tile-size sweep from 8 to 256, and the split-accumulator transposed-B variant for the zero-diff testblas_bimodal.cpp: best-of-20 SGEMM per process, run 30 times to chase the Accelerate slow modethreaded_ikj.cpp:ikjsplit acrossstd::threadat 1, 2, 4, 8, 10 threadssimd_probe.cppandsimd_probe.s: arm64 NEON/FMA assembly for dot product and SAXPYmetal_probe.mm: Objective-C++ host code with runtime-compiled Metal naive matmul and 16x16 threadgroup-tiled matmulmetal_probe_ext.mm: Metal naive, tiled at 8/16/32, and the 2x2 register-blocked kernel
References I leaned on:
- Apple, Metal overview
- Samuel Williams, Andrew Waterman, David Patterson, Roofline: An Insightful Visual Performance Model for Multicore Architectures
- Tri Dao et al., FlashAttention: Fast and Memory-Efficient Exact Attention with IO-Awareness
- Apple, MPSMatrixMultiplication