经典论文精读——《Anatomy of High-Performance matrix multiplication》

注：
（$1$）加粗部分是个人认为比较重要的地方；
（$2$）加粗且斜体部分是还没理解的部分；
（$3$）部分术语保留了原文没有翻译，因为我觉得用原文的词更便于理解，
（$4$）有的术语翻译在后面的括号里注明了原文对应的词，便于理解；
（$5$）加中划线的词是感觉翻译得不好，没能理解原文的意思。

正文开始

​我们介绍了高性能矩阵乘的基本原理，它是广泛使用的 $GotoBLAS$ 库的一部分，我们的设计和决策是通过不断细化具有多层次内存的架构模型来证明的，成果是一个 用于执行这个操作的 简单但有效的算法，在更多架构上的实现显示出可以实现接近峰值的性能。

Section1 Introduction

​实现矩阵乘法以实现接近最优性能，需要结合微观层面的高性能 $kernels$ 的周密的 $engineering(n.)$，深入理解如何在宏观层面上分层操作（$layered$），本文主要讨论宏观问题，即如何利用一个高性能 $inner-kernel$，而不是与该 $inner-kernel$ 的 $design$ 和 $engineering$ 相关的微观问题。

​在 $Gunneletal.[2001]$ 中提到了一种矩阵乘实现的分层（$layered$）方法，该方法可以最优地分摊（$amortize$） 在具有复杂的多级内存的架构上的 两个相邻内存层之间所需的数据移动（开销）。和该地区的其他工作一样 $[Agarwaletal.1994;Whaleyetal.2001]$，$Gunnel$ 等人用一个 “$inner-kernel$” 来转换（$cast$）计算，该 “$inner-kernel$” 用于计算 $C:=\tilde{A}B+C$，其中 $\tilde{A}$ 是一个规模为 $m_c\ast k_c$ 的、以某种打包（$packed$）的格式连续存储的、适应（$fitin$）于 $cachememory$ 的矩阵。不幸的是，他们使用的内存层次结构模型至少在以下两方面是不现实的：

$——$其假设该 $inner-kernel$ 用于计算的 $\tilde{A}$ 驻留在 $L_{1}cache$ 中。

$——$其忽略了与 $TranslationLook-asideBuffer(TLB)$ 相关的问题。

本文扩展了一个相关的技术报告$ [Goto_andvan_deGeijn~2002]$，该报告的观察结果如下：

$——$浮点单元执行浮点操作的速率（$flops$），浮点数从 $L_{2}cache$ 流（$stream$）到寄存器的速率，这两者的比率通常相对较小，这意味着 $\tilde{A}$ 可以从 $L_{2}cache$ 中流出。

$——$通常，$TLB$ 能够寻址（$address$）的数据量大小（$amount$）是 $\tilde{A}$ 大小（$size$）的限制因素。

此外，本文观察到以下：

$——$事实上有 $6$ 种 $inner-kernel$ 应被考虑用来构建（$building$） 用于高性能矩阵乘的 块（$blocks$），且已经证明，其中一种天然地（$inherently$）比其他几个更优。（在 $Gunneletal.[2001,2005]$中，其中三个已被鉴别（$identified$））

对这些观察结果的仔细考虑是 $DGEMMBasicLinearAlgebraSubprograms(BLAS)$ 例程（$routine$）实现的基础，该例程是广泛使用的 $GotoBLAS$ 库的一部分。

​在图 $1$ 中，我们预览了这些技术的效能。图中展示了我们的实现和供应商的实现（$Intel’sMKL(\mathrm{8.1.1})andIBM’sESSL(\mathrm{4.2.0})libraries$）分别在在因特尔 $Pentium4Prescott$ 处理器、$IBMPower5$ 处理器以及因特尔 $Itanium2$ 处理器上的性能。需要注意的是供应商实现采用了和本文所描述十分类似的技术。不能孤立地判断矩阵乘的性能，它通常是其他操作的基石，如$level-3BLAS(matrix-matrixoperations)[Dongarraetal.1990;Kagstrometal.1998]$ 和 $LAPACK[Andersonetal.1999]$。本文所描述的技术是如何影响 $level-3BLAS$ 的，在 $GotoandvandeGeijn[2006]$中讨论过。

​本文尝试在高层次描述这些问题以便有更广泛的受众，低层次的问题仅在需要时介绍。在 $Section2$ 我们介绍在本文剩余部分所使用的记号，在 $Section2$ 介绍一种用于实现矩阵乘的分层方法（$layeredapproach$），$Section4$ 讨论 $inner-kernels$ 的高性能实现，$Section5$ 给出 用于矩阵乘最常见情形的 实际算法，在 $Section6$ 给出 在实践中用于 确定 为了优化性能而必须调整的 参数的更多细节，在 $Section7$ 中给出在各种架构上使用高度调优的实现所获得的性能结果，$Section8$ 进行总结。

Section2 Notation

矩阵划分是描述矩阵乘算法的基础，给出一个 $m\times n$ 的矩阵 $X$，我们仅考虑将 $X$ 进行列分块和行分块：

若进行列分块，共分为 $N$ 个列块，表示为 $X_j(j=0,1,...,N-1)$，每一个列块又有 $n_b$ 列（除了 $X_{N-1}$，它可能不足 $n_b$ 列）；若进行行分块，共分为 $M$ 个行块，表示为 $\check{X_i}(i=0,1,...,M-1)$，每一个行块又有 $m_b$ 行（除了 ${\check{X}}_{M-1}$，它可能不足 $m_b$ 列）。

​矩阵乘法的实现将由子矩阵的乘法组成，我们对这些计算进行了命名，如图 $2$, $3$ 所示。我们注意到这些特殊的形状经常作为 用于其他线性代数操作的算法 的一部分，例如 $LAPACK$ 支持的各种操作的计算大部分都可由 $GEPP$，$GEMP$ 和 $GEPM$ 来转换（$cast$），甚至给定单个密集线性代数操作，通常存在多种算法，每一个算法都能由这些 $GEMM$ 乘法的不同情况来实现$[Bientinesietal]$。

