Design and Implementation of a Highly Efficient DGEMM for 64-bit ARMv8 Multi-Core Processors

Design and Implementation of a Highly Efficient DGEMM for 64-bit ARMv8 Multi-Core Processors 针对64位 ARMv8八核处理器，设计并实现了一种基于 OpenBLAS 的高效 DGEMM。作者首先为此架构开发性能模型，然后根据理论指导用汇编语言系统地开发高度优化的 GEBP 内核。

性能模型表明，优化 DGEMM 的峰值性能（效率）需要在内存层次结构的所有级别上最大化其计算内存访问比率。而提高 GEBP 的性能的主要挑战在于为其最内部的寄存器内核选择合适的寄存器块大小。为了实现找到的最佳比率，文章通过以下方式优化寄存器内核中的操作：

• 利用循环展开、指令调度和软件实现的寄存器旋转；
• 利用 A64 指令支持高效的 FMA 操作、数据传输和数据预取。

最终，与 ATLAS 中的 DGEMM 相比，单核峰值性能（效率）从 3.88 Gflops（80.9%）提高到 4.19 Gflops（87.2%），八核从 30.4 Gflops（79.2%）提高到 32.7 Gflops（85.3%）。此外，单核上的实现效率非常接近从微基准测试获得的理论上限 91.5%，同时并行实现在评估的一系列矩阵大小的不同线程数下实现了良好的性能和可扩展性。

II. BACKGROUND

A. The 64-bit ARMv8 Multi-Core Processor

图 1 描绘了一个64位 ARMv8八核处理器。每个核心都有一个32 KB L1指令缓存和一个32 KB L1数据缓存。两个相邻的核心构成了所谓的双核心模块。同一模块中的两个内核共享一个256 KB 的 L2缓存，四个模块（具有8个内核）共享一个8 MB 的 L3缓存。每个核心都有一个支持 FMA 的浮点计算流水线，运行频率为2.4 GHz，提供4.8 Gflops 的峰值性能。

每个核心都是乱序执行的单线程四发射超标量设计。在每个内核中，64位 ARMv8架（aarch64）定义了32个128位浮点寄存器，v0–v31，可用于浮点向量计算。如图 1 所示，四个模块、L3缓存、两个内存桥和 I/O 桥都连接到一个缓存一致性结构。

B. Open-Source BLAS implementations

几个著名的开源 BLAS 实现有 netlib BLAS、ATLAS、GotoBLAS [5][6]、OpenBLAS 和 BLIS：

• netlib BLAS 是最初的参考实现，它没有利用内存层次结构，因此在现代计算机体系结构上表现不佳。
• ATLAS 依靠自动调优来提高缓存性能，使其非常适合在各种平台上生成 BLAS 库。
• OpenBLAS 是广泛使用但已停产的 GotoBLAS 的扩展，可在一系列架构上提供具有竞争力的性能。在 OpenBLAS 中，其内部内核称为 GEBP（基本计算单元），通常以汇编方式实现。
• BLIS 是一个用于快速生成 BLAS 实例的新框架，可以看作是对 GotoBLAS 的重新实现。BLIS 还采用 GotoBLAS 中的分层方法来增加代码重用，并将 GEBP 分解为所谓微内核上的双循环，从而促进 3 级 BLAS 的优化。

C. Overview of DGEMM in OpenBLAS

DGEMM 计算 $C:=\alpha \times AB+\beta \times C$，其中 $C$、$A$ 和 $B$ 分别是大小为 $M\times N$、$M\times K$ 和 $K\times N$ 的矩阵。在不失一般性的情况下，我们假设 $\alpha =\beta =1$，从而将 DGEMM 简化为 $C+=AB$。我们进一步假设矩阵总是以列优先顺序存储。

图 2 展示了 Goto[5][6] 提出的 DGEMM 算法，包括用于阻塞（以最大化缓存性能）和打包（以使数据能够高效地移动到寄存器）的多层。在本工作中，我们通过使用系统汇编开发高度优化的 GEBP，并通过分析选择其性能关键块大小，为 ARMv8体系结构获得了高效的 DGEMM。

第1层的最外层循环将 $C$ 和 $B$ 分别划分为尺寸为 $M\times n_c$ 和 $K\times n_c$ 的列板。第2层的下一个循环将 $M\times K$ 矩阵 $A$ 和 $B$ 的 $K\times n_c$ 子矩阵分别划分为大小为 $M\times k_c$ 的列板和大小为 $k_c\times n_c$ 的行板。$C$ 更新为一系列秩 $k_c$ 更新，即 DGEMM 由若干个通用板—板乘法( GEPP )组成，然后 $A$ 的每一个 $M\times k_c$ 面板通过第3层的第3个循环划分为 $m_c\times k_c$ 块。本质上，GEPP 分解为多个块板乘法（GEBP）调用，由于64位 ARMv8 架构中存在一个 L3缓存，我们假设 $B$ 的 $k_c\times n_c$ 面板将始终完全驻留在 L3缓存中[12]。

第4层处理的内核 GEBP 将 $C$ 的一个 $m_c\times n_c$ 面板更新为 $A$ 的一个 $m_c\times k_c$ 块和 $B$ 的一个 $k_c\times n_c$ 面板的乘积。为了保证连续访问，OpenBLAS 将 $A$ 的一个块和 $B$ 的面板打包成连续缓冲区。

如图 3 所示，打包 $A$ 涉及从 $A$ 的一个 $m_c\times k_c$ 块中提取一系列大小为 $m_r\times k_c$ 的薄片(子块)，并组织在 L2缓存中。类似地，打包 $B$ 相当于从 $B$ 的 $k_c\times n_c$ 面板中提取一系列大小为 $k_c\times n_r$ 的薄片，并将这些碎片组织到L3缓存中。为了摊销包装成本，将每个 $k_c\times n_r$ 薄片逐一移入 L1缓存。

