闫昊(1996-), 男, 硕士, 主要研究领域为高性能数值计算
刘芳芳(1982-), 女, 正高级工程师, CCF专业会员, 主要研究领域为高性能扩展数学库, 超级计算机评测软件
马文静(1981-), 女, 副研究员, CCF专业会员, 主要研究领域为高性能计算, 代码生成与优化
陈道琨(1994-), 男, 博士, 主要研究领域为高性能计算, 异构并行, 稀疏矩阵相关的算法研究
稠密矩阵乘法(GEMM)是很多科学与工程计算应用中大量使用的函数, 也是很多代数函数库中的基础函数, 其性能高低对整个应用往往有决定性的影响. 另外, 因其计算密集的特点, 矩阵乘法效率往往也是体现硬件平台性能的重要指标. 针对国产申威1621处理器, 对稠密矩阵乘法进行了系统性地优化. 基于对各部分开销的分析, 以及对体系结构特点与指令集的充分利用, 对DGEMM函数从循环与分块方案, 打包方式, 核心计算函数实现, 数据预取等方面进行了深入优化. 此外, 开发了代码生成器, 为不同的输入参数生成不同版本的汇编代码和C语言代码, 配合自动调优脚本, 选取最佳参数. 经过优化和调优, 单线程DGEMM性能达到了单核浮点峰值性能的85%, 16线程DGEMM性能达到16核浮点峰值性能的80%. 对DGEMM函数的优化不仅提高了申威1621平台BLAS函数库性能, 也为国产申威系列多核处理器上稠密数据计算优化提供了重要参考.
General matrix multiply (GEMM) is one of the most used functions in scientific and engineering computation, and it is also the base function of many linear algebra libraries. Its performance usually has essential influence on the whole application. Besides, because of its intensity in computation, its efficiency is often considered as an important metric of the hardware platform. This study conducts systematic optimization to dense GEMM on the domestic SW1621 processor. Based on analysis of the baseline code and profiling of various overhead, as well as utilization of the architectural features and instruction set, optimization for DGEMM is carefully designed and performed, including blocking scheme, packing mechanism, kernel function implementation, data prefetch, etc. Besides, a code generator is developed, which can generate different assembly and C code according to the input parameters. Using the code generator, together with auto-tuning scripts, it is able to find the optimal values for the tunable parameters. After applying the optimizations and tuning, the proposed single thread DGEMM achieved 85% of the peak performance of a single core, and 80% of the performance of the entire chip of 16 cores. The optimization to DGEMM not only improves the performance of BLAS on SW1621, but also provides an important reference for optimizing dense data computation on SW series multi-core machines.
矩阵乘法(GEMM)是使用最广泛的基本矩阵操作函数之一. 它不仅在很多计算应用中被频繁调用, 而且是基础代数库BLAS中三级函数的基础. 很多其他BLAS三级函数的性能提高也依赖于GEMM的优化. 因此, 对很多科学和工程应用来说, 提高稠密矩阵乘法的性能, 对加快整体计算速度有重要作用. 而GEMM函数的优化, 又严重依赖于硬件平台, 函数实现必须充分利用处理器特性才能发挥出计算平台的计算能力. 因此, 各处理器制造厂商都在致力于提供针对自己生产的处理器进行过优化的BLAS函数库, 尤其是进行了深度优化的GEMM函数. 例如, Intel公司提供优化的数学库MKL[
申威1621处理器是我国自主研发的多核高性能处理器, 采用自主sw指令集. 由于申威处理器的架构, 存储层次及所使用的指令集与国外商用CPU有很大不同, 因此在申威处理器上无法使用专为国外商用处理器优化的BLAS库. 而对开源库进行简单适配的版本性能又远远低于处理器所提供的性能. 基于此种状况, 研制面向申威CPU的BLAS库成为一个亟待实现的任务. 刘昊等人对早期申威处理器(申威1600)上的GEMM进行了优化[
我们在申威1621平台上对GEMM从访问优化, 指令排布以及减少不必要计算和访存等方面, 进行了深度优化, 并开发了子函数代码生成器, 为不同规模的矩阵生成优化的代码并进行参数空间搜索. 相比于GotoBLAS, 单线程DGEMM实现7.39倍加速, 达到了系统峰值性能的85%, 16线程性能达到了系统峰值的80%.
稠密矩阵乘法GEMM的计算公式是:
公式(1)中, 矩阵
算法
1. for
2. for
3.
4. for l = 0 :
5.
6. end for
7. end for
8. end for
GEMM的浮点运算量是2
目前GEMM的优化主要有3个方面: (1)针对特定的平台, 在充分考虑硬件架构特征的基础上设计分块方案和单线程/多线程算法; (2)对核心代码进行精细手工优化, 例如, GotoBLAS就提供了几种不同CPU上的汇编核心代码; (3)基于自动调优(auto-tuning), 通过在不同配置下运行程序收集性能表现数据, 自动选择算法和参数. 该方法的代表是ATLAS数学库. 在库函数的开发当中, 往往需要将3种优化思路结合使用, 以期在提高开发效率的同时, 充分挖掘平台性能潜力.
各经典BLAS库中的GEMM函数及各方研究人员开发的GEMM函数中, 用到了多种体系结构相关的优化. 在广泛使用的开源软件GotoBLAS中, Goto等人考虑到硬件上缓存容量和TLB表项容量的限制, 给出了一系列分块的实际原则以及分块参数选择方法, 具有普适的指导性[
ATLAS数学库采用了自动调优的方法, 在不同配置下运行程序, 收集性能数据, 选择性能表现最优的参数组合作为调优结果[
国产申威1621平台[
本节介绍单线程DGEMM在申威1621平台上的实现方案及优化方法. 我们将文献[
本节介绍文献[
为有效利用现代处理器上的分层缓存结构, GEMM采用分块算法, 每个维度上的循环迭代步长是此维度的分块大小, 最内层循环计算一个分块的矩阵乘法. 在循环顺序上, 我们采用了与GotoBLAS2相同的N-K-M的循环顺序. 鉴于传统BLAS库中矩阵使用列优先存储, 矩阵
算法2展示了单线程DGEMM的分块算法. 矩阵
算法
1.
2. for
3. for
4. pack §
5. for
6. pack §
7. 调用计算函数KERNEL(), 计算
8. end for
9. end for
10. end for
函数KERNEL()完成
函数KERNEL图示
算法
1. for
2. for
3. for
4. load
5.
6. end for
7. load
8.
9. store
10. end for
11. end for
循环最内层共有
为了提高数据访问的效率, KERNEL中插入了针对3个矩阵的预取操作. 由于cache line长度为128字节, 每次预取操作能够读取16个双精度浮点数. 因此, 对
在现代计算机体系结构上实现高性能矩阵乘法, 分块方法的优化是必不可少的[
首先, 算法2中的分块算法无法兼顾打包开销的减少和KERNEL性能的优化. 一方面, 较大分块意味着更多的缓存缺失率, 影响KERNEL访存效率; 另一方面, 如果分块较小, 又会产生大量的重复打包操作, 增加打包操作引起的访存开销. 虽然通过调节分块大小能在一定程度上改善程序性能, 但是要想在打包和KERNEL两个环节都有性能提升, 需要对循环结构进行进一步优化以减少这两种开销. 其次, pack函数的运行时间是不可忽略的, 故对此类内存拷贝函数, 也需要进行精细的分析和优化. 此外, KERNEL的实现中也有一些可改进之处.
基于上述分析, 我们提出了兼顾减少打包开销和提高KERNEL效率的二层分块方案, 并对打包函数及KERNEL函数的实现都进行了优化. 基于优化之后的程序, 我们进一步对KERNEL的实现在减少访存量和提高数据预取质量方面进行了优化. 下面详细介绍各项优化技术.
算法2中, 函数packA()执行次数是
当输入规模固定时, 函数packB()用时只取决于机器带宽; 函数packA()用时反比于分块参数
算法
1.
2. for
3. for
4. pack §
5. for
6. pack §
7. for
8. for
9. for
10. 调用计算函数KERNEL(), 计算
11. end for
12. end for
13. end for
14. end for
15. end for
16. end for
如前所述, 在基础版本中, 对
通过简单的计时统计, 我们发现, 打包操作的时间是不可忽略的, 因此, 对packA和packB操作, 我们对其实现进行了深入分析和优化.
由于KERNEL函数处理的是
改进的函数packA示意图
类似的, 函数packB()也要将
改进的函数packB示意图
在本节及第2.2.5节中, 我们介绍对计算KERNEL的优化, 包括减少对
在KERNEL计算中, K循环结束时, 需要将矩阵
在K方向第1次迭代中, 缓冲区
算法
1. if
2. KERNEL_INIT();
3. else if
4. KERNEL_UPDATE();
5. else
6. KERNEL_GENERIC ();
7. end if
原始版本的KERNEL虽然对寄存器使用等进行了精细调优[
首先, 原有KERNEL中对
其次, 在原始程序中, 由于对
矩阵
对多线程矩阵乘法, 我们采用了与GotoBLAS相同的并行方案[
多线程并行矩阵乘法任务划分示意图
单线程GEMM的优化技术均适用于多线程DGEMM. 区别在于, 多线程GEMM在部署分块方案时, N方向需要在二层分块的基础上, 增加一层线程级分块.
在实现过程中, 我们开发了KERNEL和打包函数的代码生成器, 使用Python语言编写代码生成脚本. 通过调整参数, 该脚本可以生成不同版本的KERNEL函数和打包函数代码.
KERNEL函数的代码生成器将预取步长, 循环展开次数等作为参数传递, 可以生成使用不同优化策略和参数的KERNEL汇编代码. 对
为了选取最优参数, 我们编写了自动调优脚本, 对参数空间进行遍历搜索. 需要调优的参数包括第1层分块和第2层分块的分块大小(
经过实验遍历搜索后得出性能最佳的分块和预取参数组合如
单线程DGEMM循环分块参数表
维度 | 第1层数据分块 | 第2层数据分块 |
1024 | ||
4096 | 256 | |
256 | 128 |
DGEMM KERNEL内部预取步长参数表
矩阵 | L1预取步长 | L2预取步长 |
96 | 160 | |
64 | 160 |
对多线程GEMM, 我们采用类似的方法, 对分块参数进行遍历搜索, 得出的多线程GEMM最优分块参数如
多线程DGEMM关键分块参数表
维度 | 第1层数据分块 | 第2层数据分块 |
512 | ||
4096 | 256 | |
128 | 128 |
本节介绍我们在申威1621处理器上对双精度浮点实数矩阵乘法进行的测试, 包括单线程DGEMM使用各项优化技术后的性能, 以及优化的多线程DGEMM程序性能.
实验采用的申威1621处理器有16个计算核心, 主频1.6 GHz, 内存128 GB. 单核双精度浮点运算峰值25.6 Gflops, 16核双精度浮点运算峰值409.6 Gflops. 三级缓存的大小分别是32 KB, 512 KB, 512 MB.
对于单线程DGEMM, 我们选择了7个版本进行了测试. 第1个版本Gemm_Goto使用稍加改动的GotoBLAS代码(GotoBLAS原始代码在申威处理器上不能直接运行). Gemm_baseline采用文献[
DGEMM各版本说明
版本号 | 算法描述 |
Gemm_Goto | Goto版本 |
Gemm_baseline | 文献[ |
Gemm_opt1 | 使用二层分块 |
Gemm_opt2 | 对 |
Gemm_opt3 | 改进的packA(), packB() |
Gemm_opt4 | 对K方向上不同阶段的计算使用不同KERNEL函数 |
Gemm_opt5 | 针对预取操作的指令级优化 |
本节我们将展示在使用了上文所述的优化方法后, 单线程和多线程GEMM的性能提升.
单线程DGEMM实验结果见
1621处理器上单线程DGEMM性能测试结果
多线程DGEMM实验在两组数据规模上进行, 都采用M, N, K相同的方阵, 大小为32768, 分别在1, 2, 4, 8, 16个线程上进行测试. 结果见
多线程DGEMM扩展性测试结果(输入规模32768)
核数 | 性能 (Gflops) | 加速比 | 扩展性效率 | CPU效率 (%) |
1 | 21.84 | 1 | 1 | 85.3 |
2 | 43.3 | 1.98 | 0.99 | 84.6 |
4 | 86.02 | 3.94 | 0.98 | 84.0 |
8 | 169.89 | 7.78 | 0.97 | 83.0 |
16 | 329.86 | 15.1 | 0.94 | 80.5 |
多线程DGEMM中计算每个分块用时统计(s)
线程数 | 同步用时 | 计算函数用时 |
2 | 0.040861 | 50.4315 |
4 | 0.039397 | 25.34561 |
8 | 0.105054 | 12.76301 |
16 | 0.255258 | 6.468235 |
本文面向国产申威多核平台对稠密矩阵乘DGEMM单线程算法和多线程版的并行算法进行了重构与优化. 通过改变分块方案和打包方式等优化手段, 使程序运行最大化地减少额外开销, 同时又充分利用了数据局部性. 对KERNEL函数, 我们对函数在K方向迭代的不同阶段所读写的数据进行分析, 为每个阶段定制KERNEL, 并对迭代中每个部分进行精准的数据预取. 使用我们的优化方法之后, 实验数据显示, 单线程版性能达到浮点峰值性能的85%. 多线程版性能达到浮点峰值性能的80%. 这些优化技术对于典型分层存储架构和现代微结构上的应用优化都有参考价值. 此外, 我们开发了代码生成器和自动调优脚本, 为优化参数搜索和高效调优提供了实用工具. 在未来工作中, 我们准备将对GEMM函数的优化应用到其他BLAS函数中, 并对小规模矩阵乘法优化进行研究. 另外, 我们将对各种优化方法进行提炼, 进一步完善代码生成器, 增加可移植性, 形成较完整的优化工具链.
®-optimized math library for numerical computing. 2021. https://software.intel.com/content/www/us/en/develop/tools/oneapi/components/onemkl.html]]>
https://developer.amd.com/amd-aocl/]]>
https://rocmdocs.amd.com/en/latest/ROCm_Tools/rocblas.html]]>
Goto K, Van De Geijn R. High-performance implementation of the level-3 BLAS. ACM Transactions on Mathematical Software, 2008, 35(1): 4. [doi: 10.1145/1377603.1377607]
http://www.openblas.net/]]>
Abdelfattah A, Keyes D, Ltaief H. KBLAS: An optimized library for dense matrix-vector multiplication on GPU accelerators. ACM Trans. on Mathematical Software, 2016, 42(3): Article 18.
刘昊, 刘芳芳, 张鹏, 杨超, 蒋丽娟. 基于申威1600的3级BLAS GEMM函数优化. 计算机系统应用, 2016, 25(12): 234–239. [DOI: 10.15888/j.cnki.csa.005456]
Liu H, Liu FF, Zhang P, Yang C, Jiang LJ. Optimization of BLAS level 3 functions on SW1600. Computer Systems & Applications, 2016, 25(12): 234–239 (in Chinese with English abstract). [doi: 10.15888/j.cnki.csa.005456]
解庆春, 张云泉, 李焱, 逄仁波, 吴再龙, 鲁永泉, 高鹏东. 基于神威蓝光处理器的向量数学软件包. 软件学报, 2014, 25(S2): 70–79.
Xie QC, Zhang YQ, Li Y, Pang RB, Wu ZL, Lu YQ, Gao PD. Package of the vector math library based on the Sunway processor. Journal of Software, 2014, 25(S2): 70–79 (in Chinese with English abstract).
Goto K, van de Geijn R. Anatomy of high-performance matrix multiplication. ACM Transactions on Mathematical Software, 2008, 34(3): 12. [doi: 10.1145/1356052.1356053]
Lim R, Lee Y, Kim R, Choi J. An implementation of matrix–matrix multiplication on the Intel KNL processor with AVX-512. Cluster Computing, 2018, 21(4): 1785–1795. [doi: 10.1007/s10586-018-2810-y]
Whaley RC, Petitet A. Minimizing development and maintenance costs in supporting persistently optimized BLAS. Software: Practice and Experience, 2005, 35(2): 101–121.
http://www.swcpu.cn/show-190-254-1.html]]>
http://www.swcpu.cn/show-190-254-1.html]]>
Masliah I, Abdelfattah A, Haidar A, Tomov S, Baboulin M, Falcou J, Dongarra J. Algorithms and optimization techniques for high-performance matrix-matrix multiplications of very small matrices. Parallel Computing, 2019, 81: 1–21. [doi: 10.1016/j.parco.2018.10.003]