Cache-Aware Matrix Multiplication: Loop Order, Tiles and the Register Block That Did the Work

The roofline post ended on a mismatch. Multiplying two $N\times N$ matrices would perform $N/16$ flops per byte of DRAM traffic if every matrix crossed the memory bus only once, 64 at $N=1024$ and far to the right of the 1.1 flop/byte ridge, yet the textbook inner loop loads two doubles for every multiply-add. This post measures how far four loop nests in plain C get between those two numbers: the textbook i-j-k order, the interchanged i-k-j order, a cache-tiled i-k-j, and the same tiles with a small block of $C$ held in registers. Every kernel was built twice with Apple clang 17.0.0 (clang-1700.0.13.3), once at -O3 -mcpu=native and once with -fno-vectorize -fno-slp-vectorize added, and run on one core of the same Apple-silicon laptop (Apple M1, 8 cores).

All matrices are row-major double. Every variant produced a result bit-identical to the i-j-k version, which is less surprising than it sounds: each one accumulates the $N$ products of an element of $C$ in the same order of $k$, and clang fuses every multiply-add in all of them.

Working sets and the choice of tile

sysctl -n hw.l1dcachesize hw.l2cachesize prints 65536 and 4194304 on this machine. Those are the caches of the efficiency cores; hw.perflevel0.l1dcachesize and hw.perflevel0.l2cachesize report 131072 and 12582912 for the performance cores, where these runs execute, and the 12 MiB L2 is shared by the four of them. The working sets:

In a tiled i-k-j loop the tile that is reused most is the tile of $B$: for every row of the $A$ tile, the kernel sweeps all $T\times T$ elements of it. A $64\times 64$ tile of doubles is 32 KiB, half of the smaller L1 that hw.l1dcachesize reports, which leaves room for the rows of $A$ and $C$ that stream past it; all three tiles together (96 KiB) fit in the 128 KiB L1 of a performance core. So $T=64$, as a compile-time constant. A sweep further down checks that choice.

Four loop nests

• i-j-k
• i-k-j
• cache-tiled
• register-tiled

• __attribute__((noinline)) void mm_ijk(int N, const double *restrict A,
                                        const double *restrict B, double *restrict C) {
      for (int i = 0; i < N; i++)
          for (int j = 0; j < N; j++) {
              double s = 0.0;
              for (int k = 0; k < N; k++) s += A[i * N + k] * B[k * N + j];
              C[i * N + j] = s;
          }
  }
  

  Each element of $C$ is one chain of $N$ dependent multiply-adds, and B[k * N + j] steps a full row ($8N$ bytes) at a time.

• __attribute__((noinline)) void mm_ikj(int N, const double *restrict A,
                                        const double *restrict B, double *restrict C) {
      memset(C, 0, sizeof(double) * N * N);
      for (int i = 0; i < N; i++)
          for (int k = 0; k < N; k++) {
              const double a = A[i * N + k];
              for (int j = 0; j < N; j++) C[i * N + j] += a * B[k * N + j];
          }
  }
  

  The inner loop is now an axpy over a row of $C$ and a row of $B$, both with stride 1.

• __attribute__((noinline)) void mm_blocked(int N, const double *restrict A,
                                            const double *restrict B, double *restrict C) {
      const int T = TILE;                       /* 64, a compile-time constant */
      memset(C, 0, sizeof(double) * N * N);
      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++)
                      for (int k = kk; k < kk + T; k++) {
                          const double a = A[i * N + k];
                          for (int j = jj; j < jj + T; j++) C[i * N + j] += a * B[k * N + j];
                      }
  }
  

  The same i-k-j statement, restricted to one $T\times T$ tile of each matrix at a time. $N$ must be a multiple of $T$.

• __attribute__((noinline)) void mm_blocked_reg(int N, const double *restrict A,
                                                const double *restrict B, double *restrict C) {
      const int T = TILE;
      memset(C, 0, sizeof(double) * N * N);
      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 += 4)
                      for (int j = jj; j < jj + T; j += 8) {
                          double c[4][8];                               /* 4x8 tile of C in registers */
                          for (int r = 0; r < 4; r++)
                              for (int q = 0; q < 8; q++) c[r][q] = C[(i + r) * N + j + q];
                          for (int k = kk; k < kk + T; k++)
                              for (int r = 0; r < 4; r++) {
                                  const double a = A[(i + r) * N + k];
                                  for (int q = 0; q < 8; q++) c[r][q] += a * B[k * N + j + q];
                              }
                          for (int r = 0; r < 4; r++)
                              for (int q = 0; q < 8; q++) C[(i + r) * N + j + q] = c[r][q];
                      }
  }
  

  Inside each cache tile, a $4\times 8$ block of $C$ stays in local variables for the whole $k$ range of the tile and touches memory only at the start and the end.