让我们回到图 2 中的 GEBP。GEBP 包括由第 5 层的第五个循环和第 6 层的第六个循环形成的双循环。这个双循环负责实现上面讨论的打包过程。特别地，第五个环路将 $B$ 的一个 $k_c\times n_c$ 面板分割成 $k_c\times n_r$ 薄片，第六个循环将 $A$ 的一个 $m_c\times k_c$ 块分割成 $m_r\times k_c$ 的细条。第5层和第6层的计算分别称为 GEBS 和 GESS。GESS 又称 BLIS 中的微核，它利用 $A$ 的 $m_r\times 1$ 列子列和 $B$ 的 $1\times n_r$ 行子列对 $C$ 的 $m_r\times n_r$ 子块进行秩1更新，其中选择 $m_r$ 和 $n_r$ 作为微核的寄存器块大小。

当第7层执行每个秩1更新时，$C$ 的一个 $m_r\times n_r$ 子块、$A$ 的两个 $m_r\times 1$ 列子条和 $B$ 的两个 $1\times n_r$ 行子条将驻留在寄存器中，如图 3 所示。这里，$A$ 和 $B$ 的一个子分组是针对当前正在使用的 $A$ ($B$ )元素，另一个子分组用于下一个1级更新中所需要的即将到来的 $A$ 和 $B$ 元素。这最后一层称为寄存器内核。

III. PERFORMANCE MODELING

为64位ARMv8体系结构开发高效DGEMM的理论基础是什么？我们的出发点是关于 CPU 速度和内存速度之比的通用性能模型，即计算访存比。该模型清楚地表明，优化 DGEMM 的峰值性能(效率)需要在该体系结构的所有内存层次结构中最大化其计算内存比。此外，我们的模型还给出了 DGEMM 实现性能的下界。

考虑到图 4 所示的 ARMv8 架构中的内存层次结构，必须注意相对于数据上执行的计算量，有效地处理从内存到寄存器的数据移动量。我们假设执行单个操作需要固定成本 $\mu$，而不区分加法、减法和乘法。对于通信，我们还假设将一个单词从第 $i$ 级移动到第 $j$ 级的固定成本为 ${\nu }_{ij}$，将消息从第 $i$ 级移动到内存层次结构中的第 $j$ 级的固定成本为 ${\eta }_{ij}$。这里，一个词是一个浮点值，一个消息表示由连续几个词组成的缓存行。因此，${\nu }_{ij}$ 可以看作带宽的逆，${\eta }_{ij}$ 可以看作时延。如果我们暂时忽略计算和通信之间的重叠，那么程序的执行时间 $T$ 可以估计为：
$$T=F\mu +\sum _i\sum _j{W}_{ij}{\nu }_{ij}+\sum _i\sum _j{M}_{ij}{\eta }_{ij}$$
其中，$F$、$W_{ij}$ 和 $M_{ij}$ 分别表示操作、单词和消息的数量。例如，$W_{10}$ 表示从 L1缓存加载到寄存器的字数。

假设打包数据连续存储在慢速内存中，如图 2 中的第 4 层所示，我们假设在连续的计算中需要一个消息中的所有单词，即它们可以作为一个消息(一个缓存行)一起读或写。因此，移动消息数与移动词数之比几乎为常数：${\sum }_i{\sum }_j{M}_{ij}\backsimeq \kappa {\sum }_i{\sum }_j{W}_{ij}$。由于 ${\nu }_{ij}\ge 0$ 和 ${\eta }_{ij}\ge 0$，我们有：
$$T\le F\mu +(1+\kappa )\sum _i\sum _j{W}_{ij}\times (\sum _i\sum _j{\nu }_{ij}+\sum _i\sum _j{\eta }_{ij})$$

为了方便起见，设 $\pi ={\sum }_i{\sum }_j{\nu }_{ij}+{\sum }_i{\sum }_j{\eta }_{ij}$，$W={\sum }_i{\sum }_j{W}_{ij}$。然后，程序的计算内存访问比 $\gamma$ 可以表示为：
$$\gamma =\frac{F}{W}=\frac{F}{{\sum }_i{\sum }_j{W}_{ij}}$$

那么我们可以得到：

$$T\le F\mu +(1+\kappa )W\pi$$

由于重叠计算和通信对于提高性能而言是重要且必要的优化方法，因此我们提出了一个关于 $\gamma$ 的重叠因子函数 $\psi (\gamma )$。使用这个重叠因子，我们可以将上面公式细化为：

$$T_{\mathrm{opt}}\le F\mu +(1+\kappa )W\pi \psi (\gamma )$$

注意，若 $\gamma \to 0$ 则 $\psi (\gamma )\to 1$ 并且若 $\gamma \to \infty$ 则 $\psi (\gamma )\to 0$。此外，$\psi (\gamma )$ 通常是关于 $\gamma$ 的单调递减函数。从而：

$$T_{\mathrm{opt}}\le F(\mu +(1+\kappa )W\pi \frac{\psi (\gamma )}{\gamma })$$

最后，我们获得了 DGEMM 实现性能的下界：
$$Per{f}_{opt}=\frac{F}{T_{\mathrm{opt}}}\ge \frac{1}{(\mu +(1+\kappa )\pi \frac{\psi (\gamma )}{\gamma })}$$

这清楚地表明，较大的计算内存比 $\gamma$ 总是导致更好的峰值性能(效率)。

