[BLIS参数分析，论文个人阅读分享]Analytical Modeling Is Enough for High-Performance BLIS

Analytical Modeling Is Enough for High-Performance BLIS

* Authors: Tze Meng Low, Francisco D. Igual, Tyler M. Smith, Enrique S. Quintana-Orti

• DOI: 10.1145/2925987
  

---

阅读内容总结

1 论文摘要(Abstract)

展示了类似BLAS的库实例化软件(BLIS)框架，提供了GotoBLAS(现在维护为OpenBLAS)实现的更详细的分层，允许在数学上分析地确定矩阵乘法的高端实例化的调优参数。从多级cache的角度建立了关于GotoBLAS方法的五个参数$m_r$,$n_r$,$m_c$,$k_c$,$n_c$的数学模型，为BLIS在复杂高性能体系结构的应用提供了一种参数分析方法

2 研究动机(Motivation/Introduction)

实现更多高性能BLAS的需求在于：密集线性代数(Dense Linear Algebra, DLA)领域是最早认识到、理解并促进对基本原语的接口进行标准化(事实上的)以支持可移植的高性能的领域之一，在近四十年来，这些基本原语一直是BLAS程序，而新的体系结构的出现就需要新的BLAS实现。为此目标，Intel 维护了MKL库，AMD维护了ACML库，IBM维护了ESSL库，同时还有长久以来的OpenBLAS库和ATBLAS，问题在于对于每个待支持的目标体系结构，这些解决方案都需要大量的时间和人力。
为了减少开发高性能实现的工作量，ATLAS和其处理器PH i PAC提出了对高性能计算领域的自动调优，这两个库背后的基本原理是，优化对于人类专家来说既困难又耗时，而自动生成与自动调优相结合就是答案。但是Yotov的工作表明，被自动调优社区鼓吹为成功自调优案例的自动调优操作可能是不必要的操作，原因在于ATLAS使用的内部核是由人工手动调优，编码得到的，确切地说，这个黑盒内核的存在使得ATLAS的分析建模存在问题，因为一些调优参数隐藏在内部。
GotoBLAS也是由人持续维护的内核，这导致让不完全理解内核复杂细节的人分析变的很难，直到最近，还没有详细说明 GotoBLAS 内部内核的复杂性的论文。
BLIS (BLAS-like Library Instantiation Software) 是一个用于实例化BLAS的新框架。从实现的角度来看，BLIS模块化了GotoBLAS2底层的方法(现在也被许多其他实现所采用)，以隔离执行GEMM的微内核，在该微内核上可以构建所有3级BLAS功能。因此，对于给定的体系结构，程序员只需要为BLIS开发一个高效的微内核，并调整几个关键参数值，以优化围绕该微内核的循环，从而自动实例化所有3级BLAS。重要的是，BLIS的分层特性自然地暴露了需要调优的参数。虽然BLIS简化了BLAS(类BLAS)操作的实现，但为目标体系结构优化BLIS的关键步骤是开发人员确定微内核及其周围循环的特定参数值。传统的解决方案是应用实证研究。在本文中，我们采用另一种模型驱动的方法来分析地确定BLIS所需的参数值。特别是，我们对算法和架构进行了类似于Yotov等人的分析，以确定BLIS框架中根深蒂固的算法的数据使用模式，并基于现代处理器架构中的硬件特征构建分析模型。
本文主要工作如下：
1. 分析最为人知的BLAS算法，也就是BLIS的底层分析。不选择ATLAS的原因是手工编码的BLIS/GotoBLAS/OpenBLAS实现通常比通过ATLAS自动生成的实现产生更高的性能。
2. 建立了更现代的模型，相较于Yotov和Whale用于分析的架构。具体地说，我们的模型比以前的模型更通用，因为它允许现代硬件功能，如矢量寄存器文件和具有多级集关联数据缓存的内存层次结构。
3. 建立了更全面的模型。分析模型更加全面，因为它包含了表征微内核的参数值。这些值是开发人员必须根据经验确定的。这一点很重要，因为ATLAS提供的最佳实现通常涉及围绕手工编码(黑盒)内核的循环。由于手工编码内核的内部未知，Yotov等人[2005]中的模型无法识别全局算法的参数值。
4. 建立了更具可操作性的模型。我们的分析模型旨在提供实际值，开发人员可以使用这些值为特定的体系结构定制BLIS框架。因此，它不同于现有的努力，如Whaley等人[2001]所描述的那些“模型”被用作启发式方法来指导经验搜索引擎。
5. 与人工编码的实现进行性能比较。先前工作的一个缺点是，它不能生成与使用手工编码内核的实现相竞争的代码。我们的文章表明，通过分析获得的参数值可以产生与手动调优实现相媲美的性能，从而实现一流的性能。

