Project 02 — CPU Matmul Benchmark

Make the same matrix multiply 1,000× faster without changing the math — and explain every factor.

1. What you build#

A benchmark harness that multiplies two n × n float32 matrices six different ways, reports GFLOP/s for each, measures your machine’s memory bandwidth, and draws your CPU’s roofline with each implementation placed on it.

2. Why it matters#

GEMM is most of inference. The distance between a naive loop and a tuned BLAS — two to three orders of magnitude on the same CPU — is the clearest demonstration that performance is about how bytes move through the memory hierarchy, not about the operation count. You will re-live this exact lesson on the GPU in Sections VI and X.

3. Read first#

• I.08 — FLOPs, bandwidth, arithmetic intensity
• II.02 — Caches and the memory hierarchy
• II.04 — SIMD and vectorization
• III.02 — Matmul by hand

4. Spec#

Implement in Go, with one variant in C called through cgo:

v1  Go, loop order i-j-k
v2  Go, loop order i-k-j               (unit-stride inner loop)
v3  v2 + cache blocking (tile size as a parameter)
v4  v3 + goroutines across tiles
v5  C through cgo: v3 compiled with -O3 -march=native (auto-vectorized), then AVX2/NEON intrinsics
ref an optimized BLAS through cgo      (OpenBLAS / MKL / Accelerate `cblas_sgemm`)

Go’s compiler does not auto-vectorize, so v4 → v5 is where you see what SIMD is worth. That gap is the reason numerical libraries are written in C and assembly, and called from everything else.

Metric: GFLOP/s = 2·n³ / seconds / 1e9. Sweep n ∈ {64, 128, 256, 512, 1024, 2048, 4096}.

5. Milestones#

1. Harness. A timer that warms up, runs ≥5 repetitions, reports median. Verify every variant against v1 to within a small floating-point tolerance.
2. v0→v2. Observe that reordering two loops, with identical arithmetic, gives a multiple. Explain it with stride and cache lines.
3. Blocking. Sweep tile size 8…512. The best tile relates to your L1/L2 size — compute the prediction from lscpu before running.
4. Vectorize and thread. Inspect the assembly (objdump -d or Compiler Explorer) to confirm SIMD was actually emitted.
5. Bandwidth. Write a STREAM-style triad (a[i] = b[i] + s*c[i]) on arrays ≥4× your L3. Report GB/s single-threaded and all-threads.
6. Roofline. Peak GFLOP/s = cores × SIMD width × FMA units × 2 × GHz. Plot the compute ceiling and the bandwidth slope; place v1…v5 and GEMV on it.

6. Starter skeleton#

// v2: i-k-j — the inner loop walks B and C with unit stride
void matmul_ikj(int n, 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];
        }
}

// v3: blocked
void matmul_blocked(int n, int T, const float *A, const float *B, float *C) {
    for (int ii = 0; ii < n; ii += T)
      for (int kk = 0; kk < n; kk += T)
        for (int jj = 0; jj < n; jj += T)
          for (int i = ii; i < ii+T && i < n; i++)
            for (int k = kk; k < kk+T && k < n; k++) {
              float a = A[i*n + k];
              for (int j = jj; j < jj+T && j < n; j++)
                C[i*n + j] += a * B[k*n + j];
            }
}

func gflops(n int, d time.Duration) float64 { return 2 * math.Pow(float64(n), 3) / d.Seconds() / 1e9 }

7. What to measure#

All into numbers.md. You will reuse the bandwidth number in almost every later section.

8. Done when#

• All variants agree with v1 within tolerance.
• You have a speedup ladder with a one-line cause for each step.
• Your best hand-written version is within ~3-5× of BLAS, and you can say what BLAS does that you do not (packing, micro-kernels, prefetch).
• You have a roofline plot with a ridge point, and can say which side GEMV sits on.

9. Common pitfalls#

The compiler deleted your benchmark. If the result is unused, -O3 removes the loop. Consume the output (checksum).

Frequency scaling and turbo. Pin the governor or at least run long enough to stabilize.

Measuring allocation. calloc of a large matrix returns untouched pages; the first write page-faults. Touch memory before timing.

NUMA on big machines. First-touch decides placement. Use numactl or initialize in the same threads that compute (II.03).

Aliasing blocks vectorization. Use restrict.

10. Stretch goals#

• Add FP16/INT8 variants and see what happens to the roofline.
• Implement packing (copy tiles into contiguous buffers) as in BLIS and close the gap to BLAS.
• Reproduce the huge-page effect on TLB misses for n=4096 (II.06).
• Plot energy per multiply using RAPL counters.

11. Interview questions this project answers#

1. Why does loop order change performance when the operation count is identical?
2. What is the ridge point of a roofline and how did you find yours?
3. Why is GEMV memory-bound and GEMM compute-bound on the same machine?
4. Why does thread scaling of a memory-bound loop saturate before you run out of cores?

12. Next#

Project 03 — First GPU kernel