IV. FAST IMPLEMENTATION

基于我们的性能模型，我们通过在汇编中系统地开发高度优化的 GEBP 内核，为64位 ARMv8架构获得了高效的 DGEMM 实现。如图 2 所示，GEBP 包括第4-7层。它的开发涉及实现在第7层( 在第 II - C节中称为寄存器核 )执行的每个秩1更新，并确定四个层使用的各种块大小。我们将从第 7 层（寄存器内核）到第4层从内到外描述我们的 GEBP 实现，即从最快到最慢，跨越 ARMv8架构中内存层次结构的四个级别。因此，在某个级别确定的一些块大小稍后将用于确定较低级别的其他块大小。

在开发高度优化的 GEBP 时，主要挑战在于为其寄存器内核选择合适的寄存器块大小。我们以最大化从 L1数据缓存到寄存器的计算内存访问比为目标，分析地做出这样的选择。为了实现找到的最优比率，我们通过以下方式优化寄存器内核中的操作：

• 利用循环展开、指令调度和软件实现的寄存器旋转；
• 利用 A64指令来支持高效的 FMA 操作、数据传输和数据预取。

随后，我们通过最大化 GEBP 在所有三个级别的缓存中的计算访存比来优化 GEBP。考虑设置的关联性和替换策略，我们通过分析确定使用的其他块大小来实现这一点。

在第 IV-A 节中，我们描述了如何确定寄存器内核的寄存器块大小 $m_r\times n_r$ 以及相关的优化。在第 IV-B 节中，我们我们找到了分别对应于第6、5和4层的块大小 $k_c$、$m_c$ 和 $n_c$。如图 3 所示，$k_c$、$m_c$ 和 $n_c$ 分别由所使用的 L1、L2和 L3缓存决定。此外，我们还确定了如何插入预取指令来将数据预取到 L1 数据缓存中，以进一步加速寄存器内核中的操作。在第 IV-C 节中，由于缓存共享，我们在从串行实现转移到并行实现时调整块大小 $m_c$ 和 $n_c$。

A. Register Blocking

让我们关注第7层寄存器内核中发生的计算和数据移动（图 2 和图 3）。我们关心的是在其执行过程中最大化从 L1缓存到寄存器的计算内存访问比。当执行当前的 $2m_r{n}_r$ flops 时，$A$ 的一个 $m_r\times 1$ 列子条和 $B$ 的一个 $1\times n_r$ 行子条正在从 L1缓存加载到寄存器中。在这里，$B$ 的一个$1\times n_r$ 行子条总是驻留在 L1缓存中，$A$ 的一个 $m_r\times 1$ 列子条也通过有效预取数据的方式纳入到 L1缓存中。$C$ 子块中的 $m_r$ 元素总是驻留在寄存器中以参与连续操作。因此，计算内存访问比为：
$$\frac{2m_r{n}_r}{(m_r\times 1{)}_{L1\to R}+(n_r\times 1{)}_{L1\to R}}$$

其中下标 $L1\to R$ 表示数据从 L1缓存移动到寄存器。因此，优化问题是：

$$\max \gamma =\frac{2}{\frac{1}{n_r}+\frac{1}{m_r}}$$

受以下约束：
$$(m_r{n}_r+2m_r+2n_r)\times element\_size\le (n_f+n_{rf})\times p_f$$

这里，寄存器块大小 $m_r\times n_r$ 由可用寄存器的数量决定，元素大小是矩阵中的一个元素的字节大小，例如双精度为8个字节，$n_f$ 是可用的浮点寄存器的数量，$n_{rf}$ 是这些寄存器可重用于寄存器预加载的数量，最后，$p_f$ 是以字节为单位的浮点寄存器的大小，例如64位 ARMv8架构中为16字节。注意 $n_{rf}$ 必须满足以下约束：
$$0\le n_{rf}\times p_f\le (m_r+n_r)\times \mathrm{elementsize}$$

从这个优化公式可以看出，当 $m_r\approx n_r$ 时，可以最有效地分摊将数据加载到寄存器中的成本。然而，也有一些实际问题会影响 $m_r$ 和 $n_r$ 的选择。例如，在64位 ARMv8架构中，每个浮点寄存器可以容纳两个双精度值。因此，最好将 $m_r$ 和 $n_r$ 设置为2的倍数：
$$m_r=2i,n_r=2j,i=1,2,\cdots ,j=1,2,\cdots$$

在我们的设置中，$n_f=32$，$p_f=16$，元素大小为8。基于目标函数 (8) 及其三个相关约束 (9)–(11)，我们在图 5 中根据 $m_r$ 和 $n_{rf}$ 绘制了寄存器内核计算内存比的曲面图。为了获得最高的比率6.857，只需将重用浮点寄存器的数量设置为 $n_{rf}=6$，并且选择寄存器块大小 $m_r\times n_r$ 为 $8\times 6$ 或 $6\times 8$ 即可。由于一个缓存行有64个字节，当寄存器块大小选为 $m_r\times n_r=8\times 6$ 时，可以方便地预取 $A$ 中的元素，这将在第 IV-B 节中详细讨论。因此，我们使用24个寄存器 v8–v31 来存储 $C$ 的48个元素，使用8个寄存器 v0–v7 来存储 $A$ 的8个元素和 $B$ 的6个元素。