Section3 A layered approach for GEMM

在图 $4$ 中，我们展示了 $GEMM$ 能够怎样分解为图 $2$ 中所描述的特殊情形，通用的 $GEMM$ 能够被分解为多个 $GEPP$，$GEMP$ 或 $GEPM$ 的调用， $GEPP$，$GEMP$ 和 $GEPM$ 又可被分解为多个 $GEBP$，$GEPB$ 或 $GEPDOT$ $kernels$ 的调用。现在的想法是，如果这三个最低层的 $kernels$ 能获得高性能，则 $GEMM$ 的其它情形也能获得高性能。

在图 $5$，我们将图 $4$ 中的顶部分支与一个三层嵌套循环联系起来

其中矩阵 $A$，$B$ 和 $C$ 按如下方法分为若干子矩阵：

$$C=\left\{\begin{array}{cccc}{C}_{11}& C_{12}& \cdots & C_{1N}\\ C_{21}& C_{22}& \cdots & C_{2N}\\ ⋮& ⋮& \ddots & ⋮\\ C_{M1}& C_{M2}& \cdots & C_{MN}\end{array}\right\},A=\left\{\begin{array}{cccc}{A}_{11}& A_{12}& \cdots & A_{1K}\\ A_{21}& A_{22}& \cdots & A_{2K}\\ ⋮& ⋮& \ddots & ⋮\\ A_{M1}& A_{M2}& \cdots & A_{MK}\end{array}\right\},B=\left\{\begin{array}{cccc}{B}_{11}& B_{12}& \cdots & B_{1N}\\ B_{21}& B_{22}& \cdots & B_{2N}\\ ⋮& ⋮& \ddots & ⋮\\ B_{K1}& B_{K2}& \cdots & B_{KN}\end{array}\right\}$$
其中 $C_{ij}\in {\mathbb{R}}^{m_c\times n_r}$，$A_{ip}\in {\mathbb{R}}^{m_c\times k_c}$，$B_{pj}\in {\mathbb{R}}^{k_c\times n_r}$，块大小（$blocksizes$）$m_c$，$k_c$ 和 $n_r$ 的具体取值将在后文具体讨论。

​在 $Gunnelsetal.[2001]$ 中可以找到一个支持 本节中提到的对于一般方法的最优性主张 的理论，特别地，该论文支持这样的观察结果：如果内存层（$memorylayers$）之间的数据移动能被最优地摊销（$amortized$），那么计算就应按照图 $4$ 给出的决策树来转换（$cast$）。然而，本文是自成体系的（$self-contained$），因为我们表明该方法能很好地分摊（$amortized$）这种开销（$overhead$），因此最优性对我们的讨论并不重要。

Section4 high-performance GEBP，GEPB，GEPDOT （inner-kernels的高性能实现）

我们现在讨论用于 $GEBP$，$GEPB$ 和 $GEPDOT$ 高性能实现的技术，我们通过首先分析 使用内存层次结构（$memoryhierarchy$）的一个公认的朴素（$naive$）模型在内存层之间移动数据的开销（$cost$） 来做这件事，在 $Section4.2$ 我们为该模型添加了更多的实际细节，这为 $Section5$ 中的 $GEBP$，$GEPB$ 和 $GEPDOT$ 的算法奠定了基础。

4.1 Basics

在图 $6$（左）我们描述了一个非常简单的多级内存模型，在 $RandomAccessMemory(RAM)$ 和寄存器组（$registers$）之间插入了一个缓存层（$cachememory$），与 $GEBP$，$GEPB$ 和 $GEPDOT$ 高性能实现相关的顶级问题（$top-levelissues$）可以用这个简化的体系结构来描述。

​我们首先关注 $GEBP$，其中 $A\in {\mathbb{R}}^{m_c\times k_c}$，$B\in {\mathbb{R}}^{k_c\times n}$，$C\in {\mathbb{R}}^{m_c\times n}$，且 $B，C$ 进行列分块
$$B=(B_0|B_1|\cdots |B_{N-1}),C=(C_0|C_1|\cdots |C_{N-1})$$
并做以下假设：

​假设（$a$）：$A$ 的尺寸 $m_c$ 和 $k_c$ 足够小，使得 $A$、$B$ 的一个列分块里的 $n_r$ 列（即 $B_j$）和 $C$ 的一个列分块里的 $n_r$ 列（即 $C_j$）三者能一起放进缓存中（$fitsinthecache$）。

​假设（$b$）：如果 $A$，$C_j$ 和 $B_j$ 都在缓存中，那 $C_j:=A{B}_j+C_j$ 就能以 CPU 峰值速率计算。

​假设（$c$）：如果 $A$ 在缓存中，则其一直在缓存中，直到不需要为止。

​在这些假设下，图 $7$ 中的 $GEBP$ 方法分摊了（$amortizes$）了在主存（$mainmemory$）和缓存之间移动数据的开销：

具体计算如下：

