稀疏三角线性方程组求解(SpTRSV)是迭代方法预条件子步骤的重要组成部分^[1], 其中结构化SpTRSV问题^[2]是一类极具代表性的问题, 在计算电磁学、计算气象学、计算流体动力学等领域^[2-7]的应用十分广泛. 结构化SpTRSV问题耗时较多, 往往是迭代方法预条件子步骤的主要性能瓶颈^[8]. 因此针对结构化SpTRSV问题, 设计高性能并行算法十分必要.

以Summit^[9]及Sierra ^[10]等为代表的超级计算机平台, 采用GPU作为异构加速卡部件加速迭代方法的绝大部分运算操作. 针对GPU平台, 目前AMD的ROCSPARSE^[11]以及NVIDIA的CUSPARSE^[12]等商用高性能数学库, 以分层调度(level scheduling)方法作为并行求解SpTRSV问题的主要手段. 该方法求解SpTRSV问题之前需对输入问题的非零元分布特点进行细致的预处理分析. 然而分层调度方法预处理步骤耗时较多, 导致分层调度方法的整体加速效果出现明显下降^[13-15]. 对此, 部分研究工作^[13-15]采用简易直接的任务划分方式来避免预处理步骤引起的开销, 但由于缺少数据依赖信息的支撑, 此类方法对计算任务的调度规划方面存在一定不足, 在性能上需要做出妥协.

另外, 现有分层调度方法多采用基于逐元素处理的非零元分配策略^{[16, 17]}. 但结构化SpTRSV问题由于非零元数目稀少, 使用逐元素处理策略将导致大量GPU线程资源处于闲置状态, 进一步影响分层调度方法的加速效果. 由于结构化SpTRSV问题的非零元分布与偏微分方程组问题的离散化方法之间关系紧密, 此类问题通常呈现出较为规律的非零元分布结构. 目前已有一些针对特定异构平台的研究^{[2, 3, 18]}, 通过利用这种非零元分布结构, 有针对性地设计了结构化SpTRSV问题的并行算法.

本文在现有工作的基础上, 构造了一种面向GPU平台的高效并行结构化SpTRSV算法, 其实测性能显著优于分层调度方法, 而且该算法无需对输入问题进行任何预处理分析. 本文针对结构化SpTRSV问题非零元数目较少的特性, 对逐元素处理的分配策略进行了改进, 有效缓解了GPU线程资源的闲置问题. 此外, 本文还充分结合GPU体系架构的特点, 对算法的实现进行了深度优化.

1 背景介绍 1.1 结构化稀疏三角线性方程组问题

稀疏三角线性方程组问题要求解的是下列形式:

$Ax=b$

问题的解向量 $x$ , 其中A是一个规模为 $n\times n$ 的稀疏三角矩阵, $b$ 为 $n\times 1$ 的右端项向量. 假定A是下三角矩阵, 若其对角元素不为0, 可通过回代求解步骤:

${x}_{i}\leftarrow {A}_{ii}^{-1}\cdot \left({b}_{i}-{{\sum}}_{0\leqslant j < i}{A}_{ij}\cdot {x}_{j}\right) ,$

计算解向量分量 ${x}_{i}$ 的取值. 为节约存储空间, 稀疏矩阵多以行压缩格式(compressed row, CSR)^[19]存储. CSR格式包括3个数组 $Ai\left[\cdot \right], Aj\left[\cdot \right], Aa\left[\cdot \right]$ , 矩阵第r行包含 $\left(Ai\left[r+1\right]-Ai\left[r\right]\right)$ 个非零元素, 第r行第j个非零元素的列下标和数值分别为 $Aj\left[Ai\left[r\right]+j\right]$ 及 $Aa\left[Ai\left[r\right]+j\right]$ . 本文中等式的索引下标以0起始.