由于每个寄存器包含两个元素，因此很容易将24个寄存器分配给 $C$ 的48个元素。但是，剩下的8个寄存器怎么分配给 $A$ 的8个元素和 $B$ 的6个元素呢？一个简单的方法是只使用7个寄存器，剩下一个。本文采用了一种软件实现的寄存器轮换方案，为64位 ARMv8架构提供了一个更高效的寄存器内核，该体系结构比 x86 具有更少的物理寄存器用于寄存器重命名。由此，我们从以下优化问题入手：
$$\underset{s\in S_c}{\max }\underset{v_i\in V}{\min }Loc{(}^{\prime }{R}^{\prime }{,}^{\prime }N{F}^{\prime },v_i,s)-Loc{(}^{\prime }R{0}^{\prime }{,}^{\prime }C{L}^{\prime },v_i,s)$$

其中 V = { v 0 , v 1 , ⋯   , v 7 } V = \{v_0,v_1, \cdots, v_7\} V={v0​,v1​,⋯,v7​} 是 $A$ 和 $B$ 使用的寄存器集合，‘R’表示一个使用 $v_i$ 的读指令( fmla )，‘CL’和‘NF’表示 $v_i$ 中的当前值最后一次读取，同一个寄存器中的下一个值第一次读取，令 $S_c$ 为寄存器内核中 Read (fmla)指令所有执行顺序的集合，$Loc$ 表示某个指令在特定顺序中的位置。

在制定上式过程中，我们有目的地忽略了定义所有寄存器所使用的值的所有相应加载指令。在最佳解决方案中，$Loc{(}^{\prime }{R}^{\prime }{;}^{\prime }N{F}^{\prime },v_i,s)$ 处的指令与 $Loc{(}^{\prime }R{0}^{\prime }{,}^{\prime }C{L}^{\prime };v_i;s)$ 处的指令之间的距离尽可能大。在这两条指令之间有一条 load 指令，用于将后一个指令使用的值写入 $v_i$。当 $Loc{(}^{\prime }{R}^{\prime }{;}^{\prime }N{F}^{\prime },v_i,s)-Loc{(}^{\prime }R{0}^{\prime }{,}^{\prime }C{L}^{\prime };v_i;s)$ 较大时，有很好的机会定位加载指令，以防止管道停滞，这样加载的值在 $Loc{(}^{\prime }R{0}^{\prime }{,}^{\prime }C{L}^{\prime };v_i;s)$ 需要时需要时就已经可用。

为了解决上式，我们提出了一种软件实现的寄存器轮换方案，如表 I 所示，它展开控制寄存器内核执行的循环，即图 2 中描绘的最内层循环按8倍展开。其中，$\#i$ 表示循环体在原始寄存器内核中的第 $i$ 个副本。要执行 $\#i$，总共需要7个(128位)寄存器来存储 $A$ 的8个元素和 $B$ 的6个元素。为了预取 $\#(i+1)\%8$ 执行期间使用的 $A$ 和 $B$，还需要另外7个寄存器。然而，总共只有8个寄存器， { v 0 , v 1 , ⋯   , v 7 } \{v_0,v_1, \cdots, v_7\} {v0​,v1​,⋯,v7​} 可用。因此，如前所述，$n_{rf}=6$。这意味着在原始寄存器内核的连续两次迭代之间，即 $\#i$ 和 $\#(i+1)\%8$，重用了6个寄存器。图 6 说明了我们对 $\#i$ 和 $\#1$ 的寄存器分配方案，以及在每种情况下执行计算的顺序。求解上式得到了最佳距离7。

在为原始寄存器内核的每个副本 $\#i$ 调度指令时，我们还需要考虑它们的 WAR ( write-after-read，读后写)和RAW (write-after-read，写后读)依赖性。在每个副本 $\#i$ 中，执行24条浮点 FMA 指令(fmla)和7条加载指令(ldr)，以及一条预取指令(prfm)。由于寄存器重命名，必须仔细考虑 RAW，但 WAR 并不重要（稍后验证）。但是，我们仍然必须努力隐藏加载指令的延迟，防止它阻塞流水线，以便在需要时加载值可用。我们通过解决以下优化问题来实现：
$$\underset{s\in S_c}{\max }\underset{v_i\in V}{\min }Loc{(}^{\prime }{R}^{\prime },v_i,s)-Loc{(}^{\prime }{W}^{\prime },v_i,s)$$

其中，${}^{\prime }{R}^{\prime }$ 和 ${}^{\prime }{W}^{\prime }$ 分别表示读（fmla）和写（ldr）指令， V = { v 0 , v 1 , ⋯   , v 7 } V = \{v_0,v_1, \cdots, v_7\} V={v0​,v1​,⋯,v7​}，$S$ 是所有可能的指令执行顺序的集合，$Loc$ 表示指令在特定顺序中的位置。

在每个副本 $\#i$ 中，执行其24个浮点 FMA 指令的顺序是固定的，沿着图 6 中的之字形路径。因此，上式的优化问题可以简化为寻找适当的点以插入所有需要的加载指令。图 7 显示了这些指令是如何调度的，以行优先执行顺序展示。绿色寄存器在 $\#i$ 中加载，红色寄存器在 $\#(i-1)\%8$ 中加载(也就是说 ,在原始寄存器内核中较早的一次迭代)。我们可以观察到，上式的最优距离解为9。最后，图 8 给出了8×6大小的寄存器内核循环体展开的汇编语言代码片段。

B. Cache Blocking for Locality

让我们考虑图 2 中第 6 层的 GESS。来自 $B$ 的 $k_c\times n_r$ 子条的元素被多次重用，因此保留在 L1缓存中(图 3)。同时，$A$ 的 $m_r\times k_c$ 子条中的元素应该从 L2缓存加载到 L1缓存中。在执行 $2m_r{n}_r{k}_c$ flops 时，需要将 $C$ 子块中的 $m_r\times n_r$ 元素加载到寄存器中，并将更新后的结果从寄存器存储到内存层次结构中。那么可以估计 GESS 的计算内存访问比为：
$$\frac{2m_r{n}_r{k}_c}{(m_r{k}_c{)}_{L2\to L1}+(m_r{k}_c{)}_{L1\to R}+(k_c{n}_r{)}_{L1\to R}+(2m_r{n}_r{)}_{M\leftrightarrow R}}$$

