基于X-DSP 的GEMM 算法实现∗

王华龙 陈小文

（国防科技大学计算机学院 长沙 410073）

1 引言

随着超大规模集成电路技术的发展，处理器的体系结构也由于追求高性能，低功耗发生了重大变化。在高性能计算应用中，由于能耗的因素，促进了多核、向量处理、可编程可配置、可重构等体系结构的发展［1］。数字信号处理器（DSP）一直是嵌入式系统的核心组成部分，主要优势在于高功耗效率，价格低廉及易编程。但是多核环境下并行访问存储器的性能依然是多核系统设计和应用开发的难点之一。

为了发挥不同体系结构处理器的性能，通常需要提供高效的标准函数库。BLAS算法库是大规模科学计算最常用的核心数学算法库之一，而BALS库中最常用的核心算法是矩阵乘法（General Matrix-Matrix Multiplication，GEMM），其在高性能计算中扮演着重要角色，是一种同时具有计算密集和访存密集的运算［5］，对处理器的运算能力、访存带宽及延迟的要求非常高，甚至占用高性能标准测试程序HPL 主要运算量的90%以上。吴少刚［7］等在基于Intel 微处理器的PC 作为节点机的类Beowulf 机群系统上使用分支优化、访存优化、译码优化以及执行优化技术，将矩阵乘运算单机性能较原始BLAS 提升5 倍。Sohl 等［9］面向一种异构并行DSP体系结构研究了大规模矩阵乘法，研究表明针对处理器的并行结构和存储系统优化设计的矩阵乘法映射能够有效地隐藏数据访问开销。Gunnel［11］在基于Cache 的多级存储的体系结构下，提出了分层计算的方法实现BLAS 中的矩阵乘算法，这种方法降低了存储层次间搬运数据的平均开销。Goto［12］等发现，二级Cache 到寄存器的带宽足以满足核心循环运算速度的要求，对于经典矩阵乘法C=AB+C，他提出把矩阵A 的分块存放在二级Cache 中，B矩阵的分块存放在一级Cache 中，这样就增大了A矩阵分块的大小，同时也就提升了B 矩阵的运算访存比。

本文主要在X-DSP 平台实现GEMM 算法设计。X-DSP是一款高性能多核DSP，每个内核均拥有一级程序Cache（L1P）和一级数据Cache（L1D）各32KB，二级统一存储器512KB 的L2，并能够配成Cache（默认）、SRAM，或是两者的组合。SM（Shared Memory）负责处理系统中所有master对SM和DDR3 的访问请求。首先针对X-DSP 的硬件资源和体系结构，选择最终的分块方式以及数据搬移方案，实现单核GEMM 的设计；然后再对存储空间进行合理划分，采用分块的方式对数据进行并行分解，使块与块之间的计算相互独立，结合X-DSP 多核通信及同步机制将并行的分块算法映射至多个核中，完成多核GEMM 设计；最后通过对单核与多核GEMM的性能测试，实现高性能GEMM的设计。

2 单核GEMM的设计与实现

GOTO-BLAS是目前使用最广泛且性能最佳的BLAS 实现之一。从早期版本开始，它已经成功移植到许多可扩展的架构中［13］。其核心就是分块，将GEMM 分解为不同的层，每个层映射到存储结构中不同层次。同时分块计算也会隐藏数据搬移的时间，从而接近最佳性能。

2.1 GEMM分块方案

一般情况下，矩阵乘法可以分为小矩阵乘法和矩阵向量乘法［14］（矩阵向量乘法一般性能比较差），而小矩阵乘法又可以分为GEPP、GEMP 和GEPM，可再进一步分解为GEBP、GEPB 和GEPDOT，如图1 所示。Watts D J［15］等比较了这些分块方法以及不分块方法的性能，发现分块的方法性能远高于不分块的方法，进一步证明了分块方法的有效。

图1 GEMM分块实现方案

