高效求解稀疏线性方程组的算法

在13.1节中已经提到，在整个数值模型的实现程序中，稀疏线性方程组的求解时间在总模拟时间中所占的比重十分大。需要提及的是，对文中所提及的程序而言，需要求解的稀疏线性方程组的系数矩阵必定是对称正定矩阵，对这种线性方程组，如果矩阵各行的带宽或矩阵的外形不大，则采用第12章混凝土细观力学分析串行程序中的变带宽算法将是十分有效的，如果是从二维问题得到的稀疏线性方程组，还可以采用已有的许多行之有效的排序技术来大幅减小外形，从而减少存储需求与计算量。在本文所述的程序中，针对的是三维问题，所以矩阵的外形很大，相应的线性方程组求解时间也很长。为减少存储需求与计算量，可以采用共轭梯度法（CG法）（Saad，1996）来求解这种稀疏线性方程组。

CG法需要不断进行迭代，直到迭代得到的近似解满足给定精度要求为止。取定一个初始解向量并计算出对应的初始残向量、方向向量之后，利用上一步的残向量与方向向量确定本迭代步的方向向量，从而更新解向量。在该方法中，由于主要的计算只有向量内积、向量范数、向量更新、系数矩阵乘以向量这几类操作，所以只需要存储系数矩阵与少量几个向量，而系数矩阵是稀疏的，完全可以采用稀疏矩阵技术来进行存储。例如，对13.3.4节中所述试件进行模拟时，对系数矩阵每行只需要存储最多81个非零元，其他零元素不需要存储。

同时，在每个迭代步内，向量内积、向量范数与向量更新等操作的计算量很小，与线性方程组的阶数成正比，所以计算量主要表现在系数矩阵与向量的乘法，在采用稀疏矩阵存储技术后，计算量也可以减少到大约只要矩阵非零元个数2倍这么多个浮点操作。由此可知，每一次迭代的计算量相对较小，且与总结点数成正比。

为快速求解稀疏线性方程组，采用CG法时每次迭代所耗费的时间很短，但可能需要很多次迭代才能使得迭代得到的近似解满足给定的精度要求。CG法的收敛速度与系数矩阵的特征值分布有关，分布越集中，收敛越快。为使得特征值分布更为集中，现在较为流行且较为有效的一个技术是采用预条件，即将原线性方程组Ax＝b转化为另外一个与之等价的线性方程组MAx＝Mb，其中M是A的逆矩阵的某种近似，且计算M与任意向量的乘积时计算量必须比较小。(www.xing528.com)

针对一般的稀疏矩阵，已经有许多预条件技术（Saad，1996；Benzi，2002；吴建平等，2004），其中最为有效的是ILUT预条件（Saad，1996）与MRILUT预条件（吴建平等，2005）。但是，ILUT与MRILUT都是针对一般稀疏矩阵的，即使对一个对称矩阵，构造得到的ILUT与MRILUT预条件也可能是不对称的，这样就无法利用预条件CG法来求解对应的稀疏线性方程组。吴建平等（2005）针对对称正定矩阵，构造了一类类似于ILUT的预条件，但实际进行的是不完全Cholesky分解而不是不完全LU分解，得到的预条件具有对称性，同时探讨了当主元很小或为0时的处理方法，称为ICT预条件。在对文中所述的程序进行算法改进时，就采用了这种预条件。

此外，由不完全分解构造得到的预条件能不能有效近似原系数矩阵的逆，要看原系数矩阵的对角占优性如何。系数矩阵的对角占优性越好，不完全分解与真实分解的差异越小，不幸的是对混凝土力学性能模拟时得到的线性方程组，其系数矩阵虽然对称正定，但对角占优性比较差，甚至绝大多数行都没有对角占优性。为改善对角占优性，从而提高不完全分解预条件的有效性，在求解线性方程组之前，利用对角矩阵来适当地增强对角占优性，将原线性方程组Ax＝b转化为D－1AD－1（Dx）＝D－1b，其中D的第k个对角元为A的第k个对角元的平方根。