通过进行类似的分析，我们还可以将第5层 GEBS 的计算访存比表示为：
$$\frac{2m_c{n}_r{k}_c}{(m_c{k}_c{)}_{L2\to L1}+(m_c{k}_c{)}_{L1\to R}\lceil \frac{m_c}{m_r}\rceil +(k_c{n}_r{)}_{L1\to R}+(2m_c{n}_r{)}_{M\leftrightarrow R}}$$
在实际中，$m_c$ 是 $m_r$ 的整数倍。那么 GESS 与 GEBS 的比值非常接近，因此简化为：
$$\gamma =\frac{2}{\frac{2}{{\mathbf{n}}_{\mathbf{r}}}+\frac{1}{{\mathbf{m}}_{\mathbf{r}}}+\frac{2}{k_c}}$$

其中粗体的 ${\mathbf{m}}_{\mathbf{r}}$ 和 ${\mathbf{n}}_{\mathbf{r}}$ 已在第 IV-A 节中固定。从上式可以看出，如果 ${\mathbf{n}}_{\mathbf{r}}$ 大于 ${\mathbf{m}}_{\mathbf{r}}$，那么将 $A$ 的元素从 L2缓存加载到 L1缓存的代价就更容易分摊。然而，在我们的实现中，考虑到一个缓存行有64个字节，我们选择了寄存器块大小为 $m_r\times n_r=8\times 6$。因此，$A$ 的每个子条（8个元素）都可以在 L2到 L1的一条缓存行中预取。此外，$C$ 元素在寄存器中的加载不能与计算重叠。但是，我们可以将 $C$ 元素存储回内存的过程与计算重叠起来。因此，如果 $m_r$ 和 $n_r$ 已知，则上式中的比率 $\gamma$ 在 $k_c$ 尽可能大的情况下最大化，以分摊更新 $C$ 元素的成本。

在我们的实现中，L1缓存采用 LRU 替换策略，并且多次重用了 $B$ 的 $k_c\times n_r$ 子条。因此，很容易将 $B$ 保存在 L1缓存中。同时，为了优化性能，还应该将 $A$ 的至少两个 $m_r\times k_c$ 列子条和 $C$ 的一个 $m_r\times n_r$ 子块加载到 L1缓存中，而不会造成与 $B$ 元素的冲突未命中。为了考虑到 L1缓存的集合关联性，我们通过施加两个约束遵循[14]：
$$\begin{array}{rl}{k}_c\times n_r\times \mathrm{elemen}{\mathrm{t}}_{\mathrm{s}}\mathrm{ize}& \le \frac{(assoc1-k1)\times L1}{assoc1}\\ (m_r\times n_r+m_r\times 2)\times \mathrm{elemen}{\mathrm{t}}_{\mathrm{s}}\mathrm{ize}& \le \frac{k1\times L1}{assoc1}\end{array}$$

其中，$m_r=8$，$n_r=6$ 如前所述，$\mathrm{elemen}{\mathrm{t}}_{\mathrm{s}}\mathrm{ize}=8$，$L1=32K$ 是 L1数据缓存的大小（以字节为单位），$assoc1=4$ 为 L1缓存的路数，$k1$ 是满足 $0<k1<assoc1$ 的整数。很容易看出 $k1$ 越小，$k_c$ 越大，这使得我们可以得出 $k1=1$ 和 $k_c=512$ 的结论。这意味着 $B$ 的一个 $k_c\times n_r$ 的子条填充了 L1数据缓存的3/4。

接下来，我们最大化图 2 中第4层 GEBP 的计算内存访问比。我们假设 $A$ 的一个 $m_c\times k_c$ 块已经驻留在 L2缓存中，$B$ 的一个 $k_c\times n_c$ 面板已经驻留 在L3缓存中(图 3)。因此，将 $2m_c{k}_c{n}_c$ 除以 $(m_c{k}_c{)}_{L2\to L1}\lceil \frac{n_c}{n_r}\rceil +(m_c{k}_c{)}_{L1\to R}\lceil \frac{n_c}{n_r}\rceil +(k_c{n}_c{)}_{L1\to R}\lceil \frac{m_c}{m_r}\rceil +(k_c{n}_c{)}_{L3\to L2}+(k_c{n}_c{)}_{L2\to L1}+(2m_c{n}_c{)}_{M\leftrightarrow R}$，即可获得其计算内存访问比。与之前类似，我们得到：
$$\gamma =\frac{2}{\frac{2}{{\mathbf{n}}_{\mathbf{r}}}+\frac{1}{{\mathbf{m}}_{\mathbf{r}}}+\frac{2}{{\mathbf{k}}_{\mathbf{c}}}+\frac{2}{m_c}}$$

