Cache Blocking and Loop Tiling

Cache blocking — also called loop tiling — restructures a loop nest so that it works on submatrices small enough to stay resident in a fast cache, reusing each loaded value many times before it is evicted. It is the single technique that takes a matrix multiply from a few per cent of a machine’s capability to most of it. This note builds the transformation from scratch in C on the machine described below, measuring every rung: naive, one-level tiling, two-level tiling, register tiling, and finally the packed five-loop structure of Goto and van de Geijn’s Anatomy of High-Performance Matrix Multiplication (ACM TOMS, 2008). The measured ladder ends at 256.9 GFLOP/s against a measured peak of 326.5 — 78.7% of the machine — from 7.0 GFLOP/s naive, a 36.6× speedup. Two findings along the way contradict what tutorials usually show: cache blocking on its own made the kernel slower, not faster, and the optimal block size was a tiny fraction of the cache, not most of it.

This note is the sequel to Naive Matmul and Why It Is Slow, which owns the baseline, the six loop permutations, and the cache-miss accounting that explains them. Read it first: everything here is measured against its ikj baseline, and its central finding — that the naive loop is instruction-bound rather than memory-bound — is the reason the blocking ladder behaves the way it does. Layers as Matrix Multiplication and im2col and Convolution as GEMM own why a neural network cares.

Uncertain

Every absolute throughput figure in this note is machine-specific and will not reproduce on your hardware. What reproduces is the shape of each result: that scalar tiling loses to no tiling, that register tiling is the large win, that the un-packed register tile collapses with problem size while the packed one does not, and that the block-size optimum sits well below cache capacity. The constants — 326.5 GFLOP/s peak, 256.9 achieved, the KC = 288 knee — belong to one AMD Zen 5 laptop part measured on 2026-08-22, performance governor, one core busy, gcc 16.1.1. Reason: throughput depends on clock, boost residency, memory configuration, compiler version and thermal state. To resolve for your machine: rebuild and re-run; the programs are self-contained C.

The Machine, and Why Its Numbers Are the Design Inputs

Block sizes are not folklore. They are derived from cache geometry, and the geometry is readable at runtime. The kernel exposes it under /sys/devices/system/cpu/cpu0/cache/index*/, documenting size as “the total cache size in kB”, ways_of_associativity as “degree of freedom in placing a particular block of memory in the cache”, number_of_sets as “total number of sets in the cache, a set is a collection of cache lines with the same cache index”, and coherency_line_size as “the minimum amount of data in bytes that gets transferred from memory to cache” (Linux ABI documentation). Read them; do not assume them.

$ for d in /sys/devices/system/cpu/cpu0/cache/index*; do
    echo "$(cat $d/level) $(cat $d/type): $(cat $d/size), \
$(cat $d/ways_of_associativity)-way, $(cat $d/number_of_sets) sets, \
$(cat $d/coherency_line_size) B lines"
  done

1 Data:        48K,  12-way,    64 sets, 64 B lines
1 Instruction: 32K,   8-way,    64 sets, 64 B lines
2 Unified:   1024K,  16-way,  1024 sets, 64 B lines
3 Unified:  32768K,  16-way, 32768 sets, 64 B lines   (shared by CPUs 0-7,16-23)

Every one of those numbers is a design input below. Alongside them, the quantities that had to be measured rather than read:

Quantity	Measured	Method
Core clock under load	5.046 GHz	dependent chain of 1-cycle addq
Peak fp32, AVX-512	326.52 GFLOP/s	12 independent FMA accumulators (2.02 FMA/cycle)
Peak fp32, AVX2 / scalar	163.32 / 20.42 GFLOP/s	same, narrower
vfmadd213ps latency	3.96 cycles	fully dependent FMA chain
L1d load-to-use latency	4.96 cycles	randomised pointer chase, 32 KiB
L2 load-to-use latency	16.55 cycles	512 KiB
L3 load-to-use latency	57.76 cycles	8 MiB
DRAM load-to-use latency	815.62 cycles	256 MiB
Read bandwidth L1 / L2 / L3 / DRAM	645 / 240 / 146 / 39.9 GB/s	8-accumulator streaming read

The latency ladder is the reason blocking exists at all: an L1 hit costs 5 cycles and a DRAM miss costs 816 — a factor of 164. Every value you can arrange to reuse from L1 instead of re-fetching from DRAM buys back two orders of magnitude.

flowchart LR
    R["<b>registers</b><br/>32 × zmm<br/>lat 0 · 326.5 GFLOP/s"]
    L1["<b>L1d</b> 48 KiB<br/>12-way, 64 sets<br/>lat 5 cyc · 645 GB/s"]
    L2["<b>L2</b> 1 MiB<br/>16-way, 1024 sets<br/>lat 17 cyc · 240 GB/s"]
    L3["<b>L3</b> 32 MiB<br/>16-way, per 8 cores<br/>lat 58 cyc · 146 GB/s"]
    D["<b>DRAM</b> 125 GiB<br/>lat 816 cyc · 39.9 GB/s"]
    R --> L1 --> L2 --> L3 --> D
    R -. "holds C micro-tile<br/><b>MR × NR</b>" .-> R
    L1 -. "holds B micro-panel<br/><b>KC × NR</b>" .-> L1
    L2 -. "holds packed A block<br/><b>MC × KC</b>" .-> L2
    L3 -. "holds packed B block<br/><b>KC × NC</b>" .-> L3