加载（$load$） $A$ 只需要一次，其数据移动开销为 $m_c{k}_c$；但 $B_j$ 和 $C_j$ 各有 $N-1$ 个，故各需要加载 $N-1$ 次，计算完还要把 $C_j$ 存回（$store$）主存，$C_j$ 又需要额外的 $N-1$ 次，因此每一次（更新 $C_j$）的数据移动开销为 $k_c{n}_r+m_c{n}_r+m_c{n}_r=(2m_c+k_c)n_r$，共 $N-1$ 次，再对 $n_r$ 累加求和，得 $(2m_c+k_c)n$。因此总的数据移动开销为
$$m_c{k}_c+(2m_c+k_c)n(memops)$$
计算 $A{B}_j$ 需要 $m_c{k}_c{n}_r$ 次乘法，$m_c(k_c-1)n_r$ 次加法，计算 $C_j:=A{B}_j+C_j$ 需要 $m_c{n}_r$ 次加法，故更新一次 $C_j$ 需要的总计算次数为 $2m_c{k}_c{n}_r$，共更新 $N-1$ 次，再对 $n_r$ 累加求和，得总的计算开销为
$$2m_c{k}_{c}n(flops)$$
总计算开销/总数据移动开销的比率为
$$\frac{2m_c{k}_{c}n}{m_c{k}_c+(2m_c+k_c)n}\frac{flops}{memops}\approx \frac{2m_c{k}_{c}n}{(2m_c+k_c)n}\frac{flops}{memops}when{k}_c\ll n(1)$$
即
$$\frac{2m_c{k}_c}{2m_c+k_c}(2)$$
应在 “$m_c{k}_c$ 个浮点数（即矩阵 $A$）填充了缓存的大部分”这个约束条件和上述的假设（$a$）-（$c$）下 ，最大化（2）式比值。在 $Section6.3$ 我们将看到，实践中还有其他影响 $k_c$ 选择的因素。然而，底线（$bottomline$）是这样的：在简化的假设下，$A$ 应该占据尽可能多的缓存空间，而且应该是方阵（即 $m_c=k_c$），同时至少需要给 $B_j$ 和 $C_j$ 留出缓存空间。如果 $m_c=k_c\approx n/100$，那么即使数据移动操作（$memops$）比计算操作（$flops$）慢10倍，数据移动操作也仅给计算操作（$computation$）增加约10%的开销。

​$GEPB$ 和 $GEPDOT$ 也可进行类似的分析。

4.2 Refinements

在讨论实践中的考虑时，我们仍关注 $GEBP$ 的高性能实现，在本文的剩余部分，我们假设矩阵都是以列优先的方式存储。

4.2.1 Choosing the Cache Layer

​图 $6$（右）描述了一个更加精确的内存层次结构，该图指出通常有多级缓存。

​第一个问题：$m_c\times k_c$ 的矩阵 $A$ 应该放在哪一层缓存？（2）式表明（当然是在假设（$a$）-（$c$）下）：$m_c\times k_c$ 越大，在 $RAM$ 和缓存之间移动数据的开销就被计算（$computation$）分摊（$amortized$）得越好。这表明在（大致）满足假设（$a$）-（$c$）的前提下，矩阵 $A$ 应加载进距离寄存器组最远的缓存层，因为它能容纳最多的数据。

​$L_1$ 缓存本身就具备这样的特性：如果用它存储 $A$，$B_j$ 和 $C_j$，则满足假设（$a$）-（$c$），但 $L_1$ 缓存往往非常小。为了让 $m_c\times k_c$ 很大，能否将 $A$ 存在 $L_2$ 缓存中呢？分别记 $R_{comp}$ 和 $R_{load}$ 为CPU执行浮点操作的速率和将浮点数从 $L_2$ 缓存流（$stream$）到寄存器组的速率，同时假设 $A$ 存在 $L_2$ 缓存中，而 $B_j,C_j$ 存在 $L_1$ 缓存中，再进一步假设在 $L_1$ 缓存和寄存器组之间有足够的带宽，使得能顺利将 $B_j$ 和 $C_j$ 的元素（$elements$）从 $L_1$ 缓存加载到寄存器组。在 $Section4.1$ 中，我们计算过完成一次 $C_j:=A{B}_j+C_j$ 需要 $2m_c{k}_c{n}_r$ 次计算（$flops$），并且需要将 $A$ 的 $m_c\times k_c$ 个元素从 $L_2$ 缓存加载到寄存器组。为了将 把 $A$ 的元素加载进寄存器组的开销 和计算的开销重叠（$overlap$），也就是说前者应不大于后者，即
$$\frac{1}{R_{load}}\times m_c{k}_c⩽\frac{1}{R_{comp}}\times 2m_c{k}_c{n}_r\Rightarrow \frac{2n_r}{R_{comp}}⩾\frac{1}{R_{load}}$$
即
$$n_r\ge \frac{R_{comp}}{2R_{load}}(3)$$

4.2.2 TLB Considerations

​第二个体系结构方面的考虑和页面管理系统相关。在我们的讨论中，考虑一个典型的现代架构（$architecture$）就足够了，其使用虚拟内存，这样可用内存的大小（$size$）就不受物理内存大小的限制了：内存被划分为一些规定大小（通常是固定的）的页面（$pages$）。页表（$pagetable$）将虚拟地址映射到物理地址并跟踪（$track$）一个页面是在内存里还是在磁盘上，问题在于页表本身也是放在内存里的，当进行虚拟地址向物理地址转换时，会产生额外的访存开销。为了克服这个问题，引入了一个更小的表 $TranslationLook-asideBuffer(TLB)$ 用于存储最近使用（$mostrecentlyused$）的页面的信息，若一个虚拟地址能在 $TLB$ 中找到，那么其转换就很快，若不能找到（发生 $TLBmiss$），就会去查找页表，并将查找的条目（$entries$）从页表移到 $TLB$，换句话说，$TLB$ 就是页表的一个缓存。最近有些架构已经引入了二级 $TLB$，其动机和引入 $L_2$ 缓存类似。

​缓存不命中（$cachemiss$）和 $TLBmiss$ 之间最显著的差异在于前者不一定会暂停（$stall$）CPU运行。通过使用算法预取技术（$algorithmicprefetchingtechniques$），可以容忍少量的缓存不命中，只要数据能够从它所在的内存中读取得足够快，并在需要计算时到达CPU。相比之下，一次 $TLBmiss$ 就会造成CPU停顿，直到 $TLB$ 更新新的地址。也就是说，预取可以掩盖（$mask$）缓存不命中，但不能掩盖 $TLBmiss$。

​$TLB$ 的存在意味着还需满足以下额外的假设：

