CN104636273B - 一种带多级Cache的SIMD众核处理器上的稀疏矩阵存储方法 - Google Patents

Abstract

本发明公开了一种带多级Cache的SIMD众核处理器上的稀疏矩阵存储方法，包括：(1)获取矩阵A中行非零元个数最大值a，处理器SIMD单元可同时计算的非零元个数b，计算大于a且为b的倍数的最小值作为临时行宽度；(2)对矩阵A，数组Value和Colidx分别顺序存储每一行的非零元值和列坐标，每行非零元个数不足临时行宽度的，在该行尾分别补0和‑1；(3)对Colidx和Value按b列分块；(4)每个列块按行压缩，使列块中有非零元的行集中于列块上部；(6)列块按b行分块，得到子块；(7)去除全零子块，将子块按行存储。本发明的存储方法将稀疏矩阵分成稠密子块，在保留非零元局部性的基础上，提高了对处理器SIMD处理单元以及寄存器的利用率，提高了稀疏矩阵向量乘的性能。

Description

一种带多级Cache的SIMD众核处理器上的稀疏矩阵存储方法

技术领域

本发明涉及并行程序设计领域，特别涉及一种带多级Cache的SIMD众核处理器上的稀疏矩阵存储方法。

背景技术

稀疏矩阵向量乘(SpMV)是很多科学和工程应用的一个重要的计算核心，它的计算效率是科学和工程应用的计算性能的关键。该算法的主要功能是计算y＝y+Ax，其中A是一个二维稀疏矩阵，x和y都是一维稠密的向量。然而该算法核心在现代具有多级Cache的SIMD众核处理器上，由于稀疏矩阵非零元分布的不规则性，其SIMD利用率很低，导致SpMV性能较差。要想提高该算法的性能，我们往往需要综合考虑输入集和计算平台的特性，然后选择合适的稀疏矩阵存储格式，这些格式可以有效利用具体计算矩阵的数据特点和计算平台的结构特点。

随着多核CPU的快速发展，其单指令多数据SIMD计算单元的计算宽度也越来越大。例如，Intel MMX指令提供了64位宽度的SIMD支持，Intel SSE指令则支持128位的SIMD支持，Intel AVX指令支持256位的SIMD操作。近年来英特尔推出了Intel Xeon Phi协处理器，其提供了512位的SIMD计算单元。这些处理器的共同特点是：每个计算核心都具有能够处理长向量的宽SIMD计算单元，并具有多层高速缓存(Cache)结构。

由于这些处理器拥有宽SIMD计算部件，所以其非常适合稠密矩阵的计算。而稀疏矩阵向量乘要想在这种带多级Cache的SIMD处理器上获得很高的计算性能，需要克服由于稀疏矩阵非零元分布不规则导致的计算瓶颈：

(1)SIMD利用率较低；

(2)x向量中的数据重用率低，使得高速缓存缺失和访存延迟开销很大；

本发明充分地利用了带多级Cache的SIMD众核处理器的各种硬件特性，如每个核上的宽SIMD计算能力，多级高速缓存等，提出了一种对稀疏矩阵进行分块压缩的存储办法，称为ERB(ELLPACK Register Blocking)，该方法是一种寄存器分块优化方法。该方法有效地解决了稀疏矩阵向量乘在带多级Cache的SIMD众核处理器上的计算瓶颈，提供了一种可以高效地利用其结构特征的稀疏矩阵存储格式。

发明内容

本发明提出一种带多级Cache的SIMD众核处理器上的稀疏矩阵存储方法ERB(ELLPACK Register Blocking)，之前的工作表明，基于ELLPACK的稀疏矩阵向量乘算法可以有效地利用宽SIMD计算能力，而且ELLPACK的变种在CPU、GPU和MIC等具有宽SIMD计算能力的架构上都表现出了不俗的计算性能。所以，本发明在ELLPACK的基础上，通过对以ELLPACK格式存储的稀疏矩阵进行分块和压缩等操作，提出一种可以有效解决背景技术中所述瓶颈的寄存器分块优化的稀疏矩阵存储方法。

本发明实现的针对具有宽SIMD计算单元的处理器的稀疏矩阵存储方法ERB的具体实现步骤如图1所示，本发明的格式转换过程共有5个步骤，具体如下：

(1)特征提取及矩阵扫描，按行扫描稀疏矩阵，获取稀疏矩阵中非零元个数最多的行的非零元个数a以及处理器SIMD处理单元中可同时计算的非零元个数b，图2、图3、图4中假设b为2，计算得到大于a且为b的整数倍的最小值作为临时行宽度。按行扫描稀疏矩阵A，用数组Value顺序存储稀疏矩阵每一行的非零元，每行非零元个数不足临时行宽度的，在该行补0。用数组Colidx存储Value中每一个非零元对应的列索引值，对Value中的零元(Value中存储全部非零元，但由于每行非零元个数不同，所以在每行最后补0，使每行长度相等，补的0即此处的零元)，在Colidx中存储为-1。稀疏矩阵A扫描结果如图2所示。

(2)列分块，将步骤1得到的稀疏矩阵按列宽为b进行按列分块。由于下面的转换过程同时应用于Value矩阵和Colidx矩阵，并且操作完全一致，故以下仅以Value数组的变换作图阐明。在步骤1得到的稀疏矩阵中，将Value矩阵和Colidx矩阵按同一固定列宽度为b分块。

(3)列压缩，在按列分块的基础上，对每一个列块，将其中全零的行去除，并将其余行依次向上压缩。Value矩阵的列分块和列压缩的过程如图3所示。图中列块1中，第2行的元素都是零元，则将第三行及其后的行向上压缩。