The measured memory hierarchy of the test machine, annotated with which GEMM operand each level is asked to hold. What it shows: the five block parameters MR, NR, KC, MC, NC are not free knobs — each one exists to make a specific operand fit a specific level. The insight to take: the design method is “assign one operand per level, then solve for the block size that makes it fit.” That is the whole of Goto and van de Geijn’s macro-level argument, and everything below is its arithmetic.

Mental Model: Blocking Converts Reuse-in-Principle into Reuse-in-Fact

A matrix multiply of size n performs 2n³ floating-point operations on 3n² data. The ratio grows without bound: for n = 4096 the available arithmetic intensity is 2n³ / (16n²) = n/8 = 512 FLOP per byte, against a DRAM ridge point of 326.52 / 39.9 = 8.18. There is, mathematically, sixty-two times more reuse available than the machine needs. Naive Matmul and Why It Is Slow measured what the naive loop actually realises: 2.01 bytes moved across the L1 boundary per FLOP — an intensity of 0.498, sixteen times below even the L1 ridge point of 0.51.

Blocking is the mechanical procedure for closing that gap. The idea in one sentence: instead of computing one full row of C before moving on — which requires walking all of B and evicting it — compute a small square of C completely, touching only the slice of A and B it needs, and keep that slice resident.

flowchart TB
    subgraph U["Unblocked: reuse exists but is not collected"]
        direction LR
        U1["for i, for k, for j:<br/>C[i][j] += A[i][k]·B[k][j]"] --> U2["row of B is walked n times<br/>but B is 4n² bytes<br/>→ evicted between uses"]
    end
    subgraph B["Blocked: reuse is collected"]
        direction LR
        B1["for ii, kk, jj (steps of BS)<br/>  for i, k, j within the tile"] --> B2["tile of B is 4·BS² bytes<br/>chosen so it FITS<br/>→ reused BS times from cache"]
    end
    U --> T["The transformation is<br/><b>strip-mining + interchange</b>:<br/>split each loop into<br/>outer-tile + inner-element,<br/>then hoist all three<br/>tile loops outward"]
    T --> B

Loop tiling as a transformation. What it shows: tiling does not change the set of n³ multiply-adds performed, nor their operands — it changes only the order, by splitting each loop into a tile loop and an element loop and moving all the tile loops outside. The insight to take: because it is a pure reordering, it is always legal for C += A·B (floating-point addition is not associative, so results differ in the last bits, but the algorithm is unchanged) — which means the only question is whether it pays, and that question is empirical.

Lam, Rothberg and Wolf put the goal precisely in 1991: “Instead of operating on entire rows or columns of an array, blocked algorithms operate on submatrices or blocks, so that data loaded into the faster levels of the memory hierarchy are reused.” Their paper is also the source of the taxonomy this note leans on — misses are intrinsic (compulsory), self-interference (“interference between elements of the same array variable”), or cross-interference (“interference between two different variables”).

Rung 1 and 2: One-Level and Two-Level Tiling — and Why They Lose

The textbook first step. Split all three loops by a block size L1BS, hoist the tile loops out, keep the fast ikj order inside:

#define BS 64
for (int ii = 0; ii < n; ii += BS)                    /* tile row of C   */
  for (int kk = 0; kk < n; kk += BS)                  /* tile of the sum */
    for (int jj = 0; jj < n; jj += BS) {              /* tile col of C   */
        int iM = MIN(ii+BS,n), kM = MIN(kk+BS,n), jM = MIN(jj+BS,n);
        for (int i = ii; i < iM; i++)                 /* --- element loops, */
          for (int k = kk; k < kM; k++) {             /*     unchanged ikj  */
              float a = A[(size_t)i*n + k];
              for (int j = jj; j < jM; j++)
                  C[(size_t)i*n + j] += a * B[(size_t)k*n + j];
          }
    }

The MIN clamps handle an n that is not a multiple of BS — the “edge tiles”. With BS = 64 the three tiles together occupy 3 × 64 × 64 × 4 = 48 KiB, exactly the L1 data cache. By the textbook argument this should be excellent. Sweeping BS at n = 1024, gcc -O2 -march=native:

BS	3 tiles occupy	GFLOP/s
8	0.75 KiB	5.14
16	3 KiB	6.88
24	6.75 KiB	4.99
32	12 KiB	5.62
48	27 KiB	6.13
64	48 KiB (= L1d)	6.64
96	108 KiB	7.35
128	192 KiB	7.77
192	432 KiB	7.98
256	768 KiB	8.12
—	untiled ikj baseline	9.02

Single-level tile-size sweep against the untiled baseline. What it shows: performance rises monotonically toward not tiling at all — BS = 256 is the best tiled result and it is still 10% slower than the plain loop. The insight to take: cache blocking, applied on its own to a good loop order, is a pure loss on this kernel. Every rung of the ladder that follows exists because of this negative result.

Adding a second level of tiling — an L2BS-sized outer tile wrapping the L1BS inner tile, the classic multi-level scheme — does not rescue it either. Across four problem sizes with L1BS = 64, L2BS = 256:

n	naive ikj	one-level tile	two-level tile
256	8.83	6.84	6.95
512	9.28	6.64	6.93
1024	9.35	6.79	6.89
2048	9.33	6.72	6.83
4096	7.02	6.43	6.49

The blocking ladder’s first two rungs, GFLOP/s. What it shows: one-level tiling costs 25–29% of the baseline’s throughput; a second level recovers about one point of that and never approaches the baseline. The insight to take: the effect is uniform across a 256× range of problem sizes — this is not a “tiles were the wrong size” result, it is structural.

Why blocking alone loses

