CN102141976B - 稀疏矩阵的对角线数据存储方法及基于该方法的SpMV实现方法 - Google Patents

Abstract

本发明公开一种稀疏矩阵的对角线数据存储方法及基于该方法的SpMV实现方法，存储方法为1)按行扫描稀疏矩阵A，以对角线编号表示非零元对角线的位置；2)以非零元对角线与矩阵A侧边的交点作水平线将矩阵A切分为多个子稀疏矩阵；3)按行顺序存储每个子稀疏矩阵中的非零元对角线上的元素到val数组。SpMV实现方法为：1)遍历稀疏矩阵，计算每个子稀疏矩阵的稀疏矩阵向量乘y＝A₁*x；2)合并所有子稀疏矩阵向量乘。本发明的数据存储方法不需要存储非零元的列索引，减小了存储空间需求和访存开销；对角线和x数组索引数组占用较小存储空间，降低了访存复杂度，计算所需的数据均为连续访问，使得编译器和硬件可以充分优化。

Description

稀疏矩阵的对角线数据存储方法及基于该方法的SpMV实现方法

技术领域

本发明涉及一种针对稀疏矩阵的数据存储方法及基于该方法的SpMV实现方法，尤其涉及针对稀疏矩阵的对角线存储方法和基于该方法的SpMV实现方法。

背景技术

稀疏矩阵向量乘(SpMV)y＝A*x是一个重要的计算内核，广泛应用在信号处理、图像处理、迭代求解算法等科学计算和实际应用中.但在基于cache存储层次的计算平台上，与稠密矩阵向量乘相比较，稀疏矩阵向量乘的性能很差，主要是因为cache存储层次的复杂性和稀疏矩阵非零元素分布的不规则性。浮点计算操作和存储访问操作的比率非常低，尤其是向量x的间接访问和访问的不可重用性，使得传统的压缩行存储方法得到的稀疏矩阵向量乘性能很差，往往低于机器浮点运算峰值的10％。稀疏矩阵向量乘作为内核如果能够提高运行速度，整个工程计算的运行效率将会得到大大地改善，时间也会大大地减少，在实际应用中有着十分重要的作用.例如，在天体物理项目中，因为要多次调用SpMV代码，稀疏矩阵向量乘的优化能明显提高项目的整体速度。

实际应用中大部分稀疏矩阵具有对角线模式，为此，现有的研究已提出了DIA以及RSDIA等系数矩阵存储方法，但在应用过程中都有不足之处。

稀疏矩阵是指矩阵A的元素大部分是零，而非零元素所占比例非常小，往往小于元素总数的1％.通过只存储和操作这些非零元，达到内存和计算代价上的节省.稀疏矩阵存储时，除了存储非零元外，还要记录非零元在矩阵中所占的位置.稀疏矩阵的一般存储方法是CSR(compressed spar se row)方法。CSR方法需要存储稀疏矩阵A的每个非零元素的值，非零元所在的列和每行第1个非零元的索引，即需要3个数组(其中矩阵A是m×n矩阵，有nz个非零元)，如下所示：

·val[nz]，记录每个非零元的值；

·col[nz]，记录每个非零元所在的列；

·ptr[m+1]，记录每行的第1个非零元在数组val[nz]col[nz]中的索引，其中ptr[m]＝nz.

稀疏矩阵向量乘的形式为y＝Ax，其中A为稀疏矩阵，x，y为稠密向量.稀疏矩阵向量乘的一般实现算法运行效率很低，往往低于机器浮点运算峰值的10％，其中一个原因是，数据引用的时间局部性和空间局部性比较差，尤其是向量x的间接访问和访问的不可重用性；另一个原因就是算法中的浮点计算与存储访问的比率太低，浮点操作和load指令混杂导致把大量的时间花费在数据的存储访问中.提高稀疏矩阵向量乘的计算效率，优化稀疏矩阵向量乘算法是很困难的，这是由存储层次的复杂性和存储访问的不规则性决定的.存储系统分为寄存器、一级cache、二级cache、内存，从前到后存储空间越来越大，访问数据所需要的时间也越来越长.所以应尽量使得数据在高层存储里面得到重复使用，减少再次从低层存储访问的次数和开销，从而提高运行效率.

