一种矩阵数据架构及其基于压缩稀疏列的加速SPMV的方法

group)的简称。
7.该数据架构通过对稀疏矩阵的格式进行处理获得，该数据架构具有占用空间少的特点，而使用两种不同的矩阵块结构可以实现存储量不同程度的压缩。使用该矩阵数据架构进行加速spmv计算时，能够快速确定矩阵块结构参数。
8.优选的，所述定长稠密vxg块和所述定长稀疏vxg块均包括矩阵块头和vxg数组；所述矩阵块头包括子矩阵的行和列范围描述、vxg的个数、每个vxg包含的列数s
vxg
、每一列包含所述矩阵块元素的个数se；所述vxg数组包括所述矩阵块元素数组、每个vxg对应x分量的下标数组和对应y向量的范围；所述矩阵块中第i个vxg对应的非零元素稠密子矩阵存放在所述矩阵块元素数组的第s
vxg
*(i-1)*se至到第s
vxg
*(i)*s
e-1个元素中，并且按照行块优先的方式存放矩阵块元素。
9.优选的，同一列上的非零元素存储在多个矩阵块元素中，稠密的矩阵块元素是长度固定为s
vvec
的向量，存储在同一个稠密矩阵块元素cscve中的矩阵元素具有相同的列下标和连续行下标。
10.优选的，所述定长稀疏vxg块的vxg数组还包括每个vxg中矩阵块元素的非零元素指示位图向量。在使用定长稀疏vxg块时，需要进行额外的展开运算，填充呈稠密的矩阵块元素。
11.优选的，同一列上的非零元素存储在多个矩阵块元素中，稀疏的矩阵块元素是一个行下标范围长度不大于s
vvec
的向量，同一个稀疏的矩阵块元素中的矩阵元素具有相同的列下标。
12.优选的，所述vxg连接共享访问相同矩阵行下标的多列矩阵块元素，并确保在相邻的矩阵块元素中要访问的矩阵行下标是相同的或者是相邻的；所述vxg对应一个稠密的子矩阵且vxg上不同列的列下标并不一定是连续的。
13.一种基于压缩稀疏列的加速spmv的方法，包括如下步骤：
14.步骤一：将输入的稀疏矩阵a矩阵格式转换为上述的矩阵数据架构形式的矩阵格式，该矩阵数据架构为稀疏矩阵存储格式；
15.步骤二：对局部临时y向量重排序和完全向量化的spmv计算。
16.通过局部重排序策略和上述的矩阵数据架构，该加速spmv的方法能够高效利用硬件向量化计算指令和内存带宽，使得spmv能够以更快的的速度运行，达到加速的目的。
17.优选的，所述步骤二的具体流程为：
18.s1:输入步骤一中的格式转换后的稀疏矩阵a(大小为m*n)、长度为n的向量x、长度为m的向量y、向量化长度为s
vvec
；输入局部行向量重排映射集合{ιk}和逆映射集合{ι
k-1
}；
19.s2：进行局部临时y向量重排序和完全向量化的spmv计算，使得z＝局部临时向量重排序(ιk,y)和～z＝完全向量化的spmv(ak,x,z,s
vvec
)；
20.s3：
21.s4：以s3中y～作为s2中的y向量，对矩阵a的每一个子矩阵块ak循环执行步骤s2和s3，直至所有的子矩阵块ak完成完全向量化的spmv计算。
22.该方法在不使用汇编代码的情况下能优化spmv代码并获得低级别simd指令的优势的原因是完全向量化的spmv算法很容易被编译器自动优化并生成所需的simd指令代码。
23.优选的，在所述s2中，具体的流程为：
24.s2.1:对矩阵块ak中以vxg存储的稠密矩阵e进行向量化运算q＝q+e*x，q的初始值是一个长度为m的0向量，过程包括：
25.s2.1.1对e中的第i个矩阵块元素mi进行向量化运算得到w＝a*mi+qi，其中，w、mi和qi均为长度为s
vvec
的向量，a为一个标量是e对应x分量的下标数组里第i mod s
vxg
个元素索引的x分量数值，qi是me对应的q向量元素数组，其下标范围是从(e对应y向量的起始值*(i-1)*s
vvec
/s
vxg
)到(e对应y向量的起始值*i*s
vvec
/s
vxg-1)，并且qi＝w写入到q向量的上述下标范围中。
26.s2.1.2:循环执行步骤s2.1.1，令上一个循环中步骤s2.1.1的q作为下一个循环中的步骤s2.1.1的q，下一个循环中的q＝上一个循环中的a*mi+qi；
27.s2.1.3：完成e中所有的矩阵块元素mi的向量化运算，得到～z＝q＝q+e*x
28.s2.2:循环执行步骤s2.1、s2.1.1、s2.1.2和s2.1.3，令上一个循环中步骤s2.1.3的q作为下一个循环中的步骤s2.1的q，下一个循环中的q＝上一个循环中的q＝q+e*x；
29.s2.3：完成ak中所有以vxg存储的稠密矩阵e的向量化运算，得到～z＝完全向量化的spmv(ak,x,z,s
vvec
)＝z+ak*x＝z+q。
30.每一次浮点数运算都直接读取一个给定长度为s
vvec
的矩阵列元素向量然后数乘一个x向量元素并且对齐加和并存储到y向量上的一个下标连续且长度为s
vvec
的子向量中。因此spmv的计算的过程是一个完全向量化的过程，不存在额外的存取非连续下标的向量元素操作。而且s
vvec
是可变参数，能使计算过程不改变的情况下适合不同宽度和硬件架构的向量化指令。
31.优选的，局部行向量重排映射集合将y向量按照投影数据区域上的i、j坐标轴排序布局变成按照平行和垂直于某一个像素点在投影数据区域的轨迹排序布局；矩阵每一列的非零元素根据y向量重排映射进行分块排列和存储。
32.与现有技术相比，本发明的有益效果是：该数据架构具有占用空间少的特点，在程序运算的过程中能够加快运行速度。再结合基于压缩稀疏列的加速spmv的方法，能够高效利用硬件向量化计算指令和内存带宽，使得spmv能够以更快的的速度运行，达到加速的目的。
附图说明
33.图1是本发明的一种基于压缩稀疏列的加速spmv的方法的流程图；
34.图2是ct影像重建过程中不同像素点及其在投影区域的轨迹分布；
35.图3是图1不同像素对应的矩阵列的cscve内存布局；
36.图4是向量y和矩阵非零元素的不同内存布局下的simd效率；
37.图5是单线程基于压缩稀疏列的向量化稀疏矩阵向量乘算法的计算全过程；
38.图6是10个软件使用平台上全部的cpu核心进行spmv时,内存使用量、每秒浮点运算量和只读访存带宽利用率的统计。
具体实施方式
39.附图仅用于示例性说明，不能理解为对本专利的限制；为了更好说明本实施例，附图某些部件会有省略、放大或缩小，并不代表实际产品的尺寸；对于本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。附图中描述位置关系仅用于示例性
说明，不能理解为对本专利的限制。
40.下面通过具体实施例，并结合附图，对本发明的技术方案作进一步的具体描述：
41.实施例1
42.如图1所示为一种矩阵数据架构的实施例，包括矩阵头和矩阵块；所述矩阵头包括了矩阵行、列数，非零元个数、矩阵块的个数和向量化长度s
vvec
；每一个所述矩阵块对应的x向量像素几何方块的边长为s
imgb
；所述矩阵块的结构至少有两种，分别为定长稠密vxg块和定长稀疏vxg块。vxg是向量化执行组(vectorized execution group)的简称，是矩阵数据架构中参与执行spmv的基本单位，其行下标是连续的，因此用以索引数据的存储和访问可以被大幅度压缩。
43.其中，所述定长稠密vxg块和所述定长稀疏vxg块均包括矩阵块头和vxg数组；所述矩阵块头包括子矩阵的行和列范围描述、vxg的个数、每个vxg包含的列数s
vxg
、每一列包含所述矩阵块元素的个数se；所述vxg数组包括矩阵块元素数组、每个vxg对应x分量的下标数组和对应y向量的范围；所述矩阵块中第i个vxg对应的非零元素稠密子矩阵存放在矩阵块元素数组的第s
vxg
*(i-1)*se至到第s
vxg
*(i)*s
e-1个元素中，并且按照行块优先的方式存放矩阵块元素。
44.具体的，同一列上的非零元素存储在多个矩阵块元素中，稠密的矩阵块元素是长度固定为s
vvec
的向量，存储在同一个稠密矩阵块元素中的矩阵元素具有相同的列下标和连续行下标。所述定长稀疏vxg块的vxg数组还包括每个vxg中矩阵块元素的非零元素指示位图向量。在使用定长稀疏vxg块时，需要进行额外的展开运算，填充呈稠密的矩阵块元素。同一列上的非零元素存储在多个矩阵块元素中，稀疏的矩阵块元素是一个行下标范围长度不大于s
vvec
的向量，同一个稀疏的矩阵块元素中的矩阵元素具有相同的列下标。
45.在本实施例中，所述vxg连接共享访问相同矩阵行下标的多列矩阵块元素，并确保在相邻的矩阵块元素中要访问的矩阵行下标是相同的或者是相邻的；所述vxg对应一个稠密的子矩阵且vxg上不同列的列下标并不一定是连续的。
46.具有本实施的矩阵数据架构形式的矩阵格式称为cscv格式矩阵，矩阵块元素称为cscve。基于压缩稀疏列向量称为cscv。
47.本实施例的工作原理：该数据架构通过对稀疏矩阵的格式进行处理获得，该数据架构具有占用空间少的特点，而使用两种不同的矩阵块结构可以实现存储量不同程度的压缩。使用该矩阵数据架构进行加速spmv计算时，能够快速确定矩阵块结构参数。
48.本实施例的有益效果:该数据架构具有占用空间少的特点，在程序运算的过程中能够加快运行速度。
49.实施例2
50.一种基于压缩稀疏列的加速spmv的方法，包括如下步骤：
51.步骤一：将输入的稀疏矩阵a的矩阵格式转换为实施例1的矩阵数据架构形式的矩阵格式；
52.步骤二：对局部临时y向量重排序和完全向量化的spmv计算。
53.通过局部重排序策略和上述的矩阵数据架构，该加速spmv的方法能够高效利用硬件向量化计算指令和内存带宽，使得spmv能够以更快的的速度运行，达到加速的目的。
54.具体的，所述步骤二的具体流程为：
55.s1:输入步骤一中格式转换完成后的稀疏矩阵a(大小为m*n)、长度为n的向量x、长度为m的向量y，向量化长度为s
vvec
；输入局部行向量重排映射集合{ιk}和逆映射集合{ι
k-1
}；
56.s2：对矩阵a的子矩阵块ak进行局部临时y向量重排序和完全向量化的spmv计算，使得z＝局部临时向量重排序(ιk,y)和＝完全向量化的spmv(ak,x,z,s
vvec
)；＝完全向量化的spmv(ak,x,z,s
vvec
)的具体的流程为：
57.s2.1:对矩阵块ak中以vxg存储的稠密矩阵e进行向量化运算q＝q+e*x，q的初始值是一个长度为m的0向量，过程包括
58.s2.1.1对e中的第i个矩阵块元素mi进行向量化运算得到w＝a*mi+qi，其中，w、mi和qi均为长度为s
vvec
的向量，a为一个标量是e对应x分量的下标数组里第i mod s
vxg
个元素索引的x分量数值，qi是me对应的q向量元素数组，其下标范围是从(e对应y向量的起始值*(i-1)*s
vvec
/s
vxg
)到(e对应y向量的起始值*i*s
vvec
/s
vxg-1)，并且qi＝w写入到q向量的上述下标范围中。
59.s2.1.2:循环执行步骤s2.1.1，令上一个循环中步骤s2.1.1的q作为下一个循环中的步骤s2.1.1的q，下一个循环中的q＝上一个循环中的a*mi+qi；
60.s2.1.3：完成e中所有的矩阵块元素mi的向量化运算，得到～z＝q＝q+e*x
61.s2.2:循环执行步骤s2.1、s2.1.1、s2.1.2和s2.1.3，令上一个循环中步骤s2.1.3的q作为下一个循环中的步骤s2.1的q，下一个循环中的q＝上一个循环中的q＝q+e*x；
62.s2.3：完成ak中所有以vxg存储的稠密矩阵e的向量化运算，得到＝完全向量化的spmv(ak,x,z,s
vvec
)＝z+ak*x＝z＝q。
63.s3：令矩阵a的向量＝局部临时向量重排序
64.s4：以s3中作为s2中的y向量，对矩阵a的每一个子矩阵块ak循环执行步骤s2和s3，直至所有的子矩阵块ak完成完全向量化的spmv计算。
65.该方法每一次浮点数运算都直接读取一个给定长度为s
vvec
的矩阵列元素向量然后数乘一个x向量元素并且对齐加和并存储到y向量上的一个下标连续且长度为s
vvec
的子向量中。因此spmv的计算的过程是一个完全向量化的过程，不存在额外的存取非连续下标的向量元素操作。而且s
vvec
是可变参数，能使计算过程不改变的情况下适合不同宽度和硬件架构的向量化指令。
66.该方法在不使用汇编代码的情况下能优化spmv代码并获得低级别simd指令的优势的原因是完全向量化的spmv算法很容易被编译器自动优化并生成所需的simd指令代码。
67.本实施例的有益效果：实施例1的矩阵数据架构形式的矩阵格式结合加速基于压缩稀疏列的加速spmv的方法，能够高效利用硬件向量化计算指令和内存带宽，使得spmv能够以更快的的速度运行，达到加速的目的。
68.实施例3
69.一种基于压缩稀疏列的加速spmv的方法的实施3，基于实施例2，用于面向影像重建过程，其线积分算子方程离散化矩阵的局部向量重排映射(简称为ioblr)是将y向量按照投影数据区域上的i、j坐标轴排序布局变成按照平行和垂直于某一个像素点在投影数据区域的轨迹排序布局；矩阵每一列的非零元素根据y向量重排映射进行分块排列和存储。矩阵的排列称为ioblr优先布局。ioblr的逆映射是将y向量重新排列成为按照投影数据区域上
的i、j坐标轴平行的布局。
70.本实施例中选取的ioblr在投影数据上的显示如图3所示。图3的箭头折线表示了y向量的排列所用的ioblr，而同一条箭头折线上的方格子表示矩阵元素在ioblr下重排后的结果。图3显示了图1示例子矩阵块中三列转换为实施例1的矩阵数据架构前后的内存布局。在转换为矩阵数据架构之前，每一列的非零元素对应于一个像素投影轨迹上的点，用蓝色格表示。列中的非零元素存储在以探测单元为主的布局中。转换后，以矩阵数据架构记录的非零元素以红色平行折线上的网格表示。转换之后，新添加的填充零用黄色格表示。每个矩阵块元素(cscve)对应一个箭头折线段。矩阵块元素(cscve)中的非零值存储在ioblr为主的布局中，其中沿箭头方向的元素在内存中是连续的。平行折线的形状由参考像素的最小探测单元数量曲线决定，参考像素由像素块的中心点决定。图1中不同的颜色代表不同的像素点，黄色和灰色表示多个像素点共同的轨迹。像素点代表着spmv中x向量元素的在影像重建问题中的几何分布，而投影数据代表着spmv中y向量元素的在影像重建问题中的几何分布，同一像素点轨迹的轨迹代表着spmv中同一列矩阵非0元素在影像重建问题中的几何分布。图2中蓝色格子和黄色格子分别表示系统矩阵中的非零元素和零元素；箭头折线中的彩色格构成cscve。折线显示了ioblr中参考像素的轨迹。
71.在本实施例中，选取的ioblr在投影数据上的显示如图3所示，包括了行向量y的重排映射，以及矩阵元素在行向量重排下的内存布局(简称为ioblr优先存储布局)。在图4中，三种不同内存布局：探测单元优先、角度优先、ioblr优先，由不同颜色的箭头线表示。simd效率由箭头线和蓝色网格的交点的面积定义，其中蓝色网格表示矩阵中的非零元素。当s
vvec
为8时，探测单元优先、角度优先、ioblr优先的区域范围分别为3,2-6和7-8。
72.ioblr为基于压缩稀疏列的spmv算法提供了高效的y向量数据局部重排和矩阵元素重排方案以及快速确定矩阵块结构参数的能力。ioblr的有效性源于影像重建过程中积分算子的几何性质。该性质在不同大小的成像矩阵中都存在。另外，选取局部的ioblr可以无需后验分析就可获得一个有效的排序方案。因此cscv矩阵块结构参数的选择，无需逐案分析；只要对每一种硬件平台根据影像的成像几何选取一个成像矩阵并进行少量的参数组合spmv实测，比较运行时间就可以获得运行参数。
73.在本实施例中，基于压缩稀疏列的spmv算法的软件实现包括了两个分支：cscv-z和cscv-m。基于压缩稀疏列的spmv算法可以结合其他优化方法进行调整，例如填充0元的去除技术。这两种实现的计算过程只有一点不同，那就是是否去除cscve中的填充0元。cscv-m采用了填充0元的去除技术，可以进一步减少内存使用量，在矩阵cscve的读取过程中要进行vexpend操作；而cscv-z的内存带宽使用率更高，无需额外的矩阵元素处理。cscv-m的输入矩阵块采用定长稀疏vxg块格式，而cscv-z的输入矩阵块采用定长稠密vxg块格式。
74.在本实施例中，cscv-z和cscv-m两个实现的单线程计算过程(包括局部的临时重排序和计算过程中完全向量化的spmv计算)如图5所示。
75.为了展现cscvspmv计算过程的性能优势，本实施例选择了两个x86架构的硬件平台以及10个当前主流的spmv软件实现进行对比测试。测试使用的案例是用平台上全部的cpu核心根据平行光束几何ct重建一张1024*1024像素大小的影像。测试过程只进行前向投影计算，并且计算的矩阵大小为700800*1048576，非零元素个数约为13亿。测试硬件平台包括了：intel xeon gold 6130(简称为skl)和amd epyc 7452(简称为zen2)。两种平台的向
量指令宽度是不同的。skl平台使用的是512bit宽度的向量化指令，而zen2平台使用的是256bit宽度的向量化指令。两种平台的cpu核心数目是不同的，skl为32，zen2为64。
76.本实施例的测试使用的spmv软件实现包括：cscv的两个软件实现cscv-z和cscv-m，以及spc5、cvr、merge、csr5、vhcc、esb、intel mkl csc(简称为mkl-csc)和csr(简称为mkl-csr)。测试中还包含了两种数据精度，单精度和双精度。因为医院的临床ct成像都是使用单精度，所以该种测试案例更加重要。
77.本实施例的测试主要考察三个方面的内容：矩阵存储的内存使用量、每秒浮点运算量和只读访存带宽利用率。具体的测试数据见图6。根据测试数据可得，cscv-m的内存使用量m
rit
(gb)降低为mkl-csr的53％～68％，每秒浮点运算量performance(gflops)可达mkl-csr的3.51倍(即计算速度是3.51倍)，只读访存带宽利用率r
em
(％)大于85％且最高为98％。每秒浮点运算量的前两名均为cscv的软件实现。和实验的第三名相比，每秒浮点运算量为后者的1.6～2倍。对短字长数据(单精度数据)，加速倍数更明显。
78.为了展示cscv spmv具有适应宽simd指令以及在不使用汇编代码的情况下具可获得低级别simd指令的优势，我们比较cscv-m和spc5的每秒浮点运算量。在本实施例的测试中，cscv-m没有使用汇编代码，而每秒浮点运算量明显是使用底层汇编代码的spc5的1.6倍，而且可以跨硬件平台使用。spc5的spmv计算内核是使用intel avx-512汇编代码编写向量化计算过程并且手动调优的，只能在intel的平台运行。cscv-m只是采用高级语言c++编写，用编译器的自动优化实现向量化过程。cscv spmv计算流程具有全向量化计算内核，容易被编译器优化，因此仍然有很好的性能，而且可以跨硬件平台使用。
79.为了展现本实施例的方法的计算过程具有快速确定块结构参数的优势，本实施例的引入不同4个不同大小的ct成像矩阵进行与前面相同的性能测试。矩阵的参数见表1。cscv的软件实现采用与前面相同的参数(具体见表2)，包括了向量长度s
vvec
、块的大小s
imgb
和vxg的长度s
vxg
.各种软件的在四个矩阵测试的每秒浮点运算量统计见表3。结果显示，cscv软件实现具有明显的计算优势，每秒浮点运算量的最大数值均比对手优胜。
80.表1待测试的4个ct成像矩阵的参数
[0081][0082]
表2不同软件实现在不同平台上运行采用的块架构参数
[0083][0084][0085]
表3
[0086]
[0087]
显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。