解大型稀疏方程组的ICCG方法及其计算机实现

　　摘要　针对ICCG算法中的关键步骤，提出了快速有效的计算技巧，以利于实际应用。
　　关键词　大型稀疏方程组　不完全Cholesky共轭梯度法　解　实现
　　中国图书资料分类法分类号　O241.6＋P631.34

ICCG METHOD FOR SOLVING LARGE SPARSE EQUATIONS AND ITS COMPUTER PROGRAMMING

Wu Xiaoping　Xu Guoming　Li Shican
(University of Science and Technology of China)

　　Abstract　In this paper,an efficient method for the computation of critical step in ICCG iteration is put forward,and some necessary programs (in Fortran) are developed for sake of the practical application.
　　Keywords　Large sparse equations;Incomplete Cholesky conjugate gradient(ICCG);resolution;implement

1　引言

　　随着数值计算及计算机技术的发展，三维复杂结构的地球物理反演解释渐渐成为可能。不过研究中发现，三维正演计算的低效率是其主要的障碍。因为正演是反演的基础，且反演中的大部分计算量在正演计算上；而地球物理三维正演一般无解析解可寻，只能利用诸如有限差分、有限元等方法求得数值解，最后都是归结为解大型稀疏线性方程组问题：

Ax=b。

(1)

　　这里A一般为大型稀疏带状对称正定矩阵，正演的关键正是解此大型方程组。近年来利用不完全Cholesky共轭梯度(ICCG)迭代方法解大型稀疏线性系统取得较好效果^［1，2］，应用于地球物理三维正演计算极具潜力^［3］。吴小平和徐果明^{［4，5，6］}将ICCG方法应用于点源二、三维地电场的有限差分计算，同时引入系数矩阵按行索引一维稀疏存储模式，无论内存需求和计算速度方面均获较大成功。上述诸文对ICCG算法理论与方法的探讨详多。然而由于一维稀疏存储模式的引入使得编程颇为复杂，所以实际的计算机实现还是具有相当的难度并且也鲜有文献做这方面的讨论。正是基于这种情况，本文针对ICCG算法中的关键步骤，提出了快速有效的计算技巧，以期对ICCG算法的实现有更大帮助。

2　ICCG方法

　　解(1)式线性系统的ICCG方法在文献［4］、［5］、［6］中均有详细描述，这里直接给出其迭代过程如下：
　　令r₀=b－Ax₀,p₀=(CC^T)^－1r₀，则

(2)

i=0,1,2,…。

　　其中C为下三角矩阵，是Meijerink and Van^［1］提出的不完全Cholesky分解因子。分解如下：

A≈CC^T。

(3)

　　除对角元素不一样外，矩阵C与A下三角阵的非对角元素完全一样，当然它们也有完全一样的稀疏性。所以矩阵C的计算非常简单快速，且无需开辟另外的空间来存储矩阵C，详见文献［4］、［5］、［6］。

3　ICCG方法的计算机实现

3.1　大型稀疏矩阵按行索引的一维稀疏存储模式
　　按行索引的一维稀疏存储模式只要求存储系数矩阵A的非零元素，而无需存储带宽内的大量零元素，大大节省了存储空间。
　　为描述矩阵A_N×N，要求建立两个一维数组sa和ija。sa是一实型数组，ija是一整型数组。存储原则如下：
　　a.　sa的前N个元素按顺序存储矩阵A的对角元素，包括零对角元素。这并不增加多少存储量，且多数实际问题中系数矩阵A的对角元素是非零的。
　　b.　ija的第一个元素总等于N＋2。读出ija的第一个元素则可确定N的值。
　　c.　从元素ija(2)到ija(N＋1)的值由各行中不为零的非对角元素的个数决定；要求满足ija(i＋1)－ija(i)等于第i行中不为零的非对角元素的个数，其中i=1,2,…，N。
　　d.　ija(N＋1)是矩阵A最后一行的最后一个不为零的非对角元素在sa中的顺序值加1。因此读取ija(N＋1)可确定矩阵A中非零元素的个数或数组sa和ija的维长。sa(N＋1)是没用的，可设置成任意值。
　　e.　在sa大于等于N＋2的元素中按行顺序存储矩阵A的非对角元素值，在同一行中则按从小到大的列顺序排列。
　　f.　在ija大于等于N＋2的元素中存储对应sa中元素的列号。
　　这一规则初看象是任意的，实际上是非常巧妙的存储方法，只要求2倍非零元素个数的存储量。例如，某一矩阵A为：

按以上存储原则，矩阵A由两个一维数组sa和ija描述如表1所示。

表1　描述矩阵A的一维数组sa和ija表