由于CSR中每次计算乘加都需要前一次计算的结果yt，且浮点运算是多级流水，产生了stall，为此，出现了L-CSR结构，应用unroll-and-jam(loop jamming)方法，将拥有非零元个数一样的行存放在一起为一组，并保存一个长度索引数组指向所有不同长度的组，使用2级unroll-and-jam并对长度为2和大于2的情况分别行编写了代码。要求每行的非零元的长度在一个小范围内，可能破坏CSR算法中访问向量x的局部性，尤其是当矩阵为对角线模式时，并且，写入向量y的步长可能为1。

而CSR-DU算法减少了索引大小和访问带宽需求，只存储同一行中每个元素与前一个元素列索引的差值，并将row_ptr和col数组合并，在压缩的索引中标识新行的开始，在计算SpMV是解压缩索引，用减少的带宽需求来覆盖解压缩的时间达到提升性能的目的。但是要求矩阵不存在大量非零元较少的行。CSR-VI减少了val数组大小，只保存不同的val值，并为原始CSR中的所有非零元建立一个间接访问val数组的索引，只在不同非零元值数目较小时有用，它也减少了利用CPU进行流水计算的可能性。

对称存储方法针对对称矩阵，寄存器分块(大小为r*c)和向量分块(矩阵乘以多个向量，对向量进行分块，大小为v，以利用矩阵A的局部性)进行优化。提出了一个类似Sparsity中的模型分析r，c和v：给定机器参数，对称矩阵A和向量个数k，首先进行off-line benchmark，对稠密矩阵用稀疏存储方法，所有的1≤r，c≤8和1≤v≤10的所有组合都进行测试，记录峰值P(r，c，v)；然后进行run-time search，选择A中部分行，对所有的1≤r，c≤8计算非零元填充率f(r，c)；最后用启发式算法，选择P(r，c，v)/f(r，c)最大的r，c和v值。

改进的CSR还有DCSR和RPCSR。DCSR使用6个command code进行模式压缩，将3个连续的模式压缩为四字节，称为一组，非常适用于非零元多为连续的情况。并将row_ptr和col数组合并，在压缩的索引中标识新行的开始，在计算SpMV是解压缩索引，用减少的带宽需求来覆盖解压缩的时间达到提升性能的目的。RPCSR对四字节的组进一步压缩，选取出现多次的模式编写专门代码优化。

Belgin等提出了基于Pattern的压缩：选择分块大小r，c，共有2^(r*c)中块的非零元模式，忽略少于3个非零元的块和出现次数少于1000的块模式，这些元素用CSR存储，选择r，c使得尽量减少用CSR存储的元素，为每一个块模式编写代码，并利用SSE2/3指令对每个r*c大小的块中的2*2大小的块编写SSE代码。

Mellor-Crummey等使用寄存器分块，cache分块和多向量乘优化SpMV。BCSR减少索引存储空间和访问带宽需求，减少向量x间接访问次数。增加了访问块内零元的开销，增加了带宽需求和计算量。

Geus等对对称稀疏矩阵向量乘使用软件流水(显式预取数据，即先取下一次计算的数据，然后进行本次计算，增加指令并行度)，寄存器分块和矩阵重排三种优化方法。

Buluc等提出了CSB存储结构以及并行算法。分块大小为b，共有n^2/b^2个块，存储每个块在val数组中的起始位置的索引，在每个块内用三元组存储，这样可以减少行和列索引所需的大小。

Toledo等认为CSR算法有四个性能瓶颈：对val和col是顺序访问，可能会引起一些Cache缺失，或者是L1 Cache，或者是更低层次的Cache缺失；对x的访问缺少时空局部性，导致大量Cache缺失；每次浮点乘加操作需要加载三个浮点数，所以CPU的load-store单元的性能十分关键；利用col数组访问x需要间接寻址。

