Interactive tutorial · ACM Transactions on Mathematical Software 34(3), 2008
Anatomy of High-Performance Matrix Multiplication
Kazushige Goto and Robert A. van de Geijn, The University of Texas at Austin
This is the paper behind the GotoBLAS matrix multiply. It shows how to get within a few percent of a CPU's peak speed by deciding which memory level holds each block of the matrices, then writing one small kernel per machine.
Move with → / ←, Space, the buttons at the bottom, or just scroll. Sliders and ▶ buttons on the slides are live.
Where each block lives in the final algorithm: Ã in L2, Bj in L1, the C micro-tile in registers.
Percentages are medians for matrix sizes 1800–2000, read off the paper's Fig. 1.
The answer first
Copy the operands into packed buffers, then keep each one at the memory level its reuse needs
à sits in L2, not L1. The L2 cache can feed the floating-point units fast enough, as long as each element of à is reused across \(n_r\) columns.
The TLB, not the cache, usually caps Ã's size. Copying A into a contiguous buffer à lets it fit in the fewest TLB entries.
Of six ways to layer GEMM, one wins. With column-major storage the block-times-panel kernel, gebp, should beat the other five.
The whole algorithm, five loops deep
for p in panels of k, width kc: // gemm → gepp pack B(p,:) → B̃ // kc × n, in RAM for i in blocks of m, height mc: // gepp → gebp pack A(i,p) → Ã // mc × kc, stays in L2 for j in slivers of n, width nr: // B̃j: kc × nr, in L1 for each mr rows of Ã: c = 0 // mr × nr, in registers for q in 0 … kc−1: c += Ã[:,q] · B̃j[q,:] C(rows, cols) += c
The next slides build it up one loop at a time, and slide 14 animates it.
§1 · Fig. 1 Why it matters
GotoBLAS matched or beat each vendor's own library on three very different CPUs
\(C := AB + C\) does \(2mnk\) flops on only \(mk + kn + mn\) numbers. Each number can be reused up to n times. Whether it actually is depends on how the code moves data, not on the arithmetic.
It matters beyond GEMM: the level-3 BLAS and LAPACK do most of their work inside matrix multiplies.
Every library timed used assembly-coded inner kernels, ATLAS included. The paper notes the vendors had adopted similar techniques.
dgemm speed against matrix size, the plot's m = n. Drag or tap the chart to read values. Source: the paper's Fig. 1, digitized.
§4 · Fig. 6 The machine model
Memory is a pyramid: each level down is bigger, cheaper and slower
Tap a level. The numbers are for the 64-bit Pentium 4 Prescott in the paper's Fig. 12, the machine of its case study.
The simple model (one cache) is enough for the first cost analysis. The refined model adds the levels that decide the real design: two caches, and the slice of memory the TLB can address without a miss.
§2 · Figs. 2–3 Vocabulary
Decide which of m, n and k are small, and GEMM splits into eight named shapes
Names read ge·x·y: x is the shape of A and y the shape of B. m = matrix (both sides large), p = panel (one side small), b = block (both small). gepdot breaks the pattern: it is a generalized dot product.
C is m × n, A is m × k, B is k × n. LAPACK casts most of its work as gepp, gemp and gepm.
§3 · Fig. 4 A layered approach
Cut GEMM into gepp, gemp or gepm, cut those again, and six kernels remain
If the three bottom kernels run near peak, every shape of GEMM does too. So all the hard engineering goes into one small routine.
Tap a node. Highlighted: the path that always takes the top branch, which the paper argues is best (Fig. 8).
§3 · Fig. 5 From tree to code
Always taking the top branch gives three nested loops over blocks
for p = 1 : K // gemm_var1: k in steps of kc → gepp for i = 1 : M // gepp_var1: m in steps of mc → gebp for j = 1 : N // gebp: n in steps of nr Cij += Aip Bpj
Here \(C_{ij} \in \mathbb{R}^{m_c \times n_r}\), \(A_{ip} \in \mathbb{R}^{m_c \times k_c}\) and \(B_{pj} \in \mathbb{R}^{k_c \times n_r}\). The rest of the paper is about choosing \(m_c\), \(k_c\) and \(n_r\).
Gunnels et al. (2001) argue that this layering amortizes data movement optimally. This paper needs less: only that it amortizes well.
Tap a loop line to see what it walks over.
§4.1 · Fig. 7 · eqs. (1)–(2) One cache, three assumptions
With one cache, gebp loads A once and B and C once per sliver, so a large, square-ish A pays off
Assume (a) A, one \(n_r\)-column sliver of B and the matching sliver of C fit in the cache together, (b) with all three in cache the CPU runs at peak, and (c) A stays cached until no longer needed. Then gebp costs
\(m_c k_c\): load A once. \(k_c n\): every column of B, once. \(2 m_c n\): load and store every column of C.
§4.2.1 · eq. (3) Choosing the cache layer
à can live in the big L2 cache, as long as the \(2n_r\) flops each element feeds outlast its load
L1 meets assumptions (a)–(c) easily but holds only 16 KB on a Pentium 4: a 45 × 45 block of doubles. L2 holds 2 MB. So keep \(B_j\) and \(C_j\) in L1 and stream à from L2.
Each element of à multiplies a row of \(n_r\) elements of \(B_j\): \(2n_r\) flops, taking \(2n_r/R_\text{comp}\) cycles. Loading it from L2 takes \(1/R_\text{load}\) cycles. The load hides behind the arithmetic when
\(R_\text{comp}\): flops per cycle. \(R_\text{load}\): doubles per cycle from L2 to registers, measured ("sustained bandwidth" in Fig. 12).
§4.2.2 TLB considerations
Prefetching can hide a cache miss but not a TLB miss: a TLB miss stalls the CPU
Memory is mapped in pages, 4 KB on most machines in Fig. 12. Every access needs its page's physical address. The page table holding them lives in memory, so the CPU caches recent translations in the TLB: 64 entries on a Pentium 4, covering 256 KB.
So the paper adds two assumptions: (d) Ã, \(B_j\) and \(C_j\) are all addressable by the TLB at once, and (e) Ã stays addressable until it is no longer needed.
The demo sweeps a block over and over, the way gebp re-reads à for every sliver of B. With least-recently-used replacement, one page too many makes every access miss.
§4.2.3 Packing
Copied into one contiguous buffer, the same block of A needs far fewer TLB entries
A is a submatrix of a big matrix, so each of its \(k_c\) columns lies in a different stretch of memory, often on its own page. Packing it into a work array à fills consecutive pages instead.
The TLB budget then reads
Why the 2: while moving from \(B_j, C_j\) to \(B_{j+1}, C_{j+1}\), the old entries must coexist with the new ones. Least-recently-used replacement would otherwise evict part of Ã.
When L2 is smaller than the TLB's reach, the same argument applies to the cache. The paper says the TLB is the more common limit: on a Pentium 4 it reaches about 256 KB, against a 2 MB L2.
§4.2.3 · §5.1 What the copies cost
Each element copied into à feeds 2n flops and each element of B feeds 2m, so packing is cheap until m or n is small
Copies \(m_c k_c\) elements, then reused for \(2 m_c k_c n\) flops. The copy also leaves à in L2 with its TLB entries loaded, so it replaces a load that had to happen anyway. The paper finds it the cheaper copy in practice.
Copies \(k_c n\) elements, reused for \(2 m n k_c\) flops across all the gebp calls of one gepp. A plain memory-to-memory copy, it also disturbs the TLB.
When m is small, packing B dominates
Share of run time spent packing on a Pentium 4 Prescott, n = k = 2000, m varied. Source: the paper's Fig. 15, digitized.
§5.1 · Fig. 8 gepp via gebp
The chosen algorithm packs B once per panel and A once per block, and streams C
// C := gepp_blk_var1(A, B, C)Pack B into B̃for i = 0 … M−1 Či += Ǎi B̃ // gebp_opt1 // C := gebp_opt1(A, B, C)Pack A into Ãfor j = 0 … N−1 Caux := Ã Bj Unpack Cj := Cj + Caux
Appropriate when m and n are both large and k is not too small. Once à and B are packed, the gebp loop runs at almost the peak of the floating-point unit.
Where the TLB entries go
Packed, \(B_j\) needs one entry and \(B_{j+1}\) another. Writing \(\tilde A B_j\) to a small contiguous \(C_\text{aux}\) costs one more, and unpacking into \(C_j\) up to \(n_r\). So
Figs. 5, 8 and §6.1, animated
Watch the loops: which block sits in which memory level at every step
for p: pack B(p,:) → B̃ for i: pack A(i,p) → Ã for j: B̃j into L1 for r: c = 0 (registers) for q: c += Ã[r,q] B̃j[q,:] C(r, j) += c
§6.1 Register blocking
Inside the kernel, an mr × nr block of C stays in registers for all kc rank-1 updates
\(C_\text{aux} := \tilde A B_j\) is computed one \(m_r \times n_r\) block at a time. Each step loads one column of the à strip and one row of the \(B_j\) sliver, then does \(m_r n_r\) multiply-adds.
The block is written out once, after \(k_c\) steps: \(2 m_r n_r k_c\) flops for \(m_r n_r\) stores. So C need not stay in L1, or even L2. That is why \(k_c\) is chosen large.
§4.2.4 · §6.1 Accessing data contiguously
à is stored as mr-row strips, so the kernel reads it with stride 1
In the paper's implementation each \(m_r \times k_c\) strip of à is contiguous, and stored column by column inside. The kernel's reads then march through memory one address at a time.
Read in place, the same block jumps by the leading dimension at every step and touches many more cache lines and pages. The paper notes other implementations store à as the transpose of A, which needs a slightly more complex access pattern.
One cell = one double. Rows of 16 cells = 128 bytes. Toy sizes: mc = 8, kc = 6, mr = 4, leading dimension 22, cache line 8 doubles, page 32 doubles (pages alternate in shade).
§6.2 · Fig. 12 Choosing mr × nr
All 18 processors in the table give exactly half their registers to the C block
- Half the registers hold C. The rest prefetch elements of à and B.
- \(m_r \approx n_r\) amortizes register loads best: \(m_r + n_r\) loads buy \(2 m_r n_r\) flops.
- \(n_r \ge R_\text{comp} / (2 R_\text{load})\), so loads of à from L2 hide behind arithmetic.
The chart takes \(R_\text{load}\) from whichever cache holds Ã: L1 on the Alpha EV4, L3 on Itanium2 and PPC440, L2 elsewhere. By the table's own numbers, the SPARC IV would need \(n_r \ge 17\). The paper says nothing more about that machine.
Required nr (bar) vs chosen nr (dot)
Bars: \(R_\text{comp}/(2R_\text{load})\) from Fig. 12. Dots: the nr GotoBLAS used. Labels: mr × nr and registers in doubles.
§6.3–6.4 · Fig. 12 Choosing kc and mc
kc is set by L1 and the page size; mc by the L2 cache and the TLB's reach
\(B_j\) is reused for every strip of Ã, so it must stay in L1: \(k_c n_r\) doubles under half of L1. In the authors' experience the best \(k_c\) fills half a page: 256 doubles with 4 KB pages.
the smaller of (1) the memory the TLB can address and (2) the L2 cache. Set-associativity and replacement policy keep it to about half.
§7.2 · Fig. 12 18 processors
Pick a processor: its hardware, the block sizes chosen for it, and the paper's rules checked
§5.6 Discussion
Three rounds of argument eliminate five of the six layerings
The paper stresses that the ranking assumes blocking for L2 pays off, which held on essentially all processors of the day. It could change if blocking for L1 were preferred.
§7.3 · Fig. 13 Pentium 4 Prescott, 3.6 GHz, 64-bit
On the Pentium 4 the best à was about 1 MB: half the L2, four times the TLB's reach
- \(n_r \ge 2/(2 \times 1.03) \approx 0.97\). 16 registers hold 2 doubles each, and 8 of them hold C: \(m_r \times n_r = 4 \times 4\).
- \(k_c = 192\), not 256: the inner loop is fully unrolled, and the unrolled code has to fit in the instruction cache.
- The one exception to the TLB rule. TLB misses cost less on this chip, and filling half the L2 ran slightly faster than filling half the TLB's reach.
- Prescott's SSE3 needs no extra copying. On the earlier Northwood, SSE2 needed each element of B duplicated, which was folded into packing B.
Speed vs size of Ã, for many (mc, kc) choices
Each line is one series of runs; mc and kc were both varied from 8 to 2000 in steps of 8. The caption says 696 × 196; the table and the plot's own label say 696 × 192. Source: Fig. 13, digitized.
§7.3 · Figs. 14–15 Pentium 4 Prescott
The gap between the kernel and full dgemm is almost exactly the time spent packing
Near size 2000 the kernel runs at 6.9 GFLOPS and dgemm at 6.7: a 2.8% gap. Packing A and B takes 2.6% of the time.
Small k is the most common case in practice: rank-k updates, i.e. gepp. With m = n = 2000, dgemm is within 5% of full speed from k = 48 up.
When m is small (gepm), B is not amortized. The paper suggests skipping that copy, or the Fig. 9 algorithm. When n is small (gemp), Ã is not amortized: skip packing A, or use Fig. 11.
Top: GFLOPS. Bottom: % of run time in the packing routines. Source: Figs. 14–15, digitized.
§7.4 · Figs. 16–20 Other architectures
Retuned for each machine, the same approach ran at 89–97% of peak on four more designs
On this device A benchmark you can run
Run it: the same ideas in plain JavaScript, timed on whatever you're reading this on
Four versions of \(C \mathrel{+}= AB\) on doubles. Register blocking with a 4 × 4 micro-kernel does most of the work; packing adds stride-1 reads and fewer pages, and helps most on machines with small TLBs.
- Dot products: i-j-p loops, B read down its columns.
- Row updates: i-p-j loops, every read stride 1.
- Blocked 4 × 4: 16 accumulators in registers, no copies.
- Goto: pack B (kc × n) and à (mc × kc), then the same 4 × 4 kernel.
GFLOPS on this device
Press Run. Each version runs after a warm-up, and the faster of two runs is shown. Nothing is sent anywhere.
The 4 × 4 micro-kernel
for (let p = 0; p < kc; p++) {
const a0 = At[ai], a1 = At[ai+1], a2 = At[ai+2], a3 = At[ai+3];
const b0 = Bt[bi], b1 = Bt[bi+1], b2 = Bt[bi+2], b3 = Bt[bi+3];
c00 += a0*b0; c01 += a0*b1; c02 += a0*b2; c03 += a0*b3;
c10 += a1*b0; c11 += a1*b1; c12 += a1*b2; c13 += a1*b3;
c20 += a2*b0; c21 += a2*b1; c22 += a2*b2; c23 += a2*b3;
c30 += a3*b0; c31 += a3*b1; c32 += a3*b2; c33 += a3*b3;
ai += 4; bi += 4; // packed: both reads are stride 1
}
§8 Conclusion
Most of LAPACK and the level-3 BLAS reduce to gepp, so this one kernel carries them
Almost all LAPACK routines do the bulk of their work in gepp, gemp or gepm, and the level-3 BLAS can be cast in the same three shapes.
For almost every such routine there is an algorithmic variant that casts most of the work as gepp, the case this kernel serves best.
QR factorization is the exception: about half its work needs gemp or gepm, and its panels must stay narrow. The paper calls for more work on those shapes.
OpenBLAS grew from the GotoBLAS2 code. BLIS (Van Zee and van de Geijn, 2015) rebuilt the same loops around a small micro-kernel.
Five things to remember
- Count memops per flop at every level. GEMM wins by reuse, not arithmetic.
- Block for L2 when \(n_r \ge R_\text{comp}/(2R_\text{load})\) holds. L1 alone is too small for a large Ã.
- Size à to the TLB's reach. A TLB miss stalls; a cache miss can be prefetched.
- Pack, so every operand is contiguous, on few pages, and read with stride 1. The copy is amortized over O(m) or O(n) flops.
- Keep an \(m_r \times n_r\) block of C in half the registers, with \(m_r \approx n_r\) and \(k_c\) large.
Sources
Read the original
- Kazushige Goto and Robert A. van de Geijn. Anatomy of High-Performance Matrix Multiplication. ACM Transactions on Mathematical Software 34(3), article 12, May 2008. doi:10.1145/1356052.1356053 · PDF at UT Austin
- K. Goto and R. A. van de Geijn. On reducing TLB misses in matrix multiplication. Tech. Rep. CS-TR-02-55, UT Austin, 2002.
- K. Goto and R. A. van de Geijn. High-performance implementation of the level-3 BLAS. FLAME Working Note #20, TR-2006-23, UT Austin, 2006.
- J. A. Gunnels, G. M. Henry and R. A. van de Geijn. A family of high-performance matrix multiplication algorithms. ICCS 2001, LNCS 2073.
- R. C. Whaley and J. J. Dongarra. Automatically tuned linear algebra software. SC'98.
- F. G. Van Zee and R. A. van de Geijn. BLIS: A framework for rapidly instantiating BLAS functionality. ACM TOMS 41(3), 2015.
About this deck: every chart is redrawn from the paper's own vector figures, digitized from the PDF, and the processor data is Fig. 12 as printed. Derived numbers (percent of peak, required nr, footprints) are computed from those. The browser benchmark is new and measures your device, not the paper's.