The explanation is in the previous note’s central measurement: the naive ikj loop was never memory-bound. It was measured flat at 8.3–10.0 GFLOP/s from n = 64 (48 KiB working set, entirely L1-resident) to n = 2048 (48 MiB, DRAM-resident) — a 1,024× change in working set with no change in throughput. GCC leaves the inner loop scalar (vmovss / vfmadd213ss, one element per iteration), and it runs at 0.93 scalar FMAs per cycle against a scalar issue limit of 2.

So the loop is limited by issuing instructions, not by waiting for memory. Blocking reduces memory traffic. Reducing a resource you were not short of buys nothing — and tiling is not free: it adds three loop levels, three MIN computations per tile, and re-entry overhead into the innermost loop (n/BS)³ times instead of once. Measured cost of that overhead: 25%.

There is a sharper way to state the rule. Blocking pays exactly in proportion to how memory-bound your kernel already is. Two corollaries follow, and both are measurable.

Corollary one: blocking rescues a bad loop order dramatically. The ijk order — inner loop over k, walking a column of B with stride n — is memory-bound, catastrophically so. Tiling it at n = 1024:

BS	16	32	48	64	96	128	256	untiled
tiled ijk GFLOP/s	4.37	2.97	2.06	1.45	1.36	1.21	0.49	0.41

Tiling applied to the bad loop order. What it shows: a 10.7× speedup from tiling — the win every tutorial promises — but only because the untiled baseline was 0.41 GFLOP/s. The best tiled ijk (4.37) is still less than half the untiled ikj (9.02). The insight to take: fix the loop order before you reach for blocking. A tutorial that demonstrates tiling against an ijk baseline is measuring the cost of the baseline, not the value of tiling.

Corollary two: the optimal block size is tiny. The best BS above is 16 — three tiles of 16 × 16 × 4 bytes, 3 KiB, or 6% of the 48 KiB L1. Performance falls monotonically as the block grows toward cache capacity. This is a direct, on-2026-hardware confirmation of Lam, Rothberg and Wolf’s 1991 conclusion, which is worth quoting in full because it is still routinely ignored:

“The conventional wisdom of trying to use the entire cache, or even a fixed fraction of the cache, is incorrect. If a fixed block size is used for a given cache size, the block size that minimizes the expected number of cache misses is very small.”

The mechanism is their self-interference: a BS × BS tile of a row-major matrix is not contiguous — it is BS rows of BS floats, each 4n bytes apart. With n = 1024 that stride is 4096 bytes, which on this machine’s L1 (64 sets × 64 B = a 4096-byte set span) maps every row of the tile to the same cache set. A 64×64 tile then needs 64 ways in one set and has 12. It thrashes regardless of its nominal size. A 16×16 tile needs only 16 — still over 12, but close enough that most of it survives.

This is the observation that motivates the two remaining rungs. You cannot make blocking work by choosing a cleverer BS, because the problem is not size, it is layout. The fixes are to hold the tile in registers (which have no sets or ways) and to pack the operands into contiguous buffers (which have no stride).

Rung 3: Register Tiling

Register tiling — Goto and van de Geijn’s mr × nr micro-tile — computes a small block of C entirely in vector registers, so that across the whole k loop C is never loaded or stored. This attacks the actual bottleneck: the naive inner loop’s three memory operations per FMA become, per micro-tile, MR × NR / 16 FMAs against a handful of loads.

#define MR 6                        /* rows of C in the micro-tile      */
#define NR 32                       /* cols of C — 2 zmm registers wide */
#define NV (NR/16)                  /* vector registers per row         */
 
for (int i = 0; i + MR <= n; i += MR)
  for (int j = 0; j + NR <= n; j += NR) {
      __m512 c[MR][NV];                                  /* the micro-tile */
      #pragma GCC unroll 32
      for (int r = 0; r < MR; r++)
        #pragma GCC unroll 8
        for (int v = 0; v < NV; v++)
          c[r][v] = _mm512_loadu_ps(C + (size_t)(i+r)*n + j + 16*v);
 
      for (int k = 0; k < n; k++) {                      /* the ONLY memory loop */
          __m512 b[NV];
          #pragma GCC unroll 8
          for (int v = 0; v < NV; v++)
            b[v] = _mm512_loadu_ps(B + (size_t)k*n + j + 16*v);
          #pragma GCC unroll 32
          for (int r = 0; r < MR; r++) {
              __m512 a = _mm512_set1_ps(A[(size_t)(i+r)*n + k]);
              #pragma GCC unroll 8
              for (int v = 0; v < NV; v++)
                c[r][v] = _mm512_fmadd_ps(a, b[v], c[r][v]);
          }
      }
      /* store the micro-tile once, after the whole k loop */
      ...
  }

Line by line: __m512 is a 512-bit vector of sixteen fp32 lanes; _mm512_loadu_ps loads sixteen consecutive floats; _mm512_set1_ps broadcasts one scalar of A into all sixteen lanes; _mm512_fmadd_ps(a, b, c) computes a*b + c lane-wise in one instruction. The k loop performs MR × NV FMAs for NV vector loads and MR broadcasts — an arithmetic-to-memory ratio of roughly MR·NV / (NV + MR), versus the naive loop’s 1 FMA per 3 memory operations.

The #pragma GCC unroll lines are load-bearing, and finding that out was the largest single win

Written without them, this kernel measured 107 GFLOP/s. The disassembly shows why:

401140: vbroadcastss (%rdx),%zmm0            ; load A element
40114e: vmovaps %zmm0,%zmm1
401154: vfmadd213ps -0x40(%rax),%zmm2,%zmm0  ; accumulator FROM MEMORY
40115b: vfmadd213ps -0x80(%rax),%zmm3,%zmm1  ; accumulator FROM MEMORY
401162: vmovaps %zmm0,-0x40(%rax)            ; accumulator BACK TO MEMORY
401169: vmovaps %zmm1,-0x80(%rax)            ; accumulator BACK TO MEMORY

