TileSpTRSV：一个在GPU上的并行分块稀疏三角解算法

中文题目：TileSpTRSV：一个在GPU上的并行分块稀疏三角解算法

论文题目：TileSpTRSV: a tiled algorithm for parallel sparse triangular solve on GPUs

录用期刊/会议：CCF Transactions on High Performance Computing (CCF THPC, CCF-C类期刊)

原文链接：https://www.ssslab.cn/assets/papers/2023-lu-THPC.pdf

录用/见刊时间： 2023-03-03（录用时间）

封面图片：

封面摘要：

稀疏三角解（SpTRSV）是稀疏基本线性代数子程序（BLAS）中最重要的二级内核之一。目前，很多相关工作都集中在减少SpTRSV的层次集算法和无同步算法中的依赖性和同步性上。然而，能在GPU上充分利用稀疏结构实现SpTRSV的工作却较少。本文提出了一种叫做TileSpTRSV的分块算法，其通过利用稀疏矩阵的二维空间结构来优化GPU上的SpTRSV。本文分别基于层次集算法和无同步算法为其设计了两种实现，即TileSpTRSV_level-set和 TileSpTRSV_sync-free。在最新的NVIDIA GPU上的实验结果表明，我们的TileSpTRSV_level-set算法分别对比cuSPARSE、Sync-free和Recblock算法的平均加速比可达5.29 倍（最高 38.10 倍）、5.33 倍（最高 21.32 倍）和 2.62 倍（最高 12.87 倍）。

作者列表：

1） 芦正阳 中国石油大学（北京） 信息科学与工程学院 计算机科学与技术系 硕21

2） 刘伟峰 中国石油大学（北京） 信息科学与工程学院 计算机科学与技术系教师

背景与动机:

稀疏基本线性代数子程序（Sparse BLAS）用于对稀疏向量或稀疏矩阵进行运算，其在数值计算、机器学习、应用数学等领域都有着广泛应用。稀疏三角解（SpTRSV）是Sparse BLAS中最重要的二级子程序之一。与另一个二级稀疏BLAS子程序稀疏矩阵-向量乘（SpMV）相比，SpTRSV 具有天然的串行计算模式，其一般来说更难在多核处理器（如 GPU）上找到较高的并行度。如今，很多 SpTRSV 的优化工作都集中在减少层次集和无同步算法中的依赖关系和同步代价上。然而，在GPU 上能够很好地利用输入矩阵的稀疏结构进行 SpTRSV 计算的工作却比较少，因此这激励我们提出了一种分块的并行稀疏三角解算法。

设计与实现:

为更好的利用输入矩阵的稀疏特性，本文在GPU上实现一种名为TileSpTRSV的分块算法，通过利用稀疏矩阵的二维空间结构来优化SpTRSV 计算。

(1) 预处理阶段

TileSpTRSV 在预处理阶段首先将输入矩阵划分为多个16乘16大小的非零稀疏块并根据稀疏块的位置及其内部非零元分布特征将这些稀疏块匹配到不同的稀疏存储格式中。我们的TileSpTRSV使用两级存储结构来存储输入稀疏三角形矩阵的信息。图1显示了一个例子，输入矩阵被分成10个稀疏块，并存储在一个包含多个数组的两级存储结构中。此外，我们还为每种稀疏存储格式设计了对应的块级别的GPU核函数以进行块内的计算。

图1. 一个输入矩阵经过分块及格式匹配后的例子

(2) 基于层次集方法的算法实现

我们首先在层次集算法的基础上实现了一种名为TileSpTRSV_level-set的分块SpTRSV 算法。该算法先对输入矩阵按其层次集顺序进行排序随后将分块后的矩阵块划分为多个层次集合，每个层次集合内的块可以进行并行计算，而层次集合之间需要按照层次顺序串行计算。注意到，每个层次集合内会首先对非对角块进行计算随后才会对所有对角块进行计算。此外，我们还为TileSpTRSV_level-set算法设计了warp-level和thread-level的SpTRSV核函数和一种自适应SpTRSV内核选择策略，其会根据输入矩阵的块级别并行性以调用不同的SpTRSV 算法。

图 2. 基于层次集方法的TileSpTRSV算法计算过程

(2) 基于无同步方法的算法实现

我们还在无同步算法的基础上实现了一种名为TileSpTRSV_sync-free的分块SpTRSV 算法。该算法同样会先将输入矩阵划分为多个16乘16大小的块，并会将每个块视为一个单元以使用无同步算法进行计算。同样地，在具体计算每个块时，我们会根据其稀疏存储格式来匹配不同的GPU核函数来进行计算，图3展示了TileSpTRSV_sync-free算法的计算过程。

图 3. 基于无同步方法的TileSpTRSV算法计算过程

实验结果及分析:

通过分析在NVIDIA GPU上测试16个代表性矩阵的测试数据，结果表明与cuSPARSE、Sync-free和Recblock算法相比，TileSpTRSV-levelset算法在这16个代表性矩阵上的平均加速分别可以达到5.29倍(最高38.10倍)、5.33倍(最高21.32倍)和2.62倍(最高12.87倍)。

图4. TileSpTRSV和其他三种SpTRSV算法在16个代表性矩阵上的性能比较

结论:

在本文在NVIDIA GPU上分别实现了两个版本的TileSpTRSV算法，即在层次集算法基础上实现的TileSpTRSV_level-set，以及在无同步算法基础上实现的TileSpTRSV_sync-free。我们还为TileSpTRSV_level-set算法设计了自适应SpTRSV内核选择策略，根据输入矩阵的块级别并行性使用不同的SpTRSV算法。实验结果表明我们的分块算法实现相对于其他SpTRSV算法来说有着较为明显的性能优势。

作者简介:

芦正阳，中国石油大学（北京）硕士三年级学生。研究方向为高性能计算，稀疏矩阵计算。联系方式：zhengyang.lu@student.cup.edu.cn

通讯作者简介:

刘伟峰，博士、中国石油大学（北京）教授、博士生导师，欧盟玛丽居里学者。2002年和2006年于中国石油大学（北京）计算机系获学士与硕士学位。2006年至2012年在中国石化石油勘探开发研究院历任助理工程师、工程师和高级研究师，其间主要研究领域为石油地球物理勘探的高性能算法。2016年于丹麦哥本哈根大学获计算科学博士学位，主要研究方向为数值线性代数和并行计算，其中尤其关注稀疏矩阵的数据结构、并行算法和软件。研究工作发表于SC、PPoPP、DAC、ASPLOS、ICS、IPDPS、ICPP、TPDS、JPDC、FGCS和Parco等重要国际会议和期刊。担任TPDS、SISC和TKDE等多个重要国际期刊审稿人，以及SC、ICS、IPDPS和ICPP等多个重要国际会议的程序委员会委员。他是IEEE高级会员、CCF高级会员、ACM和SIAM会员。

联系方式：weifeng.liu@cup.edu.cn