Pinar等基于Toledo的工作，引入变长的分块，BCRS，1*m分块，用大小为分块数目的Nzptr数组存储每个1*m分块的在val数组中的起始位置。Row_ptr指向每行第一个块在Nzptr中的位置，算法使用三重循环。首先是对行循环，然后对每行中的块循环，最后是块内循环。但是BCRS效果并不好。文章进一步提出了矩阵重排的方法，增加矩阵每一行中的稠密块大小。将矩阵重排问题化为TSP问题。用TSP的启发式算法解决。稀疏矩阵中两列在相同的行有非零元的数目做为这两列的权重，可以通过计算W＝A*(A的转置)，A为稀疏矩阵，1代表非零元，则W矩阵为权重。最后，重排后的矩阵的分块数目＝矩阵非零元数目-TSP路径的权重(即所有相邻列的权重)。TSP两大类启发式：构建和改进。构建直接构建一个解，然后改进从一个初始解开始，寻找邻居节点进行改进。

Rice等提出了ELLPACK/ITPACK(ELL)格式，两个相同大小的数组，都为m*mnz，mnz是每行拥有最大非零元的个数，少于mnz个非零元的行用0来填充，另外一个数组存储每行所有非零元的列值，对于填充的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中设置对角线，比如n×n的矩阵共有2n-1条对角线，n为正整数，其编号如图1所示，上三角和下三角中对角线编号分别用D_i和D_-j表示，i≥0，j≥0，i和j为整数，i为上三角对角线中第一个元素的列号减1，j为下三角对角线中第一个元素的行号减1。则主对角线编号i＝j＝0，为D₀，

D_i对角线的大小为n-|i|，D_-j对角线的大小为n-|j|，图1的右侧说明了计算每条对角线需要的相应的向量x和y的范围。

扫描子稀疏矩阵时，若相邻行同一对角线上的元素相同，则具有相同元素的相邻行被进一步切分为一可压缩的子稀疏矩阵，在val数组中仅需存储该可压缩子稀疏矩阵一行的非零元对角线的元素。此时稀疏矩阵的数目为进一步切分为可压缩的稀疏矩阵后的矩阵个数。

所述val数组的存储方式为按行顺序存储非零元对角线上的所有元素，包括非零元和零元。

所述稀疏矩阵A为对角线模式的稀疏矩阵，即非零元主要集中在若干条对角线上。

在此过程中，还存储若干个数组及数据：对角线数目；diag_index[diag_cnt]，所有对角线的编号；row_block_cnt，划分的子稀疏矩阵的个数；row_cnt[row_block_cnt]，每个子稀疏矩阵内的行数；compressed[row_block_cnt]，每个子稀疏矩阵是否为可压缩稀疏矩阵，row_region_diag[row_block_cnt]，每个子稀疏矩阵所用到的对角线范围。由于对角线模式的diag_cnt与row_block_cnt个数相比非零元个数非常小，所以这些数组的影响可以忽略不计。

基于上述的稀疏矩阵数据存储方法，下面提供一种SpMV实现方法

1)遍历稀疏矩阵中的每一个子稀疏矩阵，计算每个子稀疏矩阵的稀疏矩阵向量乘y＝A₁*x，A₁为第l个子稀疏矩阵，l为正整数。

第l个子稀疏矩阵的SpMV实现方法为：

a)根据该子稀疏矩阵中m条非零元对角线的编号和当前计算行的行号j+1，确定第j+1行x的元素索引；m为正整数。

b)从val数组中读取第j+1行的m个元素值，与x的元素索引对应的m个x值分别相乘后累加，得到y_j+1；

c)继续计算y_j+2；

2)合并所有子稀疏矩阵计算得到的稀疏矩阵向量乘，得到稀疏矩阵A的稀疏矩阵向量乘。

从val数组中读取第j+2行的m个元素值，与x的元素索引+1对应的m个x值分别相乘后累加，计算得到y_j+2；根据row_region_diag[l]和diag_index数组，得到第l个子稀疏矩阵的非零元对角线的个数m及位置，根据这m个对角线位置和当前计算的y的行号，可以确定该子稀疏矩阵内第j+1行每个非零元所对应的x元素索引，而剩余的其他行的x位置，由于对角线是连续存储的，所以可通过前一行相同对角线上非零元所对应的x索引加1得到。