GCC did not unroll the r loop, so c[MR][NV] stayed an array — allocated on the stack (sub $0xb88,%rsp, 2,952 bytes) and loaded and stored on every single FMA. The register tile was not in registers. It ran acceptably only because the stack array sat in L1.

GCC documents #pragma GCC unroll n as controlling “how many times a loop should be unrolled”, noting it “must be placed immediately before a for, while or do loop… and applies only to the loop that follows”, and that “whether the loop is unrolled or not will be determined by target costing” (GCC manual, Loop-Specific Pragmas). Adding it forces full unrolling, which turns the array into scalar variables, which lets the register allocator keep them in zmm registers. The inner loop becomes 24 vfmadd instructions with no stack traffic, and:

MR × NR = 12 × 32	GFLOP/s
accumulators as an un-unrolled array (spilled to stack)	107.06
accumulators fully unrolled into registers	245.31

A 2.3× speedup from three pragmas, with the C semantics completely unchanged. This is the same failure mode that made a peak-FLOP microbenchmark for this note read 98 GFLOP/s instead of 326 until it was hand-unrolled. If you write a register-tiled kernel and it lands near a third of peak, disassemble it before doing anything else and look for (%rsp) in the inner loop.

Uncertain

Verify: that #pragma GCC unroll is the only reliable way to force this in portable C, and that clang behaves the same. Reason: the GCC manual states unrolling is subject to “target costing”, so the pragma is a request rather than a guarantee, and this note tested only gcc 16.1.1. clang 22.1.8 is installed on this machine but was not measured. To resolve: rebuild with clang -O2 -march=native and #pragma clang loop unroll(full), and compare disassembly. Note also that SIMD and Vectorization measured GCC’s __attribute__((vector_size(64))) generic-vector extension at 159.78 GFLOP/s, beating hand-written AVX-512 intrinsics on its kernel — so intrinsics are not automatically the fastest route to a micro-kernel.

Choosing MR × NR: the analytical model, and what it predicted

Low, Igual, Smith and Quintana-Ortí’s Analytical Modeling Is Enough for High-Performance BLIS (ACM TOMS) derives the micro-tile shape from instruction latency rather than search. Their argument: each element of the micro-tile Cr is updated once per k iteration by a dependent FMA, and Lvfma cycles must elapse between two dependent FMAs, during which Nvfma · Lvfma independent FMAs must be available to issue. Hence

  mr · nr  ≥  Nvec · Lvfma · Nvfma

walking the symbols: Nvec is lanes per vector (16 for fp32 in a 512-bit register), Lvfma is FMA latency in cycles, Nvfma is FMAs issued per cycle. On this machine all three were measured directly: Nvec = 16, Lvfma = 3.96 ≈ 4 cycles, Nvfma = 2.02 ≈ 2. So

  mr · nr  ≥  16 × 4 × 2  =  128 elements of C

The measurement, sweeping shapes at n = 1024 (all with unrolled accumulators):

MR	NR	tile elements	zmm for C	register tile	packed kernel
4	32	128	8	65.72	227.86
8	16	128	8	66.94	225.18
6	32	192	12	108.77	260.86
8	32	256	16	152.33	239.46
16	16	256	16	142.18	231.29
4	64	256	16	131.80	238.51
6	48	288	18	173.68	226.66
10	32	320	20	200.50	228.70
12	32	384	24	243.28	211.07
8	48	384	24	245.06	239.43
6	64	384	24	183.44	243.09
24	16	384	24	79.15	227.38
9	48	432	27	279.74	234.44
14	32	448	28	261.33	239.18
28	16	448	28	67.00	221.11

Micro-tile shape sweep, GFLOP/s at n = 1024. What it shows: the model’s lower bound is confirmed — both shapes at exactly 128 elements are the two worst results (65.72 and 66.94), and throughput climbs steeply above it, peaking at 9 × 48 = 432 elements / 27 zmm registers for 279.74 GFLOP/s, 85.7% of peak. It also shows the bound is necessary but not sufficient: 24 × 16 and 28 × 16 have 384 and 448 elements yet score 79 and 67, because a tile only one vector wide gives NV = 1, so every FMA needs its own broadcast — the load ports, not the FMA latency, become the limit. The insight to take: aim for mr ≈ nr above the latency bound, and keep at least two vectors of width; the model tells you the floor, the sweep finds the peak.

For reference, real libraries land in the same region. BLIS’s Skylake-X configuration uses fp32 MR = 32, NR = 12 — 384 elements, 24 registers (config/skx/bli_cntx_init_skx.c); its Zen 3 configuration uses MR = 6, NR = 16 for AVX2’s sixteen 256-bit registers (config/zen3/bli_cntx_init_zen3.c). Goto and van de Geijn state the register-budget rule directly: “Typically half the available registers are used for the storing mr × nr submatrix of C. This leaves the remaining registers for prefetching elements of Ã and B̃,” and “it can be shown that amortizing the cost of loading the registers is optimal when mr ≈ nr.”

Uncertain

Verify: that BLIS upstream ships no Zen 4 or Zen 5 configuration, so the closest published parameters for this machine are Zen 3’s (AVX2) and Skylake-X’s (AVX-512, Intel cache geometry). Reason: the BLIS config/ directory listing fetched on 2026-08-22 contains zen, zen2, zen3 but no zen4/zen5; AMD maintains a downstream fork (AOCL-BLAS) that does. To resolve: check AMD’s AOCL-BLAS release for bli_cntx_init_zen5.c and compare its MR/NR/MC/KC with the values measured here.

