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. 93% of peak, 3.6 GHz Pentium 4 (6.7 of 7.2 GFLOPS) 95% of peak, 1.9 GHz POWER5; IBM's ESSL: 84% 98% of peak, 1.5 GHz Itanium2; Intel's MKL: 97% 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 1 · §4.2.1 à 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. 2 · §4.2.2–4.2.3 The TLB, not the cache, usually caps Ã's size. Copying A into a contiguous buffer à lets it fit in the fewest TLB entries. 3 · §5.6 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 This shape is 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 \[ \begin{aligned} \frac{\text{flops}}{\text{memops}} &= \frac{2 m_c k_c n}{m_c k_c + (2 m_c + k_c)\, n} \\ &\approx \frac{2 m_c k_c}{2 m_c + k_c} \quad \text{when } k_c \ll n \end{aligned} \] \(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. — flops per memop — time added by memory traffic The paper's examplemc = kc = n/100, a memop costs 10 flops §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 \[ n_r \;\ge\; \frac{R_\text{comp}}{2\,R_\text{load}} \] \(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. ▶ SweepReset §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 \[ T_{\tilde A} + 2\,(T_{B_j} + T_{C_j}) \le T \] 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 Pack A → à 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. Pack B → B 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 \[ T_{\tilde A} \le T - (n_r + 3). \] Figs. 5, 8 and §6.1, animated Watch the loops: which block sits in which memory level at every step à · L2 Bj · L1 C tile · registers B panel · RAM ▶ Play Step Next j Next i Reset Now 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. ▶ PlayStep §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. ▶ Play readsReset 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 kc · as large as possible, but \(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. mc · à fills about half of 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 Processor §5.6 Discussion Three rounds of argument eliminate five of the six layerings Next round ▸Reset Why 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