所述子稀疏矩阵为可压缩子稀疏矩阵时，j+1行的元素值与第j+2行的元素值相同，val数组中仅存储第j+1行的元素值，计算y_j+2时读取val数组存储的该可压缩子稀疏矩阵的m个元素值，与x的元素索引加1对应的m个x值分别相乘后累加。由于可压缩子稀疏矩阵内的所有行的非零元均相同，因此只需要读取一次val数组。根据row_region_diag[l]和diag_index数组，得到第l个子稀疏矩阵的非零元对角线的个数m及位置，然后从val数组中读取下面m个元素，并用这m个元素计算该矩阵内所有行。

所述第j+1行x的元素索引由非零元对角线的编号加上j+1获得。

该稀疏矩阵的SpMV实现方法的代码可见表1：

表1基于对角线数据存储方法的SpMV实现代码表

本发明的有益效果：

本发明的数据存储方法不需要存储非零元的列索引，减小了存储空间需求和访存开销；当对角线数目较少时，本发明的SpMV实现方法所需要子矩阵，对角线和x数组索引数组占用较小存储空间，不需要col数组的访问，降低了访存复杂度，计算所需的数据均为连续访问，使得编译器和硬件可以充分优化。例如使用预取指令，而CSR算法中很难对x数组进行预取。矩阵A的访存复杂度为ne。

与现有的CSR算法相比，在不考虑本发明数据存储方法中对角线零元的情况下，CSR和本发明的SpMV实现方法计算y和访问x，A的顺序是一致的，只是通过对角线存储方式，减少了col数组的访问开销。

计算复杂性(CC：computational complexity)分析：CSR计算复杂性为O(m+nz)，DDD-SPLIT计算复杂性为O(ne)，可见当对角线数多为稠密时，两种算法的计算复杂性相当。

访存复杂度(MAC：memory access complexity)分析：CSR：将可以重用x的对角线合并为一组，设由于间隔较大而不能重用x的对角线组个数为n_dig_group。访问val和col数组的复杂度为2×nz；向量y在寄存器中重用，所以访存复杂度为n；虽然每条对角线需要访问的向量x的范围不一样，但是可以统一记为，O(n)，所以向量x的访存复杂度为O(n×n_dig_group)。总访存复杂度为2×nz+n+O(n×n_dig_group)。DDD-SPLIT：矩阵A的访存复杂度为ne，向量y在寄存器中重用，所以访存复杂度为n；类似于CSR算法中向量x的分析，向量x的访存复杂度均为O(n×n_dig_group)。总访存复杂度为ne+n+O(n×n_dig_group)。

因此根据复杂度的分析，本发明的稀疏矩阵对角线存储方法，向量y在寄存器中可重用，对于某些具有对角线模式的稀疏矩阵，在相同对角线上存在连续的相同非零元，而由于这些相同非零元属于同一对角线，压缩这些相同非零元，并进一步划分子稀疏矩阵，使得每个子稀疏矩阵，算法只需要存储第一行非零元，而节省了存储行块内剩余的行的元素的空间，降低了带宽需求，并且，由于非零元是连续的，不用单独开辟存储索引的数组。

本发明的方法对对角线稠密的稀疏矩阵具有更好的效果。对角线稠密是指非零元对角线上的非零元元素较多，即超过阈值的对角线上元素都是非零元的稀疏矩阵。阈值在不同的测试平台上有所不同。

附图说明

图1是本发明的存储方法的存储结构图

图2是本发明的存储方法中对稀疏矩阵的分割方法示意图

图3是实施例中的实验矩阵的填充率结果示意图

图4是实施例中测试平台1上的加速比结果示意图

图5是实施例中测试平台2上的加速比结果示意图

图6是实施例中可压缩子稀疏矩阵在测试平台1上的加速比结果示意图