其中粗体的 ${\mathbf{m}}_{\mathbf{r}}$、${\mathbf{n}}_{\mathbf{r}}$ 和 ${\mathbf{k}}_{\mathbf{c}}$ 之前已确定。我们观察到，数据移动量 $(k_c{n}_c{)}_{L3\to L2}$ 可以与 $2m_c{k}_c{n}_c$ 运算重叠，但由 $(k_c{n}_r{)}_{L2\to L1}$ 表示的数据移动量只能与 $B$ 中涉及的相同 $k_c\times n_r$ 薄片的最后 $2m_r{k}_c{n}_r$ 运算重叠。$m_c$ 越大，数据移动量 $(k_c{n}_r{)}_{L3\to L2}$ 就可以隐藏的越好。与建立公式 (15) 的情况类似，我们得到：
$$\begin{array}{rl}{m}_c\times k_c\times \mathrm{elemen}{\mathrm{t}}_{\mathrm{s}}\mathrm{ize}& \le \frac{(assoc2-k2)\times L2}{assoc2}\\ k_c\times n_r\times \mathrm{elemen}{\mathrm{t}}_{\mathrm{s}}\mathrm{ize}& \le \frac{k2\times L2}{assoc2}\end{array}$$

其中 $k_c=512$，$n_r=6$，$L2=256K$ 是 L2缓存的大小（以字节为单位），$assoc2=16$ 是 L2缓存中的路数，$k2$ 是满足 $0<k2<assoc2$ 的整数，且 $\mathrm{elemen}{\mathrm{t}}_{\mathrm{s}}\mathrm{ize}=8$。由公式(16)可得，$m_c$ 尽可能大。由公式(17)得到 $k2=2$，$m_c=56$。因此，$A$ 的 $m_c\times k_c$ 块填充了 L2缓存的7/8，$B$ 的 $k_c\times n_r$ 薄片占据了 L2缓存的1/8以下。

我们现在准备分析 $A$ 和 $B$ 所需的预取距离。由于 与 $A$ 的每个 $m_r\times kc$ 子条相乘时，$B$ 的 $k_c\times n_r$ 薄片总是驻留在 L1缓存中，因此后者不需要预取。只需要在 $B$ 的当前条与 $A$ 的最后一条相乘的过程中，预取 $B$ 的下一个 $k_c\times n_r$ 条到 L2缓存中即可。此时设置预取距离为 ${\mathrm{PRE}}_B=k_c\times n_r\times \mathrm{elemen}{\mathrm{t}}_{\mathrm{s}}\mathrm{ize}=24576$。为了让所有对 $A$ 的子条的访问都命中 L1 缓存，我们使用较短的距离来预取 $A$：${\mathrm{PRE}}_A={\alpha }_{prea}\times num\_unroll\times m_r\times element\_size=2\times 8\times 8\times 8=1024$。

最后讨论如何选择 $n_c$。从图 2 第3层执行的 GEPP 可以看出，$n_c$ 应该尽可能大，以分摊 $A$ 的一个 $m_c\times k_c$ 块从 L3缓存到 L2缓存的数据移动成本。因此，对于 nc，我们有公式(15)公式17)的类似物：
$$\begin{array}{rl}{k}_c\times n_c\times \mathrm{elemen}{\mathrm{t}}_{\mathrm{s}}\mathrm{ize}& \le \frac{(assoc3-k3)\times L3}{assoc3}\\ m_c\times k_c\times \mathrm{elemen}{\mathrm{t}}_{\mathrm{s}}\mathrm{ize}& \le \frac{k3\times L3}{assoc3}\end{array}$$

其中 $m_c=56$，$k_c=512$ 如前所述，$L3=8M$ 是 L3缓存的大小（以字节为单位），$assoc3=16$ 是 L3缓存的路数，$k3$ 是满足 $0<k3<assoc3$ 的整数。很容易得到 $k3=1$ 和 $n_c=1920$，意味着 $B$ 的一个 $k_c\times n_c$ 面板占据了 L3缓存的15/16，而 $A$ 的 $mc\times kc$ 薄片低于1/16。

C. Cache blocking for Parallelism

我们讨论如何根据先前发现的块大小 $m_r$、$n_r$ 和 $k_c$，在多线程环境下调整块大小 $m_c$ 和 $n_c$。如图 9 所示，并行化第3层的循环，使得共享 L3缓存具有更好的局部性，其中存储了 $B$ 的所有 $k_c\times n_c$ 行板。在这种情况下，每个线程分配到 $A$ 的一个不同的 $m_c\times k_c$块，所有线程共享同一个 $B$ 的 $k_c\times n_c$ 行板，然后每个线程将自己的 $A$ 块与 $B的共享行板相乘。文献[15]中广泛讨论了在第3层并行此特定循环的策略。在这里，64位A RMv8多核处理器也采用了它。

再次考虑图 1 中的64位 ARMv8架构。在不失一般性的情况下，我们考虑了一个 DGEMM 的8线程并行实现，每个内核一个线程。由于一个模块中的两个核心共享一个256KB 的 L2缓存，因此对应的两个线程在同一个 L2缓存中将有各自的 $A$ 块。因此，公式(17)需要修改为：
$$\begin{array}{rl}2\times m_c\times k_c\times \mathrm{elemen}{\mathrm{t}}_{\mathrm{s}}\mathrm{ize}& \le \frac{(assoc2-k2)\times L2}{assoc2}\\ 2\times k_c\times n_r\times \mathrm{elemen}{\mathrm{t}}_{\mathrm{s}}\mathrm{ize}& \le \frac{k2\times L2}{assoc2}\end{array}$$

给定 $n_r=6$ 和 $k_c=512$，求解这两个方程得到 $m_c=24$ 和 $k2=4$。类似地，公式(18)的多线程版本是：
$$\begin{array}{rl}{k}_c\times n_c\times \mathrm{elemen}{\mathrm{t}}_{\mathrm{s}}\mathrm{ize}& \le \frac{(assoc3-k3)\times L2}{assoc3}\\ 8\times m_c\times k_c\times \mathrm{elemen}{\mathrm{t}}_{\mathrm{s}}\mathrm{ize}& \le \frac{k3\times L3}{assoc3}\end{array}$$

