稀疏矩阵计算：矩阵划分和负载均衡

第二个挑战性问题是矩阵划分和负载均衡。随着多核/众核体系结构的发展，矩阵的合理划分和任务在多个核或节点之间的均衡分配对于充分利用硬件资源至关重要。

（1）SpGEMM的形式化表示：SpGEMM公式化和矩阵划分密切相关，总共有四种SpGEMM的计算方式。

①基于外积的形式化表示。该表示分别基于输入矩阵A和B的按列和按行划分，通过将A的每一列a*i与B的对应行bi*的外积相加得出结果矩阵C，即

运算符⊗表示向量外积。

②基于内积的形式化表示。该表示分别基于输入矩阵A和B的按行和按列划分，结果矩阵C由矩阵A的每一行和矩阵B的每一列的内积组成，即

i=1，2，…，p，j=1，2，…，r ，I（i，j）表示满足aik和bkj都不为零的索引k的集合。

③按行相乘的形式化表示。此表示基于两个输入矩阵的按行划分。结果矩阵C的每一行ci*是通过将ai*中的每个非零项aik和B中相应行bk*的相乘结果相加而得，即

Ii（A）表示A的第i行中非零项的列索引k的集合。

④按列相乘的形式化表示。此表示基于两个输入矩阵的按列划分，类似于按行相乘的公式化。结果矩阵C的每一列c*j是通过将b*j的每个非零项bkj和A中相应列a*k的相乘结果相加而得，即

Ij （B）表示B的第j列中非零项的行索引k的集合。(www.xing528.com)

（2）块划分：SpGEMM的四种形式化表示代表了四种不同的一维矩阵划分方案：列−行，行−列，行−行，列−列。但是考虑到稀疏矩阵中非零项的不规则分布，基于普通切片的矩阵划分通常会导致负载不均衡。因此一些研究人员致力于在特定体系结构上达成负载均衡。

刘伟峰教授和Brian Vinter提出了一种基于输出空间分解的负载均衡方法。与Dalton等人提出的方法有类似的思想，即对中间结果矩阵的行进行排序，将其分为3个子矩阵，其中包括大小范围不同的行，并对子矩阵使用差分ESC方法。而他们的方法不是严格地对中间矩阵的行进行排序，只是通过CPU上更快的线性时间遍历将行分配给固定数量的任务组，然后以更详细的方式分解了输出空间，以确保更加有效的负载均衡。为了平衡负载并有效利用GPU资源，Nagasaka等人在符号阶段将矩阵行按中间乘积的数量（即输出矩阵的矩阵行中非零元数的上限）分为几组。在数值阶段，按每行中非零元的数量进行分组。Winter等人还设计了一种适用于GPU的均衡的任务分配方案。给每个线程块分配矩阵相等数量的非零项，而忽略矩阵的行边界。Deveci等人提出了两种适用于KNL（Knight Landing）和GPU的分块方法。区别在于，对于KNL，按列对输入矩阵A进行划分，而对于GPU，因为其具有多级内存的架构，故对A采用二维划分。Deveci等人也对两个输入矩阵采用了按行的划分方案，但在高性能计算平台上使用了不同的任务分配方案。Chen等人设计了一个三级划分方案，以改善SpGEMM在神威·太湖之光超级计算机上的负载均衡。第一级表示将A划分为NP个子矩阵，其中NP是计算中使用的核组数。在第二级中，B被划分为NC个子矩阵，其中NC是每个核组中使用的计算处理单元（CPE）核心的数量。在第三级中，A和B的子矩阵进一步被划分为几组。

