稀疏矩阵的CSR格式存储

稀疏矩阵的具有多种存储格式，本文仅对压缩稀疏行存储（CSR）格式进行介绍。

1 压缩稀疏行存储（CSR）格式

如上图1所示，在存储稀疏矩阵的时候，需要有3个数组来表示所有信息。一个为value数组，长度为NNZ，用于存储非零元素；一个为colIndex数组，长度为NNZ，用于存储非零元素所对应的列索引信息，一个为rowIndex数组，长度为N+1，用于存储非零元素所对应的行索引信息。需要注意的是rowIndex数组还包含其他有用的信息：rowIndex(i)表示第i行的第一个元素在value中的索引，所以有rowIndex(1)=1和rowIndex(N+1)-rowIndex(1)=NNZ，此外m(i)=rowIndex(i+1)-rowIndex(i)表示第i行非零元素的数量。

在有限元中，常采用CSR格式来存储总体矩阵（总体刚度、质量和阻尼矩阵等），原因有两个：一方面可以节省内存开销，有限元的总体矩阵大多为稀疏矩阵且稀疏比较大，以某一具体案例为参考，求解自由度N=669330，NNZ=14049513，稀疏比为0.00612266%。考虑对称性，按照稠密矩阵以double类型存储，占用内存大小为669330×(669330+1)/2×8byte=1.63TB。按照稀疏矩阵CSR以double类型存储，占用内存大小为14049513×8byte=0.10GB。另一方面可以加快矩阵求解（运算）速度，且第三方数学运算库会提供稀疏矩阵CSR存储格式。

在结构有限元中，如何得到总体稀疏矩阵的全部信息，即得到总体稀疏矩阵的rowIndex数组，colIndex数组以及value数组是一大难点，也是后续求解方程的基础。其中rowIndex和colIndex仅与节点单元之间的连接相关，即网格的拓扑信息，因此其可以事先确定好。value与单元矩阵有关，需要先计算单元矩阵，然后将其添加到value中相应位置。

以下图2所示的简单示例来说明总体稀疏矩阵的形成过程，示例有2个四边形单元构成，总共6个节点，假设每个节点2个自由度，则系统共有12个自由度。图2（a）显示了网格拓扑关系，图2（b）显示了系统总体矩阵的稀疏情况，红色的圆圈表示非零元素，黑色的圆圈假设其也表示非零元素，实际上其是零元素。以节点1和5关系来看，其并非直接连接，但其都属于单元1上面的节点，在通过积分点计算单元矩阵并没有判断其是否非0，直接将其添加到总体矩阵的相应位置（单元矩阵当成稠密矩阵处理了），因此在总体结构中认为其位置为非零元素。此种假设对生成总体稀疏矩阵格式和判断相关节点都带来很大的方便。

图2 总体稀疏矩阵的形成过程

相关节点就是指在网格信息中直接相连接的节点，需要说明的是同一单元内的节点都属于相关节点，自身与自身是也是相关节点。这样就可以通过网格拓扑信息得到每个节点的相关节点列表（不可重复且升序排列），具体做法是：先遍历每个单元，在对每个单元遍历其每个节点，得到每个节点的相关节点列表。

如下图3所示，图3（a）所示的conNodeList数组为每个节点的相关节点列表，需要说明的是conNodeList数组的第一维度为总节点数nNode，第二维度大小事先并不知道，有两种处理方法，一是事先假设一个较大的维度，二是采用数据结构，在C++可采用set关联式容器，每次将相关节点insert即可，set属性满足了不可重复且默认升序排列的特性。图3（b）所示的conDofList数组为每个节点的相关节点自由度列表，需要说明的是conDofList数组的第一维度为总自由度数nDof（也为总系统方程数量），第二维度为每个节点自由度相关联的自由度数量，也就是每个方程相关的方程数量，即总体刚度矩阵中每行的非零元数量，且每个元素的具体数值表明了其在总体矩阵中的列位置。对于对称矩阵，只想要存储上三角情况（列>=行），则只需选择数组的列号大于行号的部分即可。这样，通过对相关节点自由度列表conDofList数组的解析，就可以得到总体稀疏矩阵的rowIndex和colIndex信息。

图3 相关节点及相关自由度列表

最后就是确定总体稀疏矩阵的value信息，需要先计算单元矩阵，然后将其添加到value中相应位置。计算完单元刚度矩阵后，可根据该单元上节点对应的自由度判断其该元素在总体矩阵中的行列信息，然后在通过rowIndex和colIndex数组即可确定其在value中的具体 位置。具体做法为：得到行列信息(iR, iC)后，通过遍历第iR行所有元素，找到colIndex(IP)=iC的IP，则(iR, iC)在总体稀疏矩阵中对应的位置为IP，最后将单元矩阵中的值添加在value(IP)即可。