3 符号系统(Notation)

4 结论(Conclusion)

• 开发了一个分析模型，该模型可以为 以BLIS 为基础以及其他手动调整的高性能密集线性代数库的 GEMM 算法生成最佳参数值。分析方法的关键是认识到 BLIS 的许多参数都受到硬件特征的限制。通过理解算法如何与处理器架构交互，可以选择更积极地利用硬件特性(如系统的多层高速缓存)的参数值。

• 将通过文章分析方法获得的值与人工优化实现(例如，由OpenBLAS背后的专家，一个没有作者参与的开源项目)的值进行比较，并报告它们即使不相同也相似。与与ATLAS相比的类似工作不同，本文是第一篇证明识别参数值的分析方法可以与同类最佳的专家优化实现相媲美的文章。这是证明可以在不依赖昂贵的经验搜索的情况下获得GEMM的高性能实现的重要一步。

• 作者关注了模型的扩展，以包含类似 LAPACK 的复杂线性代数操作。具体而言，矩阵因子分解和 GEMM 操作（通用矩阵乘法）作为块算法的一部分成为焦点，提出了寻找最佳块大小的问题，这需要利用对于 GEMM 操作性能有利的参数知识。

• 文章分析方法存在的一些问题是：

  1. 文章讨论了纳入数据移动考虑因素的必要性。目前微内核（$m_r$ 和 $n_r$）参数的确定仅考虑浮点算术指令的延迟，但这基于一个假设，即系统具有足够大的带宽。然而，在低功耗系统中，这一假设可能不成立，因此需要调整参数（$m_r$ 和 $n_r$）以确保在加载下一组元素的时间内执行足够长的计算时间，以隐藏数据加载的时间。

  2. 文章指出当前分析模型仅考虑了双精度算术，而在复杂算术中，相同大小的微分块涉及更多的浮点运算和数据量，因此需要对模型的公式进行更新。

  3. 作者提到了与专用硬件设计的关系。问题涉及到是否可以利用当前的分析模型来更好地设计未来的计算架构。此处引用了 Pedram 等人的研究，这些研究探讨了为线性代数计算设计专用硬件的方法。作者认为当前用于微内核参数值确定的模型可能被用于确定专用硬件的理想缓存大小和/或缓存替换策略。

5 主体工作

5.1 GEMM的高性能算法基础分析（GotoBLAS算法，rank-1 更新）

1. 算法简述
   

   利用外积操作计算GEMM的分块乘法示例，6重循环可以将LOOP 6简化为$C_c(i_r:i_r+m_r-1,j_r:j_r+n_r-1)+=A_c(i_r:i_r+m_r-1,0:k_c-1)\cdot B_c(0:k_c-1,j_r:j_r+n_r-1)$
   这一算法将B先分块为$k_c\times n_c$并打包到L3cache，再将A分块到$m_c\times k_c$并打包到L2cache，然后通过更小的分块$m_r\times n_r$计算C矩阵的值。