为了提高有效的并发性，Kurt等人实现了“虚拟warp”，该方法将线程块划分为多个“虚拟warp”，并且只有虚拟warp中的线程参与处理矩阵的同一行，而不是用线程块中的所有线程处理矩阵的一行。例如，对于大小为128的线程块和大小为4的虚拟warp，每个线程块将同时负责A的四行。然后，使用32（ = 128/4）个线程来循环处理A中每一行对应的B中的非零项。另外，他们首次对SpGEMM中数据移动需求的上限进行了系统的探索和分析。使用的任务分配方案类似于基于GPU子warp的行合并方法。Ballard等人对AMG中Galerkin三元乘积中的SpMM运算的通信成本提出了详细的理论预测。他们还推断出，对于第一个矩阵乘，按行相乘是最好的一维方法，而对于第二个矩阵乘，外积法是最好的一维方法。

Van de Geijn等人提出了基于分布式内存平台的二维块数据分解和映射方法，用于稠密矩阵相乘。Buluc等人将其用在了SpGEMM的计算中。针对二维块的数据分解后引入的超稀疏矩阵，他们利用各个子矩阵的超稀疏性提出了一种新的压缩格式DCSC（doubly compressed sparse column）。然后分别基于外积和按列相乘开发了两种顺序执行的SpGEMM算法，这些算法主要提升了二维划分后的超稀疏子矩阵的相乘效率。最后，他们提出了一种高效并行的可扩展通用矩阵乘算法（SUMMA）。此外，Azad等人基于Buluc等人的工作实现了三维的并行SpGEMM算法，它同时利用了第三个处理器网格维度内的节点间并行和节点内的多线程并行。

（3）超图划分：传统的基于块的矩阵划分方案很少从输入或输出矩阵的稀疏模式来考虑通信问题。Akbudak等人将一种基础的超图划分（HP）模型应用到基于外积的并行SpGEMM中。其目的在于使通信成本最小化，同时以预定义的最大允许的不平衡系数来保持平衡的计算负载。HP模型由三个部分组成：第一个部分是表示A的每一列与B的对应行的外积的顶点；第二个部分是表示C的每一个非零项的顶点，以支持输出矩阵非零元的二维划分；第三个部分是表示C的每一个非零项的超边（或网络），用于编码在中间结果累加中发生的总通信量。假设在HP模型中，输入矩阵A和B的划分为

式中，Q是由划分引起的置换矩阵，基于外积并行的SpGEMM的HP模型包括两个阶段。在第一阶段（也称为乘法阶段），每个处理器Pk拥有被置换矩阵的列和行，然后在本地完成无通信开销的SpGEMM计算：。第二阶段（也称为求和阶段）规约第一阶段产生的部分结果，从而计算出C的每一个非零项的最终值。为了在两个阶段上保持负载均衡，他们根据每一个顶点的权重提出相应的约束条件。其中，每一个顶点分配有两个权重，分别代表与该顶点相关联的两个阶段的计算开销。此外为了实现负载均衡，作者在划分中引入了不平衡系数。

Ballard等人基于HP模型，提出了一种细粒度的框架，以证明并行和顺序SpGEMM的通信下限依赖于矩阵的稀疏性。进一步，他们证明了对于给定输入矩阵，确定其最优通信算法等价于相应HP问题的求解。Akbudak等人还使用HP模型来改进基于Xeon Phi架构的外积和内积式并行稀疏矩阵乘法。在基于内积的并行SpGEMM中，HP模型尽可能将贡献于C矩阵相同项的A矩阵列和B矩阵行分派给相同组。在基于外积的并行SpMM中，HP模型将需要相同B矩阵行的A矩阵行尽可能分配到相同组，从而更好地利用时间局部性。

此外，Akbudak等人也将图模型和HP模型应用于外积、内积和逐行乘积的并行SpGEMM算法中。Demirci等人采用类似的二分图模型用于分布式内存平台上的矩阵划分，但补充了一种3−约束的边权重分配方案，以保持服务器之间的负载均衡。Selvitopi等人使用图模型和超图模型实现map和reduce任务的同步分配，以提升数据局部性和维护不同处理器在map和reduce任务上的负载均衡。然后，考虑到MapReduce中map和reduce任务与SpGEMM中间结果累加的相似性，他们应用基于图和超图模型获得的MapReduce划分方案来改善SpGEMM中的数据局部性和负载均衡。