《高性能科学与工程计算》——3.6　案例分析：稀疏矩阵向量乘

本节书摘来自华章计算机《高性能科学与工程计算》一书中的第3章，第3.6节,作者:(德）Georg Hager Gerhard Wellein 更多章节内容可以访问云栖社区“华章计算机”公众号查看。

3.6　案例分析：稀疏矩阵向量乘

上节讨论的稀疏矩阵向量乘是采用循环分块和展开优化方法的一个实际应用。该算法是大多数迭代矩阵对角化算法（Lanczos、Davidson、Jacobi-Davidson）的重要部分，而且经常会成为性能限制因素。当矩阵的非零元素数量Nnz随矩阵行数Nr的增加而线性增长时，这个矩阵称为稀疏矩阵。出于效率方面的考虑，内存中只存储稀疏矩阵的非零元素。因此，稀疏MVM（sMVM）是一个O(Nr)/O(Nr)问题，并且当Nr取值非常大时，该算法是访存受限的。尽管如此，循环嵌套使算法具有显著的优化潜力。图3-15所示的sMVM算法中，尽管矩阵的访存模式是非常有利的，但是RHS向量经常是非连续访问甚至是间接寻址的。下面章节将会进行不失一般性的讨论。

3.6.1　稀疏矩阵的存储机制
目前有多种稀疏矩阵存储机制，其中一些机制只适合存储特定类型的矩阵[N49]。内存的访问模式和sMVM的性能特征严重依赖于所使用的存储机制。CRS（Compressed Row Storage，压缩行存储）和JDS（Jagged Diagonals Storage，锯齿对角线存储）是目前两个最重要也是最常用的存储机制。下面的讨论中，我们将会看到：CRS适合基于cache的微处理器架构，而JDS则支持依赖和循环结构，该结构非常适用于向量化系统。
CRS存储机制使用长度为Nnz的数组val存储矩阵所有的非零元素（按行存储，且元素间没有间隔）。所以必须提供两个额外的整型数组来指定数组val元素在原矩阵中所属的行、列：长度为Nnz的数组col_idx和长度为Nr的数组row_ptr。前者存储了每个非零元素所属的列，后者存储了每行开始的索引（使用非零元素在val数组的下标）（参见图3-16）。使用该存储机制完成MVM的初始代码非常简单：

下面几点需要引起注意：

• 外层循环迭代次数非常多（迭代次数为Nr）。
• 相对于一般微处理器架构的流水线长度，内存循环迭代次数可能会比较少。
• 向量C的访存已经高度优化了：只从内存中加载一次。
• 数组val中非零元素的访存是连续的。
• 同期望的一样，RHS向量B的访存模式为间接地址映射。然而，依赖于矩阵结构，这可能不是一个严重的性能问题。如果非零元素主要集中在对角线周围，则甚至会有明显的空间和时间局部性。
• Bc = 5/4 W/F（如果col_idx是4字节整型）。因为cache行的部分使用忽略了可能的大规模的数据传输。

其中有些点当我们后面演示并行SMVM时（参见7.3节）将很重要。
JDS不仅要完成矩阵元素的消零，还需对矩阵元素进行重新组织。首先，移除矩阵中所有的零元素并将非零元素移到左边；其次将矩阵行按照非零元素的数目排序，使非零元素最多的行位于矩阵顶端，非零元素最少的行位于矩阵底部。矩阵行排序操作同时也生成置换映射数组并存储在长度为Nr数组perm中；最后，将矩阵的非零元素按列存储在数组val中。这些列也称为锯齿对角线（jagged diagonal，从原始稀疏矩阵的左上部遍历到右下部，参见图3-17）。每一个非零元素在原矩阵的列索引像CRS一样存储在数组col_idx中。为使RHS和LHS向量元素的顺序相同，col_idx数组根据置换映射数组重新生成。数组jd_ptr保存了

第Nj条锯齿对角线的起始索引。采用JDS存储机制的sMVM的基本实现代码只比采用CRS稍微复杂一点：

数组perm在这里没有被用到；通常情况下，所有sMVM操作都需要在置换空间内完成。这个循环有以下值得注意的属性：

• 内层循环的迭代次数非常多，且互相没有依赖。对于向量处理器来说，JDS是比CRS更好的数据存储机制。
• 外层循环迭代次数较少（锯齿对角线的数目）。
• 结果向量从内存中被加载多次（至少是部分向量），所以这里可能有比较大的优化空间。
• 数组val中的非零元素是连续访存的。
• 同CRS一样，RHS向量访存为间接地址映射。虽然一个良好的矩阵布局是使用直对角线而不是压缩行，但上述的注释依然可用。作为额外的复杂度，矩阵行和RHS向量都被置换了。
• Bc = 9/4 W/F（如果col_idx是4字节整型）。

从代码平衡值来看，sMVM似乎更倾向于CRS。
3.6.2　JDS sMVM优化
JDS sMVM要用到循环展开与合并技术，但是这个技术需要内层循环的长度独立于外层循环索引。不幸的是，锯齿对角线一般不会具有相同的长度，违反了这个条件。然而，一个称为循环剥离（loop peeling）的技术可解决这个问题：对于m路的循环展开，分割成m×x个块，剩下的少于m-1个对角线可单独处理（参见图3-18，像往常一样剩余循环被忽略）。

假定被剥离迭代的时间消耗可以忽略不计，则m路循环展开与合并将代码平衡值降低为：

如果m取值足够大，则JDS代码平衡值就会接近于CRS。然而，如之前讨论的那样，较大的m值会导致寄存器使用量的增多，因此这一做法并不总是可取。通常情况下，一个合理的循环展开结合分块方法可以用来减少内存数据传输，同时提高cache内部性能。循环分块方法对于JDS sMVM来说也适用（见图3-19）：

应用这个优化方法，当分块大小b取值不太大时，结果向量只需从内存中加载一次。尽管代码平衡值没有发生变化，但此时应该能够获得同CRS相近的性能。同之前论述的稠密矩阵转置一样，矩阵分块并没有优化寄存器的重用而是优化了cache的利用率。
图3-20显示了CRS与JDS sMVM（原始、两路循环展开、大小为400的分块）在三个不同硬件平台上的性能对比。测试矩阵来源于固态物理学（半填充情形下六角一维Holstein-Hubbard模型）。CRS似乎更适用于标准AMD和英特尔微处理器。这并不奇怪，因为它不需要手工优化就具有最低的代码平衡值，并且迭代次数较少的内层循环比较适合具有乱序执行能力的CPU。然而，采用EPIC架构的英特尔Itanium2处理器[V113]对于CRS来说性能表现平平，但在分块的JDS版本中，其性能最高。因为缺乏乱序处理能力，编译器虽然检测到内层循环所有的指令集并行，但不能重叠一行的wind-down阶段和另一行的wind-up阶段，所以这个架构不能很好应付CRS中的短循环。当工作集可以加载到cache中时，这个效果会更加明显[O56]。