Rung 4: Packing, and the Five-Loop Structure

Register tiling produced the largest single jump — from 9.0 to 279.7 GFLOP/s, 31× — but it has a defect that only appears at scale. The un-packed micro-kernel reads A[(i+r)*n + k] for MR different rows on every k, i.e. MR separate streams each striding by 4n bytes, and B[k*n + j] striding likewise. Those are exactly the conflict-prone strided accesses that ruined single-level tiling.

Measured across problem size, best-shape register tile versus packed kernel:

n	working set	9×48 register tile	14×32 register tile	6×32 register tile	packed kernel
1024	12 MiB	271.02	259.65	107.97	251.01
2048	48 MiB	149.51	140.71	69.44	253.20
4096	192 MiB	30.56	45.70	21.85	254.37

The crossover that justifies packing, GFLOP/s. What it shows: the best un-packed register tile is the fastest thing on the machine at n = 1024 and collapses by 8.9× by n = 4096 (271.02 → 30.56), while the packed kernel is essentially flat (251.01 → 254.37) across a 16× growth in working set. The insight to take: packing is not a peak-performance optimisation — at small n it loses — it is a scaling optimisation. Its entire value is that it makes performance independent of problem size, which is exactly what a library must guarantee.

xychart-beta
    title "Register tiling vs packing as the problem outgrows cache"
    x-axis "matrix dimension n" [1024, 2048, 4096]
    y-axis "GFLOP/s (single core, peak = 326.5)" 0 --> 300
    line "register tile 9x48 (no packing)" [271.02, 149.51, 30.56]
    line "register tile 14x32 (no packing)" [259.65, 140.71, 45.70]
    line "packed micro-kernel 6x32" [251.01, 253.20, 254.37]

The same data as a chart. What it shows: two lines that start at the top and fall off a cliff, and one that starts slightly lower and does not move. The insight to take: if you benchmark a GEMM at one problem size you will draw the wrong conclusion. The un-packed kernel wins at n = 1024 and is 8× slower at n = 4096.

What packing actually does

Goto and van de Geijn identify the problem exactly: “The fundamental problem now is that A is typically a submatrix of a larger matrix, and therefore is not contiguous in memory. This in turn means that addressing it requires many more than the minimal number of TLB entries. The solution is to pack A in a contiguous work array, Ã.” They stress that a TLB (Translation Look-aside Buffer) miss is worse than a cache miss because “a TLB miss… causes the CPU to stall until the TLB has been updated… prefetching can mask a cache miss but not a TLB miss.”

Packing rewrites the operands into the exact order the micro-kernel will consume them:

  Packed A block (Ac), MC × KC, split into MC/MR micro-panels:

  panel 0                        panel 1
  ┌──────────────────────────┐   ┌──────────────────────────┐
  │ p=0: a[0] a[1] ... a[MR-1]│  │ p=0: a[0] ... a[MR-1]    │  ...
  │ p=1: a[0] a[1] ... a[MR-1]│  │ p=1: ...                 │
  │  ...  (KC rows)          │   │                          │
  └──────────────────────────┘   └──────────────────────────┘
   ← MR floats contiguous →       ← next panel starts at +MR·KC floats →

  Packed B block (Bc), KC × NC, split into NC/NR micro-panels:

  panel 0                              panel 1
  ┌────────────────────────────────┐   ┌───────────────────┐
  │ p=0: b[0] b[1] ....... b[NR-1] │   │ p=0: ...          │  ...
  │ p=1: b[0] b[1] ....... b[NR-1] │   │                   │
  │  ...  (KC rows)                │   │                   │
  └────────────────────────────────┘   └───────────────────┘
   ← NR floats contiguous = NV zmm loads →

Packed panel layout. An ASCII box diagram is used rather than mermaid because this is a memory-layout figure and mermaid has no construct for byte ordering within a buffer. What it shows: after packing, the micro-kernel’s entire k loop walks both Ac and Bc strictly forward, contiguously — ap += MR and bp += NR per step. The insight to take: stride is gone, so conflict misses are gone, TLB pressure is gone, and the hardware prefetcher works perfectly. Goto and van de Geijn call this “accessing data contiguously” and note their packing stores each mr × kc submatrix contiguously precisely so that the kernel can “stride strictly contiguously through memory.”

The cost is real but well amortised. Packing B costs O(KC · NC) copies amortised over O(MC · KC · NC) FLOPs; Goto and van de Geijn observe that “O(m) computations will be performed for every copied item,” and that “the overhead caused by the packing operations accounts almost exactly for the degradation in performance from the kernel curve to the dgemm curve.” That degradation is visible in the table above: 279.74 (bare kernel, n = 1024) versus 251.01 (packed) — about 10%, paid to buy size-independence.

The five loops

The full structure, as implemented and measured here, is Goto and van de Geijn’s gepp_blk_var1 / gebp_opt1 pair:

for (int jc = 0; jc < n; jc += NC) {            /* loop 5: L3 blocking   */
  for (int pc = 0; pc < n; pc += KC) {          /* loop 4: the k blocking */
    pack_B(...);                                /*   Bc: KC × NC → L3    */
    for (int ic = 0; ic < n; ic += MC) {        /* loop 3: L2 blocking   */
      pack_A(...);                              /*   Ac: MC × KC → L2    */
      for (int jr = 0; jr < nc; jr += NR)       /* loop 2               */
        for (int ir = 0; ir < mc; ir += MR)     /* loop 1               */
          micro(kc, &Ac[ir*kc], &Bc[jr*kc],     /* loop 0: the k loop,  */
                &C[(ic+ir)*n + jc+jr], n, ...); /*   MR×NR tile in regs */
    }
  }
}