​假设（$d$）：$A$ 的尺寸 $m_c$ 和 $k_c$ 足够小，使得 $A$、$B$ 的一个列分块里的 $n_r$ 列（即 $B_j$）和 $C$ 的一个列分块里的 $n_r$ 列（即 $C_j$）三者能同时被 $TLB$ 寻址，从而在计算 $C_j:=A{B}_j+C_j$ 时没有 $TLBmiss$ 发生。

​假设（$d$）：如果 $A$ 能被 $TLB$ 寻址，则其一直保留在那，直到不需要为止。

4.2.3 Packing

​现在最重要的问题是，$A$ 通常是一个更大矩阵的子矩阵，因此在内存中不连续，这意味着对其寻址所需要的 $TLB$ 条目（$entries$）远超（实际） $TLB$ 条目的最小值，解决方案是将 $A$ 打包（$pack$）成一个连续的工作阵列（$workarray$）$\tilde{A}$，然后再选择使得 $\tilde{A},B_j$ 和 $C_j$ 三者能一起放进 $L_2$ 缓存且都能被 $TLB$ 寻址的参数 $m_c$ 和 $k_c$。

Case 1 The TLB is the limitting factor

​假设 $TLB$ 有 $T$ 个可用条目，记 $T_{\tilde{A}},T_{B_j}$ 和 $T_{C_j}$ 分别为 $\tilde{A},B_j$ 和 $C_j$ 专用的 $TLB$ 条目的数量，则有
$$T_{\tilde{A}}+2(T_{B_j}+T_{C_j})⩽T$$
​引入因子 $2$ 的原因如下：当 $B_{j+1}$ 和 $C_{j+1}$ （注：除了 $B_{N-1}$ 和 $C_{N-1}$，它们各有 $n_r$ 个列）第一次被寻址时（即访问的下一个块（$block$）属于 $B_{j+1}$ 和 $C_{j+1}$）时，用于寻址（$address$）它们的 $TLB$ 条目本应替换掉那些用于寻址 $B_j$ 和 $C_j$ 的条目，然而，完成计算 $C_j:=\tilde{A}{B}_j+C_j$ 后，一些和 $\tilde{A}$ 相关的 $TLB$ 条目会成为最近最少使用的（$LeastRecentlyUsed,LRU$），因此会被用于寻址 $B_{j+1}$ 和 $C_{j+1}$ 的条目替换掉。因子 $2$ 允许和 $B_j$ 和 $C_j$ 相关的条目能够和与 $\tilde{A}$，$B_{j+1}$ 和 $C_{j+1}$ 相关的条目共存，这样当 $B_{j+2}$ 和 $C_{j+2}$ 第一次被寻址时，和 $B_j$ 和 $C_j$ 相关的条目会成为最近最少使用的，从而被替换掉，这就避免了和 $\tilde{A}$ 相关的 $TLB$ 条目被替换。

（个人理解：整个 $Section4$ 讨论的都是 $GEBP$ 的实现，对于 $GEBP$ 操作来说，$B_j$ 和 $C_j$ 用了一次就不用了，但 $A$ 会反复用到，因此应保证用于寻址 $A$ 的 $TLB$ 条目不被替换）

​如果 $A\to \tilde{A}$ 的打包（$packing$）满足以下条件：打包完成后 $\tilde{A}$ 驻留在 $L_2$ 缓存中，且能够被 $TLB$ 寻址，以为后续计算做准备，则其无需产生比 把 $A$ 加载进 $L_2$ 缓存和 $TLB$ 的开销还大的开销（$overhead$）。访问 $A$ 来完成打包过程的开销（$cost$）无需远大于将 $A$ 移动到 $L_2$ 缓存的开销大，即使不进行 $A$ 的打包，后者也是必要的。

​$GEBP$ 操作是在 $GEPP$ 或者 $GEPM$ 的上下文（$context$）执行的，对于 $GEPP$ 来说，会多次调用 $GEBP$ 操作，且每次都会重用矩阵 $B$，因此将 $B$ 也打包（$pack$）成一个连续的工作阵列（$workarray$）$\tilde{B}$ 是值得的，从而在进行 $C:=\tilde{A}\tilde{B}+C$ 计算时能够降低 $T_{\widetilde{B_j}}$ （用于 $\widetilde{B_j}$ 的 $TLB$ 条目）的值。

Case 2 The size of the L2 cache is the limiting factor

​该情形也能进行类似的讨论，但由于当前的受限因素主要在于 $TLB$ 能够寻址的内存大小（例如，当前一代的 $Pentium4$ 上的$TLB$ 可以寻址约 $256Kbytes$，而 $L_2$ 缓存可以容纳 $2Mbytes$），我们不做具体细节的详细说明。

4.2.4 Accessing Data Contiguously

​为了最高效地将数据移动到寄存器组，组织计算使得 在内存中连续（$consecutive$）的数据能够连续的操作使用 是很重要的，实现这一点的方法不仅是将 $A$ 打包（$pack$）成一个连续的工作阵列（$workarray$）$\tilde{A}$，还需要仔细安排（$arrange$）它，详见 $Section6$。

4.2.5 Implementation of $GEPB$ and $GEPDOT$

​与 $GEBP$ 的分析类似。

Section5 Practical Algorithms

在相对粗略地分析了这些方法之后，我们现在讨论图 $4$ 中所有六个选项的实际算法，同时结合其他体系结构方面的考虑。

5.1 Implementing $GEPP$ with $GEBP$

现在总结 $Section\mathrm{4.2.1}-\mathrm{4.2.4}$ 的观察结果，即图 $8$ 中的用 $GEBP$ 来实现 $GEPP$。

打包（$packing$）和计算（$computation$）被安排（$arrange$）来最大化 $\tilde{A}$ 的尺寸（$size$）：

在算法 $GEPP\_VAR1$ 中将 $B$ 打包成 $\tilde{B}$ 使得 $B_j$ 通常只需要一个 $TLB$ 条目，引进 $B_{j+1}$ 需要第二个 $TLB$ 条目，缓冲区 $C_{aux}$ 仅需一个 $TLB$ 条目，以及 $C_j$ 所需的最多 $n_r$ 个条目（如果 $C_j$ 的主要尺寸（$leadingdimension$）很大，那就需要 $n_r$ 个），因此 $T_{\tilde{A}}$ 至少为