(4)行分块，在列压缩的基础上，对每一个列块按行分块大小为b进行按行分块得到子块。并用相同的办法处理子块的非零元的列信息矩阵，此外还要保存子块每一行的行信息。Value矩阵的行分块过程如图4所示。

(5)按行存储，将子块内的数值按行存储，并将所有的块按行存储，则得到ERB存储格式的稀疏矩阵。

在按行存储后，转换过程执行完毕，稀疏矩阵最终的存储格式如图5所示，接下来就需要将分好的子块平均分配到处理器各核上进行稀疏矩阵向量乘的计算。

与现有技术中存在的问题相比较，本发明所具有的优点和积极效果主要体现在：

(1)本发明通过对稀疏矩阵进行列分块和列压缩，使得矩阵中的非零元集中，在最后ERB存储格式的每个子块中，非零元的密度高。相比较目前在具有宽SIMD计算单元的处理器上的稀疏矩阵向量乘算法可以提高每个计算核的SIMD寄存器和SIMD计算单元的利用率，是一种寄存器分块优化的存储方法；

(2)本发明通过对稀疏矩阵进行细粒度的划分，然后对每个块进行稀疏矩阵向量乘的过程中保留了稀疏矩阵非零元之间固有的局部性，相比较目前基于ELLPACK的优化办法，有效地保留了稀疏矩阵非零元固有的局部性，从而很好地利用了具有宽SIMD计算单元的处理器结构高速缓存中的局部性；

本发明提出的方法相较于Intel MKL中已有的稀疏矩阵向量乘算法，算法运行时的实际表现好很多，而且本优化方法对稀疏矩阵的模式没有限制，是一个通用的能够有效提高稀疏矩阵向量乘算法在带多级Cache的SIMD众核处理器上计算效率的存储方法。

附图说明

图1、针对具有宽SIMD计算单元的处理器的稀疏矩阵存储格式转换过程示意图；

图2(a)、稀疏矩阵A扫描结果示意图，普通的稀疏矩阵A；

图2(b)、为扫描后的结果格式；

图3(a)、矩阵的Value矩阵列分块及列压缩示意图，按列块大小为2进行按列分块，形成3个列块；

图3(b)，在每个列块中，分别进行按列压缩，列块1的第2行，列块2的第1行和第2行由于是全零的行，则被其后的行所压缩；

图4、矩阵的Value矩阵行分块示意图，图中先按行块大小为2进行按行分块，去掉全为0的块后，则剩余4个子块；

图5、稀疏矩阵ERB存储格式示意图，转换结果由三个数组组成，其中value数组按列存储矩阵的非零元的值；colidx数组按列存储矩阵非零元的列坐标；row_perm数组存储每个块中每行的原始行坐标。

具体实施方式

本节将本发明应用于一个典型的在带多级Cache的SIMD众核处理器上的稀疏矩阵向量乘计算。从而进一步说明本发明的目的、优点和关键技术特征。这个实施只是该方案的一个典型范例，凡采取替换或者等效变换而形成的技术方案，均落在本发明要求保护的范围之内。

对于一个要计算的稀疏矩阵A：

先将矩阵A按本发明提出的优化方法经过特征提取及矩阵扫描、列分块、列压缩、行分块、按行存储。转换成如图5的ERB存储格式，或者直接将要计算的稀疏矩阵存储成如图5的存储格式，并保存到文件中。

在计算时，从文件中读入图5中存储的稀疏矩阵，然后根据小子块大小确定每个子块所在的位置从Value数组中读入各子块的非零元的值，从Colidx中读入各子块中非零元的列坐标信息，从row_perm中读入各子块的行信息。然后将各子块按固定的分配算法分配到不同的核上去计算。在每个子块内部，一次读取一行放入到SIMD计算单元中去计算，然后将结果加入到y向量中。其中x向量和y向量参与计算的元素分别由当前非零元的列坐标和行坐标决定。

Claims (3)

1.一种带多级Cache的SIMD众核处理器上的稀疏矩阵存储方法，其特征在于：将稀疏矩阵通过五次矩阵变换后进行存储，包括以下步骤：

特征提取及矩阵扫描：获取稀疏矩阵中非零元最多的行的非零元个数a以及处理器SIMD处理单元中可同时计算的非零元个数b，计算得到大于a且为b的整数倍的最小值作为临时行宽度，按行扫描稀疏矩阵A，用数组Value顺序存储稀疏矩阵每一行的非零元，每行非零元个数不足临时行宽度的，在该行后补0，用数组Colidx存储Value中每一个元素对应的列索引值，对Value中的零元，其在Colidx数组中对应的值为-1；

列分块，将特征提取得到的b的值作为列分块宽度，对Colidx数组和Value数组按列划分为多个列块；

列压缩，在按列分块的基础上，对每一个块，将其中全零的行去除，并将其余行依次向上压缩；

行分块，在列压缩的基础上，对每一个列块在Value和Colidx的值按行宽为b进行按行分块得到子块，并用row_perm数组保存每个子块的行索引信息；

按行存储，将子块内的元素按行存储，并将所有的块按行存储，则得到ERB存储格式的稀疏矩阵。

2.根据权利要求1所述的方法，其特征在于：在计算时，从文件中读入与稀疏矩阵相关的Value，Colidx和row_perm数组，然后根据子块大小确定每个子块所在的位置从Value数组中读入各子块的非零元的值，从Colidx中读入各子块中非零元的列坐标信息，从row_perm中读入各子块的行信息；然后将各子块按固定的分配算法分配到不同的核上去进行稀疏矩阵向量乘的计算。

3.根据权利要求1所述的方法，其特征在于：矩阵变换过程中参数包括：

1)稀疏矩阵：行数，列数，非零元个数，每行非零元个数，按行和按列分块大小；

2)处理器：SIMD计算单元宽度；

3)数据类型：单精度浮点类型，或双精度浮点类型。