本文主要研究基于GEMM-GEPP-GEBP 进行分块计算，单核计算C=A×B 时，首先外层循环沿着A 矩阵的列方向进行分块，采用GEPP 分块方式计算GEMM，最后内层循环沿着A 矩阵的行方向进行分块，采用GEBP 分块方式计算GEPP。单核运算中基于SM 的分块方案如图2，灰色分块便是GEBP的运算部分，虚线箭头表示GEBP 逐渐完成GEPP的分块，真正解决的是单核的分块运算和存储层次间的映射，数据的放置、移动及计算和数据移动的重叠。

图2 基于GEBP的GEPP

2.2 映射到X-DSP核心架构

矩阵A 的分块如图3，选择合适大小分块放入L2 ，基于Goto-BLAS 诸多文献的研究结果［16～17］和X-DSP 的实际情况考虑，用L2 Cache 的一半（256KB）存放该分块。矩阵B 被分成小块，选择合适大小放入L1 中，通常选择一半L1（16KB）作为缓存，实际中，L1 Cache 还要用来缓存其它的数据（如程序中的其它变量、循环的开销）。已经保存在L2中的分块被进一步细分为Mkernel*KSM的较小子块。矩阵B 的子块也被分为KSM*Nkernel。通过驻留在L2中的A 子块与L1 中的B 子块相乘得到Mkernel*Nkernel的小块矩阵C。通常情况下，为了减少L1 Cache 的冲突，B 矩阵分块KSM*Nkernel容量应该小于L1 Cache容量的一半。Mkernel与Nkernel的选择取决于内核中可用的寄存器数量。Goto BLAS 针对以往处理器研究的大量文献［18～21］的实验结果都表明Mkernel*Nkernel分块大致为方阵时效果更好，不过在实际的优化设计中，对分块参数Mkernel与Nkernel选择还会有其他的影响因素。内核的内部循环使用内部向量指令进行优化，以充分利用X-DSP的内核架构。

2.3 数据搬移方案

数据搬移要求尽可能地将数据移动与计算重叠，从而需要充分利用X-DSP 的DMA 功能。如图3，说明了GEPP 即A 矩阵的一列乘以B 矩阵的一行，A、B矩阵最初驻留在DDR存储器中并以列顺序存储。我们在X-DSP 的存储器中创建暂存缓冲区，用来存储A和B矩阵的子块。

SM 有三个数据buffer 用于接收由DMA 从DDR 中搬移的数据，一个buffer 接收A 矩阵子块m_SM * k_SM，另外两个buffer 接收B 矩阵子块k_SM*n，双buffer 的设计可以将GEBP 数据搬移和计算重叠。

L2 中也有两个buffer 分别用于接收C 矩阵子块和A矩阵子块。

L1中缓存区主要存放B矩阵子块。

图3 基于GEPP-GEBP数据搬移

最终目的是将驻留在L2 存储器中的B 矩阵子块与位于L1 中的A 矩阵的子块相乘。首先，DMA搬运A 矩阵和B 矩阵子块到SM 中相应位置，传输完成后，A矩阵子块最终存储在L1缓冲区中，SM 中缓冲区变空准备接收下一个A 矩阵子块；B 矩阵子块最终存储在L2 缓冲区。此时，启动DMA 搬移下一次乘法所需的A 矩阵子块与B 矩阵子块到SM。这次搬移与A 矩阵与B 矩阵的GEBP 同时发生。在封装B Panel 子矩阵的同时，需要更新的Cj的子矩阵由DMA 搬移至相应的L2 buffer 中，待所有的数据准备好后，核心循环的运算开始运行（图3 第四栏）。

3 多核GEMM的设计与实现

在多核的环境下，追求高性能的GEMM依然非常重要。分析多核的存储体系结构如何进行数据的并行分解，以便可以更好地使用多核进行并行的相互独立的运算处理，从而达到最佳的性能。根据文献［22～23］对多核DSP 的并行编程研究，同时考虑X-DSP 存储体系结构，以便可以更好地使用多核进行并行的相互独立的运算处理，并且要使得多核间的负载均衡，最终实现kernel设计。

3.1 多核并行性分析

多核并行性的重要性在于可以让不同的核独立的运算每个分块矩阵，而且每一个独立运算都可以按照单核的计算方式进行，这样就可以大大提高执行速度。

