Page 39 - 《软件学报》2021年第8期

水超洋 等:国产异构系统上 HPL 的优化与分析 2321

叠探索了一种粗粒度的流水线算法.
作为 HPL 的核心运算,双精度矩阵乘法的优化在 HPL 的优化中占有重要地位并且得到了全面而深入的研
究.在传统 CPU 平台上,文献[9]给出了 CPU 上的矩阵乘法的分层算法,通过对 CPU 存储层次的模拟设计出相应
的多级分块缓存策略,以尽可能地利用高速缓存中的数据.在包含加速器的异构架构上,李佳佳等提出了五阶段
流水的异构矩阵乘算法来掩盖 CPU 与加速器之间的 PCIe 数据传输 [10,11] .MAGMA [12] 通过细粒度的任务划分并
且灵活地在 CPU 和加速器上调度这些任务来实现负载的动态平衡.为充分利用加速器如 GPU 的计算能力,文献
[13]通过微基准测试来探测 GPU 的体系结构,在汇编语言层面做了多种优化,以实现接近 GPU 理论浮点峰值性
能的双精度矩阵乘效率 [13] .
关于 HPL 的性能建模,文献[14]以预测 HPL 的扩展性为目的给出了 HPL 在 CPU 系统上的性能模型.HPL
的求解时间被建模为 panel 分解,panel 广播,行交换和尾矩阵更新的时间之和,并给出了模型中一些常量系数的
经验值.王申等人在这个基础上考虑了 look-ahead 算法中部分广播开销可以被尾矩阵的更新和行交换所掩盖的
情况 [15] ,给出了更为精确的模型.文献[16]认为上述模型的计算都不够精确,因为模型中的常量值都是经验值,而
CPU 计算效率以及通信带宽等模型常量都会受数据量大小的影响.他们主张将这些常量系数视为可变的,通过
已有的测试结果去学习这些系数,用学习得到得系数建模预测大规模求解的性能.
已有文献中对 HPL 的优化主要集中在同构架构上的简单算法优化和性能建模以及异构架构上的双精度
矩阵乘的优化上,而缺少异构架构上的 HPL 的算法建模和算法流水层面的优化.本文在 CPU 的 HPL 性能模型
的基础上建立了国产处理器-国产加速器异构架构上的 HPL 性能模型,并提出了多线程细粒度的 HPL 流水线算
法,以充分发挥异构系统中国产处理器的巨大计算能力.在实现上,我们实现了一个轻量级的跨平台异构加速框
架 HPCX,并用生产者消费者模型来协调多线程和多流的协同计算.

2 HPL 算法和性能模型

2.1 HPL算法简介
HPL 算法通过迭代法求解 N 阶线性方程组 Ax=b.求解过程包含两个步骤,首先通过带行交换的高斯消元法
对系数矩 阵 进 行 LU 分解得到 [ ] [[ , ] ]A b = L U y ,然后进行三角 回 带求解 x.其中 LU 分解的计算量为
2N 3 /3 N− 2 2 / ,三角回带的计算量为 2N 2[1] .给定系统的理论浮点峰值性能和 HPL 的求解时间 T,系统的实测浮
点峰值性能 R max 可以表示为 R max = (2N 3 /3 N− 2 / )2 /T ,HPL 的效率 E=R max /R peak .
⎛ A A 12 ⎞ ⎛ L ⎞ ⎛ U U 12 ⎞ ⎛ L U L U ⎞
A = 11 ⎟ = ⎜ ⎜ 11 ⎜ 11 ⎟ = ⎟ ⎜ 11 11 11 12 ⎟ (1)
⎝ A 21 A 22 ⎠ ⎝ L 21 L 22 ⎠ ⎝ U 22 ⎠ ⎝ L U 11 L U + L U 22 ⎠
12
21
21
22
3
2
HPL 的两个步骤中 LU 分解的计算量为 O(N ),三角回带的计算量为 O(N ),相比于 LU 分解的时间,三角回
带的时间基本可以忽略 [17] ,所以我们的优化也集中在 LU 分解上.式(1)给出了 LU 分解的符号表示,详细的算法
数学证明请参考文献[1].在实现上,LU 分解以 NB 列为迭代步进行迭代求解,算法 1 给出了 LU 分解的详细算法
描述.每一轮迭代包含 4 个子过程,分别是 panel 分解、panel 广播、行交换和尾矩阵更新.其中 panel 分解
(panel_factorization)通过递归高斯消元求解得到 L 11 ,U 11 和 L 21 ;panel 广播(panel_bcast)将 L 11 和 L 21 广播给同行
的行进程;行交换(row_swap)根据 panel 广播收到的行交换信息做行交换;尾矩阵更新(update)首先执行双精度
三角矩阵求逆(DTRSM)得到 U 12 ,然后通过将 L 21 和 U 12 做双精度矩阵乘(DGEMM)更新 A 22 矩阵.
算法 1. 并行 LU 分解算法.
1: for i=0; i<N; i+=NB do
2: panel_factorization(A,i,NB); // recursive panel factorization using Gaussian elimination
3: panel_bcast(L 11 ,L 21 ); // broadcast panel in row processes, including row swap information
4: row_swap(A 12 ,A 22 ); // pivoting according to swap information
5: update(A 22 ,L 21 ,U 12 ); // first calculate U 12 by DTRSM, then update A 22 using DGEMM
6: end for