flowchart TB
    L5["<b>loop 5</b> — jc, step NC<br/>partition C and B by columns"]
    L4["<b>loop 4</b> — pc, step KC<br/>partition A by cols, B by rows<br/><i>pack B → Bc (KC×NC), targets L3</i>"]
    L3["<b>loop 3</b> — ic, step MC<br/>partition C and A by rows<br/><i>pack A → Ac (MC×KC), targets L2</i>"]
    L2["<b>loop 2</b> — jr, step NR<br/>pick a B micro-panel<br/><i>KC×NR, targets L1</i>"]
    L1["<b>loop 1</b> — ir, step MR<br/>pick an A micro-panel"]
    L0["<b>micro-kernel</b> — k = 0..KC<br/>MR×NR tile of C lives in zmm registers<br/><i>no loads or stores of C at all</i>"]
    L5 --> L4 --> L3 --> L2 --> L1 --> L0
    L0 -.->|"2·MR·NR·KC flops<br/>for MR·NR memops"| L0

The five-loop GEMM. What it shows: each loop level exists to place one operand in one level of the hierarchy, with the two packing steps sitting exactly at the loops whose blocks they produce. The insight to take: the innermost loop performs 2·MR·NR·KC FLOPs while writing C only MR·NR times — Goto and van de Geijn’s point that “during the computation of Cj it is not necessary that elements of that submatrix remain in the L1 or even the L2 cache.” C is written once per micro-tile, no matter how large KC is, which is why KC should be as large as the L1 will allow.

Deriving KC and MC from the cache geometry, then checking

Low et al. derive kc by requiring that a newly loaded micro-panel of A evict exactly the previous one, leaving the B micro-panel resident. Their result:

  kc = C_Ar · N_L1 · C_L1 / (mr · Sdata)      with   C_Ar ≤ ⌊ (W_L1 − 1) / (1 + nr/mr) ⌋

Symbol by symbol: N_L1 = 64 sets, C_L1 = 64 B lines, W_L1 = 12 ways, Sdata = 4 bytes per fp32, and C_Ar is “the number of cache lines taken up by a micro-panel Ar in each set of the L1 cache”. The −1 reserves one line per set for C, because “at least one cache line must be used when the entries of a micro-tile Cr are loaded into the registers.” In their column-major convention mr is the vectorised dimension; in this note’s row-major kernel that is NR = 32, with MR = 6 playing the role of nr. So

  C_Ar ≤ ⌊ 11 / (1 + 6/32) ⌋ = ⌊ 11 / 1.1875 ⌋ = 9
  kc   = 9 × 64 × 64 / (32 × 4) = 36,864 / 128 = 288

The model predicts KC = 288. Measured, sweeping KC at n = 4096 with MR = 6, NR = 32, MC = 288:

KC	Ac footprint	GFLOP/s
64	72 KiB	133.85
128	144 KiB	204.52
192	216 KiB	239.51
256	288 KiB	247.35
288	324 KiB	249.20
320	360 KiB	253.67
384	432 KiB	254.56
448	504 KiB	253.75
512	576 KiB	256.96
768	864 KiB	253.70

KC sweep. What it shows: a steep climb (133.85 → 247.35 as KC goes 64 → 256, an 85% gain) that flattens into a broad plateau of 249–257 from KC = 288 onward. The analytical prediction of 288 sits precisely at the knee. The insight to take: the model gets you to within a few per cent in one calculation and with no search — which is the entire claim of the paper’s title. Beyond the knee the curve is flat, so exact tuning is not worth much; below it, every doubling matters.

MC is set so the packed Ac block occupies a fair share of L2. With KC = 384, MC · 384 · 4 ≤ 512 KiB (half of the 1 MiB L2) gives MC ≤ 341. Measured:

MC	Ac footprint	% of L2	GFLOP/s
48	72 KiB	7%	257.63
96	144 KiB	14%	259.20
192	288 KiB	28%	255.90
288	432 KiB	42%	257.17
384	576 KiB	56%	252.47
456	684 KiB	66%	250.87
576	864 KiB	84%	249.18
768	1,152 KiB	112%	244.08
1152	1,728 KiB	168%	248.92
2304	3,456 KiB	337%	250.00

MC sweep. What it shows: a shallow curve with the worst point at exactly the size where Ac stops fitting in L2 (112%, 244.08 GFLOP/s) and the best at 14% of L2. The insight to take: Lam et al.’s 1991 warning holds on 2026 silicon and at the outermost blocking level too — “trying to use the entire cache… is incorrect”, and the optimum is “very small”. The model’s job is to give the ceiling past which performance degrades; the optimum lives well below it.

Goto and van de Geijn reach the same conclusion from experience rather than simulation: “In practice, mc is typically chosen so that Ã only occupies about half of the smaller of (1) [the TLB-addressable memory] and (2) the L2 cache.” Half was the ceiling here; an eighth was the optimum.

The Complete Ladder

Final configuration MR=6, NR=32, KC=384, MC=96, NC=3072, gcc -O2 -Wall -Wextra -march=native, single core, best of three trials:

Rung	n = 1024	n = 2048	n = 4096	% of 326.5 peak (n=4096)
0 · naive ikj	9.02	8.78	7.02	2.2%
1 · one-level tiling	6.66	6.68	6.43	2.0%
2 · two-level tiling	6.76	6.75	6.49	2.0%
3 · register tiling (this shape)	106.06	64.80	23.49	7.2%
3 · register tiling (best shape 9×48)	271.02	149.51	30.56	9.4%
4 · packed micro-kernel	258.77	250.11	256.94	78.7%

