利用超松弛预处理共轭梯度法求解大型稀疏方程组

研究简报
数 学
利用超松弛预处理共轭梯度法 求解大型稀疏方程组
陈春香 尹洪东
(中国石油大学 (北京 ) 机电工程学院 ,北京 102249)
摘 要 利用有限差分法构造大型稀疏方程组对井地电位成像测量非均质电阻率的三维正演进行研究 。对于线性方程组
A x = b, A 是大型稀疏的带状矩阵 , 解大型稀疏方程组的直接共轭梯度法 , 一般要求巨大的计算机内存来存储系数矩阵 A, 而且
预处理共轭梯度法可以求解该类问题 。
用 . 岩石力学与工程学报 , 2007; 26 ( 1) : 2820 —2826
(下转第 2394 页 )
© 1994-2010 China Academic Journal Electronic Publishing House. All rights reserved.
3 林绍忠 . 用预处理共轭梯度法求解有限元方程组及程序设计 . 河
海大学学报 , 1998; 26 ( 3) : 112 —115
4 J I N X Q. Fast iterative solvers for symmetric Toep litz system s—a sur2 vey and an extension. Journal of Computational and App lied Mathe2 matics, 1996; 66 ( 1 /2) : 315 —321 5 吴小平 ,徐国明 , 李时灿 . 解大型稀疏方程组的 ICCG方法及其计
(W
- T T
元素 , 包括零元素 , 这不会增加太多存储量 , 而且实 际问题中对角元素一般都不为零 。
( 2 ) 数组 b 的第一个元素为 N + 2, 知道数组 b
AW
- 1
) W x =W
- T
b
( 2)
式 ( 2 ) 中 :M 为预处理因子矩阵 , 当 M 与矩阵 A 结构 越近似时 , W
采用预处理共பைடு நூலகம்梯度法
[3 ]
, 首先要引入预处理
© 1994-2010 China Academic Journal Electronic Publishing House. All rights reserved.
10 期
陈春香 ,等 : 利用超松弛预处理共轭梯度法求解大型稀疏方程组
3
2391
- T
AW
- 1
的条件数越小 (条件数越小 , 收
的第一个元素即可知道 N 的值 。
( 3 ) 数组 b中的元素 b ( 2 ) 到 b ( N + 1 ) 的值由
敛速度就越快 , 条件数最小为 1 ) 。 关于预处理因子矩阵 M 的确定 , 应满足两个基 本条件 , 一是矩阵 M 应该与矩阵 A 相接近 ; 二是矩 阵 M 应易于求逆 。这里采用超松弛法 方法如下 :
http://www.cnki.net
2394
科 学 技 术 与 工 程
10 卷
A C la ss of Cartesian A then tica tion Codes from PG ( n , q )
L I Xiu 2li, ZHAO Pi2qing
http://www.cnki.net
2390
科 学 技 术 与 工 程
10 卷
为描述矩阵 AN × N , 首先要建立两个一维数组 a 和 b, a 为实数数组 , b为整型数组 , 存储原则如下 :
( 1 ) 数组 a 的前 N 个元素为矩阵 A 的对角线
因子矩阵 M , 目的是降低系数矩阵的条件数 , 以加快 收敛速度 。 对于方程 A x = b, 令 M =W W , 则方程变为
10 时 ,共轭梯度法迭代就非常慢了 。本文采用超
2
松弛预处理的目的就是降低系数矩阵的条件数 , 从 而减少迭代次数与计算量的目的 。
1 超松弛预处理共轭梯度法及其迭代格式
1. 1 大型稀疏矩阵按行索引的稀疏存储模式
2009 年 12 月 24 日收到
按行索引的稀疏存储模式只要求存储系数矩 阵 A 的非零元素 , 而不用存储带宽内的大量零元 素 ,大大节省了存储空间
计算速度极其慢 。因此引入按行索引的稀疏存储模式及超松弛预处理共轭梯度算法 , 充分利用系数矩阵 A 的稀疏性 , 使得需 要的内存大大减小 , 充分提高运算速度 。这种方法对井地电位成像测量非均质电阻率的三维正演具有一定的实用价值 。 关键词 共轭梯度法 ; 超松弛预处理 ; 行索引的稀疏存储模式 中图法分类号 O151. 21; 文献标志码 A
其中 , D 为一对角矩阵 , D ij为其子元素 ; C 为一上三角 矩阵 ; L 为 A 的严格上三角矩阵 (即不含对角线元 素 ) ;ω为超松弛因子 , 其取值范围为 0 <ω < 2, 默认 值为 1。 超松弛预处理共轭梯度法的迭代格式如下 :
( 1 ) 设 x 的初值为 x0 , 计算 r 0 = b - A x0 , 同时可
[2]
第一作者简介 : 陈春香 ( 1985 —) ,女 ,汉族 ,中国石油大学 (北京 ) 检 测技术与自动化装置硕士生 ,研究方向 : 井地电位成像理论 。
。
© 1994-2010 China Academic Journal Electronic Publishing House. All rights reserved.
( 5 ) 在数组 a 中从 N + 2 个元素起按行顺序存
ω
0,
- 1
, ( i = j) ( i≠ j)
( 3) ( 4) ( 5)
C =D M =
(D + L )
2
1 T C C 2 -ω
储矩阵 A 中的非对角线的非零元素值 , 在同一行中 按从大到小的列顺序排列 。
( 6 ) 在数组 b中从 N + 2 个元素开始存储数组 a 中对应元素所在的列号 。
图 1 均匀半空间模型及网格剖分
3
3 结论
探讨 。
参 考 文 献
1 孟 颖 . 基于共轭梯度的工业 CT 图像重建 . 科学技术与工程 , 2008; 8 ( 8) : 2085 —2088 2 张永杰 , 孙 秦 . 大型稀疏线性方程组的全稀疏存储策略 . 陕西
理工学院学报 , 2005; 21 ( 4) : 67 —68
分析检验 。图 1 给出了大小为 20 × 20 × 22 m 的均 匀半空间模型 , 剖分单元大小为 1 × 1 × 1 m , 电阻 率为 100 Ω ・m ,供电电流为 1 A ,井深为 10 m ,用超 松弛预处理共轭梯度法求得的数值解电位值和解 析解之间的相对误差曲线如图 2 (横坐标表示 100 个点 , 纵坐 标 表 示 数 值 解 与 解 析 解 的 相 对 误 差 ) 所示 。 用直接共轭梯度法计算对计算机的内存和速 度要求过高 ,本文采用的行索引存储可以大大减小 内存需求 ,超松弛共轭梯度法解大型方程组节省了 不少时间 ,因此 , 可以相信本文所说的方法在速度 和计算机内存方面的优势 , 在解决该问题时 , 具有 一定的实用价值 。 本文将解大型稀疏方程组的超松弛预处理共 轭梯度法与系数矩阵的稀疏存储模式结合在一起 , 可以实现井地电位成像对地下非均质电阻率的三 维正演计算 , 为基于三维电阻率反演进行了前期
由图 2 可见 , 利用超松弛预处理共轭梯度法求 得的数 值 解 与 解 析 解 之 间 的 相 对 误 差 几 乎 都 在
10%以内 ,基本能满足精度要求 ,这说明利用超松弛
然科学版 ) , 2003; 30 ( 4) : 432 —435
8 郑 超 ,张建海 . 预处理共轭梯度法在岩土工程有限元中的应
A ij D ij =
[8 ]
各行中不为零的非对角线元素的个数决定的 , 满足
b ( i + 1 ) - b ( i) 等于第 i 行中不为零的非对角线元
, 具体选取
素的个数 , i = 1, 2, …, N 。
( 4 ) b (N + 1 ) 是矩阵 A 的最后一行最后一个不
为零的非对角线元素在数组 a 中的顺序值加 1. 因此 读取 b (N + 1 )可确定矩阵 A 中的非零元素的个数或 数组 a 和 b 的维长 , a ( N + 1 ) 是没用的 , 可设为任 意值。
图 2 数值解与解析解的相对误差
算机实现 . 煤田地质与勘探 , 1999; 27 ( 6) : 54 —56
6 CaiM C, J in X Q ,W ei Y M. A generalization of T . Chanπ s p recondi2 tioner . L inear A lgebra and Its App lications, 2005; 40 ( 7) : 11 —18 7 范啸涛 ,季光明 . 预优矩阵及其构造技术 . 成都理工大学学报 (自
井地电位成像测量技术是 20 世纪 80 年代末发 展起来的一种新型电测方法 , 通过开发井的套管向 井下供入大功率的电流 , 在地表测量由套管流入地 层中的“ 漏电流 ” , 通过地下电性介质后 , 在地表形 成的电位分布 , 据此研究地下介质的电阻率分布 。 利用含油饱和度与电阻率的相关关系等方法 , 可以 得到储层中剩余油饱和度和剩余储量的空间分布 , 从而解决注水开发油藏油水分布的动态监测和高 含水开发期剩余油监测等 。 由于油田的注水开发 、 压裂 、 酸化等措施的采 用 ,地下储层电阻率是及其不均匀的 , 因此简单的 电阻率求解公式不能满足要求 , 在正演计算中 , 探 索采用地下空间模型的离散化 , 利用有限差分法构 造出大型方程组进行求解 。 迄今为止 , 求解线性方程组的方法有很多 , 如
这一规则实际上是非常巧妙的存储方法 , 只需 要 2 倍非零元素的存储量 。例如某一矩阵 A 为 :
1 0
A = 0
0 2 6 0 0
5 0 3 0 0
0 0 7 0 9
0 0 0 8 4
( 1)
求解 M h0 = r0 。
( 2 ) 令 p0 = h0 , 迭代次数为 k = 1, 2, … 进行迭代 hk , rk pk , A pk
直接法和迭代法 , 其中迭代法又有共轭梯度法
[1]
,
超松弛迭代法 ,牛顿迭代 , 高斯迭代 , 雅克比迭代等 等 。对于大型方程组问题 , 当系数矩阵的条件数不 是太大 ,右端项又较少时 , 迭代法比直接法的求解 效率更高 。迭代法的种类如此之多 , 到底该选取哪 种迭代法比较适合 , 这主要取决于系数矩阵的条件 数 。共轭梯度法用于解大型线性方程组本身是有 不足之处的 , 当系数矩阵 A 接近单位矩阵时 , 共轭 梯度迭代方法收敛快 ; 当系数矩阵 A 的条件数大于
表 1 中 X 是任意数 。按存储原则 , N 的值是 b ( 1 ) 2 等于 5, 每个数组长度为 b ( b ( 1 ) - 1 ) - 1 = 11。 1. 2 超松弛预处理共轭梯度法及其迭代格式
2 计算结果及其分析
为了验证该方法的可靠性 , 设计了一个模型来
http://www.cnki.net
第 10 卷 第 10 期 2010 年 4 月 1671 21815 (2010) 10 22389 204
科 学 技 术 与 工 程
Science Technology and Engineering
Vol110 No110 Ap r12010
Ζ 2010 Sci1 Tech1 Engng1
表 1 描述矩阵 A 的一维数组 a 和 b 的表
顺序 k 数组 a 数组 b
1 1 7 2 2 8 3 3 8 4 0 10 5 4 11 6 X 12 7 5 3 8 6 2 9 7 4 10 8 5 11 9 4
β k =
pk + 1 = hk + 1 +β k pk
其中 : xk 为每次迭代得到的方程的解 ; rk 为每次迭 代的残差 ; 其余变量均为每次迭代的过程变量 。
0 0
αk =
( 6) ( 7) ( 8) ( 9) ( 10 ) ( 11 )
xk + 1 = xk +α k pk rk + 1 = rk - α k A pk M hk + 1 = rk + 1 hk + 1 , rk + 1 hk , rk
按上述存储原则 , 矩阵 A 由两个一维数组 a 和 b 描 述如表 1 所示 。