X-DSP 的多核体系结构采用共享SM 和DDR的存储结构，分块矩阵乘法的不同子块间的计算不需要考虑数据交换，因此可以适合数据并行。考虑一个经典的矩阵乘法C=A*B+C，假设共有P 个核，标号分别为P0，P1……PP-1。本文将A 矩阵按行进行划分为P块，这些块依次记为A0，A1……AP-1，矩阵C 同时按行进行相应的划分，每块子矩阵Ci和Ai含有连续的mi=m/P 行，Ci=AiB+Ci便可以由相应的核进行独立的计算。然后再将共享矩阵B 也按照列进行相应的划分为P 块，依次记为B0，B1……BP-1，C矩阵对应的按照列划分为P 块，每块子矩阵Cj和Bj含有连续的nj=n/p 列，Cj=ABj+Ci便可以由相应的核进行独立的计算。

3.2 多核GEMM的设计与实现

从图4 中可以看出多核并行运算时矩阵B 是每个核都需要访问的，根据X-DSP 的体系结构，矩阵B 的数据可在SM 中共享，则矩阵乘法计算只需要搬移A，C 矩阵即可。但是由于SM 大小的限制，矩阵B 必须分块进入SM，其最大子块在SM 中直至所有的核不再访问为止。为防止不同的核运行速度不同而把SM 中的数据冲掉，再访问B 矩阵的下一个分块之前需要做一次同步，同步之后所有的核将全速运行，而且同步之后各个核将互不影响彼此的运行。

图4 多核GEMM设计

本次所研究的X-DSP 共有10 组DMA 通道可以同时传输数据，因此其不同核可以使用不同的DMA 资源，数据的搬移与计算在不同核中可以满足重叠条件的。SM 中Buffer 中数据到位后，每个核只需要包装搬移自己所需要处理的数据且更新的是不同的C 矩阵子块，计算完成后由DMA 将结果搬回到DDR 中即可。为了避免多个核之间重复的数据搬移，设计中对SM 做了存储段的数据分配：设计两个buffer 接收KSM*Ncore的B Panel，此外再设计8 个buffer 接收每个核所需的KSM*SM。L2 和L1是每个核私有的，无需作特殊的处理。

4 性能测试与分析

4.1 单核性能

图5 中显示了基于GEBP 的GEPP 在X-DSP中单核上的性能曲线，从曲线中可以看出本次设计所达到的单精度浮点数的最高性能为8.49GFLOPS，性能损失的主要原因大部分与体系结构相关，取决于X-DSP 中的缓存机制。通过性能测试及时间的统计，分析其原因主要在于：

1）搬移新的A 矩阵子块时，每次调用kernel 时cache未命中；

2）循环开销受限于于L1 和L2 的大小；

3）跨不同层和数据打包，移动的开销。

4.2 多核性能

图6 显 示 了 基 于GEPP-GEBP 的GEMM 在X-DSP的8个核上运行实现的性能曲线。8核上并行的BLAS 实现，单精度GEMM 的峰值性能为52.8GFLOPS。分析影响性能的主要因素在于核间的同步造成了并行时的串行等待。此外，多核并行对于共享存储器进行划分时，其实限制了运算分块的大小，虽然是多核同时计算数据，但是分块的大小会影响核心循环的性能，从而可能会影响最终的性能。

图6 X-DSP多核GEMM性能

5 结语

本文主要针对X-DSP 的体系结构特点研究GEMM 的运算性能。主要设计实现基于GEBP 的GEPP 在X-DSP 单核上的性能，对X-DSP 的存储空间进行了合理划分，得到单核的单精度浮点数的最佳性能达8.49GFLOPS。对于GEBP 的GEPP 在X-DSP多核上的性能，采用分块的方式对数据进行并行分解，使块与块之间的计算相互独立，结合X-DSP 多核通信及同步机制将并行的分块算法映射至多个核中，同时对X-DSP 的硬件传输机制和共享存储进行了优化设计。经过测试，X-DSP的多核最佳性能达52.8GFLOPS。