Tiled vs naive matrix multiplication: loop order was the bigger win

Syndicated from the original on lkforge.com, where the same numbers drive two inline charts. The matrix tools it links to are at lkforge.com/tools/math.

Tiling — splitting the matrices into cache-sized blocks — is the famous trick for making matrix multiplication fast, and it's the first optimization most guides reach for. So I wrote all three versions in C, compiled them the same way (cc -O2), and timed them against the thing you'd actually use in practice: a BLAS library. The numbers said something the tutorials usually skip: the biggest free win isn't tiling at all.

Everything below was measured on Apple Silicon (macOS), Accelerate BLAS backend. Absolute rates vary by machine — the pattern is what reproduces.

Four ways to multiply the same two matrices

Same matrices, same machine, four implementations, GFLOP/s (higher is faster):

All four produced the identical result (checksums matched); only the speed differs. Three findings jump out at n=1024:

1. Loop order alone was 6.3× faster than naive (2.64 → 16.73), with no blocking at all.
2. Tiling was situational — my best tiled version (14.83) never beat the plain reordered loop on this hardware.
3. BLAS was another 24.5× ahead of my best loop (410 vs 16.73), and 155× the naive loop — and that gap is not cache tiling.

The free lunch: loop order

The naive order is for i, for j, for k: C[i][j] += A[i][k]·B[k][j]. The problem is the inner loop over k: A[i][k] walks along a row (contiguous, good), but B[k][j] walks down a column — every step jumps a full row ahead in memory. At n=1024 that's an 8 KB stride per multiply, so almost every access misses cache. The CPU spends its time waiting on memory, which is why naive sits at just 2.6 GFLOP/s. The multiplies are cheap; the memory stalls are the cost.

Swap the last two loops to for i, for k, for j and the inner loop becomes C[i][j] += a·B[k][j] with a = A[i][k] hoisted out. Now C[i][j] and B[k][j] both march along contiguous rows — unit stride. The hardware prefetcher and the compiler's auto-vectorizer both love that pattern, and the same arithmetic runs 6.3× faster. Not one line of blocking code; just touching memory in the order it's laid out.

Where tiling helps — and where it doesn't

Tiling computes the answer one small block at a time so a T×T patch of each matrix stays hot in cache while it's reused. That's a real effect, and you can watch it work — sweeping the tile size at n=1024:

Throughput climbs steadily as the block grows — bigger tiles reuse more before eviction, right up until three tiles (3·T²·8 bytes) stop fitting in L2. But here's the honest part: every one of those bars is below 16.73 — the same loops reordered with no tiling at all.

Why the famous trick underperformed: this CPU has a large L2, and the compiler already auto-vectorizes the clean reordered inner loop into wide SIMD stores. A hand-written blocked loop adds index arithmetic and loop overhead the reordered version doesn't pay, and the cache pressure it relieves wasn't the bottleneck yet. Tiling earns its keep on smaller-cache CPUs, at much larger matrices, or — crucially — as one layer of a multi-level blocking scheme (registers → L1 → L2), which is exactly what real libraries do. As a bolt-on to a loop that already streams cache-friendly, it's a wash.

The lesson isn't "tiling is useless." It's measure on your hardware; the famous optimization isn't automatically the winning one.

The 24.5× that's left over

My best hand loop reached 16.73 GFLOP/s. numpy reached about 410 — still 24.5× more, from the same silicon. None of that remaining gap is cache tiling. It's SIMD (one instruction multiplying a whole vector at once), register blocking (keeping a tiny sub-block of the result in registers so it never touches memory mid-accumulation), multithreading across cores, and microkernels hand-tuned per CPU. Cache tiling is one rung on that ladder, not the ladder.

The practical takeaway is the oldest one in performance work: for real matmul, call the library — and spend your own effort only where no library exists.

Reproduce it

The three loops are below. Compile with cc -O2, time each, and divide 2·n³ by the seconds for GFLOP/s. Your absolute numbers will differ by machine, but the shape — naive slow, reorder ~6×, tiling situational, BLAS far ahead — reproduces.

// naive — inner loop reads B down a column (stride n): cache-hostile
void ijk(const double*A,const double*B,double*C,int n){
  for(int i=0;i<n;i++)for(int j=0;j<n;j++){
    double s=0;
    for(int k=0;k<n;k++) s+=A[i*n+k]*B[k*n+j];
    C[i*n+j]=s;
  }
}

// reordered — inner loop over j is unit-stride: the 6.3x win
void ikj(const double*A,const double*B,double*C,int n){
  for(int i=0;i<n*n;i++)C[i]=0;
  for(int i=0;i<n;i++)for(int k=0;k<n;k++){
    double a=A[i*n+k];
    for(int j=0;j<n;j++) C[i*n+j]+=a*B[k*n+j];
  }
}

// tiled — same math, blocked into T-sized patches for cache reuse
void tiled(const double*A,const double*B,double*C,int n,int T){
  for(int i=0;i<n*n;i++)C[i]=0;
  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++){
      double a=A[i*n+k];
      for(int j=jj;j<jj+T&&j<n;j++) C[i*n+j]+=a*B[k*n+j];
    }
}

The BLAS baseline is one line of Python — C = A @ B with numpy — timed the same way.

Original writeup with the charts: lkforge.com/blog/tiled-vs-naive-matrix-multiplication · try the browser matrix tools: matrix multiply · all math tools.