$T-(n_r+3)$，保守起见就取它。此外，$C_j$ 存储不连续没什么影响，因为进行 $GEPP$ 操作时不会重复使用 $C_j$。

（个人理解：要确定 $\tilde{A}$ 能始终放在 $L_2$ 中的最大尺寸，它由留给 $\tilde{A}$ 的 $TLB$ 条目来确定，因为要保证不出现 $TLBmiss$）

​一旦 $B$ 和 $A$ 分别被复制（$copy$）到 $\tilde{B}$ 和 $\tilde{A}$ 中，算法 $GEBP\_OPT1$ 中的循环几乎能以浮点单元的峰值性能执行。

$——$$B$ 的打包是一个内存到内存（$memory-to-memory$）的复制，其开销正比于 $k_c\times n$，该开销被分摊（$amortized$）到 $2m\times n\times k_c$ 的计算量（这是更新整个 $C$，也就是一次 $GEPP$ 操作所需的计算量）里，**因此对每一个复制项（$copieditem$，个人理解：一个复制项就是 $B$ 的一个元素，即一个浮点数），有 $O(m)$的计算量去摊销。**这种打包操作破坏（$disrupt$）了 $TLB$ 之前的内容。

$——$如果精心设计，从 $A$ 到 $\tilde{A}$ 的打包将 $A$ 的数据从内存里重排（$rearrange$）到一个 可能保留在 $L_2$ 缓存的 缓冲区，且能让 $TLB$ 加载有用的条目。其开销正比于 $m_c\times k_c$，该开销被分摊（$amortized$）到 $2m_c\times k_c\times n$ 的计算量（这是更新一个 $C_j$，也就是一次 $GEBP$ 操作所需的计算量）里，因此对每一个复制项（理解同上）有 $O(n)$ 的计算量去摊销。在实践中，此开销通常很小（$lessexpensive$）。

当 $m$ 和 $n$ 都很大，且 $k$ 不太小时，该方法很适合用来实现 $GEMM$

5.2 Implementing $GEPM$ with $GEBP$

在图 $9$ 提出了一个类似的方案来用 $GEBP$ 实现 $GEPM$，此时 $C$ 会被重复更新，因此在把结果加回 $C$ 前先用 $\tilde{C}=AB$ 来累加是很值得的，$\check{B_p}$ 没有重用，因此无需打包。$B_j$ 需要至多 $n_r$ 个 $TLB$ 条目，$B_{temp}$，$C_j$ 和 $C_{j+1}$ 各需要一个，因此 $T_{\tilde{A}}$ 同样至少为 $T-(n_r+3)$。

5.3 Implementing $GEPP$ with $GEPB$

图 $10$ 展示了怎样用 $GEPB$ 来实现 $GEPP$，算法 $GEPP$ 中将 $A$ 进行打包并转置以改善对其元素的连续访问，在算法 $GEPB$ 中，$B$ 也被打包了，并放在 $L_2$ 缓存，因此 $T_{\tilde{B}}$ 是我们想要最大化的值。$C_i$ 需要至多 $n_r$ 个 $TLB$ 条目，$C_{aux}$，$A_i$ 和 $A_{i+1}$ 各需要一个，因此 $T_{\tilde{B}}$ 至少为 $T-(n_r+3)$。

5.4 Implementing $GEMP$ with $GEPB$

图 $11$ 展示了怎样用 $GEPB$ 来实现 $GEMP$，引入一个临时的 $\tilde{C}$ 来累计 $\tilde{C}=(AB{)}^T$，且打包的副本 $\tilde{B}$ 填充了 $L_2$ 缓存的大部分。${\check{A}}_i$ 需要至多 $n_r$ 个 $TLB$ 条目，$C_i$，$C_{i+1}$ 和 $A_{temp}$ 各需要一个，因此 $T_{\tilde{B}}$ 至少为 $T-(n_r+3)$。

5.5 用 $GEPDOT$ 来实现 $GEPM$ 和 $GEMP$

类似地，$GEPM$ 和 $GEMP$ 也可通过 $GEPDOT$ 操作来实现。将 $C$ 的一个块（$block$）放在 $L_2$ 缓存中，并通过将 $A$ 的一些列和 $B$ 的一些行来更新它。因为我们将在下面讨论这种方法可能较差，所以我们在这里不详细讨论。

5.6 讨论

图4展示了6种不同的 $GEMM$ 实现方法，其中4种的细节分别在图8-11给出。我们现在论证以下事实：如果矩阵都按照列优先（$column-majororder$）存储，那么图8中的方法（也就是对应图4决策树的最顶部分支，同时也在图5中给出）在实践中可能获得最好的性能。

​我们首先关注的是从 $L_2$ 缓存中获得最好的带宽使用，注意基于 $GEPDOT$ 实现的方法将 $C$ 的一个块（$block$）放在 $L_2$ 缓存中，每次都要从内存流出（$streamed$）$A$ 的一些列和 $B$ 的一些行用于相乘来更新 $C$ 的块，（前提假设是 $A$ 的一些列和 $B$ 的一些行都太大了，不能放在 $L_2$ 缓存里），计算出的 $C$ 又太大，不能放在寄存器，又得写回。因此该方法要求的 $L_2$ 和寄存器组之间的带宽是基于 $GEBP$ 和 $GEPB$ 实现的方法的两倍，首先舍弃。