Compiling the four kernels with -Rpass='loop-vectorize|slp-vectorizer' -Rpass-missed=loop-vectorize gives the compiler's account (lines 11, 21, 34 and 49 are the innermost loops in the order of the tabs; line 55 is the final store of the register tile):

matmul.c:11:13: remark: the cost-model indicates that vectorization is not beneficial [-Rpass-missed=loop-vectorize]
matmul.c:21:13: remark: vectorized loop (vectorization width: 2, interleaved count: 4) [-Rpass=loop-vectorize]
matmul.c:34:25: remark: vectorized loop (vectorization width: 2, interleaved count: 4) [-Rpass=loop-vectorize]
matmul.c:49:25: remark: the cost-model indicates that vectorization is not beneficial [-Rpass-missed=loop-vectorize]
matmul.c:55:80: remark: Stores SLP vectorized with cost -1 and with tree size 7 [-Rpass=slp-vectorizer]
matmul.c:55:80: remark: Stores SLP vectorized with cost -2 and with tree size 7 [-Rpass=slp-vectorizer]
matmul.c:55:80: remark: Stores SLP vectorized with cost -4 and with tree size 7 [-Rpass=slp-vectorizer]

The i-j-k loop stays scalar: it is a floating-point reduction over a strided column. The register tile is not vectorised by the loop vectoriser at all; the SLP (superword) vectoriser starts from the eight stores of each row of the tile and packs the accumulators in pairs. The resulting $k$ loop is the kind of code one would write by hand:

	ldr	d24, [x30, x9]              ; A[i+0][k]
	ldp	q25, q26, [x26, #-32]       ; B[k][j..j+3]
	fmla.2d	v0, v25, v24[0]
	fmla.2d	v1, v26, v24[0]
	ldp	q27, q28, [x26]             ; B[k][j+4..j+7]
	fmla.2d	v2, v27, v24[0]
	fmla.2d	v3, v28, v24[0]
	ldr	d24, [x22, x9]              ; A[i+1][k]
	...                             ; 16 fmla.2d in total, one per pair of C values
	cmp	x9, #512
	b.ne	LBB3_7

Sixteen independent accumulators, the number the FMA-chain experiment in the roofline post needed to fill the pipes, fed by four scalar loads of $A$ and two paired vector loads (eight doubles) of $B$ for 64 flops.

Results

GFLOP/s ($2N^3$ divided by the time), median of three passes of nine timings each:

i-j-k runs at the speed of one dependency chain. 1.3 to 2.1 GFLOP/s, against 1.52 GFLOP/s for one dependent scalar FMA chain in the roofline post. Being a little above that at $N=256$ suggests the core overlaps the end of one element's chain with the start of the next. The decline with $N$ most likely comes from the column walk: every multiply-add touches a different 128-byte line of $B$, and those lines have to survive in cache until the next few values of $j$ reuse them, which gets harder as the column grows.

i-k-j is a triad in disguise. Swapping two loops gives 7 GFLOP/s scalar and about 12.5 vectorised, and the number barely moves between $N=256$ and $768$ although the working set grows from 1.5 to 13.5 MiB, so neither DRAM nor L2 is the limit. The inner loop loads a double of $B$, loads a double of $C$ and stores it back for every multiply-add: 24 bytes of L1 traffic per 2 flops, exactly the triad's intensity of 1/12. At the 196 GB/s that the triad reached from L1 that is 16.3 GFLOP/s, and i-k-j gets about three quarters of it.

Cache tiling helps when the sizes cooperate. The vectorised tiled kernel is 1.4 to 1.7 times faster than i-k-j at $N=256$, 384 and 768, but only 1.17 times at 512 and 1.11 times at 1024; the next section is about those two. In the scalar build tiling is a loss: clang keeps the plain axpy body, which now runs 64 iterations at a time instead of $N$, and the result is slower than scalar i-k-j. I have not established which part of the extra loop overhead costs the most.

The register tile is the step that matters. 32.5 to 46.1 GFLOP/s, 65 to 93 % of the 49.8 GFLOP/s NEON peak and roughly three times the vectorised i-k-j. Each $k$ step moves 96 bytes for 64 flops, eight times the axpy's intensity, so L1 bandwidth is no longer the limit. Built without vectorisation the same source manages only 11 to 13 GFLOP/s, and the assembly shows why: 32 scalar accumulators do not fit beside their operands in 32 floating-point registers, and the loop is full of register moves and stack spills.

Power-of-two sizes

The dips at 512 and 1024 have a classic cause. With $N=1024$, consecutive rows of a matrix are exactly 8 KiB apart, so the 64 rows of a tile map onto the same few sets of a set-associative cache and evict each other long before the cache is full: conflict misses rather than capacity misses. The standard remedy is to pad the leading dimension. The same kernels with rows of $N+8$ doubles (code otherwise unchanged):

Padding by 64 bytes makes the tiled kernel 1.5 to 1.6 times faster and the register-tiled one about 1.2 times, while i-k-j, which streams whole rows and reuses nothing across them, does not notice. That pattern is what conflict misses predict; I did not read the cache's miss counters to confirm it.

TIP

If a blocked kernel is mysteriously slow at exactly the sizes people like to benchmark (512, 1024, 2048), pad the leading dimension by a cache line before touching anything else. BLAS interfaces take a separate leading-dimension argument for this reason.

Choosing the tile, measured

The tile sweep below is at $N=768$, which every tile size divides and which has no power-of-two stride:

The tile derived from the cache sizes is not the fastest. For the register-tiled kernel it is 11 % behind $T=32$ and the curve is shallow. For the cache-tiled kernel the best tiles are much smaller than the cache would allow, and the assembly explains it: with $T=16$ or 32 known at compile time, clang unrolls the $j$ loop completely and keeps the whole row segment of $C$ in 8 or 16 vector registers for the entire $k$ loop, which has no stores left in it. The compiler built a $1\times T$ register tile on its own. At $T=64$ it attempts the same thing, needs 32 vector registers for $C$ alone, and the loop fills with spills and reloads. The scalar build does it too at $T=16$ (18.6 GFLOP/s, sixteen scalar accumulators) and not beyond. The cache sizes bounded the tile from above; the register file decided the answer.

WARNING

Make the tile a compile-time constant. An earlier version of this harness read $T$ from a global variable; the cache-tiled kernel then ran at 7.0 to 14.9 GFLOP/s instead of 11.3 to 20.9, because clang could neither unroll nor specialise the short inner loop. The register-tiled kernel, whose inner loop length does not depend on $T$, was within 13 % either way.

On the roofline

On this plot every variant sits at the same horizontal position, the ideal intensity $I_{\text{ideal}}=N/16=48$ flop/byte for $N=768$, forty-four times the ridge. That position is an assumption, not a measurement: I did not count DRAM traffic, and the seven variants certainly do not all move the same bytes, so the plot alone cannot rule out a DRAM limit. Two measurements make one unlikely. A 46 GB/s bus would hold the register-tiled kernel to 37 GFLOP/s only if it moved about sixty times the ideal traffic, and i-k-j runs at nearly the same speed at $N=256$, where all three matrices fit in L2, as at $N=768$, where they do not. What separates the variants by more than twenty times is then inside the core: the single-chain latency for i-j-k, the L1 traffic of the axpy for i-k-j and the cache-tiled loop, and, once a block of $C$ lives in registers, the FMA units themselves. The roofline still did its job: it pointed away from DRAM bandwidth and at the intensity of the innermost loop, measured against the cache that feeds it, as the number to raise. Production libraries (BLIS, OpenBLAS and the like) take the same steps further, with packed copies of the tiles so that the register kernel reads contiguous memory, several levels of tiling for L1, L2 and beyond, and kernels written in assembly.

Further reading

• Kazushige Goto and Robert van de Geijn, Anatomy of High-Performance Matrix Multiplication, ACM TOMS 2008: the layered tiling and packing behind production GEMM.
• Monica Lam, Edward Rothberg and Michael Wolf, The Cache Performance and Optimizations of Blocked Algorithms, ASPLOS 1991, including self-interference in tiles with unlucky strides.
• BLIS, a framework that exposes the micro-kernel and tiling structure explicitly.
• Ulrich Drepper, What Every Programmer Should Know About Memory, section 6.2, which works through a blocked matrix multiplication.
• Auto-Vectorization in LLVM, on the loop and SLP vectorisers whose remarks appear above.
• Sergey Slotin, Matrix Multiplication in Algorithms for Modern Hardware, the same progression on x86.