2. Packing
   
   该图是算法中对A，B元素的访问顺序。
   通过引入数据封装可以：

   1. 确保这两个矩阵的元素在微内核内部以跨行访问的方式读取；
   2. 确保指令访问的是对齐的数据；
   3. 数据可以预加载进对应的cache层级中

3. 数据重用率展示
   

4. 需要确定的关键参数
   主要是5个：$m_c$,$k_c$,$n_c$,$m_r$,$n_r$；每层循环中利用的参数不尽相同，循环6由五个参数中的三个定义($m_r$,$n_r$和$k_c$)，而循环3,4和5由五个参数中的四个定义。此外，循环3、4和5的四个参数中的三个也在下一个内部循环中。现在，当循环5的参数已经确定时，循环4的唯一未知参数是$n_c$。这一观察结果表明，参数应该从最内层的循环向外识别。在下一节中，我们将利用这一观察结果对参数集合进行分析优化。

5.2 参数的分析与确定

该图是分析的关键，就是各级存储层次中需要驻留的数据内容

1. 首先进行计算机结构的一些假设和符号设置
   • 假设处理器中有$N_{REG}$个向量寄存器，每个向量寄存器中有$N_{VEC}$个单元，每个单元是一个$S_{DATA}$位的浮点数；
   • 处理器的浮点部件一个时钟周期可以吞吐$N_{VFMA}$个VFMA指令，每个VFMA指令的延迟是$L_{VFMA}$，这个值是两个连续的VFMA指令之间的最小周期数；
   • Li表示是Level-i级cache，$C_{Li}$是cache行的大小，$W_{Li}$是cache分组参数，$N_{Li}$是cache组数，$S_{Li}$是cache的总大小，有等式$S_{Li}=N_{Li}{C}_{Li}{W}_{Li}$,所有级别的cache的替换策略都是LRU；
   • 假设cache与寄存器之间的带宽无限大，不是影响计算的瓶颈
2. $m_r$和$n_r$的确定
   • 主要思考最底层循环的优化策略：首先是需要足够大，防止在重复更新$C_r$的循环中由于指令延迟与依赖在浮点管道中引入时钟挂起延迟，同时要尽量小，来让更多的$A_r$和$B_r$驻留于内存中；
   • 通过延迟确定$m_r$和$n_r$的下限值，为了不让计算部件发生等待，则每次更新中必须要至少$N_{VEC}\times L_{VFMA}\times N_{VFMA}$个数值计算，则有： $$m_r{n}_r\ge N_{VEC}{L}_{VFMA}{N}_{VFMA}$$
   • 尽量让寄存器中剩余空间大，根据[[@Goto2008]]的研究，$m_r$和$n_r$应该相等，但是这很难取到，因此使用如下原则：$$m_r=\lceil \frac{\sqrt{N_{VEC}{L}_{VFMA}{N}_{VFMA}}}{N_{VEC}}\rceil N_{VEC}$$$$n_r=\lceil \frac{\sqrt{N_{VEC}{L}_{VFMA}{N}_{VFMA}}}{m_r}\rceil$$
