面向GPU和DCU架构的SpGEMM算法策略选择及优化方法

面向gpu和dcu架构的spgemm算法策略选择及优化方法
技术领域
1.本发明涉及高性能计算技术领域，具体涉及面向gpu和dcu架构的spgemm算法策略选择及优化方法。


背景技术：

2.spgemm(general sparse matrix-matrix multiplication)是形如c＝αa*b+βd的稀疏矩阵-矩阵乘法运算(其中，a、b、d是稀疏矩阵，结果矩阵c也是稀疏矩阵，α、β是两个标量)，它在代数多重网格(algebraic multigrid，amg)求解器、多源广度优先搜索、最短路径问题、填色问题、三角形计数问题和子图问题等都有着广泛的应用。
3.现有的spgemm算法大多数分为符号阶段和计算阶段两个阶段。符号阶段主要用于计算结果矩阵的行偏移，为结果矩阵的列索引和值分配内存空间。计算阶段则是根据分配的内存空间，对结果矩阵进行计算。常见的计算方法如esc法、hash法、breakpoint法，三种算法的思想如下：
4.1)esc法：即e(expansion,扩展)、s(sorting,排序)、c(compression,压缩)。由于esc方法先扩展再压缩的特性，使其需要额外的空间，因此该方法通常与上界法一起使用。扩展操作得到中间结果矩阵，然而中间结果矩阵是无序的，为了方便后续的压缩过程计算，需要再对中间结果矩阵以列索引进行排序，最后的压缩则是对已经排序过的中间结果矩阵进行压缩得到最终的结果矩阵。
5.2)hash法：将矩阵a的一行和矩阵b的若干行分配给一个vector(线程组，组内线程数小于64)、wavefront或block。计算时分别遍历矩阵a一行每一个非零元素以及对应在矩阵b中的行，之后通过hash表对相同列的值进行合并，hash冲突的处理主要是采用开放地址法和链地址法。
6.3)breakpoint法：一个block处理矩阵a的一行与对应矩阵b的若干行相乘，一个vector处理矩阵a的一个非零元素与矩阵b对应的一行相乘。同时对矩阵b列分块计算，每块大小为chunk_size，将每个块计算完的中间结果存入lds中，利用lds的原子操作，合并结果矩阵属于同一列的结果。
7.上面提及的spgemm算法分别适用于一些特定条件的矩阵的计算，对于不同的矩阵特征，则可能出现计算效率严重低下，甚至无法求解的情况，算法的通用性较差。因此在选择计算策略时需要综合考虑矩阵的特征，选用最适合的计算方法，从而最大化的利用资源。鉴于此，本发明提出了一种面向gpu和dcu架构的spgemm算法策略选择及优化方法。


技术实现要素：