From 7.02 to 256.94 GFLOP/s at n = 4096: 36.6×. Correctness was checked at every rung against the naive reference: rungs 0–3 are bit-identical, and the packed kernel’s maximum absolute error is 3.815 × 10⁻⁶ against a maximum reference magnitude of 8.977 — a relative error of 4.25 × 10⁻⁷, single-precision rounding, caused purely by summing k in a different order.

How far from peak, honestly

78.7% of the measured single-core peak. The remaining 21% is not mysterious; it is the packing copies (about 10%, matching Goto and van de Geijn’s observation), the edge-tile handling for n not divisible by MR/NR (this implementation falls back to a scalar path), the absence of software prefetching, and the read-modify-write of C once per KC block rather than accumulating in place.

For calibration, OpenBLAS on the same machine, same problem, single thread. The only OpenBLAS on this system lives inside the Steam Linux Runtime — version 0.3.5, built March 2019 — because there is no system OpenBLAS package installed. Its DYNAMIC_ARCH dispatcher does not recognise Zen 5’s CPUID and falls back to the oldest AMD kernel it has:

Dispatched kernel	n = 1024	n = 2048	n = 4096	% of peak
Barcelona (auto-selected; AMD K10, 2007, SSE-only)	60.65	62.93	63.24	19.4%
Zen (AVX2 + FMA)	127.94	131.73	135.27	41.4%
Haswell (AVX2 + FMA)	135.10	134.23	136.16	41.7%
SkylakeX (AVX-512)	224.95	225.08	227.74	69.8%
this note’s packed kernel	258.77	250.11	256.94	78.7%

OpenBLAS 0.3.5 by dispatched micro-kernel, forced with OPENBLAS_CORETYPE. What it shows: a 3.6× spread (63.24 → 227.74 GFLOP/s) from CPU dispatch alone, on identical code and identical data. Left to itself the library picked a 2007 SSE kernel and ran at 19% of the machine. The insight to take: always print openblas_get_corename() before trusting a BLAS baseline.

Uncertain

Verify: whether a current OpenBLAS or BLIS with a Zen 4/Zen 5 micro-kernel would beat this note’s 256.94 GFLOP/s. Reason: the only BLAS available on this machine is OpenBLAS 0.3.5 from 2019, which predates Zen 4 entirely; its best available kernel (SkylakeX) is tuned for Intel’s 32 KiB / 8-way L1 and 1 MiB L2, not this part’s 48 KiB / 12-way L1. Beating a seven-year-old mis-targeted library is a weak claim and must not be read as “hand-written C beats OpenBLAS”. To resolve: build OpenBLAS ≥ 0.3.26 or AMD AOCL-BLAS from source with the native target and re-run; expect it to win.

Failure Modes and How to Diagnose Them

The register tile is not in registers. By far the most common and most expensive failure — worth 2.3× here. Symptom: the kernel lands at 30–35% of peak and does not respond to block-size tuning. Diagnosis: objdump -d the micro-kernel and look for (%rsp) operands or sub $0x...,%rsp with a large constant in the prologue. Fix: force full unrolling of the accumulator loops, or macro-expand them into named variables.

Too many accumulators. MR · NR / 16 + NV + 1 must fit in 32 zmm registers. Exceeding it spills, silently. 14 × 32 needs 28 for C plus 3 = 31, and it worked; 9 × 48 needs 27 plus 4 = 31, and it was the fastest. Push past that and throughput falls off a cliff.

Tiling a kernel that was not memory-bound. Symptom: tiling makes it slower, monotonically, at every block size. Diagnosis: measure the untiled kernel at a working set that fits entirely in L1 (for a matmul, n such that 12n² ≤ L1 size). If it is no faster there than at n = 2048, the kernel is instruction-bound and blocking cannot help it. Fix the instruction stream first — vectorise, then register-tile — and only then block.

Benchmarking at one problem size. The crossover table above is the cautionary tale: the un-packed kernel wins by 8% at n = 1024 and loses by 8.3× at n = 4096. Always sweep at least three sizes spanning L2-resident, L3-resident and DRAM-resident.

Power-of-two leading dimensions. Naive Matmul and Why It Is Slow measured a 12.1× swing at n = 1024 from padding rows by eight floats. Packing makes this moot for the packed operands, but C is still written in place at the user’s leading dimension, so it can still bite the store path.

Edge tiles. n is rarely a multiple of MR and NR. A scalar fallback for edges is correct and simple; it is also a throughput cliff for small n. Libraries handle this by packing edge panels with explicit zero padding so the full-width micro-kernel always runs — this implementation zero-pads inside pack_A/pack_B but still branches on the store, which is why n = 4096 (a clean multiple) reads slightly better than nearby sizes.

Trusting the timer. Boost-clock ramp produced a 4.3× measurement error in an early sweep for these notes. Always run a discard pass, always take best-of-k.

Alternatives and When to Choose Them