具体实施方式

针对对角线模式的稀疏矩阵，采用如下步骤进行数据的存储：

1)按行扫描稀疏矩阵A，确定具有非零元的对角线的位置，以对角线编号来表示，主对角线编号为D₀，上三角和下三角中对角线依次编号，分别用D_i和D_-j表示，i为上三角对角线中第一个元素的列号-1，j为下三角对角线中第一个元素的行号-1。如图1所示。存储非零元对角线的编号到diag_index[diag_cnt]数组；计算非零元对角线的数目并存储；

2)参见图2，图2(a)为一对角线稀疏矩阵示例，以非零元对角线与矩阵A侧边的交点作水平线水平切分矩阵A，得到多个子稀疏矩阵，如图2(b)；

3)扫描子稀疏矩阵，按行顺序存储每个子稀疏矩阵中的非零元对角线上的元素到val数组，计算子稀疏矩阵的个数和每个子稀疏矩阵的行数并存储；计算每个稀疏矩阵中对角线的范围。即每个子稀疏矩阵块内存在非零元的对角线的范围并存储。

在扫描子稀疏矩阵时，判断子稀疏矩阵中相邻行中同一对角线上的元素是否相同，若相同则将该子稀疏矩阵进一步切分为可压缩子稀疏矩阵，由于图2中的矩阵的第1、2行对角线元素相同，第3、4行对角线元素相同，可以将1、2行划分为一个可压缩子稀疏矩阵，第3、4行划分为一个可压缩子稀疏矩阵，每个可压缩子稀疏矩阵仅存储第一行的元素，如图2(c)所示。val数组中也仅存储第一行的元素。val数组的存储方式为按行顺序存储非零元对角线上的元素，包括非零元和零元。

根据上述的对角线存储方法，对稀疏矩阵A进行SpMV实现：

1)遍历稀疏矩阵中的每一个子稀疏矩阵，计算每个子稀疏矩阵的稀疏矩阵向量乘y＝A₁*x，A₁为第l个子稀疏矩阵，

第l个子稀疏矩阵的SpMV实现方法为：

a)根据该子稀疏矩阵中m条非零元对角线的编号和当前计算行的行号j+1，确定第j+1行x的元素索引；第j+1行x的元素索引可以由非零元对角线的编号加上计算行的行号j+1获得。

b)从val数组中读取第j+1行的m个元素值，与x的元素索引对应的m个x值分别相乘后累加，得到y_j+1；

c)若子稀疏矩阵不是可压缩子稀疏矩阵，则从val数组中读取第j+2行的m个元素值，与x的元素索引+1对应的m个x值分别相乘后累加，计算得到y_j+2；

若子稀疏矩阵为可压缩子稀疏矩阵，读取第j+1行的m个元素值，与x的元素索引+1对应的m个x值分别相乘后累加，计算得到y_j+2。

2)将所有子稀疏矩阵的稀疏矩阵向量乘顺序合并，得到稀疏矩阵A的稀疏矩阵向量乘。

实施例

采用以上介绍的技术，我们用两个测试平台进行验证，测试平台1为AMD Opteron 8378，2.4GHz。测试平台2为Intel Xeon X5472 3.00G。选用的实验矩阵的相关数据见表2，一共选取23个矩阵，具有较强代表性。这些实验数据基于天文台项目，天文台既有稠密对角线，也有可用CSD存储的稀疏对角线。表2中的前10个稀疏矩阵存在非零元的对角线数目较少，使用本发明的对角线数据存储方法(简称DDD-SPLIT方法)进行存储；后面13个稀疏矩阵存在较多对角线，并且大部分对角线上的非零元个数较少，如果都用DDD-SPLIT方法储存，则浪费较多存储空间，对于这些矩阵，我们采用和DIAG一样的方法进行处理，即对于较稀疏的对角线(存在较少非零元)采用DIAG方法储存，对于其他对角线，采用DDD-SPLIT方法储存。表3则给出了后面13个稀疏矩阵在两个平台上获得最优性能时用DDD-SPLIT方法存储的稀疏对角线数目及其他相关信息，左半部分的数据是在Intel平台上的测试结果，右半部分的数据是在AMD平台上的测试结果(由于测试平台的硬件不同，其测试结果有所区别)。从表3可以看出，大部分的对角线采用的是DDD-SPLIT储存方法。而对于4，5，6，7，9这几个稀疏矩阵，可以进一步划分为可压缩子稀疏矩阵，表4显示了这几个矩阵的压缩率。