​比较图8和9的算法，主要区别在于：前者将 $B$ 打包，并从主存中来回流动 $C$ 的元素，后者从主存中流出 $B$，在一个中间的（$temporary$）缓冲区里计算 $C$，最后通过将此中间结果加回 $C$ 来实现解包（$unpack$）。在实践中，前者能够用计算来隐藏（$hide$）在主存来回移动 $C$ 的元素的开销，但暴露了对 $B$ 打包的开销，后者能够隐藏从内存取出 $B$ 的开销，但暴露了对 $C$ 解包的开销。$C$ 的解包是一个更复杂的操作，因此其开销比 $B$ 的打包更大，这使得图8的算法比图9的更可取。图10和11的算法也可进行类似的比较。

​现在我们要从图8和10中的算法选出最好的一个，二者表面上看似乎是对称的，因为 $A$ 和 $B$ 的作用刚好相反。注意前者每次访问 $C$ 的一些列，后者每次访问 $C$ 的一些行，如果矩阵按列优先存储，那按列访问矩阵块（$block$）是更可取的。因此图8的算法最优。

​我们在本文后续关注图8的算法。

​本小节（$Section5.6$）的结论是基于 “对 $L_2$ 缓存 $blocking$ 更具优势” 这一前提。

Section6 More Details Yet

​现在，我们就“像 $GEBPOPT1$ 这样的 $kernel$ 是如何使用寄存器组的”这一问题提出一些最终的见解，然后再讨论在实践中如何选择参数。前文已说图8算法可能获得最好的性能，因此我们关注它。

6.1 Register Blocking

​考虑图8中的 $C_{aux}:=\tilde{A}{B}_j$，$\tilde{A},B_j$ 分别在 $L_2$ 和 $L_1$ 缓存，该操作可通过在寄存器组里完成 $C_{aux}$ 的尺寸为 $m_r\times n_r$ 的子矩阵的 计算来实现。

​这意味着在 $C_j$ 的计算过程中，没必要把 $C_{aux}$ 的这些子矩阵存在 $L_1$ 甚至 $L_2$ 缓存中：每一次计算完 $C_{aux}$ 的子矩阵，需要 $m_r{n}_r$ 次$memops$ 将其从寄存器组（$store$）存到任一内存层，而计算一次子矩阵需要 $2m_r{n}_r{k}_c$ 次 $flops$ ，因此访存时间能被计算时间隐藏。后面我们将看到 $k_c$ 一般取得相对较大。

​此时我们更仔细地讨论一下将 $A$ 打包成 $\tilde{A}$ 是怎么实现的。在我们的实现中，$\tilde{A}$ 的存储应满足以下条件：它的每一个尺寸为 $m_r\times n_r$ 的子矩阵连续地存在内存里，每一个子矩阵本身又是按列优先存储的，这使得计算 $C_{aux}:=\tilde{A}{B}_j$ 时访问 $\tilde{A}$ 的元素是连续的。其他人的实现通常把 $\tilde{A}$ 存储为 $A$ 的转置，这在访问 $A$ 时需要一个稍微复杂一些的方法。

6.2 Choosing $m_r\times n_r$

以下考虑影响到 $m_r\times n_r$ 的选择：

$——$通常用一半的可用寄存器来存储 $C$ 的尺寸为 $m_r\times n_r$ 的子矩阵，另一半用于 $\tilde{A},\tilde{B}$ 元素的预取。

$——$当 $m_r\approx n_r$ 时，加载寄存器（$loadingtheregisters$）的开销的分摊方案是最优的。

$——$如 $Section\mathrm{4.2.1}$​ 提到的，从 $L_2$ 缓存中取出 $\tilde{A}$ 的一个元素到寄存器的时间应不超过前一个元素的计算时间，因此有 $n_r\ge \frac{R_{comp}}{2R_{load}}$，$R_{comp}$ 和 $R_{load}$ 对应图12中的 “$flops/cycles$” 列和 “$SustainedBandwidth$” 列（后者仅针对于 $L_2$ 缓存）。

寄存器短缺会限制算法 $GEBP\_OPT1$ 的性能，因为会削弱 隐藏（$hide$）与到 $L_2$ 缓存带宽相关的约束的 能力

6.3 Choosing $k_c$

为了分摊 更新 $C_j$ 的 $m_r\times n_r$个元素 的开销，$k_c$ 应越大越好。以下考虑影响到 $k_c$ 的选择：

$——$$B_j$的元素被重用很多次，因此应保持在 $L_1$ 缓存中。此外，缓存的相联度和替换策略进一步限制了 $B_j$ 能占用 $L_1$ 缓存多少空间。实践中，$B_j$ 的 $k_c{n}_r$ 个浮点数应只占据$L_1$ 缓存空间的一半以下，以保证 $B_j$ 元素不会被 $\tilde{A}$ 和 $C_{aux}$ 的元素驱逐。

$——$$\tilde{A}$ 空间（$m_c{k}_c$个浮点数）应占据 $L_2$ 的相当大一部分。

在我们的实验中，最优选择是使得 $k_c$ 个双精度浮点数占用半个页面，此选择通常也满足其他约束以及本文之外的其他体系结构的约束。

6.4 Choosing $m_c$

前面已经讨论过，$m_c\times k_c$ 规模的矩阵 $\tilde{A}$ 应占用以下二者中较小者 的相当大一部分：（1）能被 $TLB$ 寻址到的内存；（2）$L_2$ 缓存。

实际上，它还会被 $L_2$ 缓存的相联度和替换策略进一步影响。在实践中，$m_c$ 通常选择为使得 $\tilde{A}$ 仅占用上述二者中较小者一半的 值。

Section7 Experiments

本节展示用前几节描述的技术来实现的 $DGEMMBLAS$ 例程（$routine$）获得的性能。

7.1 Algorithm Chosen

在典型的架构上实现 $Section5$ 中讨论的所有算法是一项艰难的工作，因为之前讨论了图8中的算法可能获得最优性能，因此我们实现了它。$GEBP\_OPT1$ 算法针对不同的体系结构使用不同的汇编代码实现，$A\to \tilde{A}$ 和 $B\to \tilde{B}$ 的例程采用 C 语言实现，因为编译器似乎能优化这些操作。

7.2 Parameters

