专利名称：稀疏矩阵的对角线数据存储方法及基于该方法的SpMV实现方法
技术领域：
本发明涉及一种针对稀疏矩阵的数据存储方法及基于该方法的SpMV实现方法， 尤其涉及针对稀疏矩阵的对角线存储方法和基于该方法的SpMV实现方法。
背景技术：
稀疏矩阵向量乘(SpMV)y = A*x是一个重要的计算内核，广泛应用在信号处理、图 像处理、迭代求解算法等科学计算和实际应用中.但在基于cache存储层次的计算平台上， 与稠密矩阵向量乘相比较，稀疏矩阵向量乘的性能很差，主要是因为cache存储层次的复 杂性和稀疏矩阵非零元素分布的不规则性。浮点计算操作和存储访问操作的比率非常低， 尤其是向量χ的间接访问和访问的不可重用性，使得传统的压缩行存储方法得到的稀疏矩 阵向量乘性能很差，往往低于机器浮点运算峰值的10%。稀疏矩阵向量乘作为内核如果能 够提高运行速度，整个工程计算的运行效率将会得到大大地改善，时间也会大大地减少，在 实际应用中有着十分重要的作用.例如，在天体物理项目中，因为要多次调用SpMV代码，稀 疏矩阵向量乘的优化能明显提高项目的整体速度。实际应用中大部分稀疏矩阵具有对角线模式，为此，现有的研究已提出了 DIA以 及RSDIA等系数矩阵存储方法，但在应用过程中都有不足之处。稀疏矩阵是指矩阵A的元素大部分是零，而非零元素所占比例非常小，往往小于 元素总数的1%.通过只存储和操作这些非零元，达到内存和计算代价上的节省.稀疏矩阵 存储时，除了存储非零元外，还要记录非零元在矩阵中所占的位置.稀疏矩阵的一般存储 方法是CSR(compressed spar se row)方法。CSR方法需要存储稀疏矩阵A的每个非零元 素的值，非零元所在的列和每行第1个非零元的索引，即需要3个数组(其中矩阵A是mXn 矩阵，有nz个非零元)，如下所示· val [nz]，记录每个非零元的值；· col [nz]，记录每个非零元所在的列；·ρ ι·[πι+1]，记录每行的第1个非零元在数组Val[nZ]C0l[nZ]中的索引，其中 ptr[m] = nz.稀疏矩阵向量乘的形式为y = Ax，其中A为稀疏矩阵，x，y为稠密向量.稀疏矩 阵向量乘的一般实现算法运行效率很低，往往低于机器浮点运算峰值的10%，其中一个原 因是，数据引用的时间局部性和空间局部性比较差，尤其是向量χ的间接访问和访问的不 可重用性；另一个原因就是算法中的浮点计算与存储访问的比率太低，浮点操作和load指 令混杂导致把大量的时间花费在数据的存储访问中.提高稀疏矩阵向量乘的计算效率，优 化稀疏矩阵向量乘算法是很困难的，这是由存储层次的复杂性和存储访问的不规则性决定 的.存储系统分为寄存器、一级cache、二级cache、内存，从前到后存储空间越来越大，访问 数据所需要的时间也越来越长.所以应尽量使得数据在高层存储里面得到重复使用，减少 再次从低层存储访问的次数和开销，从而提高运行效率.
由于CSR中每次计算乘加都需要前一次计算的结果yt，且浮点运算是多级流水， 产生了 stall，为此，出现了 L-CSR 结构，应用 unroll-and-jam (loop jamming)方法，将拥 有非零元个数一样的行存放在一起为一组，并保存一个长度索引数组指向所有不同长度的 组，使用2级imroll-and-jam并对长度为2和大于2的情况分别行编写了代码。要求每行 的非零元的长度在一个小范围内，可能破坏CSR算法中访问向量χ的局部性，尤其是当矩阵 为对角线模式时，并且，写入向量y的步长可能为1。而CSR-DU算法减少了索引大小和访问带宽需求，只存储同一行中每个元素与前 一个元素列索引的差值，并将roW_ptr和col数组合并，在压缩的索引中标识新行的开始， 在计算SpMV是解压缩索引，用减少的带宽需求来覆盖解压缩的时间达到提升性能的目的。 但是要求矩阵不存在大量非零元较少的行。CSR-VI减少了 val数组大小，只保存不同的val 值，并为原始CSR中的所有非零元建立一个间接访问val数组的索引，只在不同非零元值数 目较小时有用，它也减少了利用CPU进行流水计算的可能性。对称存储方法针对对称矩阵，寄存器分块(大小为r*c)和向量分块(矩阵乘以 多个向量，对向量进行分块，大小为v，以利用矩阵A的局部性)进行优化。提出了一个类 似Sparsity中的模型分析r，c和ν 给定机器参数，对称矩阵A和向量个数k，首先进行 off-line benchmark，对稠密矩阵用稀疏存储方法，所有的1彡r，c彡8和1彡ν彡10的 所有组合都进行测试，记录峰值P (r，c, ν)；然后进行run-time search,选择A中部分行，对 所有的1 < r，c < 8计算非零元填充率f (r，c)；最后用启发式算法，选择P (r，c, v)/f(r, c)最大的r，c和ν值。改进的CSR还有DCSR和RPCSR。DCSR使用6个command code进行模式压缩，将3 个连续的模式压缩为四字节，称为一组，非常适用于非零元多为连续的情况。并将roW_ptr 和col数组合并，在压缩的索引中标识新行的开始，在计算SpMV是解压缩索引，用减少的带 宽需求来覆盖解压缩的时间达到提升性能的目的。RPCSR对四字节的组进一步压缩，选取出 现多次的模式编写专门代码优化。Belgin等提出了基于I^ttern的压缩选择分块大小r，c，共有2~(r*c)中块的非 零元模式，忽略少于3个非零元的块和出现次数少于1000的块模式，这些元素用CSR存储， 选择r，c使得尽量减少用CSR存储的元素，为每一个块模式编写代码，并利用SSE2/3指令 对每个r*c大小的块中的2 大小的块编写SSE代码。Mellor-Crummey等使用寄存器分块，cache分块和多向量乘优化SpMV。BCSR减少 索引存储空间和访问带宽需求，减少向量χ间接访问次数。增加了访问块内零元的开销，增 加了带宽需求和计算量。Geus等对对称稀疏矩阵向量乘使用软件流水(显式预取数据，即先取下一次计算 的数据，然后进行本次计算，增加指令并行度)，寄存器分块和矩阵重排三种优化方法。Buluc等提出了 CSB存储结构以及并行算法。分块大小为b，共有rT2/V2个块， 存储每个块在val数组中的起始位置的索引，在每个块内用三元组存储，这样可以减少行 和列索引所需的大小。Toledo等认为CSR算法有四个性能瓶颈对val和col是顺序访问，可能会引起一 些Cache缺失，或者是Ll Cache,或者是更低层次的Cache缺失；对χ的访问缺少时空局部 性，导致大量Cache缺失；每次浮点乘加操作需要加载三个浮点数，所以CPU的load-store单元的性能十分关键；利用col数组访问Χ需要间接寻址。Pinar等基于Toledo的工作，引入变长的分块，BCRS，l*m分块，用大小为分块数目 的Nzptr数组存储每个l*m分块的在val数组中的起始位置。Row_ptr指向每行第一个块 在阪Ptr中的位置，算法使用三重循环。首先是对行循环，然后对每行中的块循环，最后是 块内循环。但是BCRS效果并不好。文章进一步提出了矩阵重排的方法，增加矩阵每一行中 的稠密块大小。将矩阵重排问题化为TSP问题。用TSP的启发式算法解决。稀疏矩阵中两 列在相同的行有非零元的数目做为这两列的权重，可以通过计算W = A*(A的转置)，A为稀 疏矩阵，1代表非零元，则W矩阵为权重。最后，重排后的矩阵的分块数目=矩阵非零元数 目-TSP路径的权重(即所有相邻列的权重)。TSP两大类启发式构建和改进。构建直接 构建一个解，然后改进从一个初始解开始，寻找邻居节点进行改进。Rice等提出了 ELLPACK/ITPACK (ELL)格式，两个相同大小的数组，都为m*mnz，mnz 是每行拥有最大非零元的个数，少于mnz个非零元的行用O来填充，另外一个数组存储每行 所有非零元的列值，对于填充的0，存储-1为列值。Bell等在GPU上对SpMV进行优化，一个warp处理一行，每个线程处理一个乘法， 结果存储在local memory中，然后进行规约操作。并采用ELLPACK和COO混合格式进行存储。Kourtis设计了新的SpMV算法，压缩相同非零元，并设置索引数组来检索非零元， 由于某些矩阵只存在较小的不同非零元个数，因此，这种方法能够降低算法的带宽需求。
发明内容
本发明欲解决的技术问题是提高具有对角线模式的稀疏矩阵SpMV实现的性能， 提供一种针对稀疏矩阵的对角线存储方法，该方法明显减低了存储空间需求和计算时间。原始的针对对角线模式的稀疏矩阵存储方法最主要的为DIAG，它是将所有非零元 对角线上的元素经过零元填充与主对角线对齐后，按照编号从小到大连续存储到val数组 中。进行SpMV实现时一条对角线算一次，在对角线一级进行循环。但是这种方法会多存储 矩阵范围外的非零元，对于远离主对角线的对角线，增加了存储矩阵外虚位的负担。本存储 方法只存储矩阵范围内的非零元及相应对角线上的零元，为了类似DIAG —样存储在寄存 器中重用y，需要获得计算每个y需要的对角线数目。根据过对角线与矩阵两侧交点的水平 线对矩阵进行切分，使得每个子稀疏矩阵行块内所需要的对角线相同，因此计算模式一致， 可以重用y。这就是本发明的稀疏矩阵的对角线数据存储方法，简称DDD-SPLIT算法。由于 这种改进算法是按行来计算，所以对稀疏矩阵中存在非零元的对角线，首先填充零元为非 空元，设置非空元的值为0，然后将所有非空元，按照类似CSR的存储方法，按行主序进行存 放，每行内的非空元按照列值进行顺序存放。此外，通过大量观察发现，很多稀疏矩阵对角线上存在相同非零元，因此，通过对 子稀疏矩阵进一步切分，可以分为更细的可压缩子稀疏矩阵，对每个可压缩子稀疏矩阵，只 需要储存第一行的非零元即可，进一步降低了储存空间需求。本发明的技术方案为，一种对角线数据存储方法，包括如下步骤1)按行扫描稀疏矩阵A，确定具有非零元的对角线(即非零元对角线)的位置，以 对角线编号来表示；
2)除了主对角线，每条对角线都与矩阵的两侧交于一点，上三角的对角线交矩阵 的右侧，下三角的对角线交矩阵的左侧，用非零元对角线与矩阵侧边的交点水平切分矩阵 (可能会有左右两侧同一点同时切分矩阵)，主对角线不对矩阵进行切分，最多有两条对角 线在同一行对矩阵进行切分。将矩阵切分为若干子稀疏矩阵行块，将每个行块视为一个独 立的子稀疏矩阵，每个子稀疏矩阵内计算均需要相同的矩阵对角线，即在每个子稀疏矩阵 中，每个非零元的位置都在相同的对角线上。3)扫描子稀疏矩阵，按行顺序存储每个子稀疏矩阵中的非零元对角线上的元素到 val数组。计算并存储每个子稀疏矩阵中非零元对角线的编号和数目。稀疏矩阵A中设置对角线，比如nXn的矩阵共有2η_1条对角线，η为正整数，其 编号如
图1所示，上三角和下三角中对角线编号分别用Di和D_j表示，i彡0，j彡0，i和j 为整数，i为上三角对角线中第一个元素的列号减1，j为下三角对角线中第一个元素的行 号减1。则主对角线编号i = j = 0，为Dtl，Di对角线的大小为n-l i ,D_j对角线的大小为n_| j |，图1的右侧说明了计算每条 对角线需要的相应的向量χ和y的范围。扫描子稀疏矩阵时，若相邻行同一对角线上的元素相同，则具有相同元素的相邻 行被进一步切分为一可压缩的子稀疏矩阵，在val数组中仅需存储该可压缩子稀疏矩阵一 行的非零元对角线的元素。此时稀疏矩阵的数目为进一步切分为可压缩的稀疏矩阵后的矩 阵个数。所述val数组的存储方式为按行顺序存储非零元对角线上的所有元素，包括非零 元和零元。所述稀疏矩阵A为对角线模式的稀疏矩阵，S卩非零元主要集中在若干条对角线 上。在此过程中，还存储若干个数组及数据对角线数目；diag_indeX[diag_Cnt]，所 有对角线的编号；row_b 1 ock_cnt，划分的子稀疏矩阵的个数；row_cnt [row_b 1 ock_cnt]， 每个子稀疏矩阵内的行数；compressed [row_block_cnt]，每个子稀疏矩阵是否为可压缩稀 疏矩阵，row_region_diag[row_block_cnt]，每个子稀疏矩阵所用到的对角线范围。由于对 角线模式的diag_Cnt与r0W_bl0Ck_Cnt个数相比非零元个数非常小，所以这些数组的影响 可以忽略不计。基于上述的稀疏矩阵数据存储方法，下面提供一种SpMV实现方法1)遍历稀疏矩阵中的每一个子稀疏矩阵，计算每个子稀疏矩阵的稀疏矩阵向量乘 y = A1^x, A1为第1个子稀疏矩阵，1为正整数。第1个子稀疏矩阵的SpMV实现方法为a)根据该子稀疏矩阵中m条非零元对角线的编号和当前计算行的行号j+Ι，确定 第j+Ι行χ的元素索引；m为正整数。b)从val数组中读取第j+Ι行的m个元素值，与χ的元素索引对应的m个χ值分 别相乘后累加，得到c)继续计算y押；2)合并所有子稀疏矩阵计算得到的稀疏矩阵向量乘，得到稀疏矩阵A的稀疏矩阵 向量乘。
从val数组中读取第j+2行的m个元素值，与χ的元素索引+1对应的m个χ值分 别相乘后累加，计算得到；根据r0W_regi0n_diag[l]和diagjndex数组，得到第1个 子稀疏矩阵的非零元对角线的个数m及位置，根据这m个对角线位置和当前计算的y的行 号，可以确定该子稀疏矩阵内第j+Ι行每个非零元所对应的χ元素索引，而剩余的其他行的 X位置，由于对角线是连续存储的，所以可通过前一行相同对角线上非零元所对应的X索引 加1得到。所述子稀疏矩阵为可压缩子稀疏矩阵时，j+Ι行的元素值与第j+2行的元素值相 同，val数组中仅存储第j+Ι行的元素值，计算时读取val数组存储的该可压缩子稀疏 矩阵的m个元素值，与χ的元素索引加1对应的m个χ值分别相乘后累加。由于可压缩子 稀疏矩阵内的所有行的非零元均相同，因此只需要读取一次val数组。根据r0W_regi0n_ diag[l]和diagjndex数组，得到第1个子稀疏矩阵的非零元对角线的个数m及位置，然后 从val数组中读取下面m个元素，并用这m个元素计算该矩阵内所有行。所述第j+Ι行χ的元素索引由非零元对角线的编号加上j+Ι获得。该稀疏矩阵的SpMV实现方法的代码可见表1 表1基于对角线数据存储方法的SpMV实现代码表
权利要求
1.一种稀疏矩阵的对角线数据存储方法，包括下列步骤1)按行扫描稀疏矩阵A，确定具有非零元的对角线的位置，以对角线编号表示；2)以非零元对角线与矩阵A侧边的交点作水平线将矩阵A水平切分为多个子稀疏矩阵；3)扫描子稀疏矩阵，按行顺序存储每个子稀疏矩阵中的非零元对角线上的元素到val 数组；计算并存储每个子稀疏矩阵中非零元对角线的编号和数目。
2.根据权利要求1所述的稀疏矩阵的对角线数据存储方法，其特征在于稀疏矩阵A的 上三角和下三角中对角线编号分别用Di和Dy表示，i彡0，j彡0，i和j为整数，i为上三 角对角线中第一个元素的列号减1，j为下三角对角线中第一个元素的行号减1。
3.根据权利要求1所述的稀疏矩阵的对角线数据存储方法，其特征在于扫描子稀疏矩 阵时，若相邻行同一对角线上的元素相同，则具有相同元素的相邻行被切分为一可压缩子 稀疏矩阵，在val数组中仅存储该可压缩子稀疏矩阵一行的元素。
4.根据权利要求1或3所述的稀疏矩阵的对角线数据存储方法，其特征在于所述val 数组的存储方式为按行顺序存储非零元对角线上的所有元素，包括非零元和零元。
5.根据权利要求1所述的稀疏矩阵的对角线数据存储方法，其特征在于所述稀疏矩阵 A^nXn的稀疏矩阵，η为正整数。
6.根据权利要求1或5所述的稀疏矩阵的对角线数据存储方法，其特征在于所述稀疏 矩阵A为对角线稠密的稀疏矩阵。
7.一种基于稀疏矩阵的对角线数据存储方法的SpMV实现方法，其特征在于包括如下 步骤1)遍历稀疏矩阵中的每一个子稀疏矩阵，计算每个子稀疏矩阵的稀疏矩阵向量乘y= A1^x, A1为第1个子稀疏矩阵，1为正整数，计算方法如下a)根据该子稀疏矩阵中m条非零元对角线的编号和当前计算行的行号j+1，确定计算 yJ+1所需要的χ的元素索引，m为正整数；b)从val数组中读取第j+Ι行的m个元素值，与χ的元素索引对应的m个χ值分别相 乘后累加，得到ym ；c)继续计算j+2行的稀疏矩阵向量乘yj+2；2)合并所有子稀疏矩阵计算得到的稀疏矩阵向量乘，得到稀疏矩阵A的稀疏矩阵向量乘。
8.根据权利要求7所述的SpMV实现方法，其特征在于采用如下方法计算从val数组中读取第j+2行的m个元素值，与χ的元素索引加1对应的m个χ值分别 相乘后累加，计算得到yj+2。
9.根据权利要求7所述的SpMV实现方法，其特征在于若所述子稀疏矩阵为可压缩子稀 疏矩阵，采用下述方法计算读取val数组存储的该可压缩子稀疏矩阵的m个元素值，与χ的元素索引加1对应的 m个χ值分别相乘后累加。
10.根据权利要求7所述的SpMV实现方法，其特征在于所述矩阵第j+Ι行对应的χ的 元素索引由非零元对角线的编号加上j+Ι获得。
全文摘要
本发明公开一种稀疏矩阵的对角线数据存储方法及基于该方法的SpMV实现方法，存储方法为1)按行扫描稀疏矩阵A，以对角线编号表示非零元对角线的位置；2)以非零元对角线与矩阵A侧边的交点作水平线将矩阵A切分为多个子稀疏矩阵；3)按行顺序存储每个子稀疏矩阵中的非零元对角线上的元素到val数组。SpMV实现方法为1)遍历稀疏矩阵，计算每个子稀疏矩阵的稀疏矩阵向量乘y＝A1*x；2)合并所有子稀疏矩阵向量乘。本发明的数据存储方法不需要存储非零元的列索引，减小了存储空间需求和访存开销；对角线和x数组索引数组占用较小存储空间，降低了访存复杂度，计算所需的数据均为连续访问，使得编译器和硬件可以充分优化。
文档编号G06F17/16GK102141976SQ20111000407
公开日2011年8月3日 申请日期2011年1月10日 优先权日2011年1月10日
发明者刘芳芳, 孙相征, 张云泉, 王婷, 袁良 申请人:中国科学院软件研究所