8.本发明的目的在于解决普通spgemm算法适用范围较小，通用性较差的问题而提出了面向gpu和dcu架构的spgemm算法策略选择及优化方法。
9.为了实现上述目的，本发明采用了如下技术方案：
10.面向gpu和dcu架构的spgemm算法策略选择及优化方法，结合稀疏矩阵特征和gpu、
dcu架构的特点，对spgemm算法进行计算和访存优化，所述方法具体包括以下内容：
11.s1、符号阶段：
12.s1.1、矩阵b的稀疏度小于0.005：
13.采用上界法进行求解，计算结果矩阵中非零元素上限；
14.s1.2、矩阵b的稀疏度大于等于0.005：
15.首先对矩阵b进行压缩，然后根据矩阵规模选择不同策略：
16.若矩阵b规模小于等于65536，则使用位图法；否则采用hash法；
17.将breakpoint方法作为超过hash方法限制的补充方法；
18.s2、计算阶段：
19.s2.1、结果矩阵行分组：
20.对于计算阶段的esc方法和hash方法，当结果矩阵行最大非零元素数小于等于2048时采取行分组策略，根据结果矩阵的每一行非零元素的数量，将结果矩阵按行划分为不同的组，对不同的组设置不同的核函数配置，从而调度更多的block块，使计算与访存重叠；
21.当结果矩阵行最大非零元素数大于2048时，则选择breakpoint方法计算；
22.s2.2、计算结果矩阵：
23.s2.2.1、当矩阵行最大非零元素数大于2048时，计算阶段选择breakpoint方法，直接对矩阵进行计算；
24.s2.2.2、当矩阵行最大非零元素数小于等于2048时，策略划分如下：
25.①
矩阵b的稀疏度小于0.005：
26.执行结果矩阵行分组策略，并采用esc方法进行求解；
27.②
矩阵b的稀疏度大于等于0.005：
28.若结果矩阵c的稀疏度小于等于0.05，则执行结果矩阵行分组策略，并采用hash方法求解；否则采用breakpoint方法，对矩阵b进行分块计算。
29.优选地，所述方法在符号阶段所使用的策略选择算法包括有上界法、hash法、位图法以及breakpoint法；在计算阶段使用的策略选择算法包括有esc法、hash法和breakpoint法。
30.与现有技术相比，本发明提供了面向gpu和dcu架构的spgemm算法策略选择及优化方法，具备以下有益效果：
31.(1)本发明使用了大量的矩阵特征，对于不同的稀疏矩阵可以进行再归类；
32.(2)本发明能够实现自适应的策略选择，在保证计算结果正确性及算法求解效率的前提下，对原求解算法进行优化，根据不同的矩阵特征选择不同的spgemm求解算法，解决了单一spgemm算法无法适应不同稀疏矩阵特征的高效求解问题；
33.(3)本发明实现了自适应求解参数的设置，根据矩阵元素的特征，对矩阵行进行分组，每组采用不同的计算策略(vector-row、block-row)，最大化利用block资源，从而提高spgemm算法的求解效率。
附图说明
34.图1为本发明提出的面向gpu和dcu架构的spgemm算法策略选择及优化方法的流程
图；
35.图2为本发明提出的面向gpu和dcu架构的spgemm算法策略选择及优化方法的结果矩阵行分组流程图。
具体实施方式
36.下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。
37.spgemm算法的符号阶段、计算阶段各个算法的适用场景都不相同，为了最大限度发挥各种spgemm算法的效率以达到更好的应用效果，本文提出了一种面向gpu和dcu架构的spgemm算法策略选择及优化方法，具体实例如下。
38.实施例1：
39.本发明提出的面向gpu和dcu架构的spgemm算法策略选择及优化方法在符号阶段包含上界法、hash法、位图法以及breakpoint法；在计算阶段包含esc法、hash法和breakpoint法。根据矩阵行数、行非零元素数以及矩阵的稀疏度等特征进行矩阵行分组并选取最适合的计算策略。
40.假设稀疏矩阵a、b的规模均为n
×
n，存储格式均为csr格式。csr格式包含三个数组，其中row_ptr[]存储矩阵每一行非零元素的行偏移量，col_idx_ptr[]存储矩阵中每个非零元素的列索引，values_ptr[]存储每个非零元素值。请参阅图1，具体包括如下内容：
[0041]
1、符号阶段：
[0042]
首先根据矩阵b的行数、列数和非零元素数计算出矩阵b的稀疏度值：
[0043][0044]
1)若sparseness-b＜0.005
[0045]
采用上界法，根据矩阵a和矩阵b的row_ptr和col_idx_ptr数组，计算出矩阵a每个非零元素对应的矩阵b非零元素的数量，并对所有行求前缀和，得到总内存大小。
[0046]
2)若sparseness-b≥0.005
[0047]
首先对矩阵b进行压缩，然后根据矩阵规模选择不同策略：
[0048]
a、矩阵b的规模小于等于65536
[0049]
由于较为稠密的矩阵精确计算非零元素个数需要频繁访问矩阵b的列索引，因此使用位图法对矩阵b的列索引进一步压缩以提高访存效率，用bit位标记某个元素是否存在。假设位图长度为bitslen，则压缩后的存储空间大小为：matrix.rows/bitslen。将列索引按照位图思想进行压缩后，得到ids数组和bits_cols数组。ids数组保存压缩后的id号，bits_cols保存具体的bit位。根据压缩矩阵每行非零元素数的大小，分别采取vector和block的策略，调整不同的核函数参数，在lds中计算出中间结果的位图，然后再对非零元素通过规约等操作，得到总内存的大小。
[0050]
b、矩阵b的规模大于65536
[0051]
当矩阵规模变大时，由于lds空间最大为64kb，位图法只能调度少量的甚至单一的block进行运算，计算与访存不能重叠，导致计算效率大幅下降。因此当矩阵规模大于65536时，采用hash法来压缩b矩阵，精确计算内存空间。符号阶段产生冲突时采用的是开发地址
法解决，即继续向下寻求空间hash_val＝(hash_val+1&hashmask)(hashmask为当前分配的lds大小)。不断执行循环，直到完成插入操作。循环结束后，通过对非零元素进行规约操作统计出总内存大小。
[0052]
2、计算阶段：
[0053]
2.1、结果矩阵行分组：
[0054]
请参阅图2，通过符号阶段得到的结果矩阵的行偏移数组row_ptr_c[]，对不同非零元素数量的行进行分组计算的策略。
[0055]
当矩阵行最大非零元素数量小于等于2048的时候，此时计算阶段的方法为esc方法或者hash方法。设置group_num＝{16,32,128,256,512,1024,2048}，根据结果矩阵的每一行非零元素的数量，选择合适的group计算。例如第i行非零元素数nnz_c＝row_ptr_c[i+1]
–
row_ptr_c[i]《group_num[j](0《j《7)，则对应的esc方法或hash方法采用vector-row或者block-row的策略，每个vector数组计算一行或者每个block计算一行，同时设置对应于group_num的核函数配置。结果矩阵行分组策略主要分为三个步骤：
[0056]
(1)在划分完每块的起始位置block_row_begin和结束位置block_row_end后，每个block块分别在lds上计算出上述范围内属于每个group_num的非零元素数量，然后通过规约、求和等操作将每个block块的数据写回至全局内存中。
[0057]
(2)划分lds内存，将相同大小的group_num的不同block进行求和，求出矩阵每一部分group_num的大小，并写回至全局内存。
[0058]
(3)再次遍历所有的非零矩阵行，将处于不同group_num范围的行号写入到lds的block_group_num数组中，最终将结果写回至全局内存group_row_ptr数组中，以便计算阶段对不同的行采取不同的策略选择与对应的核函数配置。
[0059]
2.2、计算结果矩阵：
[0060]
当矩阵行最大非零元素数量大于2048的时候，则计算阶段选择breakpoint方法，直接对矩阵进行计算；否则采用如下的策略划分：
[0061]
1)若sparseness_b＜0.005
[0062]
先执行结果矩阵行分组策略，求出每一行最适合的策略以及核函数配置，然后应用传统的esc方法进行求解，计算出最终的结果。
[0063]
2)若sparseness_b≥0.005
[0064]
a、矩阵c的稀疏度小于等于0.05
[0065]
先执行结果矩阵行分组策略，再采用hash法进行计算。与符号阶段不同的是，计算阶段采用链地址法处理冲突，将当前被占用的空间移动到next链中，同时增加指向next相应位置的标记。
[0066]
b、矩阵c的稀疏度大于0.05
[0067]
将矩阵b按列进行分块，每块大小为chunk_size列，每次处理一个块。将划分后的矩阵保存至分块数组bp中，bp保存划分后每一行每个chunk_size的起始位置和结束位置，每一行分配的大小为：matrix.cols/chunk_size+1。然后选用breakpoint方法计算结果矩阵的值，并写回至全局内存。
[0068]
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其
发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。