通过求解这些约束，我们得到 $n_c=1792$ 和 $k3=2$。因此，$A$ 的八个块将与 $B$ 的共享行面板一起驻留在 L3缓存中。

V. EXPERIMENTAL RESULTS

我们在表 II 中描述的64位 ARMv8八核计算平台上进行了评估。我们展示了我们在 OpenBLAS 中通过高度优化的 GEBP 实现的 DGEMM 是该体系结构中可用的最快的。我们的串行实现的峰值性能（效率），接近于从微观基准测试中获得的理论上限。我们的并行实现在一系列不同的矩阵大小上实现了良好的可扩展性。此外，我们（在 OpenBLAS 中）的 DGEMM 实现在单核心上比 ATLAS 中的 DGEMM 实现高 7.79%，在八核上高 7.70%。请注意，ATLAS 中的 DGEMM 实现也通过汇编中高度优化的 GEBP 加速（其寄存器块大小仅为5×5），这是以前64位 ARMv8架构的唯一可用的。

在性能模型的基础上，我们通过生成高度优化的 $8\times 6$ GEBP 来优化 DGEMM，使其计算内存访问比最大化。为了验证我们的性能模型，我们还开发了另外两个 GEBP 实现，如图 10 所示：(a) 4×4 GEBP，其寄存器块大小为 $4\times 4$；(b) $8\times 4$ GEBP，其寄存器块大小为 $8\times 4$。在 $4\times 4$ 内核中，执行 $4\times 2$ 矩阵 $A$ 和2×4矩阵 $B$ 的乘积。$8\times 4$ 内核可以看作是 $8\times 6$ 内核的简化版本。按照第 IV-B 和 IV-C 节所述继续进行，可以找到单线程和多线程设置的块大小，如表 III 中所示。

A. Microbenchmarking for Register Block Sizes

我们做了一些微基准测试来证明我们的8×6寄存器块大小带来了良好的性能，并且我们的指令调度优化是有效的。第7层的寄存器内核，即图2所示的最内层循环，影响了整个执行时间。在表 I 中用 $\#i$ 表示的每个迭代中，有 $(m_r+n_r)/2$ 条128位内存指令将数据加载到寄存器中，以及 $m_r{n}_r/2$条128位 NEON 浮点 FMA 指令，如图 8 所示。算术指令占总数的百分比为 $(m_r{n}_r/2)/(m_r{n}_r/2+(m_r+n_r)/2)$。

有趣的是，分析不同负载比（记为 $\mathrm{LDR}:\mathrm{FMLA}$）下 FMA 指令的效率。我们编写了一个微基准测试，其中的指令是独立且均匀分布的，以避免指令延迟对我们的实验造成任何影响。这个微基准可以始终将数据保存在 L1缓存中。我们在表 IV 中报告了我们的发现。请注意，$1:2$、$46:164$ 和 $7:24$ 分别是大致对应于 $4\times 4$、$8\times 4$ 和 $8\times 6$ GEBP 实现的 $\mathrm{LDR}:\mathrm{FMLA}$ 比率。一个关键观察是增加算术指令占总数的百分比，即 $(m_r{n}_r/2)/(m_r{n}_r/2+(m_r+n_r)/2)$，可以提高性能。$4\times 4$、$8\times 4$ 和 $8\times 6$ GEBP 实现的百分比为66.7%、72.7%和77.4%。
因此，预计 $8\times 6$ GEBP 将是表现最好的。

在 DGEMM 实现中，有两种类型的内存延迟：WAR 和 RAW。我们修改了我们的微基准测试，让将数据加载到寄存器的内存指令跟在从同一寄存器读取的 FMA 指令后面。仍是相同的效率。因此，可能由于使用了寄存器重命名机制，似乎不需要考虑 WAR 延迟。

我们还尝试了不同的指令执行顺序，以观察 RAW 延迟对效率的影响。我们发现在执行至少 4 条 fmla 指令后，加载指令中使用的寄存器可以被释放。我们在图 7 中的指令调度满足这个约束。因此，在存在数据依赖性的情况下，可以通过应用指令调度优化来隐藏内存延迟，并获得表 IV 中相同的效率。

此外，在这些微基准测试中获得的结果表明，寄存器块大小越大，可以实现的效率就越高。由于我们的实验中没有 L1 缓存缺失，表 IV 中所示的效率可以看作是我们 DGEMM 实现的上限。

B. Performance Analysis

我们考虑4种 DGEMM 实现：OpenBLAS-8×6、OpenBLAS-8×4、OpenBLAS-4×4和ATLAS-5×5。前三个是我们在 OpenBLAS 中分别使用 $8\times 6$、$8\times 4$ 和 $4\times 4$ GEBP 内核开发的。ATLAS-5×5是在 ATLAS 中用 $5\times 5$ GEBP 内核开发的。在我们的实验中，使用大小从256到6400不等的方阵，步长为128。对于每个矩阵大小，执行时间以5次运行的平均值来衡量。

我们在图 11（单线程）和图 12（八线程）中展示了我们的性能结果。OpenBLAS-8×6在几乎所有测试的输入尺寸中表现最佳。特别是，OpenBLAS-8×6在所有输入大小上均优于 ATLAS-5×5，表明 OpenBLAS-8×6是64位 ARMv8架构中最快的 DGEMM。

