1、在 CUDA 上实现有效的稀疏矩阵乘法摘要图像处理单元（GPU）的大规模并行性在许多高性能计算应用中提供了巨大的业绩。当稠密线性代数绘制到这样的平台上，利用这种可能性来进行稀疏矩阵的计算就呈现出额外的许多挑战。其在用迭代方法解决稀疏线性系统和特征值问题，以及稀疏矩阵向量乘法方面的地位在稀疏线性代数中具有独特的重要性。在本文中，我们讨论用于 SPMV（稀疏矩阵向量乘法）的数据结构和算法，使 SPMV在适合 GPU 的细粒度并行结构的 CUDA 技术平台上可以有效的执行。基于 SPMV 的有限存储特性，我们强调有效的内存带宽和紧凑的存储格式。我们认为广泛的稀疏矩阵来自那些结构严谨的规则矩阵和那些经
2、常在矩阵的每行带有不平衡的非零值分布的高度不规则矩阵。我们开发一些方法来探究许多常见的矩阵结构形式，同时也提供一些备选方案来适应更多更大的不规则矩阵。从结构上说，在一个 GeForce GTX 280 GPU 平台上，基于网格的矩阵，在单精度模式中我们可以获得每秒 36 亿次的性能，而在双精度模式中我们可以获得每秒 16 亿次的性能。对于没有结构的有限元素矩阵，我们通过分别进行超过每秒 15 亿次和 10 亿次的单双精度测试来观察它们的性能。这些结果比起之前用传统的多核处理器进行的关于 SPMV 的state-of-the-art 研究方法要好的多。我们的双精度 SPMV 性能大体上是一个 C
3、ell BE with 8 SPEs 系统的 2.5 倍，是一个四核的 INTELClovertown 系统的 10 倍还多。1. 介绍稀疏矩阵结构出现在许多计算学科中，因此评价一种方法能否有效地操作这些矩阵往往取决于他们在许多应用中的性能如何，稀疏矩阵乘法（SPMV）运算已经被证明在计算科学领域具有非常重要的意义。在许多用于科学工程应用中用来解决大规模线性系统和特征值问题的方法中，SPMV 运算呈现出突出的价值。这些迭代方法剩下的部分可以典型地概括为稠密线性代数操作，而这些操作已经被 BLAS 和 LAPACK 尽可能有效地解决。最新的 NVIDIA GPUS 适应高吞吐量的多核处理能够提供
4、非常高的计算吞吐量峰值，实现这种可能性需要执行大量的细粒度并行性和结构化计算来实现执行通络和内存存取方式的足够规律性。最近，Volkov and Demmel and Barrachina 已经证明了如何通过线性矩阵操作来获得吞吐量和带宽的峰值点。稠密矩阵操作是非常普遍的，因此也经常被限制于吞吐量浮点（floating point throughput）。相比之下，稀疏矩阵操作在存取方式中更普遍，因此其速度一般限于纯粹的带宽。在本文中，我们研究用于 GPU 之类的多核处理器上的有效的 SPMV 内核的设计。我们用 CUDA 执行这些内核程序并在 GeForce GTX 280GPU 上分析他们
5、的性能。出现在不同问题中的稀疏矩阵呈现出广泛的规律性。我们认为，从高度规则的对角线矩阵到高度不同行长度的无规则矩阵，数据交涉的实现技术都贯穿于其中。即使 SPMV 操作具有不规律性，我们可以证明他们可以被成功的映射为细粒度并行结构（fine-granted paraller archietecture）从而被 GPU 使用。对于非结构化矩阵，我们使用大约每秒 10 亿次的双精度和每秒 15 亿次的单精度来衡量他们的性能。此外，这些内核程序达到大约每秒 90Gb 带宽或者 65%的峰值时，表明在限制带宽的计算中存在现对较高的有效性。2. 稀疏矩阵格式这里有大量稀疏矩阵的代表格式，每一种都有不同的
6、存储要求，计算特征，以及访问和进行矩阵操作的方法。在稀疏矩阵矢量乘法（SPMV ）中，我们并不关心修改矩阵元素，我们仅考虑静态的稀疏矩阵格式，而不是那些适合元素快速插入和删除的特征。各种稀疏矩阵实例之间最主要的差别在于他们的稀疏模式，或者说非零值的结构。虽然稀疏格式适合高度具体的各种矩阵类具有很大的计算吸引力，同样重要的是考虑一般存储策略针对任意稀疏模式的矩阵存储。在这一章节的接下来部分我们总结了一些稀疏矩阵的范例并且讨论了他们之间的相关权衡。第四部分将讨论采用 CUDA 技术对每一种稀疏格式执行 SPMV 操作。2.1 对角线格式当非零值被限制在少量的矩阵对角线上时，对角线格式就是这样的一个
7、非常恰当的范例，尽管并不是一个通用的格式，对角线存储策略仍然可以有效地对矩阵（从实际应用的模型到规则的网格）进行编码，并且是一种普遍的离散化方法。图 4 说明了一个带有 5 条对角线的 25*25 的稀疏矩阵格式。图 4：一个带有 5 条对角线的 25*25 的稀疏矩阵格式。该矩阵近似于 5*5 的拉普拉斯算子。对角线格式由两个矩阵组成；data 矩阵，用来存储非零值， offset 数组，用来存储每一条对角线到主对角线的偏移量。规定，主对角线偏移量相当于 0，当 i0 时表示在主对角线之上的第 i 条对角线，当 i0 时表示在主对角线之下的第 i 条对角线。图 5 说明了一个带有 3 条对角
8、线的对角线格式矩阵的范例。在我们的执行过程中，data 矩阵是按列顺序存储的，因此它的连续存储元素的线性关系在相邻的对角线上。用*标记的元素是用来填充空位的，他们将被赋予一个随意值存储。图 5：data 矩阵和 offset 数组组成了矩阵 A 的对角线格式实例对角线格式的优势是双重的。首先，每一个非零值的行和列指针暗含地被他们所在对角线的位置和对应的偏移量所定义。这种暗含的索引值不仅减少了矩阵占用的内存，同时也减少了进行 SPMV 操作时的大量数据转移。第二，所有对数据的内存访问，x 和 y 变量是同时进行的，这将提高内存处理的效率。对角线格式的下降趋势也是非常明显的。DIA 在矩阵之外分配
9、变量的存储空间，并且明确地存储对角线上的零变量值，在预设的应用领域，例如适用于规则网格的模版，这都是不成问题的。然而，值得注意的是，许多矩阵的稀疏特性并不适合 DIA，就像图 6 中所描述的矩阵。图 6：不适合线性对角线格式的线性模型2.2 ELLPACK 格式另一种适合变量结构的存储策略是 ELLPACK (ELL)格式，对于一个 M*N 的矩阵，并且每行最多有 K 个非零值元素， ELLPACK 格式将非零值存储于一个 M*K 的 data 矩阵中，该矩阵中少于 K 个非零值得行被称为 zero-padded，同样的，相应的列指针被存储在indices 矩阵中，再次用 0 或者其他的 哨兵
10、值（sentinel values）来填补空缺。因为非零值不需要遵循任何特殊格式，所以 ELL 格式比 DIA 格式更普遍。图 7 说明了一个每行有 3 个非零值的矩阵应用 ELLPACK 格式的实例。与在 DIA 格式中一样，data 和 indices 矩阵都是按列顺序来存储的。图 7：data 和 indices 矩阵组成了 A 矩阵的 ELLPACK 格式实例当矩阵每行的非零值的最大数目不完全不同于平均值时，ELL 格式就是一种非常有效地范例。注意到当列指针被显式存储的同时行指针也被隐式定义了，就像 DIA 格式。除了被用于规则网格的模版操作定义的矩阵，由半机构化的网格和完全非结构化的
11、网格得到的矩阵均满足这一格式标准。图 8 说明了适合 ELL 存储策略的半结构化和非机构化网格。图 8：半结构化（左）和结构化（右）矩阵的例子，其顶点和边的连通性可以用 ELL格式有效地编码。在每种情况下，每种变量值不高于平均值。在实际中，非机构化网格并不总是满足这种要求。的确，像图 9 中说明的，每一行非零值得最大数目和平均数目之间的比率可能会任意大。明显的，单独的 ELL 格式并不适合从这些网格中得到的代表矩阵。然而，一种联合了 ELL 格式和其他的格式组成的混合存储策略就可以提供一种有用的替代策略。我们在 3.5 部分讨论这种混合格式。图 9：车轮子顶点的联通性不适合用 ELL 格式编码
12、，应为中心顶点于平均值有高度的关联性。那么 data 矩阵和 indicesj 矩阵中的绝大多数值将会被浪费掉。2.3 并列格式（COO）并列格式是一种特别简单的存储格式。数组：row,col 和 data 分别存储行指针，列指针和矩阵中的非零值。尽管 COO 格式是一种普遍的稀疏矩阵范例，对于一个任意的稀疏格式，所需要的存储总是与非零值得数量成比例的。不用于 ELL 和 DIA 格式，在 COO 格式中，行和列指针都被显示地存储。图 10 说明了一个矩阵的 COO 范例。在我们的执行过程中，row 矩阵被存储以确保具有相同行指针的变量将被一起存储。图 10：数组 row,col,dara 组
13、成了 A 矩阵的 COO 格式范例2.4 压缩稀疏行格式（CSR）压缩矩阵格式是一种流行的，具有一般用途的稀疏矩阵范例。和 COO 格式一样，CSR格式分别将列指针和非零值显式的存储在 indices 和 data 数组中。第三个数组 ptr（存储行指针）也采取 CSR 格式。对于一个 M*N 稀疏矩阵，ptr 数组长度为 M+1,并将第 i 行的偏移量存储在 ptri中。Ptr 数组的最后一个值将另外对应于第 M+1 行，存储 NNZ（矩阵中非零值的数目）。图 11 描述了一个矩阵的 CSR 格式实例。图 11：数组 ptr indices data 组成了 A 矩阵的 CSR 格式范例CS
14、R 格式可以看成是有一种 COO 存储实例（存在经常反复出现的行指针，因而采取一种简单的压缩策略）的扩展。因此 COO 格式和 CSR 格式之间的转化时直接且高效的。除了减少存储单元外，行指针还使矩阵值的查询更加方便快捷，并提供了其他一些好处，比如在一个特殊的行中（ptri+1 - ptri）非零值得数目可以随时计算得到。由于这些原因，压缩行格式经常被用于稀疏矩阵计算（例如稀疏矩阵乘法），也包括稀疏矩阵向量乘法(本文的重点) 。2.5 混合格式（HYB）虽然 ELL 格式非常适合矢量结构，但是当每个矩阵行的非零值数目变化特别大时，ELL 格式的高效性就很快消失了。 相比之下，COO 模式的存储
15、效率于每一行中非零值得分布式没有关系的。此外，就像我们将要在第 4 部分说明的，这种不变性会延伸到一个使用 COO 格式存储的稀疏矩阵向量产品的开销，尽管在最好的情况下采用 ELL 格式执行SPMV 要好于 COO 格式，但是将 ELL 和 COO 两种格式混合起来，用 ELL 格式存储矩阵中的大部分元素，而用 COO 格式存储剩余的部分，这样会产生更大的效用。考虑用一个矩阵来编码一个三角形网格的顶点连通，同时如图 9 所示，这样一个网格的最大顶点度会任意大，顶点度的平均值由拓扑量限制。因此在一个典型的平面三角形网格中，不考虑特殊异常的顶点的存在，绝大多数的顶点都有大约 6 个相邻的顶点。因为
16、相同的属性在更高层面上搁置非结构化网格，所以有必要深入研究稀疏矩阵存储策略来利用这一知识。混合格式的目的在于用 ELL 的 data 数组存储每行的非零值的典型数目，将剩余的特殊行的值用 COO 格式存储。经常每行的非零值的典型数目被称为先验，就像平面网格一样，而矩阵的 ELL 部分是很容易被提取出来的。然而。在一般情况下，这个数目必须由输入矩阵确定。我们的执行过程将计算出一个行尺寸的柱状图，并确定一个最大数目 K，这样就可以在 HYB 矩阵的 ELL 部分使用每行 K 列，满足一定得客观规律。基于试验证明结果，我们可以预测，完全化的 ELL 格式大体上要比 COO 格式快 3 倍。根据这一模
17、型假设，如果至少有三分之一得矩阵行包含 K 个或者更多的非零值，那么添加第 K 列到 ELL 结构中是有好处的。如果包含 K 个或者更多非零值的行的数目少于这个阈值，那么就将剩下的非零值存于 COO 中。2.6 其他格式在这一部分所讨论的格式仅代表了完整的稀疏矩阵范例空间的一小部分。在这些可供选择的格式中，许多是直接来自于其中一种现行的格式（比如修正的 CSR 格式）。其他的，像带有序列或者锯齿形对角线的 CSR 格式是基本格式的自然泛化。结构化的和非结构化的HYB 组合格式，例如 DIA 和 CSR 格式，也是很有用的。3稀疏矩阵向量乘法 SPMV3.1DIA kernel并行 SPMV 的
18、对角线格式是非常简单明了的：一个线程被分配给矩阵的一行。在图 15中，每个线程首先计算出一个顺序的线程索引，row，然后计算对应的矩阵行和 X 向量的点乘。图 15：SPMV 内核用于 DIA 稀疏矩阵格式回忆 3.1 节我们可知，DIA 格式的 data 数组时按列顺序存储的。正如图 16 中所示，DIA 格式中 data 数组的线性化可以确保在同一个 warp 中的线程能连续地访问内存。因为这个矩阵有三条对角线（num_diags=3 ）,data 数组会被每一个执行线程（共四个）访问三次。因为行是连续的，所以线程相当于连续的矩阵列，在同一个 warp 中的线程也可以连续的访问 X 向量。
19、正如上文所说，内存访问 data 数组和 x 向量式连续的，但是却不能与任意具体的多字边界达到普遍地一致。幸运的是，对 data 数组的一致访问已经通过数组中的填充列（图 5）轻松地解决了。我们在执行过程中使用修改的标准 DIA 数据结构。再者，DIA 内核在每个 block 中仅访问一次 offsets 矩阵中的变量，不同于图 15 中每个线程都要对上述变量进行访问。图 16：DIA data 矩阵的线性化和 DIA 格式下 SPMV 内核的内存访问模式3.2ELL 内核图 17 所示，ELL 内核，跟 DIA 内核一样，每个矩阵行使用一个线程来并行计算。ELL 内核的结构差不多与 DIA
20、内核结构是一致的，不同的是，列指针在 ELL 中是明确的，而在 DIA 中是暗含的。因此，对于一 DIA 格式存储的矩阵来时 DIA 内核的执行一般要比ELL 内核性能好。图 17：SPMV 内核用于 ELL 稀疏矩阵格式ELL 内存访问模式，图 18 中所示，类似于 DIA 方式。因为 ELL 格式的 data 矩阵和indices 矩阵是用相同的方式存储的，因此我们只表示一种访问模式。和 DIA 格式一样，我们的执行过程也使用填充值来适当的统一 ELL 数组。然而，与 DIA 不同的是，ELL 内核不一定连续的访问 X 向量。图 18：ELL data 矩阵的线性化和 ELL 格式下 SP
21、MV 内核的内存访问模式3.3 CSR 内核图 19 说明了用于 CSR 格式的串行 CPU 内核。因为矩阵的一行与 X 向量的点乘计算可能跟其他的行不能独立进行，而 CSR 模式的矩阵乘法却可以很容易的实现并行化。图20 说明了一个简单的 CUDA 执行过程，仍然在每个矩阵行使用一个线程。这些方法的一些变种已经被写入文档13，我们用来作为标量内核的参考。图 19：用于 CSR 格式的串行 CPU 内核图 20：用于 CSR 稀疏矩阵格式的 SPMV 内核，每个矩阵行对应一个线程当标量内核显示细粒度并行性时，它在性能上就会出现一些不足。在这样问题中最为重要的是，位于同一个 warp 中的线程访问 CSR indices 和 data 矩阵的方法。当列向量和一个给定行的非零值被连续的存储在 CSR 的 data 结构中时，这些值将不能被同时访问。相反，每一个线程顺序的读取其对于行的元素时，产生的模型如图 21.。