技术特征：
1.面向gpu和dcu架构的spgemm算法策略选择及优化方法，其特征在于，结合稀疏矩阵特征和gpu、dcu架构的特点，对spgemm算法进行计算和访存优化，所述方法具体包括以下内容：s1、符号阶段：s1.1、矩阵b的稀疏度小于0.005：采用上界法进行求解，计算结果矩阵中非零元素上限；s1.2、矩阵b的稀疏度大于等于0.005：首先对矩阵b进行压缩，然后根据矩阵规模选择不同策略：若矩阵b规模小于等于65536，则使用位图法；否则采用hash法；将breakpoint方法作为超过hash方法限制的补充方法；s2、计算阶段：s2.1、结果矩阵行分组：对于计算阶段的esc方法和hash方法，当结果矩阵行最大非零元素数小于等于2048时采取行分组策略，根据结果矩阵的每一行非零元素的数量，将结果矩阵按行划分为不同的组，对不同的组设置不同的核函数配置，从而调度更多的block块，使计算与访存重叠；当结果矩阵行最大非零元素数大于2048时，则选择breakpoint方法计算；s2.2、计算结果矩阵：s2.2.1、当矩阵行最大非零元素数大于2048时，计算阶段选择breakpoint方法，直接对矩阵进行计算；s2.2.2、当矩阵行最大非零元素数小于等于2048时，策略划分如下：
①
矩阵b的稀疏度小于0.005：执行结果矩阵行分组策略，并采用esc方法进行求解；
②
矩阵b的稀疏度大于等于0.005：若结果矩阵c的稀疏度小于等于0.05，则执行结果矩阵行分组策略，并采用hash方法求解；否则采用breakpoint方法，对矩阵b进行分块计算。2.根据权利要求1所述的面向gpu和dcu架构的spgemm算法策略选择及优化方法，其特征在于，所述方法在符号阶段所使用的策略选择算法包括有上界法、hash法、位图法以及breakpoint法；在计算阶段使用的策略选择算法包权利要求书括有esc法、hash法和breakpoint法。

技术总结
本发明公开了面向GPU和DCU架构的SpGEMM算法策略选择及优化方法，属于高性能计算技术领域；本发明在大量试验的基础上，对SpGEMM算法进行计算和访存优化；结合稀疏矩阵特征和GPU、DCU架构的特点，实现了SpGEMM自适应策略选择算法，并自适应的为不同算法设置相应求解参数，以达到更高的SpGEMM求解效率。以达到更高的SpGEMM求解效率。以达到更高的SpGEMM求解效率。


技术研发人员：胡长军 蒋子涵 卢旭 何远杰 储根深
受保护的技术使用者：北京科技大学
技术研发日：2023.06.02
技术公布日：2023/8/31