表 V 总结了这些方法在串行和多线程设置中的峰值和平均效率。OpenBLAS-8×6在所有评估的四个指标中表现最好。特别是，OpenBLAS-8×6的单线程峰值效率为87.2%，非常接近我们的微基准测试（表 IV）获得的理论上限效率91.5%。与 ATLAS-5×5相比，OpenBLAS-8×6的单核峰值性能（效率）从3.88Gflops（80.9%）提高到4.19Gflops（87.2%），八核上从30.4 Gflops（79.2%）提高到32.7Gflops（85.3%）。这使得峰值性能（效率）的方面，单核提高了7.79%，八核提高了7.70%。此外，OpenBLAS-8×6的平均效率较之 ATLAS-5×5在单核上提高了8.55%，在八核上提高了10.79%。

对于 OpenBLAS-8×6、OpenBLAS-8×4、OpenBLAS-4×4和 ATLAS-5×5，其寄存器内核的计算内存比分别由公式(8)估计为6.86、5.33、4和5。对于64位 ARMv8架构，用于寄存器重命名的物理寄存器少于 x86，图 13 清楚地表明了寄存器旋转方案的有效性。

图 14 展示了 OpenBLAS-8×4在四种线程配置（1、2、4和8个线程）下的性能和可扩展性，还显示了它们的块大小 $m_r\times n_r\times k_c\times m_c\times n_c$。在所有情况下，$m_r$、$n_r$ 和 $k_c$ 相同，分别在第 IV-A 和 IV-B 节中得到。在1和8线程情况下，$m_c$ 和 $n_c$ 分别在第 IV-B 和 IV-C 节中找到。对于两个线程，我们通过用2替换8并求解公式(17)得到 $m_c$，求解公式(20)找到 $n_c$。对于4个线程，我们找到 $m_c$ 和 $n_c$ 的方法与两个线程类似。在2和4线程的情况下，不同的线程总是在不同的模块上运行（图 1），这样一个模块中运行的线程就可以独自使用整个 L2缓存。对于每个线程配置，所有线程共享相同的 L3缓存。这些结果表明 OpenBLAS-8×6不仅高效而且可扩展。

表 VI 比较了 OpenBLAS-8×6 在选择 $k_c\times m_c\times n_c$ 的几种不同块大小时的性能结果。根据[5]，这些块大小通常用于使 $A$ 的一个 $m_c\times k_c$ 块占用大约一半的 L2缓存，而 $B$ 的一个 $k_c\times n_r$ 块占用大约一半的 L1缓存。在串行设置中，我们使用了根据[5]获得的 $k_c\times m_c\times n_c=320\times 96\times 1536$。当我们在选择块大小时考虑缓存的集合关联性和替换策略时，我们在两个设置中的选择会使得效率更高。

最后，图 15 比较了 OpenBLAS-8×6、OpenBLAS-8×4和 OpenBLAS-4×4的 L1-dcache-loads 数量。在串行和并行设置中，OpenBLAS-8×6的 L1-dcache-loads 数量最少。可以合理地假设所有这些 DGEMM 实现都表现出相同数量的浮点运算。因此，OpenBLAS-8×6优于其他两种，因为它可以更有效地将计算与数据移动（尤其是数据从 L1数据缓存加载到寄存器）重叠。根据表 V，OpenBLAS-8×6串行实现的峰值和平均效率略好于并行，因为并行执行受到更多 L1-dcache-loads 的影响。

表 VII 比较了 OpenBLAS-8×6、OpenBLAS-8×4和 OpenBLAS-4×4 在串行和并行设置下的 L1-dcache-load-miss 率。请注意，OpenBLAS-8×6在两种情况下都不会产生最低的L1数据缓存未命中率。然而，由于执行的 L1-dcache-loads 的数量最少，它最终获得了最高的效率。这说明在 ARMv8架构中，L1缓存缺失率并不像其他架构[7]中那样对性能至关重要。

参考资料：

• Register rotation
• CS252 Spring 2017 Graduate Computer Architecture Lecture 10: VLIW Architectures
• 717-高等计算机系统结构
• Rotating Register Allocation for Enhanced Pipeline Scheduling
• GOTOBLAS一般矩阵乘法高效实现机制的研究
• 开源矩阵计算库OpenBLAS与深度学习
• 异构HPL算法中CPU端高性能BLAS库优化
• GPU程序优化（三）——矩阵转置程序优化实例（进阶版）
• 高性能计算初步：矩阵乘法
• 面向多核向量处理器的矩阵乘法向量化方
• 基于MLIR实现GEMM编译优化
• 优化CPU矩阵乘法
• [LAFF]5.3.4 矩阵-矩阵乘法的算法_rank-1 update 线性代数从基础到尖端.德克萨斯州立大学at奥斯汀
• 使用报表向导创建报表
• CPU缓存中的路数
• 阅读ARM Memory(L1/L2/MMU)笔记
• ARMv8官方手册学习笔记(八)：Cache和内存层级
• Cache Line操作和Cache相关概念介绍
• 如何测试各级 cache 的访问时延
• CPU缓存行与伪共享
• 计算机系统基础（七）高速缓存Cache
• CPU 缓存一致性与内存屏障
• 存储器 - 高速缓存(CPU Cache)：为什么要使用高速缓存？
• CPU缓存 - CPU cache
• 细说Cache-L1/L2/L3/TLB
• google/XNNPACK
• Best Practice Guide - ARM64
• ComputeLibrary/src/core/NEON/kernels/arm_gemm/kernels/a64_hgemm_8x24
• EMLL/src/neon_armv8a/extension/HgemmKernel.c
• 数据依赖RAW WAW WAR
• 大金哥的超标量处理器学习笔记之7——寄存器重命名
• 指令cache 为什么比数据cache失效率低？
• Lecture 13: Cache V
• 计算机组成与设计 硬件/软件接口（二）