GMRES算法的收敛分析与实现
ILU预处理Newton—Krylov方法的潮流计算
ILU预处理Newton—Krylov方法的潮流计算作者：廖小兵王文超李奔来源：《计算技术与自动化》2015年第04期摘要：由于电力系统修正方程组具有高维、稀疏的特点，本文提出将预处理Krylov子空间方法应用于潮流修正方程组的求解，形成预处理NewtonKrylov的潮流计算方法。
结合ILU 预处理方法，比较了最常用的3类NewtonKrylov方法求解潮流方程的计算效果。
通过对IEEE30、IEEE118、IEEE300 和3个Poland大规模电力系统进行潮流计算，结果表明：3类NewtonKrylov方法是电力系统潮流计算的有效方法，呈现出良好的收敛特性和计算效率。
关键词：潮流计算；修正方程组；ILU预处理；NewtonKrylov方法中图分类号：TM744 文献标识码：A1引言潮流计算是电力系统分析中最古老的经典课题之一。
传统的电力系统潮流计算通常以牛顿法为主[1]。
牛顿法是求解非线性代数方程组的有效方法之一，它将非线性方程组的求解转化为线性代数方程组的求解，但由于每次迭代后都需更新雅可比矩阵的元素，导致每次都需求解高维的线性代数方程组。
传统的直接法，如Gauss消去法，LU分解等，计算量和存储量较大，且固有的前推回代过程难以并行[2-3]。
迄今为止，越来越多的国内外研究人员在电力系统潮流计算中采用NewtonKrylov方法求解潮流方程[4-8]。
NewtonKrylov方法是在不精确牛顿法的基础上，结合Krylov子空间迭代法，形成的一类新的求解非线性方程组的数值方法。
这类方法结合了Newton方法的良好收敛特性，以及Krylov子空间方法的存储量少、计算量小、易于并行等优点[9]，非常适合并行求解大规模的非线性方程组问题[10]。
文献[4]首次将Krylov子空间法中的GMRES方法应用于潮流计算中。
文献[5]将此类迭代法与不精确牛顿法相结合（NewtonGMRES），同时采用不同的预处理方法，对两个大规模电力系统进行了对比分析计算。
Krylov子空间方法II
max {q (λi )2 }
2 yi λi
min
max {q (λi ) } y Λy
5/76
= = 即
q ∈Pk , q (0)=1 1≤i≤n q ∈Pk , q (0)=1 1≤i≤n
min
max {q (λi )2 } ϵ0 Aϵ0 max {q (λi )2 } ∥ϵ0 ∥2 A,
⊺
min
11/76
(4.6)
又 Λ 是对角矩阵, 所以 ∥q (Λ)∥2 = max |q (λi )|.
1≤i≤n
设 x(k) 是由 GMRES 方法得到的近似解. 由 GMRES 方法的最优性可 知, x(k) 极小化残量的 2 范数. 因此, ∥b − Ax(k) ∥2 = = ≤
x∈x(0) +Kk (A,r0 ) q ∈Pk , q (0)=1 q ∈Pk , q (0)=1
1+ε δ
10/76
5.2 GMRES 方法的收敛性
正规矩阵情形 设 A 是正规矩阵, 即 A = U ΛU ∗ , 其中 Λ = diag(λ1 , λ2 , . . . , λn ) 的对角线元素 λi ∈ C 为 A 的特征值. 设 x ∈ x(0) + Kk (A, r0 ), 则存在多项式 p(t) ∈ Pk−1 使得 x = x(0) + p(A)r0 . 于是 b − Ax = b − Ax(0) − Ap(A)r0 = (I − Ap(A))r0 ≜ q (A)r0 , 其中 q (t) = 1 − t p(t) ∈ Pk 满足 q (0) = 1. 直接计算可知 ∥b − Ax∥2 = ∥q (A)r0 ∥2 = ∥U q (Λ)U ∗ r0 ∥2 ≤ ∥U ∥2 ∥U ∗ ∥2 ∥q (Λ)∥2 ∥r0 ∥2 = ∥q (Λ)∥2 ∥r0 ∥2 .
预处理子空间迭代法的一些基本概念
CG算法的预处理技术：、为什么要对A进行预处理：其收敛速度依赖于对称正定阵A的特征值分布特征值如何影响收敛性：特征值分布在较小的范围内，从而加速CG的收敛性特征值和特征向量的定义是什么？（见笔记本以及收藏的网页）求解特征值和特征向量的方法：Davidson方法：Davidson 方法是用矩阵( D - θI)- 1( A - θI) 产生子空间，这里D 是A 的对角元所组成的对角矩阵。
θ是由Rayleigh-Ritz 过程所得到的A的近似特征值。
什么是子空间法：Krylov子空间叠代法是用来求解形如Ax=b 的方程，A是一个n*n 的矩阵，当n充分大时，直接计算变得非常困难，而Krylov方法则巧妙地将其变为Kxi+1=Kxi+b-Axi 的迭代形式来求解。
这里的K(来源于作者俄国人Nikolai Krylov姓氏的首字母)是一个构造出来的接近于A的矩阵，而迭代形式的算法的妙处在于，它将复杂问题化简为阶段性的易于计算的子步骤。
如何取正定矩阵Mk为：Span是什么？:设x_(1,)...,x_m∈V ,称它们的线性组合∑_(i=1)^m?〖k_i x_i \|k_i∈K,i=1,2...m〗为向量x_(1,)...,x_m的生成子空间，也称为由x_(1,)...,x_m张成的子空间。
记为L（x_(1,)...,x_m），也可以记为Span（x_(1,)...,x_m）什么是Jacobi迭代法：什么是G_S迭代法：请见PPT《迭代法求解线性方程组》什么是SOR迭代法：什么是收敛速度：什么是可约矩阵与不可约矩阵？：不可约矩阵（irreducible matrix）和可约矩阵（reducible matrix）两个相对的概念。
定义1：对于n 阶方阵A 而言，如果存在一个排列阵P 使得P'AP 为一个分块上三角阵，我们就称矩阵A 是可约的；否则称矩阵A 是不可约的。
定义2：对于n 阶方阵A=(aij) 而言，如果指标集{1，2，...，n} 能够被划分成两个不相交的非空指标集J 和K，使得对任意的j∈J 和任意的k∈K 都有ajk=0, 则称矩阵 A 是可约的；否则称矩阵A 是不可约的。
【国家自然科学基金】_gmres算法_基金支持热词逐年推荐_【万方软件创新助手】_20140802
2009年 序号 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
科研热词 非线性原对偶内点法 近似牛顿方向 边界条件 谱方法 解耦 规整化 组合算法 离散惩罚 电阻抗层析成像 电力系统 波前编码 无功优化 数值模拟 彩色图像恢复 广义极小残差法 广义极小化残余法 图像处理 周期性非定常流动 tikhonov正则化 krylov子空间 gmres
推荐指数 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
2010年 序号 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
科研热词 推荐指数 流固耦合 2 预条件 1 非连续边界元 1 近似牛顿方向 1 边界元法 1 解耦 1 膜结构 1 流激振动 1 数值方法 1 广义极小残值算法(gmres) 1 广义最小化残差 1 广义变分原理 1 工频电场 1 小扰动位势理论 1 变电站 1 动态无功优化 1 动作次数限制 1 n-s方程 1 krylov子空间 1 gmres迭代算法 1 gmres算法 1 gmres 1 bicgstab算法 1
2013年 序号 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
科研热词 层析成像 预处理算法 非结构pebi网格 重启动gmres 迭代线性反演 迭代求解算法 迭代收敛准则 边界元分析 跨孔雷达 超声法 程函方程 稀疏矩阵 气固两相流 核反应堆堆芯计算 有限差分 大型浮体 多重迭代优化 变分节块法 参数化信赖域法 传感器 不完全lu分解预条件 symm积分方程 hafmm gmres迭代算法 gmres算法
求解非对称线性方程组的积多项式预处理GMRES算法
１ ５＋ ．
２０ ． ２０ ．
取 一０ ０ ， 始迭代 向量 Ｘ 一［ ，， ，］ ． ．５初 ｏ ０Ｏ ０Ｏ
记 Ａ 的 谱 口 （ ） 一 （１ ２ ３ ４ Ａ ， ， ， ）一
Ｇ Ｅ （ ）ｅ ｅｖｌ ｓ ・ ＭＲ Ｓ３ ；ｉ ｎ ａｕ ： ｇ ｅ
（． ，． ，． ，． ）ＧＭＲ Ｓ（ ） 法 前 两 次 循 环 ０ ５ １０ １５２０ ． Ｅ ３算
所对 应 的双纽线 如 图 １所示．
１ ２ 理论 基础 ．
引理 １。 设方 程组 Ａ Ｅ ］ ｘ＝６的系数 矩阵 Ａ 是 可 对 角化 的且其谱 分解 为
条件数 ， 从而加 快残量 的收敛速度． 数值试验表明， 新算 法在残量 收敛方面具有 明显 的优势． 关键词 ：多项 式预处理 ； ｙ ｖ子空间；迭代法 ； Ｋｒｌ ｏ ＧＭＲ Ｓ Ｅ
中 图分 类 号 ： ４ ． ０２ １ ６ 文 献 标 识 码 ：Ａ
Ｐｒ ｄｕ ｔｐ ｌ ｎ ｍ ｉ ｌｐｒ ｔ ｅ ｔ ｅ ｏ ｃ ｏ ｙ ｏ ａ ｅ ｒ ａ ｍ ｎｔＧＭ ＲＥＳ ａ ｇ ｒ ｔ ｌ ｏ ｉｈｍ ｏ ｏ ｕ ｉ ｎ ｆ ｆ ｒ ｓ ｌ ｔｏ ｏ
Ａｂ ｔａ ｔ ｓｒ ｃ ：Ｗ ｈ ｎ ｔ ｅ ｎ ｍ ｂ ｒｏ ｏ ｄｔｏ ｆｃ ｅｆ ｉｎ ａｒｘ ｗａ ｏ ａ ｇ ， ｔ ｅ ｐ ｅ ｒ ａｍ ｅ ｔｍ ｅｈ ｄ ｅ ｈ ｕ ｅ ｆｃ ｎ ｉ ｎ ｏ ｏ ｆｉ ｅ ｔｍ ｔｉ ｓｔ ｏ ｌｒ ｅ ｈ ｒ ｔｅ ｔ ｎ ｔ ｏ ｉ ｃ
第３ ８卷 第 ５期
２１ 年 １ ０２ Ｏ月
高等数值分析
高等数值分析历年考题高值试题一，填空题(4*7)1.求矩阵的二范数，条件数。
矩阵A=[0,2;3,0]2.把一个矩阵SVD，然后求个广义逆算个结果3.CG法，求余量的方向二、简答1.三类矩阵哪些能被Householder变换变换出来2.GMRES的性质，每次余量是否变小3.本质上是证明牛顿迭代的超线性收敛三、计算1.strum算法的使用，比书上习题还简单2.Arnoldi过程中断时的若干性质证明3.第八章Galerkin方法的求解，边值条件改成一阶导数。
定理8.4.3的变形。
答疑的时候一定要去听听，李老师在讲牛顿迭代的超线性收敛和Galerkin方法时眉飞色舞，果断就考了。
2012高等数值分析1. 用householder和givens变换做QR分解，由于矩阵特殊，非常简单。
如果拿不准，不妨用GS方法做一遍验证一下，因为不同的QR分解只是符号有差异而已，GS还是比householder简单很多的。
2.1证明rayleigh商最大值等于A的最大特征值，将x拆成各个特征向量之和就容易证。
2.2幂法求一个特征值，一步收敛4.1 去年考题。
注意到Ak*R（k-1）=R（k-1）*A（k-1），那么就类似冒泡算法把Ak移到最右边变成一个A4.2 有点小恶心，去年这个问号问的是Qk第一列是特征向量x1，只需要两边同乘以e1即可，但今年问的是最后一行是特征向量xn，顿时就不会证了，当时打眼一看觉得貌似A也不能说明是可逆的，就没往反幂法这个方向去，但是后来想想其实最多也就是lamdaN=0，其他不为0，也许可以分情况讨论下？后来有同学说假设可逆用了反幂法也没得到什么结论...不知道...还好就是5分，丢了就丢了吧5 考过多年的经典背诵题，默写rayleigh ritz方法和贾哥定理，以及Arnolid和精化Arnolid算法6.1 lagrange差值，三个点而已6.2 最佳平方逼近，解一个2*2的法方程组就完事2003高等数值分析1 证明不动点定理（存在唯一性）2 第三章习题83 共扼剃度法ak的选取，以及正交的证明4 梯形法（迭代，相容，稳定区间）具体为dy/dt + y=0 y(0)=1?5 求正交阵使H*（2/3 1/3 2/3）'=e1求I－2ww' （w的二范数为1）的特征值已知H，问计算Ha的运算量6 摄动原理误差分析7 拉各朗日插值（这里实际考的是代数基本定理的应用）8 忘了2005高等数值考题下面是B卷内容，总共六道题1.用Givens变换QR分解一个3*2的矩阵，并求解一个最小二乘2.证明：对于Minres和Gmres(1)A有k个特征值时，至多k步收敛(2)A有n个不同的特征值，r0由k个属于不同特征值的特征向量构成时，k步收敛（这里没有“至多”）3.A为m*n矩阵，m>n(1)用完全QR分解，不完全QR分解以及SVD表示A+(2)用完全QR分解以及SVD得到min||Ax-b||问题的xls和rls，并加以证明4.(1)证明Arnoldi过程中断时找到准确解(2)证明Arnoldi过程中断时不会发生方法中断(3)当A为正定对称阵时，证明Lanczos方法不会发生方法中断（即W'AV非奇异，讲义上有的）5.A=uv'。
基于CFD-GMRES方法的UAV在线航迹规划算法
第３８卷第６期 计算机应用与软件Ｖｏｌ ３８Ｎｏ．６２０２１年６月 ＣｏｍｐｕｔｅｒＡｐｐｌｉｃａｔｉｏｎｓａｎｄＳｏｆｔｗａｒｅＪｕｎ．２０２１基于Ｃ／ＦＤ ＧＭＲＥＳ方法的ＵＡＶ在线航迹规划算法方 斌１ 冯晓锋１ 李 杰２１（湖南警察学院交通管理系 湖南长沙４１０１３８）２（国防科技大学智能科学学院 湖南长沙４１００７３）收稿日期：２０１９－１０－１２。
湖南省教育厅优秀青年项目（１７Ｂ０８７）；湖南警察学院博士专项项目（２０１６ＺＸ０３）。
方斌，讲师，主研领域：任务规划，机器学习，智能交通。
冯晓锋，副教授。
李杰，讲师。
摘 要 为保证ＵＡＶ在线航迹规划结果满足执行任务的终端约束条件，提出一种基于Ｃ／ＦＤ ＧＭＲＥＳ的ＵＡＶ在线航迹规划滚动优化求解算法。
将同伦连续方法与前向差分增广最小残差法结合，在每个采样时刻计算一次残差向量线性方程，采用前向差分近似替代雅可比矩阵与向量的乘积，再基于广义最小残差法快速求解大规模线性方程。
仿真实验结果表明，该算法能够有效规避空间障碍，并规划出稳定收敛至固定终端状态的运动轨迹，规划时间为百毫秒级，满足在线规划实时性的要求。
关键词 最优控制 在线航迹规划 Ｃ／ＦＤ ＧＭＲＥＳ中图分类号 ＴＰ２７３ 文献标志码 Ａ ＤＯＩ：１０．３９６９／ｊ．ｉｓｓｎ．１０００ ３８６ｘ．２０２１．０６．０３８ＵＡＶＯＮＬＩＮＥＲＯＵＴＥＰＬＡＮＮＩＮＧＡＬＧＯＲＩＴＨＭＢＡＳＥＤＯＮＣ／ＦＤ ＧＭＲＥＳＭＥＴＨＯＤＦａｎｇＢｉｎ１ ＦｅｎｇＸｉａｏｆｅｎｇ１ ＬｉＪｉｅ２１（ＤｅｐａｒｔｍｅｎｔｏｆＴｒａｆｆｉｃＡｄｍｉｎｉｓｔｒａｔｉｏｎａｎｄＥｎｇｉｎｅｅｒｉｎｇ，ＨｕｎａｎＰｏｌｉｃｅＡｃａｄｅｍｙ，Ｃｈａｎｇｓｈａ４１０１３８，Ｈｕｎａｎ，Ｃｈｉｎａ）２（ＣｏｌｌｅｇｅｏｆＡｒｔｉｆｉｃｉａｌＩｎｔｅｌｌｉｇｅｎｃｅｓ，ＮａｔｉｏｎａｌＵｎｉｖｅｒｓｉｔｙｏｆＤｅｆｅｎｓｅＴｅｃｈｎｏｌｏｇｙ，Ｃｈａｎｇｓｈａ４１００７３，Ｈｕｎａｎ，Ｃｈｉｎａ）Ａｂｓｔｒａｃｔ ＩｎｏｒｄｅｒｔｏｅｎｓｕｒｅｔｈａｔｔｈｅＵＡＶｏｎｌｉｎｅｒｏｕｔｅｐｌａｎｎｉｎｇｒｅｓｕｌｔｓｍｅｅｔｔｈｅｔｅｒｍｉｎａｌｃｏｎｓｔｒａｉｎｔｓｏｆｔｈｅｔａｓｋ，ａｒｏｌｌｉｎｇｏｐｔｉｍｉｚａｔｉｏｎａｌｇｏｒｉｔｈｍｂａｓｅｄｏｎＣ／ＦＤ ＧＭＲＥＳｆｏｒＵＡＶｏｎｌｉｎｅｒｏｕｔｅｐｌａｎｎｉｎｇｉｓｄｅｓｉｇｎｅｄ．Ｉｔｃｏｍｂｉｎｅｄｔｈｅｈｏｍｏｔｏｐｙｃｏｎｔｉｎｕａｔｉｏｎｍｅｔｈｏｄｗｉｔｈｄｉｆｆｅｒｅｎｔｉａｌａｐｐｒｏｘｉｍａｔｉｏｎｇｅｎｅｒａｌｉｚｅｄｍｉｎｉｍｕｍｒｅｓｉｄｕａｌｍｅｔｈｏｄ．Ａｔｅａｃｈｓａｍｐｌｉｎｇｔｉｍｅ，ｔｈｅＣ／ＦＤ ＧＭＲＥＳｃａｌｃｕｌａｔｅｄａｒｅｓｉｄｕａｌｖｅｃｔｏｒｌｉｎｅａｒｅｑｕａｔｉｏｎａｎｄｔｈｅｐｒｏｄｕｃｔｏｆＪａｃｏｂｉｍａｔｒｉｘ，ａｎｄｖｅｃｔｏｒｐｒｏｄｕｃｔｗａｓａｐｐｒｏｘｉｍａｔｅｄｔｏｔｈｅｆｏｒｗａｒｄｄｉｆｆｅｒｅｎｃｅ．ＴｈｅＧＭＲＥＳｆａｓｔａｌｇｏｒｉｔｈｍｗａｓｕｓｅｄｔｏｓｏｌｖｅｌａｒｇｅｓｃａｌｅｌｉｎｅａｒｅｑｕａｔｉｏｎｓ．Ｔｈｅｓｉｍｕｌａｔｉｏｎｒｅｓｕｌｔｓｓｈｏｗｔｈａｔｔｈｅａｌｇｏｒｉｔｈｍｃａｎｅｆｆｅｃｔｉｖｅｌｙａｖｏｉｄｂａｒｒｉｅｒｓｏｆｓｐａｃｅ，ａｎｄｔｈｅｐｌａｎｎｅｄｒｏｕｔｅｃａｎｃｏｎｖｅｒｇｅｔｏａｓｔａｂｌｅｔｅｒｍｉｎａｌｓｔａｔｅ，ｍｅｅｔｉｎｇｔｈｅｒｅａｌ ｔｉｍｅｒｅｑｕｉｒｅｍｅｎｔｓ．Ｋｅｙｗｏｒｄｓ Ｏｐｔｉｍａｌｃｏｎｔｒｏｌ Ｏｎｌｉｎｅｒｏｕｔｅｐｌａｎｎｉｎｇ Ｃ／ＦＤ ＧＭＲＥＳ０ 引 言在复杂环境下执行任务的过程中，无人机（ＵＡＶ）必须具备在线航迹规划的能力以规避任务区域中的桥梁、楼宇以及突发威胁等空间障碍，确保飞行安全。
第四讲Krylov子空间方法
如果没有特别注明, 本章内容都是在实数域中讨论.
4.1 投影方法
设 K 是 Rn 的一个子空间, 维数为 dim(K) = m ≪ n. 我们需要在 K 中寻找精确解的一 个 “最佳” 近似. 由于 K 的维数是 m, 为了能够唯一确定这个近似解, 我们需要设置 m 个约 束. 在通常情况下, 我们要求残量满足 m 个正交性条件:
x˜ = x(0) + V y.
· 4-2 ·
由正交性条件 (4.5) 可知 r0 − AV y ⊥ wi, i = 1, 2, . . . , m,
即 W AV y = W r0.
如果 W AV 是非奇异的, 则可解得 y = (W AV )−1W r0. 因此, 近似解 x˜ 可表示为 x˜ = x(0) + V (W AV )−1W r0.
vj+1 = wj /hj+1,j
14: end for
如果计算到第 k (k < m) 步时有 hk+1,k = 0, 则方法会提前终止. 此时 Avk 必定可以由 v1, v2, . . . , vk 线性表出 (这里不考虑浮点运算的舍入误差).
算法 4.1 中的向量 vi 称为 Arnoldi 向量. 需要注意的是, 在该算法中, 我们是用 A 乘以 vj, 然后与之前的 Arnoldi 向量正交化, 而不是计算 Ajr. 事实上, 它们是等价的.
r = b − Ax˜ ⊥ L,
(4.2)
其中 x˜ 是我们所要寻找的近似解, L 是另一个 m 维子空间. 这就是数值计算中常用的 PetrovGalerkin 条件. 如果 L = K, 则称为 Galerkin 条件. 子空间 L 也称为 约束空间 (constraint subspace). 相应地, K 通常称为 搜索空间.
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
２００６年第２期 福建 电脑 １１９ ＧＭＲＥＳ算法的收敛分析与实现 杨利华，艾金花，程昔恩 （景德镇陶瓷学院，江西景德镇３３３０００）
【摘要】一般情况下，ＧＭＲ．ＥＳ算法收敛速度较慢，为了提高ＧＭＲ＿ＥＳ算法的加速收敛速度，使用预处理技术以 加快算法收敛，本文分析了ＧＭＲ．ＥＳ算法的加速收敛现象并实现之。 【关键词】广义最小残量；Ｋｒｙｌｏｖ子空间；Ｒｉｔｚ值；加速收敛；正交投影方法；非对称线性方程组
引言 选择一个具备明显特征值的数值例子并在计算过程中将ｍ 设定为２０．误差条件设为１０ｅ一１０。其中系数矩阵Ａ∈尺 形 如： Ａ 并取ｂ＝［１．１，１．１。１．１，１…１～ｒ，ｘ０＝［０，０。０，０ － 化的，其特征值｛Ａ１，２，３＝１；Ａｉ＝ｉ（ｉ－４，５，．．．，２ｏ）１￣部落在正半 实轴上。根据【文献１的定理５］，容易证明在这种情况下ＧＭＲＥＳ 方法有如下收敛估计： ≤２ ２（ ）（ １ｌｒＩ ｌｌ ２＝ 、／ ＋１ 其中ｋｍａｘ和ｋｍｉｎ分别表示参与迭代过程的最大特征值 和最小特征值。 显然，Ａ的谱中极小特征值是否参与迭代过程将对ＧＭＲＥＳ 收敛率变化起关键作用。为此，考察迭代过程中极小Ｒｉｔｚ值 １ （１（）的收敛情况。当ｋ＝８时，ｋｌ（ｋ）对 １已有很好的逼近，使得 成为一个适当的值。根据『５］，次后的ＧＭＲＥＳ迁昆将和 １，２，３ 不参与迭代的对照ＧＭＲＥＳ过程有着同阶的收敛率。则ｋｌｋ充 ／＿ —— 一１ 分逼近 １的前后残量的收敛因子分别是： 兰 ＝０．６３ 、／Ａ∞，Ａｌ＋１ 皿二 ’、／ ＋１ ＝０－３８ ＧＭＲＥＳ算法在第八次迭代后表现出加速收敛行为．并在第 十七次迭代收敛。其中第ｋ步残量收敛的快慢由量 的 ｌｌ ｌｌ ２ 大小表示。 从上面的数值示例可以看出，一般情况下．ＧＭＲＥＳ算法收 敛速度较慢．为了提高ＧＭＲＥＳ算法的加速收敛速度，使用一种 预处理技术以加快算法收敛。 预处理子的选择 理想的情况是．应该能够找到一个Ｐ．使得Ｐ 的条件数小 于Ａ的，且方程Ｐｚ＝ｒ很容易求解。一些简单的预处理子如下 ： １．对角化预处理子。令Ｐ为一对角矩阵，并且Ｐ的主对角 线与Ａ的主对角线相同；就是说，Ｐ－－－ｄ／ａｇ（ａ￣ａ¨… Ｉ）。如果Ａ 的主对角线上有零元素，令Ｐ－－－ｄ／ａｇ（ｃ。，ｃ１．－－ｃ 。）其中 ＾一ｌ ｃ产（∑矗） 。 ｊ：０ ２．ＳＳＯＲ（对称连续超松弛１预处理子。令ｔＬ－（Ｄ＋ｃ．ｏＥ）Ｄ （Ｄ＋ Ｄ。其中， 为松弛系数，Ｅ为一个三对角矩阵，对角线元素为 零，Ｄ为一对角化矩阵，Ｆ为一上三角矩阵，对角线元素为零，这 三个矩阵满足Ａ＝Ｅ＋Ｄ＋，。注意，如果Ａ是Ｈｅｒｍｉｔ共轭的，则 Ｓｓ０Ｒ预处理子Ｐ也是 ｒｍ／ｔ共轭的。 ３．ＬＵ（ＬＵ分解）预处理子。假设Ａ的ＬＵ分解形式Ａ＝ＬＵ， 其中Ｌ为一单位下三角矩阵。Ｕ为一上三角矩阵。那么，令Ｐ＝ ＬＵ。不完整的ＬＵ预测器为 ，其中Ｌ 和 是将Ａ中的零 元素对应的Ｌ和Ｕ中元素置为零而得到的。因此，Ａ的稀疏结构 可以用来保存Ｌ 和 。 下面举个数值例子来证明预处理技术加速ＧＭＲＥＳ算法收 敛：
对于以上６阶方程组，选取迭代步长ｍ＝２，误差ｅ＝１．０ｅ一６， 在没有选择任何预处理子时，它迭代１９次后。算法收敛。 而在保持其它条件不变的情况下，选择ＳＳＯＲ预处理子时． ＧＭＲＥＳ算法迭代１１次后收敛。 二、算法实现 选择ｃ＋＋作为编程语言（其编程环境是ＷｉｎｄｏＷｓ ＸＰ＋ ＶＣ＋＋６．０）。 首先定义一个矩阵类ＳｐａｒｓｅＭｔｘ．其定义如下： Ｃｌａｓｓ ＳｐａｒｓｅＭｔｘ｛
Ｖｅｃｔｏｒ ｐｒｅｃｏｎｄｉｎｇ（ｃｏｎｓｔ Ｖｅｃｔｏｒ＆．ｉｎｔ ｉ＝０）：，／预处理函数 ｐｕｂｌｉｃ： ｄｏｕｂｌｅ ｅｐｓ ，／输出时的绝对残差 ＳｐａｒｓｅＭｔｘ（ｉｎｔ ｎ，ｉｎｔ ｒｅ，ｄｏｕｂｌｅ’４ｉｎｔ’ｃ．ｉｎｔ Ｄ；／／存储稀疏矩阵
Ｊ； 为了减少计算存储量．采用了适用于任意编号节点生成的
大型稀疏矩阵的按行存储非零元存储格式（压缩稀疏行格式）。 ＧＭＲＥＳ算法编程： ｉｎｔ ＳｐａｒｓｅＭｔｘ：：ＧＭＲＥＳ（Ｖｅｅｔｏｒ＆ｘ．ｃｏｎｓｔ Ｖｅｃｔｏｒ＆ｂ，ｄｏｕｂｌｅ＆ｅｐｓ，ｉｎｔ＆ｉｔｅｒ， ｉｎｔ ｐｃｎ，ｉｎｔｍ）Ｉ ／／ｐｅｎ：预处理因子．ｐｅｎ－－－Ｏ时无预处理；：ｌ时，对角预处理；：２时．ＳＳＯＲ预 处理
ｃｏｎｓｔ ｍｔ ｍａｘｌｔｅｒ＝ｌｔｅ￣ ｉｎｔ ｋ： ＣＯＯＳｔ ｄｏｕｂｌｅ ｓｔｐ＝ｅｐｓ； ｄｏｕｂｌｅ ｅｐｒ； Ｖｅｃｔｏｒ ｒ＝ｐｒｅｃｏｎｄｉｎｇ（ｂ一（．ｔｈｉｓ）’ｘ＇ｐｅｎ）； ｄｏｕｂｌｅ ｂｅｔａ＝ｒ．ｔｗｏｎｏｒｍ０； ｂｏｏｌ ｃｏｎｙ－－ｆａｌｓｅ； ｉｆ（ｍ＞ｎｒｏｗｓ） ｅｎＤｒ（ ｍ ｉｓ ｂｉｇｇｅｒ ｔｈａｎ ｎｕｍｂｅｒ ｏｆ ｒＯＷＳ”）： ｉｆ（ｂｅｔａ＜￣ｔｐ）Ｉ ｅｐｒ＝ｂｅｔａ； ｉｔｅｒ＝ｌ； ｒｅｔｕｒｎ ０： Ｊ Ｖｅｃｔｏｒ＊‘ｖ＝ｎｅｗ Ｖｅｃｔｏｒ＊【ｍ＋ｌ ｆｏｒ（ｉｎｔ ｉ＝Ｏ；ｉ＜－－￣＋ｎ ｖ【ｉ Ｉ＝ｍｗ Ｖｅｃｌｏｌ ｍ
＝ １●●● ●●●●●●Ｊ ０ ０ ０ ２ ３ ０ ０ ２ ３ ０ ０ ２ ３ ０ ０ ２ ３ ０ ０ ０ ３ ０ ０ ０ ０ 维普资讯 http://www.cqvip.com １２Ｏ 福建电脑 ２００６年第２期 ｄｏｕｂｌｅ＊＊ｈ－＝ｎｅｗ ｄｏｕｂｌｅ．【ｎ ｆｏｒ（ｉ ｉｔ抽 ｉ＋＋） Ｉ．【ｉ ｎ＿ｄｏｅｂｌｅ【ｉ＋２ｋ ｉｔｅｒ＝ｌ； ｗｈｉｌｅ（ｉｔｅｓ＜－－ｎｍｈｅｄＩ ＊，［ｏＦ－ａｂｅ￣ Ｖｅｃｔｏｒ ｇ（Ｉｎ＋１）； ｂｅｌａ： Ｖｅｃｔｏｒ ｅｓ（ｍ），ｓｎ㈤； ｆｏｒ（ｋ＝０ ·锄＆＆ｉｔｅｒ＜＝ｒａａｘｉｔｅ￣＇，ｋ＋＋ｊｌｅｒ＋＋】Ｉ Ｖｅｃｔｏｒ ｗ＝ｐｒｅ￣ｎａｌｎ８（（’ｔｈｉｓ）’（＊ｖｔｋＤ．ｐｅｎ）； ｄｏａｂｈ￣．ｔｗｏｎｏｒｍＯ； ｆｏｒ（ｉ ｉ＜－ｋ；ｉ¨】Ｉ ｈ［ｋＩｉ］＝ｄｅｔ（ｗ，＊ｄｉＤ； ｗ一＝ｈ［ｋＩｉ］’ｆ．ｖ［ｉＤ； ｝ ｈ【ｋＩｋ＋ｌ｝ ．ｔｗｏｎｏｒｍＯ； ｉｆ（ｎａｍ＋１．Ｏｅ－４．ｈ ＋ｌ】＝ Ⅻ＿ｒ）ｌ ｆ Ｉ（ｉ ｉ《＝ｋｊ＋＋）Ｉ ｄｏｕｂｌｅ ｈｆｉ＝ｄｅｔ（ｗ，＊ｖ【ｉＤ； ｈ［ｋＩｉ］＋＝ｈｒｌ； ｗ－＝ｈｒｉ’ｆ．ｖｉｉｉ）； ｝ ｈｉｋＢ＋１］＝ｗ．ｔｗｏｎｏｒｍＯ； ｝ ｉｆ ｏＩ Ｉｋ＋ｌ】＝＝ ｅｎ＇ｏ￣＇ｕｎｅｘｐｅｅｔｅｄ ｔｅｍ－ｆｌｉｖｉｓｏｒ＇）； ‘ｖ【ｋ＋１］－－￣ｄｌｑｋＩｋ＋１］； ｆｏｒ（ｉ＝ｏ 吐ｊ＋＋）Ｉ ｄｏｕｂｌｅ ｔｍｐ￣ｉ］’ｈ［ｋＩｉ］＋ｍ［ｉ］’ｈ Ｉｉ＋ｌ ．】．ｈｌｔｌｉ］－￣ｉ］’ｈｉｋＩＩ＋ｌＩ ｈ［ｋＩｉｌ＝ｔｍｐ； 】 ｉｆ《ＩＩ【ｋＩｋ＋ｌ】＝＝ Ｉ ｅｓＩｋｋ－ｌ： ∞哟＝０ ｝ｄ叫 ｄｏｕｂｌｅ ｔｐＩ 甲ｔｏＩ Ｉｋ】’ｈＩｋＤ：］＋ｈｌｋＩｋ＋１］．ｈ 】【ｋ＋ｌ Ｄ； ｅｓ［ｋ］＝ｈ［ｋＩｋ］／ｔｐｍ； ５ｎ【ｋ】－ｈ【ｋＩｋ＋Ｈ／ｔｐｍ； ｌ ｈｌｋＩｋｌ＝￣ＰｈＯ，Ｉｋ ｍ ’ｈＯ，Ｉｋ＋ｌｋ ｄｏｕｂｌｅ ｔＩｎｐ ｓ【ｋ】．Ｒｋｌ＋ｓ４ｋ］’ｄｋ＋ｌ Ｒｋ＋Ｈ＝－￣ｋ］‘Ｒｋｌ＋ｅｓＩｋ］’Ｒｋ＋ｌｇ Ｒｋ］－ｔｍｐ； ｉｆ（ｆａｂ￣ｋ＋１］）＜＝ｓｔｐ）Ｉ ｋ＋＋： ｂｒｅａｋ； ｝ ｝／／ｅｎｄｆｏｒ ｆｏｒ（ｉ－ｋ—ｌｊ＞＝ —＿）ｌ ｆｏｒｉｎｔ ｊ－－ｉ＋１０＜ｋ０＋＋】ｄ 一＝ｈＤＩ ’ｄｊｋ ＝ｈ【ｉｎ ｝ ｆｏｒ Ｏ＝Ｏ ＜ｋ；ｉ＋＋】ｘ＋＝：ｄ ’ 哇ｉＤ； ｒ＝－－ｐｒｅｃｏｎｄｉａｇ（ｂ－（＊ｔｈｉｓ）’ｘ＇ｐｃｎ】； ｈｅｔａ－－－＇ｒ．ｔｗｏｎｏｒｍＯ； ｉｆ（ｂｅｔａ＜＝ｓｔｐ）Ｉ ｃｏｎｙ－－ｔｒｅｅ４ ｂｒｅａｋ； ｝ ｝／ｅｎｄ ｗｈｎｅ ｅｐＲ＝ ｆｏｒ ６＝０；ｉ＜－－ｍ；ｉ＋＋）ｄｅｌｅｔｅ ｖｉｉｉ； ｄｅＩｅｔｅ口ｖ； ｆｏｒ（ｉ＝Ｏ；ｉ＜ｍ；ｉ＋＋Ｒｌｄｅｔｅ［］ｈ【．】； ｄｅｌｅｔｅ［］ｈ： ｉｆ（ｅｏｎ＊）￣ｔｕｍ ｅｌｓｅ 咖ｍ ｌ： ｌ／／ｇｍ￣（ｍ）结柬
ｈ【ｋ】【ｋ Ｉ】 三总结 本文举了几个具有明显特征的数值例子，根据例子表明， ＧＭＲＥＳ算法的收敛效果取决于方程组系数矩阵的实际分布情 况。如谱的分布和 ｔｚ值与Ａ的特征值的逼近程度等。在实际 应用中，为了减少存储量和计算量，通常使用ＧＭＲＥＳ算法的重 新开始版本，在重新开始版本中，对重新开始迭代步长的选择也 是至关重要的，合适的ｍ能够加速ＧＭ砌 方法的收敛，过小则 可能导致迭代的收敛失败．过大则导致计算量的成倍增加及收 敛速度减缓。另外。对于预处理子的选择也对ＧＭＲＥＳ算法收敛 有重要的影响。
参考文献 ［１】Ｙ．Ｓａａｄ ａｎｄ Ｍ．Ｈ．Ｓｃｈｕｋｚ．ＧＭＲ￣：Ａ ｇｅｎｅｒａｌｉｚｅｄ ｍｉｎｉｍａｌ ｒｅｓｉｄｕａｌ ａｌｇｏｒｉｔｈｍ ｆｏｒ ｓ０ｌＶｉ ｇ ｎｏｎｓｙｎｍｌｅｕｉｃ ｌｉｎｅａｒ ｓｙｓｍｍｓ．ＳＩＡＭＪ．Ｓｃｉ．Ｓｔａｔｉｓｔ．Ｃｏｍｐｕｔ．， ７（１９８６）。８５６—８６９． 『２１钟宝江．一种灵活的混合ＧＭＲ．Ｅ５算法．高校计算数学学报，２００１。２３：２６１－２７２． ［３】 Ｙ．ａｎｄ Ｓｃｈ，ａ：．Ｍ．Ｈ．．ＧＭＲＥＳ：Ａ ｇｅｎｅｒａｌｉｚｅｄ ｍｉｎｉｍａｌ ｒｅｓｉｄｕａｌ ａｌｇｏｒｉｔｈｍ ｆｏｒ ｓｏｌｖｉｎｇ ｎｏｎｓｙｎｍｌｅｔｒｉｃ ｌｉｎｅａｒ ｓｙｓｔｅｍｓ．ＳＩＡＭＪ．Ｓｃｉ．Ｓｔａｔｉｓｔ． Ｃｏｍｐｕｔ．，１９８６，７：８５６—８６９． ［４】Ｓａａｄ，Ｙ．．Ｋｒｙｌｏｖ ｓｕｂｓｐａｃｅ ｍｅｔｈｏｄｓ ｆｏｒ ｓｏｌｖｉｎｇ ｈ ｕｍｙｍｍｃｔｒｉｃ ｌｉｎｅａｒ ｓｙｓｔｅｍｓ．Ｍａｔｈ．Ｃｏｍｐ．１９８１，１０５—２７． ｆ５１钟宝江．ＧＭＫＥＳ方法收敛率．南京学报。１９９８，５． ［６】ＮａｃｈｄｇａｌＮ…Ｍ ｇ．ｅｉｃｈｅｌＬ ａｎｄＴｒｃｆｃｔｈｅｎＬ．Ｎ．．ＡｈｙｂｒｉｄＧＭＫＥＳ ａｌｇｏｒｉｔｈｍｆｏｒ ｎｏｎｓｙｍｍｅｔｒｉｃｌｉｎｅａｒ ｓｙｓｔｅｍｓ ＳｍＮＪ．ＭａｔｒｉｘＡｎａｌＡｐｐｌ，１９９２，１３（３） ７９６—８２５． ｒ７】插道奇，王小鸽．ｃ＋＋和面向对象数值计算【ｃ＋＋。 ｄＯｂｊｅｃｔ—ＯｒｉｅｎｔｅｄＮｕｍｅｒｉｃＣｏｍｐｕｔｉｎｇ）ｔＭ］．人民邮电出版社，２００３．ＩｌＳＢＮ ７－１１５—１０６２６－６．
维普资讯 http://www.cqvip.com