图12展示了一些典型架构的物理和算法参数，并不是所有参数都在本文分析的考虑范围之内，全部列出只是为了完整性起见。

以下几个参数需要额外解释以下：

Duplicate

该参数指出矩阵 $B$ 的元素是否被复制（$duplicated$）为 $B$ 的打包的一部分。为了利用好 $Petium4(Northwood)$ 和 $Opteron$ 处理器上的 $SSE2$ 指令，这是很有必要的。尽管 $Core2Woodcrest$ 拥有 $SSE3$ 指令集，但复制指令（$instructionsforduplication$）由乘法单元发出（$issue$）且必须使用和 $Northwood$ 架构相同的技术。

Sustained Bandwidth

观察到的（$observed$）从指定（$indicated$）内存层（如果是 $L_2$ 缓存，那就是 $R_{load}$ ）存器组的持续带宽，单位 $doubles/cycle$。

Covered Area

能被 $TLB$ 寻址的内存大小；有些架构有一个（慢得多的）$level2TLB$，其相对于 $L_{1}TLB$ 的作用就和 $L_2$ 缓存相对于 $L_1$ 缓存的作用一样。是否通过 $L_{1}TLB$ 或 $L_{2}TLB$ 的条目数量来限制 $\tilde{A}$ 的大小（$size$），这取决于 $A\to \tilde{A}$ 和 $B\to \tilde{B}$ 的开销。

$\tilde{A}(Kbytes)$

为 $\tilde{A}$ 留出的内存空间。

7.3 Focus on the Intel Pentium 4 Prescott Processor (3.6 GHz, 64bit)

​等式（$3$）指出：为了用计算 $\tilde{A}$ 元素预取的开销，需满足 $n_r\ge \frac{R_{comp}}{2R_{load}}$。对于该架构，$R_{comp}=2$，$R_{load}=1.03$，故 $n_r\ge 2/(2\times 1.03)\approx 0.97$。同样，对于 $EM64T$ 架构，它与之前的唯一区别是寄存器数量由 $8$ 个增加到 $16$ 个。每个寄存器能够存储 $2$ 个双精度浮点数，取 $m_r\times n_r=4\times 4$，因此其中 $8$ 个寄存器用于存储这 $16$ 个双精度浮点。

​$k_c$的选择很复杂，因为在该架构中，对于计算 $\tilde{A}$ 的列和 $\tilde{B}$ 的列的内积的循环，其索引的更新应尽可能避免（为什么？）。因此该循环被完全展开（循环展开可以减少循环索引的更新），但这会导致循环的代码量增大，对应生成的指令就增多，使得将指令存储在指令缓存中成为一个问题。这把 $k_c$ 限制在了 $196$（$k_c$ 再增大，指令缓存就放不下了），略低于 $Section6.3$ 中论述的 $256$（ $k_c$ 个双精度浮点数占用半个页面）。

​图 $13$ 展示了不同的 $m_c$ 和 $k_c$ 组合对 $DGEMM$ 性能的影响，该架构是 “$\tilde{A}$ 应被 $TLB$ 寻址” 这条规则的例外，因为在该架构中，$TLBmiss$ 的成本比其他架构更低。当 $\tilde{A}$ 被选择为填充 $L_2$ 缓存一半大小时，性能略高于选择其为 $TLB$ 能寻址内存大小的一半。

​$Pentium4$ 的 $Northwood$ 版本依靠 $SSE2$ 指令来每周期计算 $2$ 个 $flop$，该指令要求复制 $B$ 中的条目，这种复制是被并入$B\to \tilde{B}$ 过程的一种数据移动，而 $Prescott$ 子架构所支持的 $SSE3$ 指令在 $B\to \tilde{B}$ 时无需这种复制。

​图 $14$ 展示了在 $Pentium4$ 架构上用该方法获得的性能（注：$m$，$n$，$k$ 是大矩阵的尺寸）：

$——$顶部的标签为 $Kernel$ 的曲线对应 $kernel$ 例程（即 $GEBP\_OPT1$）的性能。

$——$稍低的标签为 $dgemm$ 的曲线对应用 $GEPP$ 操作序列实现的 $DGEMM$ 例程的性能，其中 $GEPP$ 操作使用图8的算法实现的。

$——$底部的两条曲线对应 $A\to \tilde{A}$ 和 $B\to \tilde{B}$ 例程所耗时间的占比。

打包（$packing$）操作带来的开销几乎解释了从 $Kernel$ 曲线到 $dgemm$ 曲线下降的原因

​图 $15$ 的两幅图像（第二张好像缺失了）研究了 $m$ 和 $n$ 变化对实现性能的影响：

​上面一张是 $n=k=2000$ 时，$m$ 变化引起的四个指标的变化曲线，当 $m$ 较小时，就像 $GEPM$ 操作一样，$B\to \tilde{B}$ 没有被足够的计算分摊（$amortized$），从而导致相对较低的性能，一种解决方案可能是跳过 $B\to \tilde{B}$ 的过程，另一种解决方案可能是实现图 $9$ 中的算法；类似地，下面一张是 $m=k=2000$ 时，$n$ 变化引起的四个指标的变化曲线，当 $n$ 较小时，就像 $GEMP$ 操作一样，$A\to \tilde{A}$ 没有被足够的计算分摊，从而导致相对较低的性能，一种解决方案可能是跳过 $A\to \tilde{A}$ 的过程（这将要求 $GEMP$ 操作由 $AXPY$ 操作实现，而不是内积），另一种解决方案可能是实现图 $11$ 中的算法。

7.4 其他架构

对于剩余的架构，我们简要地讨论如何选择参数，并在图 $16-20$ 展示其性能

AMD Opteron processor (2.2GHz, 64bit)

​$n_r\ge 2/(2\times 0.71)\approx 1.4$；经观察，在寄存器中存储 $C$ 的条目的最佳选择是 $m_r\times n_r=4\times 4$。