表2实验矩阵相关信息表

表3.实验矩阵11-23在不同测试平台上采用DDD-SPLIT方法存储的相关数值表

表4实验矩阵的压缩率

对上述矩阵调用不同的SpMV实现方法，并重复运行1000次，然后平均其时间，进行最终各个方法的性能比较。

图3显示了23个实验矩阵的填充率，即需要填充的非空元个数与需要存储的所有元素空间个数之比，其值fr介于0和1之间，填充率fr越小，说明稀疏矩阵的非零元都集中在几条对角线上，并且对角线上的零元较少。对于一些填充率非常小的矩阵，在图3中用1-fr示意。图3左边的图显示的是矩阵1-10的填充率，右边的图显示的是在不同的测试平台上，矩阵11-23的填充率。

基于不同的存储方法对实验矩阵1-23进行SpMV实现，其加速比参见图4、5、6.加速比的基础是CSR方法，即调用DIAG或者是DDD-SPLIT存储方法运行的时间与CSR存储方法运行时间的比值。由于矩阵11-23的填充率非常大，如果填充非空元则需要巨大存储空间，可以确定如果完全采用DIAG存储方法，则对CSR存储方法的加速比会非常低，因此没有进行实验。

该专利技术包括：基于非零元压缩的对角线存储及其相应SpMV算法。设计了改进的直接对角线存储(DDD-SPLIT)，并基于新的方法提出相应的SpMV算法。实验结果表明，对于具有对角线模式的矩阵，基于本发明存储方法的SpMV算法明显减低了存储空间需求和计算时间。

Claims (6)

1.一种稀疏矩阵的对角线数据存储方法，包括下列步骤：

1）按行扫描稀疏矩阵，对稀疏矩阵中存在非零元的对角线，首先填充零元为非空元，设置非空元的值为0，再将每行内的非空元按照列值进行顺序存放；确定具有非零元的对角线的位置，以对角线编号表示；

2）以非零元对角线与矩阵侧边的交点作水平线将矩阵水平切分为多个子稀疏矩阵；

3）扫描子稀疏矩阵，按行顺序存储每个子稀疏矩阵中的非零元对角线上的元素到数组；计算并存储每个子稀疏矩阵中非零元对角线的编号和数目。

2.根据权利要求1所述的稀疏矩阵的对角线数据存储方法，其特征在于稀疏矩阵的上三角和下三角中对角线编号分别用D_i和D_-j表示，i≥0，j≥0，i和j为整数，i为上三角对角线中第一个元素的列号减1，j为下三角对角线中第一个元素的行号减1。

3.根据权利要求1所述的稀疏矩阵的对角线数据存储方法，其特征在于扫描子稀疏矩阵时，若相邻行同一对角线上的元素相同，则具有相同元素的相邻行被切分为一可压缩子稀疏矩阵，在数组中仅存储该可压缩子稀疏矩阵中第一行的非零元素。

4.根据权利要求1或3所述的稀疏矩阵的对角线数据存储方法，其特征在于所述数组的存储方式为按行顺序存储非零元对角线上的所有元素，包括非零元和零元。

5.根据权利要求1所述的稀疏矩阵的对角线数据存储方法，其特征在于所述稀疏矩阵为n×n的稀疏矩阵，n为正整数。

6.根据权利要求1或5所述的稀疏矩阵的对角线数据存储方法，其特征在于所述稀疏矩阵为对角线稠密的稀疏矩阵，所述对角线稠密稀疏矩阵是指非零元对角线上的非零元元素超过阈值的对角线上元素的稀疏矩阵。

Images