2021, Vol. 32
BLAS(basic linear algebra subprograms)是基本线性代数子程序的缩写, 是目前应用广泛的核心线性代数数学库.随着不断的发展完善, BLAS在高性能计算、科学和工程领域都得到了广泛的应用.由于不同厂商不同体系结构下的CPU差异较大, 通用BLAS函数库不能达到很好的性能.因此在大规模矩阵运算时, 针对不同CPU的架构特点来优化BLAS, 可以充分发挥处理器的计算性能, 从而使计算效率大大提升.目前主流的BLAS分为专用和通用两种.专用型BLAS是由特定的CPU公司开发, 针对特定平台芯片进行代码优化, 极大地提高了性能, Intel的MKL和AMD的ACML就是此种类型.通用BLAS由非盈利机构开发, 不同平台都可以对开源代码进行优化, 最具代表性的是ATLAS[1]、GotoBLAS[2, 3]和OpenBLAS[4, 5].
BLIS(BLAS-like library instantiation software)数学库是美国Texas大学超算中心高性能组开发的开源架构, 是在GotoBLAS的基础上进行了改写和优化, 具有可移植性, 易用性和模块化设计易于开发的特点, 支持混合类型矩阵存储和多线程[6].文献[7]提出了BLIS架构, 并做了全面的介绍; 文献[8]研究了如何使用BLIS在通用、低功耗和多核架构中实现高效的leve-3级BLAS函数; 文献[9]系统的探讨了矩阵乘法算法在BLIS五层循环中并行化的机会.BLIS具有的这些特性使其更容易针对具体平台进行BLAS优化开发.BLAS包含三级函数, 分别为第1级向量与向量计算, 第2级向量与矩阵计算, 第3级矩阵与矩阵的计算, 其中第3级函数DGEMM计算量最大, 可以达到理论计算量的95%以上[10], 因此三级函数的优化是性能提升的关键.
高性能计算系统可以分为同构或异构两类, 同构系统是由大量CPU构建, 计算和数据处理都在CPU上进行; 异构系统是CPU和加速部件组成, 加速部件一般由Cell(cell broadband engine architecture), FPGA(Field programmable gate array)和GPU(graphic processing unit) 等为代表, 集成大量的浮点计算单元的同时舍弃了一些CPU上复杂的控制单元.异构系统中CPU负责任务调度和简单计算, 主要计算在加速部件上进行, 计算能力更强.国际上对高性能系统的评测标准程序是HPL(high performance Linpack), 由美国田纳西大学教授Dongarra设计提出[11], 是在高性能计算系统中使用高斯消元法求解一个随机的N元一次线性方程组问题, 用以评价高性能计算机的浮点性能.
异构HPL浮点计算在CPU端是通过BLAS函数完成, BLAS库的性能会影响整体系统性能.本文异构系统计算单元由CPU+协处理器(co-processor)组成(如图 1所示), 其中协处理器有64个计算单元, 共有4 096个流处理器, 16GB显存, 主频最高可达1.8GHz.CPU是国产X86架构的海光Dhyana处理器, 每个DIE通过PCI-E连接一个加速部件, CPU计算核心之间和协处理器通过高速互联进行数据交换和通信.协处理器的计算能力一般都比通用CPU强大很多, 往往造成计算负载严重不均衡, 损失性能.因此在CPU端利用更高效的数学库, 能够充分地发挥CPU性能, 使负载更均衡, 是提高计算能力的可行方式之一.文献[12-15]分别针对ARM架构处理器以及异构平台的GEMM优化做了深入研究.文献[16, 17]研究了在申威众核处理器上1, 2和3级BLAS的优化方法, 文献[18]对异构多核平台上的数学库寄存器分配优化方法进行了研究.
|
Fig. 1 Heterogeneous computing node architecture 图 1 异构计算节点架构 |
本文在BLIS算法库基础上, 分析异构HPL算法中调用的各级BLAS函数, 针对体系结构特点, 利用访存优化、指令集优化和多线程并行等技术手段形成了HygonBLIS, 提高了异构系统中CPU端的计算效率, 后文详细分析了优化策略和方法.本文第1节首先对异构HPL的BLIS库进行分析, 特别是异构环境的HPL特点以及各个数学库函数调用关系.第2节~第4节对高性能BLIS库的优化策略做详细的分析和介绍, 优化HPL中涉及的各级BLAS函数, 实现高性能的HygonBLIS算法库.第5节在异构平台上做详细的性能测试与分析, 通过性能对比, 验证优化后的HygonBLIS库性能有大幅度的提高.第6节对全文进行总结和展望.
1 算法分析和CPU微体系结构 1.1 HPL算法与BLAS函数HPL算法核心计算公式如式(1)所示, 其求解过程实际是对系数矩阵[A, b]使用高斯消元法进行LU分解, 变成[A, b]=[[L, U], Y]其中L是下三角矩阵, U为上三角矩阵, 随着因式分解的进行, 方程等价于求解Ux=Y.LU分解过程中每一次迭代计算都要经过panel分解、panel广播、行交换和尾矩阵更新.HPL程序在进程通信和矩阵计算上需要MPI[19]库和BLAS或VSIPL[20]库的支持, 除基本算法不可以改变外, 可以通过配置文件调节问题规模N(矩阵大小)、进程数等测试参数, 以及使用各种优化方法来执行HPL程序, 以获取最佳的性能.当求解问题规模为N时, 浮点运算次数为2/3×N3–2×N2.因此, 只要给出问题规模N, 测得系统计算时间T, 计算系统的浮点计算能力=计算量(2/3×N3–2×N2)/计算时间(T), 测试结果以每秒浮点运算次数(FLOPS)单位.
| $ Ax = b(A \in {R^{n \times n}}, x \in {R^n}, b \in {R^n}) $ | (1) |
HPL算法过程分为LU分解和回代过程, 异构HPL环境下, LU分解过程中大量计算要在CPU和协处理器之间进行任务分配, CPU负责调度和适当的辅助计算, 而主要计算要在加速部件上进行.当问题规模N非常大时, HPL的最大计算量是更新尾矩阵的计算, 包括双精度的三角矩阵求解和稠密矩阵乘DGEMM, 因此将DGEMM放在协处理器上计算, CPU负责计算量不大的panel分解, panel广播和行交换.异构高性能计算系统优化的难点主要是在主从芯片对计算任务的划分和通信开销这两方面, 一些学者从任务静态划分和动态划分实现负载均衡[12, 21-23], 还有学者在数据重用、存储优化、软件流水和应用双缓冲等方面深入研究[24-26]达到实现更高效的通信隐藏目的.异构HPL主程序算法中任务分配调用加速部件专用编程接口进行计算, CPU端则调用BLAS库进行选主元、panel分解计算.
本文针对的异构HPL算法由中国科学院软件研究所开发[27], 其LU分解过程会将panel的一部分矩阵(N×NBMIN)分解为单位下三角L1、上三角U和L2部分, 这3部分运算在CPU上进行, 主要涉及一级BLAS函数中的DCOPY、IDMAX、DSCAL; 二级BLAS函数中的DTRSM、DTRSV、DGEMV; 三级BLAS函数中的DGEMM.各个BLAS函数作用如下:
(1) DCOPY执行数据拷贝操作, 将一个向量拷贝到另一个向量.
(2) IDMAX用于选主元操作, 在一列元素中选择最大的元素.
(3) DSCAL用于向量更新, 用一个标量乘以一个向量.
(4) DTRSM用于三角矩阵求解操作.
(5) DTRSV用于求解三角方程操作.
(6) DGEMV用于计算向量矩阵乘法.
(7) DGEMM用于每个进程上剩余矩阵数据更新.
同时panel分解的递归部分也会调用三级BLAS函数DGEMM加速panel分解过程.因此DGEMM的效率极大地影响了panel分解的效率.在异构HPL算法中, 参数NBMIN表示递归分解算法可以分解到的最小方阵大小, 随着它的增大, 一级BLAS在整体过程中的操作次数与规模不变, 二级BLAS操作次数保持不变, 而运算规模将增加, 导致panel分解非递归部分的整体时间延长.当NBMIN减小时, HPL_pdpanllN[28]函数的耗时将减少, 但是也会增加外层非递归部分对三级BLAS的调用次数.因此需要通过权衡三级与二级BLAS的调用次数与效率来调整NBMIN的大小, 获得更好的性能.另一方面, 通过DGEMV的优化可以增大NBMIN来减少三级BLAS的调用次数.针对这些函数, 本文应用多种优化策略对BLIS库中的各级BLAS函数进行了优化, 获得了性能的大幅提升, 提高了异构HPL的效率, 更好地发挥了异构系统的性能.
1.2 Dhyana处理器微体系结构Hygon Dhyana CPU是国产64位X86架构的通用服务器芯片, 采用14nm工艺, 最高主频可达3GHz.从图 1中可以简单了解CPU架构, 以CCX(core complex)微结构为基础, 单独一个CCX由4个core组成, 每两个CCX组成一个DIE, 4个DIE封装成一个Package, 共32个计算核心.其微体系结构主要特点如下:
(1) 支持所有标准X86指令集AVX、AVX2、SMEP、SSE和SSE2等.
(2) 一个CCX内部由4个核心组成一组, 每个核心独享L1和L2级缓存, 共享L3级缓存, 支持SMT技术, 每个物理核心支持同步2线程.
(3) 三级缓存结构如图 2所示, L1 cache由64K指令缓存和32K数据缓存组成, 其中L2 cache大小为512K, 共享L3 cache大小为8M.所有Cache line大小为64字节.
|
Fig. 2 Dhyana CPU cache hierarchy 图 2 Dhyana CPU缓存层次结构 |
(4) 每个cycle可以解码4条X86指令, 并存放到大小为72 entry的微操作队列, 同时可以发射6条微操作, 微操作缓存最大支持2K条指令.共有4个定点部件, 2个浮点部件和2个存储部件, 最多同时发射2条浮点指令和4条定点指令.
(5) Dhyana CPU有三级指令TLB, 两级数据TLB.访存页面大小为4KB, 映射的物理地址空间最大可以达到1GB.
2 访存优化HPL算法中DGEMM运算占据了BLAS运算中最大的计算量和访存耗时.HygonBLIS通过实现互补的两类GEMM kernels优化了DGEMM性能, 即大kernel和小kernel.大kernel用于处理矩阵A和B的规模都相对较大的情况, 而小kernel用于处理矩阵A更加“高窄”和矩阵B规模较小的情况.在DGEMM函数接口内根据矩阵规模切换两类kernel.同时引入Auto-tunning自适应调优技术对矩阵分块参数进行优化, 优化后的分块参数可以更充分地利用CPU缓存, 提高cache命中率, 从而提高DGEMM访存性能.
2.1 矩阵分块GEMM大kernel采用经典高效的GEMM分块算法[8], 如图 3左侧的6层循环伪代码所示; GEMM小kernel的算法如图 3右侧的伪代码所示.
|
Fig. 3 High performance GEMM algorithm of big and little kernels 图 3 高性能GEMM的大kernel算法和小kernel算法 |
大kernel算法描述如下:
● 最外层Loop1遍历n维度(以jc索引), 以宽度nc为单位, 对矩阵C和矩阵B进行分块, 形成列向panels.
● 第2层Loop2遍历k维度(以pc索引), 以高度kc为单位, 将矩阵A划分成列向panels, 同时将Loop1中划分出的矩阵B的列向panel进一步划分成block; Loop2中, 还会将B block(B(pc: pc+kc–1, jc+nc–1)表示)进行打包(packing).关于大kernel中的packing操作见2.2节.
● 第3层Loop3遍历m维度(以ic索引), 以高度mc为单位, 将Loop1中划分出的矩阵C的列向panel进一步划分出block, 以Cc表示(未进行packing操作), 将Loop2中划分出的矩阵A的列向panel也进一步划分成block.同样, Loop3中, 还会将A block进行打包, 以Ac表示打包后的矩阵块.
● 第4层Loop4再次遍历n维度(以jr索引), 以宽度nr为单位, 将packed block Bc和block Cc再次进一步划分成micro-panels.
● 第5层Loop5再次遍历m维度(以ir索引), 以高度mr为单位, 将packed block Ac划分成micro-panels, 将Loop4中Cc中的micro-panel进一步划分成micro-tiles.
● 最内层Loop6, 再次遍历k维度(以pr索引), 以1位单位, 计算Ac中的micro-panel与Bc中的micro-panel的点积(dot product), 并将结果累加到Cc中的micro-panel, 即所谓的“rank-1 update”.
小kernel算法描述如下:
● 前两次loops的作用在于将小矩阵B进行打包, 其目的是便于inner-loops中最内层kernel对B矩阵的连续、对齐访问, 关于小kernel中的packing B操作见2.2节.
● inner-loops中的三层循环类似于大kernel中的实现, 只是此时循环操作的上限不再是大kernel中分块后的nc、mc和kc, 而直接是整个矩阵的大小.
与大kernel算法相比, 小kernel算法有如下优点:
(1) 省去了对矩阵A的packing操作.
异构HPL中的矩阵A往往是“高窄阵”, 其规模比矩阵B大非常多, 对矩阵A进行打包操作的耗时占整个GEMM操作的比例非常高.
(2) 极大地避免了大kernel中out-loops的循环次数.
大kernel算法与小kernel算法间的动态切换, 主要依据输入矩阵的规模, 即计算量的大小.但并不存在完美的切换阈值(threshold), 因为3个矩阵A、B和C的行、列维度千变万化, 某一个特定的切换阈值无法做到全面适用于所有的使用场景, 需要在具体的应用中实际测试来找到比较合适的阈值设置.根据本文系统平台实测, 切换进入小kernel的阈值条件满足下面两者之一即可:
(1) C矩阵的面积, 即M×N≤700×700时.
(2) 维数M≤160同时维数K≤128时.
2.2 Auto-tuning自适应分块参数调优DGEMM计算量大, 计算访存比高, 除上述kernel算法层面的优化外, 对BLIS框架内矩阵分块参数进行优化可以提高cache命中率, 提高性能.矩阵乘法循环分块的大小是通过BLIS库中DGEMM的分块参数来决定的, 分块参数的选择可以影响cache命中率, 不合适的参数可能导致流水线停顿, 从而使DGEMM性能下降.因此DGEMM分块参数的选择至关重要.BLIS中DGEMM分块是通过开发人员的经验和测试选择的一组数值, 手动测试不能保证测试覆盖率, 也不能保证得到的是最优解.因此本文采用了Auto-tuning遗传算法自适应调优技术, 更高效的优化DGEMM分块参数.
遗传算法基本算法思想是从可能具有潜在解集的一个种群开始的, 种群是基因编码后形成的一组个体集合, 在这个种群基础上通过选择(selection)、交叉(crossover)和变异(mutation)操作, 优胜劣汰, 不断的进化出最优的个体.Auto-tuning自适应调优就是应用遗传算法对DGEMM分块参数进行了智能优化选择, 将分块参数作为个体, 一定范围内的多个个体形成一个种群进行自适应迭代进化, 从而针对具体平台体系结构产生最优的DGEMM矩阵分块参数.文献[29]在Intel平台上对遗传算法在DGEMM的分块参数智能优化选择上的实现做了详细介绍.如图 4所示是遗传算法的基本流程.
|
Fig. 4 Generic algorithm flow chart 图 4 遗传算法流程图 |
● 首先初始化种群, 以分块参数[mc, kc, nc]为一个个体, 随机生成一组个体的集合.
● 然后计算种群每个个体的适应度, 适应度大小决定遗传概率.
● 接着判断迭代次数, 如果没有达到预设的最大值, 则进行selection、crossover和mutation的遗传算法进化步骤.
● 最后经过selection、crossover和mutation这3个步骤后生成一个新的种群, 与原始种群大小一致, 但个体都经过优胜劣汰, 进化到了更好的水平.重复上诉步骤, 直到满足终止条件退出.
根据上述遗传算法的介绍, 自适应程序的实现需要经过编码个体、建立种群、计算适应度、选择(selection)、交叉(crossover)、变异(mutation)和终止条件判断这几个步骤.
(1) 个体编码
编码方法有二进制编码、格雷编码、浮点数编码和符号编码等多种方式, 根据DGEMM分块参数特点, 采取浮点数编码更适合.也就是将分块参数[mc, kc, nc]这3个浮点数直接作为参数传入.
(2) 建立种群(population)
个体集合的建立过程就是要在一个确定范围内去随机生成一定数量的个体.种群范围的确定可以根据cache的大小进行理论计算得出[7], HygonBLIS分块参数原始值是[1080, 120, 8600], 在原始参数值的基础上扩大范围至[512 < mc < 2400, 60 < kc < 2400, 4800 < nc < 12000], 随机生成的分块参数要保证能够被内核参数整除, 也就是被分配的寄存器数[4, 6, 6]整除.
(3) 计算适应度(fitness)
适应度函数就是DGEMM测试程序, 其性能作为判断适应度高低的标准.性能高的分块参数会遗传给下一代种群, 性能低分块参数就会被淘汰.
(4) 选择(selection)
基本遗传算法仅仅使用selection、crossover和mutation这3种算子是没有办法收敛于全局最优解的, 因为简单的进行杂交, 可能反而会把较好的组合破坏掉.本文采用精英保留策略, 也就是把目前最优个体不进行任何操作直接复制到下一代种群中.选择方法有很多种, 包括比例选择法、排序选择等.本文采用比例选择法和精英保留策略结合, 种群中每一个个体的适应度占总适应度的比例越大, 其被选择的概率就越大, 同时这一代种群最优秀个体无条件复制到下一代.
(5) 交叉(crossover)
单点交叉、双点交叉、均匀交叉等都适合二进制编码和浮点数编码, 根据GEMM分块参数只有3个, 选择单点交叉, 按某一概率将两个个体的某一位置参数互换.交叉率设置越大, 新个体产生越快, 其全局搜索能力就会越强, 但是越大可能破坏优秀个体, 综合考虑设置交叉概率设置0.8.
(6) 变异(mutation)
变异的主要目的是改善局部搜索能力以及维持种群个体的多样性, 变异方法采用均匀变异, 按某一概率在在一个范围内随机生成一个数替换原来个体中某一位置的参数.变异概率不宜设置过大, 导致丢失优秀个体而变成随机搜索.考虑设置参数范围已经在合理范围内, 本文按0.1的概率对每一个个体的某一参数进行变异.
(7) 终止条件判断(iterations)
种群中所有的个体适应度计算完成后要判断设置的终止条件是否满足来判断程序是否结束退出.如果不满足终止条件就进行selection、crossover和mutation操作, 形成下一个同等规模的个体集合(种群).
Auto-tuning自适应调优算法完成后, 需要进行DGEMM测试, BLIS库中config目录是存放不同平台kernel的配置信息, 其中DGEMM分块参数通过bli_kernel.h头文件中的BLIS_DEFAULT_MC_D、BLIS_DEFAULT_ KC_D和BLIS_DEFAULT_NC_D这3个宏定义配置.
Auto-tuning的具体测试流程如下:
(1) 在包含初始分块的更大的范围内随机生成一组满足遗传算法个体要求的分块参数[mc, kc, nc], 并用其替换掉blis_kernel.h中的分块参数.
(2) 自动编译BLIS库, 生成此分块参数下的动态链接库libblis.so.
(3) 调用libblis.so执行测试程序DGEMM, 主程序将性能值记录并返回.
通过Auto-tuning的不断迭代, 根据算法模型, 我们会得到一个全局最优解, 这组参数就是性能更高的DGEMM分块参数, 实测数据分块参数由[1080, 120, 8400]优化为[792, 822, 8628], 将这组参数替换进BLIS中, 并使用高性能计算标准测试程序HPL和单独的DGEMM测试程序进行测试, 具体测试对比数据参考第5节性能分析.
2.3 矩阵打包连续的内存访问比非连续的内存访问要快.异构HPL中传给DGEMM的矩阵A、B和C的索引是初始化矩阵的某个位置指针, 矩阵元素通过leading-dimension索引, leading-dimension表示列优先存储方式的矩阵同一行(或行优先存储方式的矩阵同一列)中两相邻元素的距离, 也就是实际存储矩阵的二维数组的第1维(或第2维)的大小.矩阵A、B和C的行或列间跨度通过LDA、LDB和LDC表示, 这3个数值可能会非常大.比如异构HPL求解问题Ax=b, 其矩阵规模N的最大取值可以使矩阵A的存储空间占用内存80%左右, 则LDA不能小于N.在算法实现中, 大kernel算法会不断地对A和B的不同分块进行打包, 形成Ac和Bc; 而小kernel算法只将原始矩阵B进行一次打包, 形成行优先存储(row-major)的Bpacked.如图 5所示, 其中假设nr等于6, 矩阵B是以列优先存储(column-major)方式存储数据, 且行数k和列数n都是nr的整数倍.
|
Fig. 5 Packing matrix B in little GEMM kernel 图 5 GEMM小kernel中矩阵B的打包操作 |
Packing数据打包操作的好处有3个:
(1) 创造连续的内存访问方式.
BLAS接口传入的矩阵存储方式是不确定的, 可能是row-major也可能是column-major, 且矩阵的行或列间的跨度会非常大, 最内层的kernel为了适应不同的数据存储方式以及不连续的数据分布, 需要通过packing制造统一的、连续的内存布局视图.
(2) 数据对齐.
对于大多数CPU微架构实现而言, 对齐到寄存器宽度, cache边界和页边界的数据访问往往比非对齐的数据访问方式更快.
(3) 提前load数据到cache.
通过Packing操作, 最内层kernel可以在尽可能接近CPU寄存器的cache上获取到需要进行运算的操作数, 这比直接从内存中获取数据快很多.
2.4 数据预取Dhyana CPU提供了硬件取功能和符合X86 ISA的软件预取指令.其中硬件预取对开发者是透明的, 软件预取指令PREFETCHn指示处理器将后续需要访问的数据提前读入到指定的某级缓存中, 当程序使用这些数据时, 可以直接读取缓存.PREFETCHn指令格式见表 1.
|
Table 1 X86 software prefetch instructions 表 1 X86软件预取指令 |
软件预取指令不要求数据行对齐, 如果预取数据行已经存在于比指定的n更低的某层cache中, 指令将被当作NOP空指令处理.
在最内层的核心汇编代码中, PREFETCHn指令的使用主要考虑两点:
(1) 何处插入PREFETCHn指令.
(2) 预取的地址跨度设置.
Dhyana CPU的cache line大小为64B, 一个缓存行可以存放8个双精度浮点数据.如图 6所示的汇编代码第1行的软件预取操作, 表示把Ac的第9个双精度数据开始的一个缓存行预取到L1缓存中.此处的预取操作是为第22行的load指令准备数据.第5行的load操作会在L1缓存发生cache miss, 但是因为对A进行了packing操作, 理论上能够在L2 cache命中.同理, 第9行、第14行和第18行处的load操作会在L1中命中, 第11行预取操作将Ac的第17个数据预取到内存.
|
Fig. 6 Code segment of inner-most DGEMM kernel 图 6 DGEMM核心汇编代码片段 |
软件预取指令的插入位置需要综合考虑体系结构特点, 合理安排汇编指令.Dhyana CPU的L2 load-to-use的latency至少是12个core cycles, 而每时钟周期只可以执行一条FMA指令, 在CPU流水线能够连续retire FMA指令的理想状态下, 预取指令和后续load指令之间至少插入12条FMA指令.
3 指令集优化高效BLAS kernel的实现方法具有通用性, 不同于其他数学库大部分代码通过汇编实现优化, BLIS只在最内层的kernel核心代码处使用了汇编代码, 框架内的其他代码使用C实现.这样给可移植性带来了极大的便利, 在适配不同架构的CPU时, 开发者可以将主要精力集中在最内层kernel的汇编代码的优化上.
3.1 X86向量指令不同架构的CPU大多会提供高效的SIMD(single instruction multiple data)指令, Dhyana CPU支持AVX2指令集, 可以使用128位的XMM寄存器和256位的YMM寄存器, 但由于浮点运算部件的数据路径宽度是128位, 每个时钟周期最多可以同时处理两个双精度浮点数, 单核理论双精度浮点性能是(128/64)×2=4 DP-FLOPs/cycle, 实际测试结果表明, 使用XMM寄存器的效率要优于YMM寄存器.AVX2指令集中主要有两类指令: 数据传输类指令和数据操作类指令.其中, VMOVDDUP、VMOVAPD、VMOVUPD、VMULPD和VFMADD231PD是在Dhyana CPU平台上高效BLAS kernels的实现中使用频率非常高的几条指令.
HPL算法中传递给BLAS函数的参数会有固定值的情况, 可以通过适当的指令替换来精简指令, 举例如下.
(1) 对于Level-3级函数DGEMM
系数Alpha固定为–1.0, Beta固定为1.0.此时的DGEMM计算公式由C: =α·AB+β·C变成了C: =C-AB可以在packing B(见2.3节)的过程中使用xorpd指令进行符号位翻转, 精简了对系数Alpha的乘法操作.矩阵C不变, 精简了对系数Beta的乘法操作.精简指令后, 原来需要4个core cycle的乘法操作和5个core cycle的乘加指令指令只需3个core cycle的加法指令即可完成.在矩阵规模较大, 最内层kernel被循环多次的情况下, 指令时钟周期的减少可以优化整体性能.
(2) 对于Level-2级函数DGEMV
转置参数TransA固定为NoTrans, 可以省略对矩阵A进行转置操作的代码, 提高CPU流水线分支预测的准确率, 避免分支预测失败的惩罚.Alpha固定为–1.0, Beta固定为1.0, 则可以在内层kernel中使用一条VFNMADD231PD指令取代VMULPD和VFMADD231PD两条指令.
(3) 对于Level-1级函数DSCAL
IncX固定为1, 可以省去对数据成员间隔较大时所需的软件预取操作的指令(见第2.4节).
异构HPL传递给BLAS接口的参数还有其他的固定值情况存在, 这些固定值的存在可以起到精简指令的作用, 为BLAS库提供了优化空间.
3.2 循环展开循环展开是编译器常用的优化方式之一, 可以减少分支预测失败带来的开销, 提高指令级的并行度.但高级语言通过编译器循环展开后的汇编代码在BLAS这种访存密集型的程序特征下, 往往性能表现较差.因此BLAS的核心循环代码, 通常采用手写汇编的方式进行展开.
手写BLAS核心汇编代码, 不仅需要考虑底层CPU架构层面的架构寄存器的数量和复用方式, 还需要考虑底层CPU的微架构特点, 比如浮点部件的宽度、存储层级的大小和层级间的延迟等.这些CPU特性共同决定了手写汇编循环展开的程度, 也就是“循环展开因子”.本文第2.1节中提到的5个参数nc、kc、mc、nr和mr, 其中内核参数nr和mr与汇编代码循环展开直接相关, 而矩阵分块参数nc, kc和mc则与CPU存储层级特性有直接相关.针对具体平台体系结构可以理论计算出这5个参数的大致取值[30], 但实际代码测试中发现, 理论数据往往并不是最优设置, Dhyana CPU有16个128位的向量寄存器XMM0-XMM15, 内核参数nr, mr取值的首要策略是最大化利用架构寄存器, 因此综合考虑nr和mr分别取值为6和4.同时矩阵参数[nc, kc, mc]依据2.2节中的自动调优技术同样获得了比理论值更好的参数设置.
如图 6所示, 通用寄存器RAX用于索引Ac的micro-panel中的数据, 不断获取一个双精度数据然后复制成两份填充到XMM0中; RBX用于索引Bc的micro-panel中的数据, 一次获取并存储6个双精度数据到XMM1~ XMM3中.XMM4~XMM15用于存储不断累加的中间结果.在此过程中, 向量寄存器XMM0实际共存储过4个Ac数据, XMM1~XMM3共存储过6个Bc数据.
4 多线程并行常用的多线程并行手段有OpenMP[31]和POSIX Threads(pthread)[32].BLIS同时支持两种并行方式, 且当前只在Level-3实现了并行.与OpenMP相比, Pthreads技术在实际编程中需要考虑临界区、线程同步原语等非常底层的操作, 故我们在GEMM以及TRSM等小kernel和Level-1与Level-2的并行化实现中, 选择使用更加便利的OpenMP进行并行化.
4.1 control-tree优化针对Level-3级BLAS的并行实现方式, BLIS提出了新颖的“control-tree”结构[7], 该结构以统一的方式实现了如图 3所示的GEMM大kernel算法的并行以及HERK、TRMM和TRSM运算的并行.Control-tree结构最突出的优点是可以在任意层loop内分别使用环境变量BLIS_JC_NT、BLIS_IC_NT、BLIS_JR_NT和BLIS_IR_NT进行多线程并行的控制, 甚至是多个维度同时并行.这种统一结构的便利性的代价就是在所有能够发生并行的地方设置同步屏障(barrier), 以便当前loop的下层所有操作全部返回后才能继续回溯.虽然文献[7]中的性能测试结果表明control-tree结构的开销是非常小的, 但在异构HPL环境中, control-tree的开销变得尤为突出.BLIS原生control-tree的结构大致如图 7所示.
|
Fig. 7 Original BLIS GEMM control-tree 图 7 BLIS原生GEMM control-tree |
对于异构HPL中的GEMM运算, 我们选择只在loop3处(Mc维度)实现并行, 当各个线程的Ac packing打包操作和后续操作都完成后, 需要在loop3的Ic barrier处进行第1次同步, 随后的Pc barrier和global barrier处再次进行同步, 才可以保证最终GEMM运算的正确性.可以看到, Pc barrier和global barrier的同步开销其实是可以避免的, 只需要将loop3放到算法的最前端, 然后使用OpenMP自身的同步操作实现各个线程的同步即可.基于此分析, 对GEMM运算的control-tree结构做出了如图 8所示的调整.
|
Fig. 8 Optimized HygonBLIS GEMM control-tree 图 8 优化后的HygonBLIS GEMM control-tree |
通过异构HPL的实际测试, 未优化control-tree时DGEMM和DTRSM运算的同步barrier等待的时间开销为2.29s, 约占CPU端计算总时间的20%.优化control-tree结构后, GEMM函数同步barrier的等待时间减少为1.49s, TRSM运算应用类似优化思路, 同步barrier的等待时间减少到基本可以忽略不计.
4.2 Level-1和Level-2级BLAS并行优化原BLIS库没有实现Level-1级和Level-2级BLAS函数运算的并行化.本文使用如下相似的策略实现了DGEMV、DTRSV、DSCAL、DAXPY以及IDAMAX运算的并行化:
(1) 通过运行脚本中设置的环境变量“OMP_NUM_THREADS”获取当前MPI进程中的线程数.
(2) 将矩阵或向量中的某个维度按设置的线程数进行切分.
(3) 最后使用OpenMP的Pragma(编译指示)启动多线程, 在不同的线程空间中复用同一段汇编kernel来处理不同的数据.在Pragma包裹的并行域(parallel region)中需要合理而准确地为每个线程划分任务空间.可以通过omp_get_thread_num()函数获得当前线程的线程编号(TID)划分任务, 避免计算冲突.
多线程并行对于BLAS运算的性能提升是非常显著的, 矩阵或向量计算达到一定规模后, 性能往往随着线程数的增加而提升, 但提升并不会完全是线性的, 通常是计算越密集的运算提升的线性程度越高, 因此计算量最大的Level-3级BLAS运算通过多线程并行会大大提升性能.如图 9所示, 在N=43776时的Level-1运算IDMAX随着线程数的增加表现出的性能提升.
|
Fig. 9 Performance improvement by multi-threading of IDMAX 图 9 多线程对IDMAX的性能提升 |
4.3 并行化问题
多线程并行化并不是总会对计算带来帮助, 为了进一步提高性能, 需要综合考虑以下问题:
(1) 动态多线程
创建和启动多线程的开销在计算量很小的情况会变的突出, 从而导致程序性能下降.所以需要设置一个或多个特定的阈值来进行线程数的切换.具体是设置一个阈值将多线程直接切换到单线程, 还是设置多个阈值将线程数分段递减, 需要具体的BLAS函数结合特定的输入数据进行大量的测试归纳(profiling).在DGEMM的多线程优化中, 我们设计了精确的运算量模型.DGEMM需要计算一个m×k阶矩阵A和k×n阶矩阵B的乘积, 其主要计算量和乘积mkn成正比.考虑到矩阵A的规模m·k, 矩阵B的规模k·n等都对最终运算性能造成影响, 我们使用最小二乘法拟合矩阵计算的性能公式, 假设不同线程数目下运算时间根据式(2)计算.
| $ T = {a_1}mkn + {a_2}mk + {a_3}mn + {a_4}nk + {a_5}m + {a_6}n + {a_7}k + {a_8} $ | (2) |
于是对于输入参数为mi、ki和ni, 运行时间为ti的测试数据, 对应拟合数据中一行如式(3)所示.
| $ \left( {{m_i}{k_i}{n_i}, {m_i}{k_i}, {m_i}{n_i}, {n_i}{k_i}, {m_i}, {n_i}, {k_i}, 1} \right) $ | (3) |
所有这些拟合数据组成N×8阶矩阵P, 其中N是采样数, 而所有的采集时间构成对应的列向量b, 于是构成式(4)的拟合方程.
| $ {P^T}Px = {P^T}b $ | (4) |
解方程得出的8阶变量x的各分量就是a1, a2, …, a8的最小二乘拟合结果.分别对单线程、4线程、8线程和32线程测试拟合公式, 在运行中根据拟合结果动态切换线程数据, 得到如图 10所示的性能结果.
|
Fig. 10 Performance speedup of DGEMM in 32 threads 图 10 DGEMM在32个线程中的性能加速比 |
实践中还发现, 不同的kernel实现方式以及硬件设置等因素也会影响到阈值的选取.比如, 在异构HPL的DGEMV函数的优化过程中, 根据不同的“OMP_NUM_THREADS”、不同的CPU主频以及矩阵A的列数N(M方向做并行)设置不同的多个阈值.
(2) 均衡余量(remainder)处理
当使用多线程将求解问题切分成多块后, 往往会遇到不能均分的情况.所谓不能“均分”, 是指整体的运算量不能均匀的分配到每个线程中去.那么此时有两个选择: 所有多出的运算量全部放到某一个线程中, 比如最后一个线程; 或者将余量再次均分到所有线程中.Remainder处理的问题在线程数较少时(比如6个线程)不会那么关键, 但在线程数较多时(比如32个线程)会变得突出.因为此时的木桶效应更加明显, 整体运算的同步结束取决于拥有最大余量的那个线程.根据本文HPL算法测试分析后, 采取了将remainder分配给每个线程处理, 得到的性能结果最优.
(3) 循环展开与多线程的均衡
循环展开的程度和线程的数量需要均衡考虑.在使用动态多线程去加速运算的过程中, 往往会有多种求解方法去完成同一个运算.考虑以下场景: 异构HPL中的DGEMM运算在M较小时, 比如M=28, 那么此时可以使用4个线程各处理6个M维度的行, 然后使用1个线程处理剩下的4个M行; 或者使用5个线程各处理5个M维度的行然后使用3个线程处理各处理一个M行, 如式(5)所示的等式.
| $ M{\rm{ = }}28{\rm{ = }}4T(thread) \times 6 + 1T \times 4 = 5T \times 5 + 3T \times 1 $ | (5) |
不同的求解方法涉及同一个函数的不同kernel函数编写以及动态多线程threshold的选取.通过实际测试发现: DGEMM运算在M较小时, 更多的循环展开比更多的线程性能更好, 如图 11所示.
|
Fig. 11 In the case that B matrix is square matrix (dimension N=K), more unrolling is better performance than more threading in small DGEMM 图 11 B矩阵是方阵(维数N=K)的情况下, 在小规模的DGEMM中, 更多的循环展开比更多的线程性能更优 |
5 性能测试与分析
异构HPL算法中综合使用了以上介绍的访存优化、指令集优化和多线程并行等优化技术对算法中调用的7个BLAS函数进行了针对性优化, 并与OpenBLAS 0.3.6和MKL(Intel Math Kernel Library 2018.0.128)进行了对比测试.实验测试平台为高性能计算集群, 单个计算节点配置如引言所述, 32核Dhyana CPU和4个协处理器组成的异构计算平台.运行单核DGEMM测试时, 单核心的主频锁定在3.0GHz; 运行多核心DGEMM测试和异构HPL测试时, 所有核心的主频锁定在2.55GHz.内存主频2667MHz, 协处理器主频为1.55GHz.
5.1 Cache性能分析PMU全称为performance monitor unit, 即性能监控单元.通过系统性能分析工具读取CPU中各种计数器的数值, 可以分析应用程序的各项指标.Dhyana处理器支持非常多的性能监控计数器, 包括浮点、访存、指令等CPU各种资源的监控.如表 2所示为Auto-tuning优化后的矩阵分块参数与初始分块参数进行DGEMM计算时的cache miss对比情况, 分析得出分块参数[mc, kc, nc]从[1080, 120, 8400]优化为[792, 822, 8628], kc和nc的增大导致矩阵打包到L1缓存时, 因L1容量小而发生更多的cache miss.但是L2和L3缓存容量比L1大得多, 其命中率提高可以掩盖L1的性能损失并提高整体DGEMM的性能.
|
|
Table 2 Cache miss rate of HPL counted by PMU 表 2 通过PMU统计HPL的缓存失效率 |
5.2 DGEMM效率
如图 12和图 13所示分别为单核和32核情况下不同数学库的DGEMM性能对比.因为Intel MKL只有在Intel平台上才能发挥最佳效果, 在Hygon平台上效率不是最好, 仅作为参考; 而Netlib为未经过任何优化的Fortran语言实现的通用标准数学库.可以看到HygonBLIS的DGEMM实现无论是在单核心还是在32核心下性能均超过了其他BLAS库, 单核DGEMM最高效率达到了98%, 32核最高效率达到了96.9%.与未优化前的BLIS和MKL相比, 单核性能分别提升19.2%和33.3%, 32核性能分别提升42.1%和33.6%.
|
Fig. 12 DGEMM efficiency comparison of single-core 图 12 单核DGEMM效率对比 |
|
Fig. 13 DGEMM efficiency comparison of 32 cores 图 13 32核DGEMM效率对比 |
5.3 异构HPL测试
HPL实测性能值是最终评价标准, 分别统计每个BLAS函数的时间作为BLAS优化的评价标准.为避免网络通信的干扰, 本文基于单节点进行异构HPL测试, 对比包括panel分解中pfact时间以及HPL算法中使用的7个BLAS函数DGEMM、DGEMV、DTRSM、DTRSV、IDMAX、DSCAL和DCOPY的计算时间.
异构HPL算法panel分解部分主要在CPU端进行计算.如图 14左侧为HPL程序panel分解中的pfact总时间统计, 调用HygonBLIS的pfact与调用MKL的相比, 计算时间节省了53%.如图 14右侧图为DGEMM函数计算时间统计, 可以看到, HygonBLIS的DGEMM性能比MKL的DGEMM性能提升了约25%.
|
Fig. 14 Performance comparison of pfact and DGEMM calling different math libraries in heterogeneous HPL 图 14 异构HPL中调用不同数学库的pfact和DGEMM的性能对比 |
如图 15所示为DGEMV和DTRSV函数性能对比.DGEMV在panel分解中时间占比仅次于DGEMM, 因此DGEMV性能对整体性能影响也仅次于DGEMM.图 15左侧为DGEMV性能, HygonBLIS比MKL提升约62%, 比原BLIS提升约65%.图 15右侧为DTRSV函数时间统计, 原BLIS性能与MKL差距较大, 经过优化后的HygonbBLIS性能超过了MKL.因为DTRSV计算量较小, 在总时间中占比也相对很小.
|
Fig. 15 Performance comparison of DGEMV and DTRSV calling different math libraries in heterogeneous HPL 图 15 异构HPL中调用不同数学库的DGEMV和DTRSV的性能对比 |
如图 16所示为DTRSM和IDMAX函数性能对比.原生BLIS算法库中的三角矩阵求解DTRSM性能较差, 优化后的DTRSM性能同MKL基本一致.在pfact第106步~第141步的性能是超过MKL的, 而在panel分解最后段MKL运算更快, 原因是DTRSM在切换到小kernel的threshold和动态多线程的threshold需要做进一步的调整, 以适应运算量的变化.IDMAX函数主要是用于选主元操作, 计算量可以忽略不计.
|
Fig. 16 Performance comparison of DTRSM and IDAMAX calling different math libraries in heterogeneous HPL 图 16 异构HPL中调用不同数学库的DTRSM和IDAMAX的性能对比 |
如图 17所示, Level-1级BLAS函数DSCAL和DCOPY优化后的性能基本与Intel MKL一致.其中DSCAL优化后的性能比原生BLIS有大幅提升, 但在panel分解后段的性能还是略低于MKL, 原因是MKL使用的多线程库IOMP的比开源的GOMP在多线程启动、调度和动态多线程的切换上更加的高效.而DCOPY代码实现简单, 优化空间非常有限, 硬件平台性能直接决定了其时间开销, 故可以看到三者DCOPY的性能基本一致.DAXPY运算只在回代过程(back substitution)过程中调用一次, 性能影响可以忽略, 故未列出详细对比.
|
Fig. 17 Performance comparison of DSCAL and DCOPY calling different math libraries in heterogeneous HPL 图 17 异构HPL中调用不同数学库的DSCAL和DCOPY的性能对比 |
6 总结与展望
本文针对Hygon CPU平台系统架构特点, 在开源BLIS算法库的基础上开发了HygonBLIS数学库, 详细介绍了HPL调用的各级BLAS函数优化方法.通过访存优化、指令集优化和多线程并行优化, HPL调用的各级BLAS函数的计算性能与MKL的相应函数相比有着普遍和显著的提升, 个别计算量小的函数性能基本与MKL相应函数性能一致, 其中最重要的DGEMM函数提升约25%, DGEMV函数性能提升约62%.在异构单节点HPL测试中, 与调用MKL库的HPL相比, 采用HygonBLIS的HPL整体性能提升2.8TFLOPs, 效率值提升11.8%, 更好地发挥了异构高性能平台的计算能力.同时本文也有不足之处, 虽然HygonBLIS库性能测试发现DGEMM、DGEMV、DTRSV和IDMAX有着很好的优化效果, 其性能与MKL和开源BLIS库相应函数相比都有大幅提升, 但是DTRSM和DSCAL优化后性能并没有优于Intel MKL, 这可能与计算量、多线程库支持和编译器优化支持都有密切关联.底层子例程库优化是一个系统工程, 不仅需要硬件架构的优化设计, 还需要编译器和线程库等底层软件的优化.本文未涉及的其他BLAS函数的优化将是HygonBLIS库下一步优化工作的方向.
| [1] |
Whaley RC, Dongarra JJ. Automatically tuned linear algebra software. In: Proc. of the 1998 ACM/IEEE Conf. on Supercomputing (SC'98). San Jose, 1998. 1-27.
|
| [2] |
Goto K, van de Geijn RA. Anatomy of high-performance matrix Multiplication. ACM Trans. on Mathematical Software, 2008, 34(3): 1-25.
http://www.ams.org/mathscinet-getitem?mr=2444070 |
| [3] |
Goto K, van de Geijn RA. High-performance implementation of the Level-3 BLAS. ACM Trans. on Mathematical Software, 2008, 35(1): 1-14.
http://dl.acm.org/doi/10.1145/1377603.1377607 |
| [4] |
Wang Q, Zhang X, Zhang Y, Qing Y. AUGEM: Automatically generate high performance dense linear algebra kernels on X86 CPUs. In: Proc. of the Int'l Conf. on High Performance Computing, Networking, Storage and Analysis (SC 2013). Denver, 2013. 1-12.
|
| [5] | |
| [6] | |
| [7] |
Van Zee FG, van de Geijn RA. BLIS: A framework for rapidly instantiating BLAS functionality.. ACM Trans. on Mathematical Software, 2015, 41(3): 1-33.
http://doi.acm.org/10.1145/2738033 |
| [8] |
Van Zee FG, Smith TM, Marker B, Low TM, van de Geijn RA, Igual FD, Smelyanskiy M, Zhang XY, Kistler M, Austel V, Gunnels JA, Killough L. The BLIS Framework: Experiments in Portability. ACM Trans. on Mathematical Software, 2016, 42(2): 1-19.
http://citeseerx.ist.psu.edu/viewdoc/download;jsessionid=7ABD5D33DA6BC7727323E53708D54277?doi=10.1.1.361.6664&rep=rep1&type=pdf |
| [9] |
Smith TM, van de Geijn RA, Smelyanskiy M, Hammond JR, Van Zee FG. Anatomy of high-performance many-threaded matrix multiplication. In: Proc. of the IEEE 28th Int'l Parallel and Distributed Processing Symp. 2014. 1049-1059.
|
| [10] |
Gu N, Li K, Cheng G, Wu C. Optimization of BLAS based on Loongson 2F architecture. Journal of Univercity of Science and Technology of China, 2008, 38(7): 854-859(in Chinese with English abstract).
https://www.cnki.com.cn/Article/CJFDTOTAL-ZKJD200807022.htm |
| [11] |
Dongarra JJ, Luszczek P, Petitet A. The LINPACK benchmark: Past, present, and future. Concurrency and Computation: Practice and Experience, 2003, 15(9): 803-820.
[doi:10.1002/cpe.728] |
| [12] |
Tan G, Li L, Triechle S, Phillips E, Bao Y, Sun N. Fast implementation of DGEMM on Fermi GPU. In: Proc. of the 2011 Int'l Conf. on High Performance Computing, Networking, Storage and Analysis (SC 2011). Seattle, 2011. 1-11.
|
| [13] |
Jiang H, Wang F, Zuo K, Su X, Xue L, Yang C. Design and implementation of a highly efficient DGEMM for 64-bit ARMv8 multi-core processors. In: Proc. of the 44th Int'l Conf. on Parallel Processing. Beijing, 2015. 200-209.
|
| [14] |
Jiang H, Wang F, Li K, Yang C, Zhao K, Huang C. Implementation of an accurate and efficient compensated DGEMM for 64-bit ARMv8 multi-core processors. In: Proc. of the IEEE 21st Int'l Conf. on Parallel and Distributed Systems (ICPADS). Melbourne, 2015. 491-498.
|
| [15] |
Wang L, Wu W, Xu Z, Xiao J, Yang Y. BLASX: A high performance level-3 BLAS library for heterogeneous multi-GPU computing. In: Proc. of the 2016 Int'l Conf. on Supercomputing (ICS 2016). Istanbul, 2016. 1-11.
|
| [16] |
Sun J, Sun Q, Deng P, Yang C. Research on the optimization of BLAS level 1 and 2 functions on Shenwei many-core processor. Computer Systems & Applications, 2017, 26(11): 101-108(in Chinese with English abstract).
https://www.cnki.com.cn/Article/CJFDTOTAL-XTYY201711014.htm |
| [17] |
Liu H, Liu F, Zhang P, Yang C, Jiang L. Optimization of BLAS level 3 functions on SW1600. Computer Systems & Applications, 2016, 25(12): 234-239(in Chinese with English abstract).
https://www.cnki.com.cn/Article/CJFDTOTAL-XTYY201612037.htm |
| [18] |
Guo Z, Guo S, Xu J, Zhang Z. Register allocation in base mathematics library for platform of heterogenerous multi-core. Journal of Computer Applications, 2014, 34(S1): 86-89(in Chinese with English abstract).
https://www.cnki.com.cn/Article/CJFDTOTAL-JSJY2014S1026.htm |
| [19] | |
| [20] | |
| [21] |
Fatica M. Accelerating Linpack with CUDA on heterogenous clusters. In: Proc. of the 2nd Workshop on General Purpose Processing on Graphics Processing Units (GPGPU 2009). Washington, 2009. 46-51.
|
| [22] |
Yang C, Wang F, Du Y, Chen J, Liu J, Yi H. Adaptive optimization for petascale heterogeneous CPU/GPU computing. In: Proc. of the 2010 IEEE Int'l Conf. on Cluster Computing. Heraklion, 2010. 19-28.
|
| [23] |
Yamazaki I, Tomov S, Dongarra J. One-sided dense matrix factorizations on a multicore with multiple GPU accelerators. Procedia Computer Science, 2012, 9(11): 37-46.
http://www.sciencedirect.com/science/article/pii/S1877050912001263 |
| [24] |
Yang C, Chen C, Tang T, Chen X, Fang J, Xue J. An energy-efficient implementation of LU factorization on heterogeneous systems. In: Proc. of the IEEE 22nd Int'l Conf. on Parallel and Distributed Systems (ICPADS). Wuhan, 2016. 971-979.
|
| [25] |
Jo G, Nah J, Lee J, Kim J, Lee J. Accelerating LINPACK with MPI-OpenCL on clusters of multi-GPU nodes. IEEE Trans. on Parallel & Distributed Systems, 2015, 26(7): 1814-1825.
|
| [26] |
Li J, Li X, Tan G, Chen M, Sun N. An optimized large-scale hybrid DGEMM design for CPUs and ATI GPUs. In: Proc. of the 26th ACM Int'l Conf. on Supercomputing (ICS 2012). New York, 2012. 377-386.
|
| [27] |
Li L, Yang W, Ma W, Zhang Y, Zhao H, Zhao H, Li H, Sun J. Optimization of HPL on complex heterogeneous computing system. Ruan Jian Xue Bao/Journal of Software, 2021,32(8):2307-2318(in Chinese with English abstract). http://www.jos.org.cn/1000-9825/6003.htm [doi:10.13328/j.cnki.jos.006003] The optimization of HPL on a complex heterogeneous computing system.
|
| [28] | |
| [29] |
Sun C, Lan J, Jiang H. Genetic algorithm for deciding blocking size parameters of GEMM in BLAS. Computer Engineering & Science, 2018, 40(5): 798-804(in Chinese with English abstract).
[doi:10.3969/j.issn.1007-130X.2018.05.006] |
| [30] |
Low T, Igual F, Smith T, Quintana-Orti E. Analytical modeling is enough for high-performance BLIS. ACM Trans. on Mathematical Software, 2016, 43(2): 1-18.
|
| [31] |
Dagum L, Menon R. OpenMP: An industry standard API for shared-memory programming. IEEE Computational Science and Engineering, 1998, 5(1): 46-55.
[doi:10.1109/99.660313] |
| [32] | |
| [10] |
顾乃杰, 李凯, 陈国良, 吴超. 基于龙芯2F体系结构的BLAS库优化. 中国科学技术大学学报, 2008, 38(7): 854-859.
https://www.cnki.com.cn/Article/CJFDTOTAL-ZKJD200807022.htm |
| [16] |
孙家栋, 孙乔, 邓攀, 杨超. 基于申威众核处理器的1、2级BLAS函数优化研究. 计算机系统应用, 2017, 26(11): 101-108.
https://www.cnki.com.cn/Article/CJFDTOTAL-XTYY201711014.htm |
| [17] |
刘昊, 刘芳芳, 张鹏, 杨超, 蒋丽娟. 基于申威1600的3级BLAS GEMM函数优化. 计算机系统应用, 2016, 25(12): 234-239.
https://www.cnki.com.cn/Article/CJFDTOTAL-XTYY201612037.htm |
| [18] |
郭正红, 郭绍忠, 许瑾晨, 张兆天. 异构多核平台下基础数学库寄存器分配方法. 计算机应用, 2014, 34(S1): 86-89.
https://www.cnki.com.cn/Article/CJFDTOTAL-JSJY2014S1026.htm |
| [27] |
黎雷生,杨文浩,马文静,张娅,赵慧,赵海涛,李会元,孙家昶.复杂异构计算系统HPL的优化.软件学报,2021,32(8):2307-2318(in Chinese with English abstract). http://www.jos.org.cn/1000-9825/6003.htm [doi:10.13328/j.cnki.jos.006003]
|
| [29] |
孙成国, 兰静, 姜浩. 一种基于遗传算法的BLAS库优化方法. 计算机工程与科学, 2018, 40(5): 798-804.
[doi:10.3969/j.issn.1007-130X.2018.05.006] |