结构化稀疏三角线性方程组主要来源于采用规则网格离散方法(如有限差分方法)求解的偏微分方程组(PDE). 此类问题若采用ILU(k)^[1]等方法构造迭代方法的预条件子矩阵, 在迭代方法预条件阶段计算该矩阵的逆与向量乘积的过程就是典型的结构化SpTRSV问题. 结构化三角线性方程组问题^[2]指问题区域为规则网格的PDE问题, 当采用ILU(k)等矩阵不完全分解方法^[1]生成的近似矩阵, 作为迭代方法的预条件子矩阵时, 所产生的稀疏三角线性方程组问题. 由于每个网格点对应SpTRSV问题中的一个等式, 假设规则网格的维度为 $X\times Y\times Z$ , 记网格点坐标为 $(x, y, {\textit{z}})$ , 方程组等式的行索引为r, 网格点与等式之间存在如下对应关系成立:

$\left\{\begin{array}{l} f\left(r\right)=\left({f}_{x}\left(r\right), {f}_{y}\left(r\right), {f}_{{\textit{z}}}\left(r\right)\right)\\ g\left(x, y, {\textit{z}}\right)=x+y\cdot X+{\textit{z}}\cdot XY \end{array}\right. ,$

其中, ${f}_{x}\left(r\right)=r \mod X$ , ${f}_{y}\left(r\right)=\left\lfloor {\dfrac{r}{X}} \right\rfloor \mod Y$ , ${f}_{\textit{z}}\left(r\right)=\left\lfloor {\dfrac{r}{XY}} \right\rfloor$ .

上式中函数 $f(\cdot )$ 是从等式索引(r)至网格点空间坐标( $x, y, {\textit{z}}$ )之间的映射关系, 函数 $g\left(\cdot \right)$ 是函数 $f(\cdot )$ 的逆映射.

结构化SpTRSV问题的矩阵非零元分布与PDE的离散化格式密切相关, 通常呈现较为规律的结构. 图1展示了与若干结构化SpTRSV问题的非零元分布规律, 图中带有不同形状标记的网格点(等式)均依赖于网格点 $g(x, y, {\textit{z}})$ 所对应等式的解向量分量 ${x}_{g\left(x, y, {\textit{z}}\right)}$ . 对结构化SpTRSV问题而言, PDE的离散化格式以及网格尺寸可视为已知条件.

1.2 GPU体系架构

GPU在采用异构架构的超级计算机平台上是一种中常见的加速卡硬件, 可提供数量庞大的硬件线程资源. GPU线程资源以流处理器(stream multiprocessor, SM)为单元分组, 流处理器内部按32个线程为一个线程束(Warp)^[20], 由硬件负责其执行调度. 流处理器通过总线互联. GPU线程共享全局内存资源, 由于GPU全局内存带宽有限^{[21, 22]}, GPU采用L1+L2的缓存结构, 用于缓冲GPU线程的读写操作, 减少实际产生的全局访存请求^{[23, 24]}. L1级缓存不提供数据一致性保证. GPU配备单独的纹理缓存, 用于缓冲只读数据, 其数据吞吐效率略高于L1级缓存^[20].

目前比较新颖的GPU微架构^{[21, 22]}, 已对访存弱一致性模型^[25]的各项特性提供了充分支持. 该模型允许GPU线程对全局内存读取操作附加“获取(acquiring)”语义^[26]. 经修饰的内存读取操作能够观察到其余GPU线程对访存目标地址的最新改动^[27]. 附带获取语义的读取操作可将目标操作数的最新副本, 推送至发起该操作的GPU线程所属L1级缓存^[26].

在同组线程束内, GPU线程在单位时钟周期内仅能执行一条公共指令. 条件跳转等分支语句可能导致组内线程跳转至不同目标地址, 导致线程束发散(warp divergence), 降低指令流水线的执行效率^[20]. GPU架构的谓词执行(predication)功能可等价实现条件分支指令的功能, 适用于优化包含少量语句的条件体或循环体代码段^[28].

1.3 采用逐元素处理策略的分层调度方法

分层调度算法^[29]是一类基于有向无环图(directed acyclic graph, DAG)拓扑排序操作的并行SpTRSV算法. 该方法首先通过DAG描述方程组等式之间的数据依赖关系, 通过拓扑排序操作筛选出相互独立的方程等式, 将它们分配至GPU线程并行处理. 面向GPU架构的特性, 现有研究^{[16, 17, 30, 31]}已对分层调度算法进行了优化, 在其基础上加入了逐元素策略, 以充分发掘GPU硬件线程层面的并行性(算法1).

算法1. 采用逐元素处理策略的分层调度算法.

Level-scheduling $\left(Ai\left[\cdot \right], Aj\left[\cdot \right], A\left[\cdot \right], b\left[\cdot \right], x\left[\cdot \right]\right)$

1 Generate levels ${S}_{i}\in \left\{S\right\}$ by topological sorting.

2 $P\left[\cdot \right]\leftarrow 0\rhd$ Clear the state variable.

3 foreach level ${S}_{i}\in S$ do

4　 parallel-foreach equation ${s}_{ij}\in {S}_{i}\;\mathbf{d}\mathbf{o}\;\rhd$ at warp level.

5　　 ${\sigma }_{t}\leftarrow 0\; \rhd$ partial sum for each warp thread t.

6　　 parallel-foreach off-diagonal $Aj\left[{e}_{k}\right]\;\text{of}\;{s}_{ij}\;\mathbf{d}\mathbf{o}\;\rhd$ at warp thread level.

7　　　 wait for dependencies $x\text{[}Aj\left[{e}_{k}\right].\rhd$ spinning while-loop on $P\left[Aj\left[{e}_{k}\right]\right].$

8　　　 ${\sigma }_{t}\leftarrow {\sigma }_{t}+Aa\left[{e}_{k}\right]\cdot x\left[Aj\left[{e}_{k}\right]\right]$

9　　 end

10 　 $x\left[{s}_{ij}\right]\leftarrow \mathrm{DIAG}{\left[{s}_{ij}\right]}^{-1}\cdot \left(b\left[{s}_{ij}\right]-\mathrm{WARP}\text{-}\mathrm{REDUCTION}\left[{\sigma }_{t}\right]\right).\rhd$ solve equation s_ij after reducing ${\sigma }_{t}$ .

11　 $P\left[{s}_{ij}\right]\leftarrow 1.\rhd$ update the state variable.

12　 end

13 end

2 基于多级逐元素处理策略的并行结构化SpTRSV算法 2.1 基于结构信息的任务划分方案

本文在结构化SpTRSV问题的网格上进行粗粒度任务划分. 首先, 本文观察到结构化问题中, 行索引落入区间范围 ${Q}_{y{\textit{z}}}:=\left[g\left(0, y, {\textit{z}}\right), g\left(X, y, {\textit{z}}\right)\right)$ 内的等式存在严格的串行数据依赖关系. 因此本文将行索引位于该区间范围内的等式视作一个整体, 构成子任务 ${Q}_{y{\textit{z}}}$ . 本文按 $\textit{z}$ 轴坐标, 将子任务 ${Q}_{y{\textit{z}}}$ 归为 $Z$ 个任务队列 $Q\left({\textit{z}}\right):=\left\{{Q}_{0{\textit{z}}}, {Q}_{1{\textit{z}}}, \ldots , {Q}_{\left(Y-1\right){\textit{z}}}\right\}$ .

本文按任务队列 $Q(\cdot )$ 的 ${\textit{z}}$ 轴坐标增序, 依次向各队列分配线程束资源. 在任务队列 $Q(\cdot )$ 内部, 线程束以动态调度的形式, 从队列中获取要处理的子任务. $y$ 轴坐标较小的子任务具有更高的调度优先级. 本文的并行结构化SpTRSV算法如算法2所示, 其中 $D$ 代表算法向每个任务队列分配的GPU线程束数量.

2.2 多级逐元素处理策略

算法1中GPU线程束每处理完一个等式, 后续非零元素与当前线程束内各线程之间的映射, 都将从标号为0的线程重新开始. 由于结构化SpTRSV问题每个等式所包含的非零元数目较少, 标号较大的线程, 在计算过程中无法获得任何非零元素, 始终处于闲置状态, 造成GPU资源的浪费(图2). 为解决现有逐元素策略的线程闲置问题, 本文修改了上述线程和非零元素的映射关系. 算法2将当前等式的非零元素映射至组内线程之后, 将继续把后续等式的非零元素, 映射至余下处于闲置状态的线程, 这样每个线程都能够获得一定数量非零元素, 从而避免GPU线程闲置(图3).

算法2. 本文提出的并行结构化SpTRSV算法.

Para-Structured-SpTRSV $\left(Ai\left[\cdot \right], Aj\left[\cdot \right], A\left[\cdot \right], b\left[\cdot \right], x\left[\cdot \right], X, Y, Z\right)$

1 $P'\left[\cdot \right]\leftarrow -1\rhd$ clear the state variable.

2 $\textit{z}\leftarrow 0$

3 while $\lfloor \textit{z}/D \rfloor < Z$ do

4　 parallel-foreach task ${Q}_{y{\textit{z}}}\in Q\left({\textit{z}}\right)\;\mathbf{d}\mathbf{o}\;\rhd$ at warp level.

5　　 foreach equation $s\in {Q}_{y{\textit{z}}}$ do

6　　　 wait for dependencies $x\left[Aj\left[{e}_{k}\right]\right]\text{}\rhd$ spinning while-loop testing $P\left[{f}_{y}\left(Aj\left[{e}_{k}\right]\right), {f}_{{\textit{z}}}\left(Aj\left[{e}_{k}\right]\right)\right]\geqslant {f}_{x}\left(Aj\left[{e}_{k}\right]\right).$

7　　　 solve equation $s.\rhd$ with multistage element-based strategy (Section 2.2).

8　　　 $P'\left[{f}_{y}\left(s\right), {f}_{{\textit{z}}}\left(s\right)\right]\leftarrow s.\rhd$ for cache footprint reduction (Section 2.3).

9　　 end

10　 end

11　 $\textit{z} \leftarrow \textit{z}+1$

12 end

GPU线程的部分和 ${\sigma }_{t}$ 的计算:

${\sigma }_{t}\leftarrow {\sigma }_{t}+Aa\left[{e}_{k}\right]\cdot x\left[Aj\left[{e}_{k}\right]\right] .$

可拆分为2个阶段. 在第1阶段, 线程将会读取非零元信息 $Aj\left[\cdot \right]$ 和 $Aa[\cdot ]$ , 在第2阶段才会进行右端项 $x\left[\cdot \right]$ 相关的计算和内存操作. 前一阶段的全局内存访问的耗时较长. 本文采用多级处理的方式来隐藏访存开销. 多级处理在现有逐元素处理策略的基础上, 加入了非零元信息的预取步骤. 本文利用少量GPU寄存器资源, 构造了一个微型数据缓冲队列, 用于容纳预取步骤的操作数并记录操作数的下标. 由线程束内部线程根据下标决定, 是否从缓冲队列中提取非零元信息, 并进行第2阶段的计算. GPU线程从队列中取出非零元素后, 将发起下一批次非零元信息的预取操作. 因为非零元信息的预取操作与问题等式间数据依赖无关, 所以多级处理的非零元素访存和计算可重叠进行, 能够隐藏非零元信息的访存延迟.

2.3 状态变量压缩技术

SpTRSV算法使用若干变量 $P\left[\cdot \right]$ 来标记等式的状态. SpTRSV计算开始前, $P\left[\cdot \right]$ 的初始值为0, 算法得到分量 $x\left[s\right]$ 之后, 会将 $P\left[s\right]$ 的值设置成1. GPU线程在访问分量 $x\left[\cdot \right]$ 之前, 会通过检测 $P\left[\cdot \right]$ 是否为1 (算法1第2行), 来推断 $x\left[\cdot \right]$ 当前是否处于可用状态, 以维持SpTRSV问题的数据依赖关系. 当大量GPU线程以循环形式(算法1第7行), 反复访问状态变量时, 将引起频繁的缓存缺失, 进而产生大量全局内存操作请求. 而这些访存操作将与矩阵非零元等关键访存操作, 争抢有限的GPU访存带宽资源.

为解决访存带宽争用的问题, 本文采用一个维度为 $Y\times Z$ 的二维数组 $P\left[\cdot , \cdot \right]$ 来记录 $x\left[\cdot \right]$ 的状态(算法2第6行), $P\left[\cdot , \cdot \right]$ 初始值置为−1 (算法2第8行). 算法解出分量 $x\left[s\right]$ 之后, 会将状态变量 $P\left[{f}_{y}\left(s\right), {f}_{\textit{z}}\left(s\right)\right]$ 的值设置为 $s$ . GPU线程在读取 $x\left[s\right]$ 之前, 将通过测试如下条件:

$P\left[{f}_{y}\left(s\right), {f}_{\textit{z}}\left(s\right)\right]\geqslant s$

判断 $x\left[s\right]$ 当前的可用状态. 本文提出的状态标记方法将状态变量的总数据量, 压缩至原有标记方法的 $1/X$ , 减少了状态变量在算法工作集数据中的占比, 有效提高了状态变量操作的缓存命中率(第3.4节). 因此算法得以大幅减少和状态变量 ${P}'\left[\cdot \right]$ 有关的全局内存请求, 节省GPU的访存带宽资源.

2.4 针对GPU体系架构的进一步优化

处理子任务 ${Q}_{y{\textit{z}}}$ 的for-each循环体(算法2第5–9行)是热点代码段. 在现有方法^[13-17]中线程束每次获取的计算任务较少, 对于循环内部非零元素计算进行优化的预期收益较少, 而且非零元素计算效率上的不足可以通过GPU硬件高并发的特性进行弥补, 因此循环体的优化对整体性能的影响十分有限. 本文提出的算法计算任务 ${Q}_{y{\textit{z}}}$ 的粒度较大. 受限于GPU的硬件设计, 线程束的指令/访存吞吐效率不高, 所以处理非零元素的循环体就会成为算法性能的主要瓶颈. 本文结合NVIDIA V100 GPU的硬件特性, 对该循环体进行了优化.

首先, GPU线程通过前文构造的寄存器缓冲队列, 获取矩阵非零元的信息( $Aj\left[\cdot \right], Aa[\cdot ]$ ). 而且由于矩阵非零元信息在计算过程中保持不变, 本文将缓冲队列预取操作中出现的全局内存取数指令( $\mathrm{l}\mathrm{d}.\mathrm{g}\mathrm{l}\mathrm{o}\mathrm{b}\mathrm{a}\mathrm{l}.\mathrm{*}$ ), 替换成了相应无数据一致性保证的版本( $\mathrm{l}\mathrm{d}.\mathrm{g}\mathrm{l}\mathrm{o}\mathrm{b}\mathrm{a}\mathrm{l}.\mathrm{n}\mathrm{c}.\mathrm{*}$ ), 把与 $Aj\left[\cdot \right]$ 以及 $Aa[\cdot ]$ 有关的访存操作, 交由GPU纹理缓存进行缓冲, 进一步提高非零元信息的访存吞吐效率.

其次, GPU线程需要根据状态变量的数值推断其余线程的计算进度, 因此GPU线程需共享压缩后的状态变量 $P\left[\cdot , \cdot \right]$ , 并且要能够观察到任意线程对 $P\left[\cdot , \cdot \right]$ 的最新改动. 由于L1级缓存不提供数据一致性保证, 在GPU的全局内存访存指令当中, 针对L2级缓存的访存指令(ld.global.cg.*)是满足上述需求的一个选项^[27], 但ld.global.cg.*访存指令的操作延迟较高, GPU线程从修改到读取 $P\left[\cdot , \cdot \right]$ 的时间间隔至少为250个GPU时钟周期(图4实线路径). 附带“获取”语义的访存操作(ld.require.gpu.*)也可用于实现与ld.global.cg.*指令类似的功能^[27], 且使用d.require.gpu.*指令(图4虚线路径)对 $P\left[\cdot , \cdot \right]$ 进行操作的延迟(140个GPU时钟周期)显著少于ld.global.cg.*指令, 因此本文采用该指令对 $P\left[\cdot , \cdot \right]$ 进行访问. 状态变量操作的高缓存命中率也让ld.require.gpu.*指令的执行效率能够维持在较为稳定的状态. 需要说明的是结合ld. require. gpu.*指令的操作延迟数据以及NVIDIA V100 GPU的体系架构特点分析, 本文推断流处理器的总线在数据推送过程中极有可能承担了数据转发的角色.

最后, 在多级逐元素处理以及寄存器缓冲队列的预取操作中, GPU线程需要根据当前非零元的索引以及列下标 $Aj\left[\cdot \right]$ , 决定是否发起新一轮数据预取, 或者判断是否需要进行一些与对角元素相关的额外操作. 实现上述功能的条件判断语句, 将在循环体内部引入众多分支指令, 导致线程在执行过程中经常出现取值地址发散的情况. 本文采用谓词执行功能消除了这些分支指令, 以缓解取值地址发散导致的性能损耗问题. 本文将条件判断语句的真值结果, 存入GPU的谓词寄存器^[27]当中, 然后给分支语句的每一条指令添加相应的标记信息, 其中包含用于控制指令执行的谓词寄存器的索引. 被标记的指令总会进入GPU指令流水线, 因此不会造成取值地址发散的问题. 如果相应谓词寄存器存储的状态不为真, 那么GPU指令流水线不会提交被标记指令的执行结果, 该指令将不会对算法2的状态产生任何影响.

3 实验验证 3.1 实验设置

本文在NVIDIA V100 GPU平台上, 对算法性能进行测试评估, 实验平台的软硬件环境配置如表1所示. 由于子任务数据依赖关系会随结构化SpTRSV问题的类型发生变化, 本文发现算法2的实测性能与任务划分阶段, 算法给任务队列所分配的线程束数量之间存在关联. 为保证算法2能够达到最佳的实验效果, 本文依据结构化SpTRSV问题的非零元分布结构, 调整了算法向任务队列分配的线程束数量, 表2列举了算法参数D在各类测试问题中的取值.

3.2 性能对比实验结果

本文按照图1所示的非零元分布结构, 生成D3N7等4种类型的结构化SpTRSV问题, 以实际验证算法2的加速效果. 参照文献[2], 本文采用有效访存带宽作为算法性能的评判标准, 有效访存带宽定义如下:

$\mathrm{有}\mathrm{效}\mathrm{访}\mathrm{存}\mathrm{带}\mathrm{宽}:=\frac{\left(\mathrm{稀}\mathrm{疏}\mathrm{矩}\mathrm{阵}\mathrm{数}\mathrm{据}\mathrm{总}\mathrm{量}+\mathrm{右}\mathrm{端}\mathrm{项}\mathrm{数}\mathrm{据}\mathrm{总}\mathrm{量}\times 2\right)}{\mathrm{S}\mathrm{p}\mathrm{T}\mathrm{R}\mathrm{S}\mathrm{V}\mathrm{计}\mathrm{算}\mathrm{时}\mathrm{间}}.$

测试结构化SpTRSV问题的规模从64³逐渐增长至256³. 本文同时还测试了CUSPARSE数学库 (csrsv2函数)^[12]以及现有同步避免(synchronization free)方法的性能, 作为算法2加速效果的对照. 测试中所采用的同步避免方法, 包括Liu等人提出的原始同步避免算法^{[13, 14]} , 以及最近由Su等人改进的Capellini框架^[15]. Capellini框架融合了一系列SpTRSV并行算法, 本文呈现的有效带宽数据为该框架的最优性能结果. 实验测量数据取10次运行的平均值. 性能实验的结果如后文图5所示, 算法2的有效访存带宽为40.1–225.2 GB/s, 与问题规模以及等式所包含的非零元数目呈正相关. 算法2相比CUSPASRE的平均加速比为2.71倍, 对比同步避免方法的加速可达15.1倍, 性能较已有方法有大幅度提升.

3.3 性能优化效果分析

本文通过对比实验评估前文所述各项优化手段的技术效果. 实验的基线(baseline)版本仅使用第2.1节所描述的任务划分策略, 不包含任何优化措施. 多级处理版本(multistage)在基线版本的基础上, 加入了多级逐元素处理策略, 完整版本(compression & others)包含所有优化手段. 实验测试了以上各版本, 处理D3N33问题时的性能数据. 图6是实验中收集的有效访存带宽数据, 图中标注了各项优化手段对有效访存带宽的提升效果. 多级逐元素处理策略以及谓词执行等后续优化手段, 在小规模问题中的性能增益有限. 对于大规模结构化SpTRSV问题, 多级逐元素处理策略的有效访存带宽, 相较基线版本有68%左右的提升, 若使用状态变量压缩技术以及第2.4节所述的优化手段, 算法性能还可进一步提升约25%. 两种优化手段叠加之后, 基线版本有效访存带宽的提升幅度约为115%.

3.4 状态变量压缩效果评估

由于算法2采用动态任务调度, 因此本文无法直接通过性能监测工具测量算法的缓存命中率数据. 本文通过测量算法2, 在使用或未使用状态变量压缩技术的情况下, GPU线程状态变量操作所消耗的时间来间接推测其缓存命中情况. 图7为在D3N33类型的结构化SpTRSV问题中, GPU线程在状态变量操作上消耗的时间. 采用状态变量压缩技术并启用GPU缓存之后(compressed, cached), 状态变量操作的用时相较禁用GPU缓存的情况(compressed, no caching)有显著下降, 显示绝大部分与 $P\left[\cdot , \cdot \right]$ 相关的访存操作, 均有效命中GPU的缓存. 如不采用状态变量压缩技术(uncompressed, no caching/cached), 那么GPU缓存对状态变量操作的时间开销无明显影响, 说明相关操作的GPU缓存命中率较低. 本实验说明采用状态变量压缩技术, 可有效提升GPU缓存命中率, 显著减少GPU线程状态变量操作的耗时, 减少状态变量操作占用的GPU访存带宽.

4 相关工作总结

分层调度方法是目前运用范围最为广泛的并行SpTRSV算法, 但该方法预处理开销是算法的主要性能瓶颈^{[16, 17]}. 分层调度方法的早期版本需通过全局线程同步来维持等式间的数据依赖关系, 同步操作严重影响算法性能^[32], 现有研究也提出了一系列方法^[13-15], 来消除线程同步带来的额外开销.

除分层调度算法等基于CSR格式的并行SpTRSV算法以外, 现有研究也提供了许多针对SpTRSV问题的专用稀疏数据结构, 采用多种存储方式混合形式, 记录稀疏矩阵的非零元分布结构^[33,34], 为并行算法的设计提供便利条件. 在特定场景中, 预条件子步骤可采用SpTRSV问题的近似解替换精确解, 着色算法(coloring)^[35-37]是一种比较常用的近似方法. 该方法通过对SpTRSV问题的DAG进行着色, 将等式分为若干集合. 着色算法将忽略集合间的等式数据依赖关系, 以计算得到SpTRSV问题的近似解. 从理论的角度, SpTRSV问题也存在其他求解方法^[38-40].

5 结　论

本文面向GPU平台, 针对结构化SpTRSV问题, 提出了一种无需预处理步骤的结构化SpTRSV问题并行算法. 该算法采用多级逐元素处理策略, 在隐藏矩阵元素访存延迟的同时, 有效缓解了结构化SpTRSV问题导致的GPU线程闲置问题. 此外, 本文还引入了状态压缩技术, 并结合GPU的硬件特性进一步改善算法的性能.

为摆脱现有单片封装工艺的技术限制, 目前最新一代^{[22, 41]}的GPU已开始采用多片封装技术(multi-chip module, MCM)^{[42, 43]}, 将多个GPU模块整合到一起来提高GPU集成的功能单元数量. 虽然MCM模式提供了更多访存及计算资源, 但GPU模块间存在与非均衡内存访问(non-uniform memory access, NUMA)类似的性能损耗现象^[44], 这给并行SpTRSV算法的任务划分以及访存优化提出了新的挑战. 如何克服MCM封装模式的不均衡效应, 找到行之有效的并行SpTRSV算法, 是本文未来需要关注的一个研究方向.