FILIPPONE S, CARDELLINI V, BARBIERI D, 等. Sparse Matrix-Vector Multiplication on GPGPUs[J/OL]. ACM Transactions on Mathematical Software, 2017, 43(4): 1-49. DOI:10.1145/3017994.

本篇文章为2017年的工作，主要讨论SpMV在GPGPU上的各类实现，包括存储格式和计算Kernel的书写、对各种格式的分类比较、多种格式的融合使用、性能分析以及自动调优的框架等等。

SpMV应用广泛的原因在于，稀疏的线性问题Ax=b求解，使用LRU分解的方式会引入额外的非零元，使得最终的计算量增加，而Krylov 提出了一种subspace projection methods，可以通过迭代进行Ax乘法计算解，而不破坏稀疏的结构。

许多的存储格式都没有直接的可以使用的代码，也就是说缺乏实际可以使用的生态，更别说和目前已有的各类框架和工具进行整合，从而投入真正的生产了。

对于针对GPGPU所提出的新表示格式和新的算法，主要考虑的点有：

  • memory footprint
  • atomic data access
  • coalesced access
  • thread mapping
  • dara reuse
  • zero padding

一般的稀疏矩阵计算格式

COO

三元组，行列索引加上对应的值

COO最为直观，同时存储开销只由非零元个数决定而与稀疏结构无关。

但是性能表现不佳，因为无法从COO的数据格式中得到矩阵的结构，也就无法将计算任务拆分给不同的线程。由于最终结果需要进行reduction，最终的结果也就无法知道要从哪些计算合并得来。

文章介绍了一些在COO基础上的改进方法，大概方法包括压缩和分块以及重排。

CSR

行坐标压缩，以row_pointer指定每行起始非零元的位置

CSR是一种非常灵活且同样的格式。这种灵活性的代价就是在将计算分散到不同的线程时会引入相应的开销。包括lack of coalescing、inefficient memory access pattern、load imbalance、thread divergence。

处于不破坏通用性考虑，，大多数工作都专注在对CSR的算法进行优化，这一部分工作可以归为针对GPGPU硬件的软件调优。优化的工作主要是根据上面会引入开销的那几个因素，进行性能调优。

最初使用CSR格式的算法包括两个，scalar和vector：

  • scalar : 每行分配一个线程

  • vector : 每行分配多个线程

    • 通过开发coalescing可以实现更好的并行度
    • 一类工作是研究每行分配的最优线程个数（依据矩阵的平均非零元个数）

除了针对那几个方面的性能调优外，还有一个工作是先做元素级别的乘法，然后再合并成最终的结果向量。

基于CSR的新格式也是针对某几个方面提出涉及到存储表达的改动和优化。

CSC

列坐标压缩，以col_pointer指定每列起始非零元的位置

对于SpMV计算来说，CSC相比CSR拥有更加不规则的存储访问。因此大多数时候不使用，但是常用于稀疏矩阵分解方面。

ELL

按行压缩，每行的长度固定，填0补齐

ELL对于类似GPGPU这样多计算核心的架构来说，在规则化方面具有天生的优势。但是对于稀疏分布不规则的矩阵，需要补0会变得很多，增加存储和计算开销。

相关的研究集中在减少存储开销，提高对于通用稀疏矩阵计算的效率上。

  • ELLPACK-R: 在ELL的基础上额外存储每行非零元，在load数据时ELL格式能够利用coalesced access模式，额外的非零元的长度可以用来减少无效的计算。

  • ELLPACK-T : 在ELLPACK-R的基础上对每行的非零元以及列号进行重排(可能会换行)，保证每行非零元素个数为16的倍数，以获得更好的coalesced粒度和对global memory的对齐的访问。

    • ELL-RT : 可以在ELL-R和ELL-T之间进行auto-tuning
    • Sliced ELL-RT : 在ELL-RT的基础上加上行方向上的slice

  • Sliced ELLPACK (SELL-C) : 对非零元进行重排和分割（slice 是指对行进行切分），再按块组织成ELL

    • SELL-C-sigma : 相近非零元个数的行被组织在一块slice内， sigma表示重排的范围（可以避免全局重排的开销）
  • BiELL : 在SELL-C使用BiSection对columns进行分割重组，目的是提升load balance
  • BELLPACK : 探索利用稀疏矩阵中的块结构，作为ELLPACK的一种显式存储矩阵稠密部分的格式扩展）

  • AdELL : 采用一种启发式的方法，每一行可以依据计算任务的大小而分配不同数量的线程

    • CoAdELL : 在AdELL的基础上，对连续的非零元列索引进行无损压缩
    • AdELL+ : 优化AdELL上用的启发式方法，同时提出了一种在线自动优化的方法
  • JAD(JDS) : 按行内非零元数量进行行排序，之后再分块存成ELL，这样每个块内需要填充0的个数不定（数量更少），但是可能会存在硬件利用率问题
  • ELLPACK-RP : ELLPACK-R + JAD

  • pJDS : JAD + ELLPACK-R，适合unstructure matrix

    • ELL-WARP : pJDS变种，针对有限元素的unstructure grids
  • BRC: Blocked Row-Column 说是结合了JDS的行变换和ELLPACK的补齐机制（没太明白）。已经在PETSc库里使用，可以作为竞品对象
  • BiJAD : 类似BiELL的Bisection技术，用在JAD上
  • TJAD : Transpose Jagged Diagonal 对列做排序

  • 文章自己提出了两种格式：类ELLPACK-R和HLL

    • 在PSBLAS库中可用

DIA

按照对角线将非零元排成一列或者一行

DIA适合天生具有自然的对角线结构的矩阵

HYB

COO和ELL的组合，适合没有一个一般性统一结构的矩阵，主要应用在GPGPU上。

此外还包括CSR和ELL的混合（HEC），根据占据主导的不同，可以适应不同的应用矩阵；
DIA和ELL的混合，DIA和CSR的混合（HDC），BCOO和BCSR的混合

新型的格式

  • BLSI : 减少数据组织和CPU到GPU传输的时间
  • CRSD: 解决DIA的zero fill-in问题

自动调优方法

主要可以分为三类

  • reorganizing for efficient parallelization
  • reducing memory traffic
  • orchestrating data movement
  • 其他： include, among the others, fine-grain parallelism, segmented scan, index compression, and register blocking

但是从结果上看，并没有一个通用最优的方案，因此需要根据实际运行情况和输入矩阵的参数选择最合适的计算方案。