Approach	Effort	Measured here	When it is the right call
Fix the loop order	one line	0.41 → 9.02 GFLOP/s (22×)	Always, first, before anything else.
Pad the leading dimension	a few lines	up to 12.1× on ijk	Whenever n is a power of two and you are not packing.
Transpose B once	n² copy	~4× (Drepper §6.2.1: 100% → 23.4%)	B reused across many multiplies, e.g. a weight matrix over an epoch.
Cache tiling alone	an hour	−26%	Only when the kernel is genuinely memory-bound — i.e. a bad loop order you cannot change.
Register tiling	half a day	9.02 → 279.74 (31×)	The highest-value rung by a wide margin. Do this before packing.
Packing + five loops	a day or two	flat 251–257 across 16× size range	When you need size-independent throughput, i.e. a library.
Link a tuned BLAS	one line	227.74 (mis-targeted)	Production. See The BLAS Interface.
Compiler generic vectors	minutes	see SIMD and Vectorization (159.78 GFLOP/s, beat intrinsics on its kernel)	When portability across ISAs matters more than the last 20%.
Parallelise over threads	a day	not measured here	After single-core is right. See Parallelizing with OpenMP, Cache-Conscious and NUMA-Aware Parallelism.

Drepper’s ladder on a 1000×1000 double-precision matmul is the same shape and worth quoting for its independence: original 100%, transposed 23.4%, sub-matrix blocked 17.3%, vectorised 9.47% of the original cycle count. Note that his blocking step bought 6.1 points and his vectorisation step bought 7.8 — the same ordering of value found here, on hardware two decades older.

The FLAME group’s how-to-optimize-gemm tutorial walks the identical sequence — “computing four elements of C at a time”, then “a 4 × 4 block of C at a time”, then “blocking to maintain performance”, then “packing into contiguous memory”. That the order is the same in a teaching tutorial, in Drepper’s 2007 paper, and in these 2026 measurements is not coincidence: it is the order in which the wins are cheapest.

Production Notes

Nobody hand-writes this in production, and that is the correct conclusion. The value of building the ladder is that afterwards you can read a performance number and know what it means. A layer running at 9 GFLOP/s is un-vectorised; at 60, it is on a mis-dispatched BLAS; at 130, it is on an AVX2 kernel when the machine has AVX-512; at 250+, it is doing what the silicon can do.

Check what kernel your BLAS chose. The 3.6× dispatch spread measured above is not exotic — DYNAMIC_ARCH builds shipped by distributions routinely predate the CPU they run on, and a CPUID the dispatcher does not recognise falls back, silently, to something ancient. openblas_get_corename() costs one line. OpenBLAS’s own param.h is where the per-architecture GEMM_DEFAULT_P/Q/R blocking constants live if you want to see how a real library parameterises the same five loops.

The parameters are per-machine, and libraries treat them that way. BLIS keeps a config/ directory with one file per microarchitecture precisely because MC, KC, NC, MR, NR do not transfer; its KernelsHowTo.md documents the micro-kernel contract an implementer must satisfy. If you ship a hand-blocked kernel, read the cache geometry at startup from sysfs and pick parameters at runtime rather than baking in constants.

Blocking interacts with threading, and the interaction is not benign. The NC loop targets L3, which on this machine is shared by eight cores (shared_cpu_list reads 0-7,16-23). Eight threads each packing a KC × NC block will evict each other. Parallelising the ic (L2) loop keeps each thread’s Ac in its private 1 MiB L2; parallelising the jc loop shares Bc in L3. Which is right depends on the shapes. Cache-Conscious and NUMA-Aware Parallelism and Memory Alignment and False Sharing own that ground, and the C micro-tiles must be laid out so two threads never write the same cache line.

Layer shape decides whether any of this is worth doing. These measurements are square GEMMs. A neural-network layer is usually far from square, and the vault’s earlier measurements showed batching worth 2.2× per example on a 784×128 layer (weights 401 KiB, L2-resident) but 32× on a 2048×2048 layer (weights 16 MiB, DRAM-resident). Blocking a layer whose weights already sit in L2 buys close to nothing. Measure your layer, not a benchmark’s.

Convolutions behave differently. For a convolution lowered through im2col the column matrix is already tall at batch size 1, so it arrives with the reuse a dense layer only gets from batching — and larger batches actively hurt (a CIFAR conv2 measured 257 → 154 GFLOP/s from N = 1 to N = 256). im2col and Convolution as GEMM owns that measurement; do not carry the dense-layer intuition across.

Keep the naive kernel forever. Every rung here was validated against mm_ikj. A blocked kernel that disagrees with the reference by more than single-precision reassociation error — 4.25 × 10⁻⁷ relative, in this implementation — is buggy, and the failures are usually in the edge tiles, which no square benchmark exercises.

See Also

• Naive Matmul and Why It Is Slow — the baseline this note is measured against; the six loop orders, the cache-miss counts, and the proof that the naive loop is instruction-bound
• The BLAS Interface — the library that already contains all of this, and how to call it
• SIMD and Vectorization — the micro-kernel’s other half; measured generic vector extensions beating hand-written intrinsics
• Layers as Matrix Multiplication — why a layer is a GEMM; available versus realised arithmetic intensity
• im2col and Convolution as GEMM — the other route into GEMM, with inverted batching behaviour
• Row-Major versus Column-Major Layout · Tensors as Flat Memory Plus Strides — the layout decisions packing works around
• Parallelizing with OpenMP · Cache-Conscious and NUMA-Aware Parallelism · Memory Alignment and False Sharing — what changes when the blocked kernel is threaded
• Roofline Analysis for LLM Inference — the bound model applied to transformer inference
• perf stat and perf top · The Performance Monitoring Unit · Cycles Per Instruction — the measurement apparatus this note had to work without
• Compiler Intrinsics · Compiler Optimization Diagnostics — writing the micro-kernel, and finding out why it did not vectorise
• Amdahl’s Law and Gustafson’s Law — the ceiling on what threading adds after this
• Building a Neural Network in C MOC — the parent map; this note is the second rung of §P7