3. $k_c$的确定
   主要思路：由于重用的块被映射到内存层次结构的不同层，这意味着缓存的大小施加了 $k_c$、$m_c$和 $n_c$ 的自然上限。此外，由于重用的数据理想情况下应该在迭代之间的缓存中保留，因此缓存替换策略和缓存组织对$k_c$、$m_c$ 和 $n_c$ 的值施加了进一步的限制。侧重于确认 $k_c$ 的解析值，$m_c$ 和 $n_c$ 的值可以以类似的方式导出。

   1. 将$B_r$驻留在L1 cache中。
      $B_r$ 在循环5的每次迭代中被重用。因此，$B_r$ 应该驻留在 L1 缓存中，执行循环 5。由于 L1 cache实现了 LRU 替换策略，传统方法是选择矩阵不同部分的大小，以便循环两个连续迭代所需的数据保留在缓存中。这意味着cache应该保留一个$B_r$、两个$A_c$和两个$C_c$(也就是传统的GotoBLAS分析方法） 。这种方法的问题是保持两个$A_c$和两个$C_c$，这减少L1 cache中$B_r$的大小，使$k_c$更小，可能重用的数据量也减少了。此外，$k_c$ 是 Loop 6 的迭代次数，减少这个值也意味着在微内核中摊销具有足够计算访存操作的机会更少。
      由此，可以得出两种限制：

      • 在循环 5 的每次迭代中，$A_c$ 中的微面板 $A_r$ 和重用的微面板 $B_r$ 被完全访问一次。这意味着首先加载的微面板将包含最近使用的元素。因此，由于 LRU 缓存替换策略，只要 $B_r$ 在 Ar 之后加载，来自 $A_r$ 的条目将在 $B_r$ 的条目之前从缓存中驱逐。
      • 由于每个微面板$A_r$在每次迭代中只使用一次，因此用将在当前迭代中使用的新$A_r$的相应条目覆盖在前一个迭代中使用的微面板$A_r$的条目是有利的。这样做可能会允许更大的$B_r$放入缓存，从而增加被重用的数据量。

   2. 将前一个$A_r$驱逐出cache
      从前面的观察中可以看出，在缓存中保持大型微面板$B_r$的策略是将在循环5的前一次迭代中加载的旧微面板$A_r$从缓存中驱逐，方法是将其条目替换为将在当前循环迭代中使用的新微面板$A_r$的条目。
      为了实现这一目标就需要保证$A_c$的连续微面板的相应元素必须位于$N_{L1}\times C_{L1}$字节间隔的整数倍的内存中,由此可得：$$m_r{k}_c{S}_{DATA}=C_{A_r}{N}_{L1}{C}_{L1}$$
      因此：$$k_c=\frac{C_{A_r}{N}_{L1}{C}_{L1}}{m_r{S}_{DATA}}$$
      这一值可以确保新读取的$A_c$微面板$A_r$映射到与现有$A_c$微面板$A_r$相同的组。

   3. 确定$C_{A_r}$
      $C_{A_r}$实际上是$A_r$所占的cache行个数，因此在$W_{L1}$路组相联下有：$$C_{A_r}+C_{B_r}\le W_{L1}$$
      在实际应用中，需要为不打包的$C_r$预留至少一个cache行，因此：$$C_{A_r}+C_{B_r}\le W_{L1}-1$$
      因为$B_c$这一分块需要打包进连续内存，所以，$C_{B_r}$满足：$$C_{B_r}=\lceil \frac{n_r{k}_c{S}_{DATA}}{N_{L1}{C}_{L1}}\rceil =\lceil \frac{n_r}{m_r}{C}_{A_r}\rceil$$
      带入上式中有：$$C_{A_r}\ge \lfloor \frac{W_{L1}}{1+\frac{n_r}{m_r}}\rfloor$$
      这一估计需要满足的前提条件是：

      • Cache中的每个组都必须均匀地填充$A_r$的微面板;也就是说，每个组中包含来自$A_r$的元素的cache行数在所有集合中是相同的。这一假设确保了不同微面板的相应元素将被分配到相同的cache组。
      • 在每组缓存中有一条cache行是为Cr元素保留的，这样加载它们就不会驱逐已经驻留在留在cache中的微面板$B_r$。

6 实验结果(Evaluation)

在固定的矩阵规模和$m_c$下观察$k_c$对矩阵乘法的影响：绿线表示分析和经验确定的值$m_c$值，蓝线表示小于经验值的$m_c$，红线表示大于经验值的$m_c$。

在固定的矩阵规模和$k_c$值下，观测mc值对矩阵乘法效率的影响：蓝圈表示的是经验最优$m_c$，绿圈表示的是分析得出的最优$m_c$。