​不像 $Pentium4$ 那样，计算 $\tilde{A}$ 的列和 $\tilde{B}$ 的列的内积的内层循环没有必要展开，$L_1$ 的大小也不是一个问题。因此，$k_c$ 可以取为使得 $\tilde{B}$ 的一列（$\tilde{B}$ 的一列就是 $k_c$ 个双精度浮点数）占据半页的值：$k_c=256$，取 $m_c\times k_c=384\times 256$ 可以使得 $\tilde{A}$ 填充 $TLB$ 可寻址空间的约三分之一。

​最新的 $Opteron$ 架构支持 $SSE3$ 指令，我们注意到复制 $\tilde{B}$ 的元素仍是有益的，这增加了$B\to \tilde{B}$ 的开销，使性能降低了约 3%。

​该架构的性能展示如图 $16$ 所示。

Intel Itanium2 processor (1.5GHz)

​对于浮点数计算，该架构天然忽略了 $L_1$数据缓存和 $L_{1}TLB$，从而 $Itanium$ 的 $L_2$ 缓存和 $L_3$ 缓存扮演了其他架构的 $L_1$ 缓存和 $L_2$ 缓存的角色，并且只有 $L_{2}TLB$ 是有意义的。因此 $n_r\ge 4/(2\times 2.0)=1.0$。由于有足够的可用寄存器（$128$ 个），取 $m_r\times n_r=8\times 8$（占 $32$ 个寄存器）。尽管 $k_c$ 的（理论）最佳选择是 $1K$（$1K$ 的双精度浮点刚好占据半页），但在实践中 $k_c=128$ 的性能几乎一样好。

​该架构有许多使得优化变得容易的特点：大量的寄存器，缓存和寄存器之间非常好的带宽，以及没有乱序执行。

​该架构的性能展示如图 $17$ 所示。

IBM POWRE5 processor (1.9GHz)

​对于该架构，$n_r\ge 4/(2\times 0.93)\approx 2.15$，$m_r\times n_r=4\times 4$。该架构有一个$D-ERAT$，相当于（其他架构的）$L_{1}TLB$，它的 $TLB$ 相当于（其他架构的） $L_{2}TLB$。$k_c=256$ 使得 $\tilde{B}$ 的一列刚好占据一个页的一半，选择 $m_c\times k_c=256\times 256$ 可以使得 $\tilde{A}$ 填充$TLB$ 可寻址空间的约四分之一。这是一种折中方案：$TLB$ 的速度相对较慢。通过将 $\tilde{A}$ 的内存占用空间控制在 $D-ERAT$ 可寻址的大小，可以获得更好的性能。

​该架构的性能展示如图 $18$ 所示。

PowerPC440 FP2 processor (700MHz)

​对于该架构，$n_r\ge 4/(2\ast 0.75)\approx 2.7$，$m_r\times n_r=8\times 4$。该架构有一个额外的复杂点：将 $\tilde{B}$ 和 $\tilde{A}$ 的元素从 $L_1$ 缓存和 $L_2$ 缓存移动到寄存器所需的组合带宽使得总带宽饱和，这意味着将 $C$ 的元素加载进寄存器不能被计算重叠（$overlapped$），从而 $k_c$ 应该取得非常大，以便在尽可能多的计算中分摊掉这个暴露的开销。选择 $m_c\times k_c=128\times 3K$，这使得 $\tilde{A}$ 填充了 $L_2$ 缓存的 $3/4$。

​到缓存的带宽缺乏、采用 $FIFO$ 策略的 $L_1$ 缓存和乱序执行的指令，使得对该架构的优化变得困难。由于页面大小较大，因此 $TLB$ 的可寻址空间很大。

​需要注意的是，$IBM$ 使用了类似于本文中讨论的技术来实现它们的矩阵乘法。

​该架构的性能展示如图 $19$ 所示。

Core 2 Woodcrest (2.66GHz) processor

​在本文的最终版修订时，$Core2Woodcrest$ 最近发布，因此该架构的性能指标特别有趣。

​对于该架构，$n_r\ge 4/(2\times 1.0)=2.7$，$m_r\times n_r=4\times 4$。对于 $Prescott$ 架构（为什么又提 $Prescott$ 架构？），有 $16$ 个可用寄存器，每个能存储 $2$ 个双精度浮点数，其中一半（$8$ 个）用于存储 $C$ 的 $m_r\times n_r$ 个条目（$entries$）。$\tilde{A}$ 的占用空间等于 $TLB$ 能寻址（$covered$）的内存空间。

​该架构的性能展示如图 $20$ 所示。

Section8 Conclusion

​我们对影响高性能矩阵乘性能设计的高层次问题（宏观层面，$highlevelissues$）给出了系统的分析，这些见解被整合到一个能在大量不同架构上取得极高性能的实现中。

​几乎所有目前 属于$LAPACK$一部分 的例程 $[Andersonetal.1999]$ 都会执行大量的 $GEPP$、$GEMP$ 或 $GEPM$ 操作。类似地，重要的 $BLASkernels$ 也能用这三种 $GEMM$ 的特殊情形来描述 $[Kagstrometal.1998]$。我们最近与 $FLAME$ 项目相关的研究表明：几乎所有这些例程都是 用 大量 $GEPP$ 来进行计算的 算法变体。当与 基于本文中的见解实现的 矩阵乘法例程进行接口时，这些可选择的算法变体将获得非常好的性能。

​其中一个不能 用大量的 $GEPP$ 来重定义 的操作是 $QR$ 分解。在 $QR$ 分解中，大约一半的计算可用 $GEPP$ 来描述，而另一半固有地需要 $GEMP$ 或 $GEPM$ 操作。而且 $panel$ 必须天然地取得很窄，因为它越宽，就必须执行更多的额外计算。这表明有必要进一步研究 这些 $GEMM$ 的特殊情形的 高性能实现。

​代码：http://www.tacc.utexas.edu/resources/software。

正文结束

以上均为个人理解，难免有理解错误或不到位的地方，欢迎指正。