Scalar matmul baseline

One kernel, one shape

kernels/matmul_scalar.c holds a single function:

void matmul_f32(const float *A, const float *B, float *C, int M, int N, int K);

Row-major, C[i*N+j] += A[i*K+k] * B[k*N+j], three loops in i,j,k order, -O2, and nothing else. That is the job: this kernel exists to be beaten. The 512 size lives in one header constant, so the timer, the reference and the checker cannot drift onto different shapes.

The FLOP count is arithmetic rather than measurement: 2*M*N*K, which is 268,435,456 at 512 cubed. Every speedup you claim later divides by that.

A reference that can fail

Write ref/matmul_ref.c from the definition, in i,k,j order, in double precision, and never call the kernel from it. A checker that calls the kernel is the kernel compared against itself.

Two invariants catch most mistakes. The total sum identity sum(C) == sum(A) * sum(B) follows from distributing the double sum over the product. A 4x4 case with hand-computed entries catches what the sum misses.

Then mutate on purpose: change B[k*N+j] to B[j*N+k] and confirm the check fails. A check you have never watched fail is a check you cannot trust. Note the blind spot while you are there: the checksum is invariant under transposing B, so it passes that mutation.

Things that break

Comparing f32 output against an f32 reference at a tight tolerance fails on accumulation order, not on bugs. Keep the reference in double, pick a relative tolerance from the magnitude of your data, and leave -ffast-math off both sides. Results that differ by a few ULPs, and only at large K, are reassociation. Results that differ in the fourth digit are an index bug.