[发明专利]针对稀疏矩阵的数据存储方法及基于该方法的SpMV实现方法

说明书

技术领域

本发明涉及一种改善稀疏矩阵向量乘性能的数据存储方法及SpMV实现方法。

背景技术

稀疏矩阵向量乘(SpMV，y＝Ax)是一个十分重要的，且经常被大量调用的计算内核，广泛应用在科学计算、信息检索、气象、航天、油藏模拟、天体物理、数据挖掘等科学计算和实际应用中，而且被重复调用的次数常常会达到成千上万次。但是现代基于Cache存储层次的计算平台上，传统CSR存储的稀疏矩阵向量乘，运行效率往往低于硬件浮点运算峰值的10％。对稀疏矩阵向量乘性能优化，可以大大改善整个工程计算的运行效率。因此优化稀疏矩阵向量乘的性能成为提高工程效率的关键，在实际应用中有着十分重要的意义。

针对稀疏矩阵向量乘的常用存储方法是CSR，其使用3个数组表示一个m×n，含有nnz个非零元的稀疏矩阵：val和colval分别保存每个非零元的取值和列值，rowptr保存每行第一个非零元在val或colval中的索引位置，见图1，其SpMV实现表1：

表1基于CSR存储方法的SpMV操作的代码表

基于CSR的SpMV实现运行效率很低，主要原因是浮点计算操作和存储访问操作比率非常低，稀疏矩阵非零元分布的不规则性导致对x访问的不规则，且每次对x访问都需间接访问colval数组。对x的访问是由非零元在稀疏矩阵的位置决定的，不同应用间使用稀疏矩阵的多样性增大了高效通用SpMV实现的难度。

针对应用矩阵非零元分布的特点，合理优化对x访问，人们提出了很多优化策略。我们分下面4类进行论述。

(1)对矩阵进行分块

Vuduc 2005年设计开发OSKI，利用稠密线性代数库中常用的分块和循环展开的方法(blocking and tiling)，使用便于寄存器分块的BCSR存储模式。BCSR与CSR类似，CSR以非零元为存储单位，而BCSR以位置相近的非零元组成的r×c非零元块为存储单位，块中的非零元共享同一个索引信息，见图2。

BCSR为了维持非零元块的结构，需要显示填充零元。填充零元有两个缺点：存取零元需要占用额外的带宽资源；零元的浮点计算需要占用额外的计算资源。为减少非零元的填充，OSKI可以在运行时根据矩阵的信息，动态的选取最优分块大小：在安装时，对具体运行环境不同大小的分块进行运行速率测试，记作PerfDenser×c；同时用非零元填充率衡量填充非零元的数目，即转换成BCSR后非零元的数目与原CSR非零元数目的比值，记作FillRatior×c，该值与运行的具体稀疏矩阵相关，在运行时为节省时间，对矩阵进行抽样计算得到。OSKI对分块大小为r×c的评估函数为：

并选取使评估值最大的分块大小作为最终的分块大小。同时，[Vuduc et al 2002]提出了寄存器分块进行SpMV计算的速率模型。以上过程就是OSKI对任意给定稀疏矩阵，在运行时自适应地实现SpMV的过程。

寄存器分块有减少索引数据的优点，同时在一定程度上重复利用存入寄存器中的x。因此可以将寄存器分块的优点归纳为3点：是一种对索引信息进行的压缩策略，减少了cache与主存的通讯量；在分块的行大小大于1时，同列的非零元与同一个x相乘，使得x重用，这一特点不与分块的列大小相关；分块的列大小仅与稠密线性代数中运用到的循环展开类似，决定在一次循环内部操作的操作次数。但BCSR在存储非零元子块时存在以下限制：

第一，需要索引坐标的行列值必须是分块大小的行列大小的倍数。会导致在原矩阵稠密分块在坐标不对齐，需要填充整列或整行的零元，如图3(a)。

第二，整个稀疏矩阵，使用单一分块大小。原始矩阵小的稠密块必须填充更多的零元去维持分隔块的结构，如图3(b)。