VC毕业论文GMRES算法的加速收敛现象分析毕业论文

《VC毕业论文GMRES算法的加速收敛现象分析毕业论文》由会员分享，可在线阅读，更多相关《VC毕业论文GMRES算法的加速收敛现象分析毕业论文（48页珍藏版）》请在装配图网上搜索。

1、学院理学学士论文 摘要摘要随着科学和工程技术的发展，越来越多的问题需要求解大规模的线性方程组，对这类方程的快速求解已成为数值代数研究的热点之一，特别是具有稀疏结构的大型方程组的求解。基于Galerkin原理的Arnoldi算法是求解这种线性代数方程组的近似算法，以下称这种方法为广义极小残余算法(GMRES算法)。GMRES方法是目前求解大型稀疏非对称线性方程组最为流行的一种迭代方法。GMRES算法在迭代过程中通常表现出一种加速收敛行为，随着迭代次数的增加，这种加速收敛现象越明显，即残量收敛会随着迭代步数的增加而逐渐得到改善。在CG方法中，这种加速收敛与Ritz值有密切关系。通过分析，我们发现G

2、MRES的加速收敛与其斜投影过程中产生的Ritz值对特征值的逼近程度有关系。在实际应用中，为了减少存储量和计算量，我们通常使用GMRES算法的重新开始版本来求解大型非对称线性方程组。本文描绘了GMRES和GMRES(m)的加速收敛现象，并通过实验给予解释。关键字： 广义最小残量； Krylov子空间； Ritz值； 加速收敛； 正交投影方法； 非对称线性方程组 I学院理学学士论文 AbstractOn The Superlinear Convergence of GMRESAbstractWith the development of science and project technolog

3、y, more and more questions need the solution of big linear systems. This solution is one of the fastest ways for researching numerical algebra, especially for the big sparse matrix. The way of Arnoldi is based upon the principle of Galerkin, which is closed to the solution of the linear numerical sy

4、stem. Here, we call the solution as Generalized Minimum Residual (GMRES). GMRES is one of the most popular iterative methods for the solution of big nonsingular nonsymmetric linear systems. It usually has a so-called superlinear convergence behavior. The rate of convergence seems to improve as the i

5、teration proceeds. For another say, the rate of residual variable will be improved as we increase its iteration. For the conjugate gradients method, this method has been related to a degree of convergence of the Ritz value. Through some analysis, we found that for GMRES too, changes in convergence b

6、ehavior seem to be related to the convergence of Ritz value. In our practical application, we also usually use GMRES(m) for reducing storage and counter solving big linear systems. This paper studies the superlinear convergence behavior of GMRES and GMRES(m), and supplies explain through experiment.

7、Keyword: GMRES; Krylov subspace; Ritz value; superlinear convergence; orthogonalization method; nonsymmetric linear system III学院理学学士论文 目录目录摘要IABSTRACTII第一章 引言1第二章 GMRES算法基础知识32.1 向量范数32.2 线性方程组最小二乘问题42.2.1 Gram-Schmidt正交化方法42.2.2 Givens变换4第三章 GMRES算法理论63.1 Krylov子空间方法的基本理论63.2 Arnoldi算法73.3 GMRES算法结

8、构8第四章 GMRES算法的加速收敛现象分析9第五章 数值示例与算法实现195.1 数值实验195.2 算法改进与实现225.2.1 预处理技术225.2.2 算法实现245.3 实验总结34致谢35参考文献36REPORT OF LITERATURE37文献报告41学院理学学士论文 第一章 引言第一章 引言关于线性方程组的数值解法一般分为两大类:直接法和迭代法。直接法是在没有舍入误差的情况下，通过有限步四则运算可以求得方程组精确解的方法。但是，在实际计算时，由于初始数据变为机器数而产生的误差以及计算过程所产生的舍入误差等都对解的精确度产生影响，因此直接法实际上也只能算出方程组真解的近似值。直

9、接法的基本思想是将结构上比较复杂的原始方程组，通过等价变换化成结构简单的方程组，使之变成易于求解的形式，然后再通过求解结构简单的方程组来得到原始方程组的解。即G通常是对角矩阵、三角矩阵或者是一些结构简单的矩阵。目前较实用的直接法是Gauss消去法的一些变形，例如选主元的Gauss消去法和矩阵的三角分解法，它们都是目前计算机上常用的有效方法。迭代法就是对任意给定的初始近似解向量，按照某种方法逐步生成近似解序列,使极限为方程组的解，即。因此迭代法是用某种极限过程去逐步逼近真解的方法，从而也可以用有限步运算出具有指定精确度的近似解。迭代法主要有Jacobi迭代法、Gauss-Seidel迭代法、逐次

10、超松弛法以及共轭斜量法。 直接法的优点是计算量小，并且可以事先估计，缺点是所需存储单元较多，编写程序复杂；迭代法的优点是原始系数矩阵始终不变，因而算法简单，编写程序也比较方便，且所需存储单元也较少，缺点是只有近似解序列收敛时才能被采用，而且存在收敛性和收敛速度的问题。对于中等规模的n阶(n100)线性方程组，由于直接法的准确性和可靠性，所以它们是经常被采用的方法。对于较高阶的方程组，特别是对某些偏微分方程离散化后得到的大型稀疏方程组(系数矩阵中绝大多数为零元素)，由于直接解法的计算代价较高，使得迭代更具竞争力。随着大规模并行计算机的快速发展，现在可以计算的规模越来越大，这些计算多数来源于偏微分

11、方程离散后得到的大型稀疏线性方程组。因此，大型稀疏线性代数方程组的求解已成为数值算法研究的热点问题。在许多应用学科和工程领域中，如流体力学、结构力学、航空航天工程、电子工程等等，经常会遇到大型稀疏非对称线性方程组问题。目前求解这类问题最流行的算法便是GMRES算法。本文主要分析GMRES算法的加速收敛现象和重新开始版本的GMRES(m)算法的加速收敛现象。2学院理学学士论文 第二章 GMRES算法基础知识第二章 GMRES算法基础知识2.1 向量范数为了研究迭代法的收敛性，我们需要对n维向量和n阶方阵引进某种度量-向量范数的概念。向量范数是三维Euclid空间向量长度概念的推广，向量范数在数值

12、分析中起着极其重要的作用。设为实n维向量空间，为复n维向量空间。我们记为或，为实数域或复数域。若上任一向量，对应一个非负实数，对于及，满足如下条件:非负性: ，且的充要条件是；齐次性: ；三角不等式:。则称为向量的范数(上述三个条件称为向量范数公理)。常用的向量范数有(1)1范数；(2)2范数；(3)范数 ；(4)一般的范数 。2.2 线性方程组最小二乘问题考虑确定一个向量，使函数 达到最小，这个问题称为线性最小二乘问题，因为它要求剩余向量 分量的平方和达到最小，而且线性地依赖于。其解亦称为方程组 的最小二乘解。现在介绍两种在GMRES算法中用到的正交变换方法，Gram-Schmidt正交化方

13、法、Givens变换。2.2.1 Gram-Schmidt正交化方法设是n维空间中n个线性无关的向量，则 这时是n个正交基向量 这时是正交矩阵，是上三角矩阵。Schmidt正交化过程就是将矩阵A分解为正交矩阵与上三角矩阵的乘积 即A的分解。2.2.2 Givens变换欲把一个向量中许多分量化为零，可以用Householder变换；如果只将其中一个分量化为零，则采用Givens旋转得到的分解。一个的Givens旋转是一个如下形式的归一矩阵 和满足。它将单位向量旋转到向量，第二项为0。 一个的Givens旋转将单位阵的对角线上的块替换为的Givens旋转 表示一个的带有在和行与列的的Givens旋

14、转。令 和则5学院理学学士论文 第三章 GMRES算法理论第三章 GMRES算法理论3.1 Krylov子空间方法的基本理论Krylov子空间投影法是一类非常有效的大型线性方程组解法，已被广泛应用于各种领域，随着左右空间的不同取法可以得到许多求解大型方程组的求解算法。考虑线性方程组 其中，是阶的大型稀疏实矩阵， 1.当时，称为Galerkin方法又称正交投影法。A为一般矩阵时，有FOM算法；A为对称正定矩阵时，有CG算法。2.当时，称为斜投影方法。当时，称为最小二乘法，当A为一般矩阵时，有GMRES方法，GCR方法等；当A为对称正定阵时，有CR方法。所以GMRES算法又属于最小二乘法。下面是迭

15、代方法的流程图:考虑系数矩阵是否对称使用CG方法是否正定是是否MINRES方法 否应用GMRES方法是否有足够的存储空间否图3-1 迭代方法流程图Fig3-1 Chart of iteration methods对于以上大型线性方程组，令右空间和左空间，其中各自线性无关。假设为初始迭代值，投影方法寻求这样的一种近以解: 通过求解一个最小值问题: 来解方程组,这就是Krylov子空间投影方法。3.2 Arnoldi算法Arnoldi算法过程:(1) 输入 和初始向量 。(2) 。(3) ， ， 。(4) 如果 ，则转步(5) ；否则 ，转步(3) 。(5) 求解得 ， ，输出 ，结束 。这一算法

16、就是所谓求解线性方程组的Arnoldi方法。容易证明，由Arnoldi方法得到的满足：，即在没有误差出现的条件下，使用Arnoldi方法将在有限步内给出方程组的精确解。在误差出现时，当然也可用 是否很小来判断迭代是否终止。但是，当j很大时，由于整个计算过程需要保存所有的向量，因而需要占用大量的内存，以至于使得所处理的问题的规模受到很大的限制。为了克服Arnoldi方法的这一缺点，通常使用时，先选择一个恰当的正整数m(一般不宜太大)，计算到 时，就去求，然后再以作为初始向量重新开始，这样周而复始直到求出所需精度的近似解为止。这就是所谓的循环Arnoldi方法。这样做，虽然解决了存储问题，但又会遇

17、到方程组不相容的情形，致使迭代半途而废。针对循环Arnoldi方法的缺点，我们讨论基于Galerkin原理的Arnoldi算法也即广义极小剩余法(GMRES)。3.3 GMRES算法结构GMRES算法步骤:45学院理学学士论文 第四章 GMRES算法的加速收敛现象分析第四章 GMRES算法的加速收敛现象分析在这一章中，我们主要分析用GMRES方法1求解大型稀疏非对称线性方程组 (4.1)中表现出的加速收敛现象。设是(4.1)解的初始估计，是初始残量，为由和产生的Krylov子空间。GMRES方法的第k步迭代解和迭代残量满足 (4.2) (4.3)记为次数不超过k且常数项等于1的所有代数多项式全

18、体，根据(4.2)容易得到，从而存在多项式使得由此(4.3)可以等价地表示为 (4.4)为GMRES残量多项式。GMRES残量收敛率会随着迭代步数的增加逐渐得到改善。其实这一现象在CG方法的应用中也出现过。A.Van der Vorst2证明了Arnoldi过程中Ritz值对特征值的逼近和CG方法残量加速收敛现象的因果关系。但未能确定GMRES残量多项式和Ritz值之间的关系而类似CG残量多项式却可以推广到FOM残量多项式。Van der Vorst3等利用这一性质解释了GMRES方法的加速收敛行为，但是这一解释只有当GMRES残量多项式和FOM残量多项式按范数大致相等时才能成立，否则，缺乏一

19、定的说服力。针对这一点，钟宝江4通过GMRES与对照GMRES过程利用GMRES残量多项式实质上是以一个伴随GMRES迭代的斜投影过程产生Ritz值为零点这一性质合理的分析了GMRES方法的收敛加速现象。但在实际的应用中执行整体的GMRES所需要的计算量和存储量会随着迭代步数的增加而变得让人难以接受。为了克服这一困难，Y.Sadd等5给出了GMRES算法的重新开始版本。重新开始的GMRES算法的收敛性在理论上得不到保证，至今还不能像直接法那样形成一种通用的数学软件，在具体的计算时需要密切观察其收敛状态已做出相应的调整，因此对研究迭代过程的理论尤为重要。在许多学科和工程领域中如：流体力学，结构力

20、学，航空航天工程，电子工程等等经常遇到求解大型稀疏非对称线性方程组问题，而目前求解此类问题的最为流行的方法便是GMRES方法。由于大型方程组的阶数可达10以上，执行算法一般为其重新开始版本。因此，对GMRES算法的迭代过程的研究具有尤为重要的现实意义和应用价值。我们来看一下求解系数的值和向量 的斜投影过程 (4.5)记为GMRES迭代过程中由Arnoldi算法生成的标准正交向量系列，则矩阵和的列分别组成子空间和的一组基。根据,设，这里。(4.5)等价于广义特征值问题 (4.6)记，根据Arnoldi算法有 其中形如记，其中为阶单位阵，则由此(4.6)可以具体地写作 (4.7)4证明了下面的结论

21、。引理1. 设,则GMRES残量多项式满足 (4.8)其中为问题(4.7)的个非平凡特征值。根据6，当GMRES方法的第k步收敛时即可等价地保证。以下首先分析当有若干个值逼近矩阵的某个特征值，比如时，对GMRES收敛过程将产生的影响。这里所采用的分析方法和2,3是一致的。设n阶非奇异矩阵约化为Jordan标准形，即 其中阶子矩阵形如约定标准形中同一特征值所对应的Jordan块按阶数从大到小相邻顺序放置。该约定可以通过右乘一个交换阵来实现。当时，定义有理分式 若注意到存在使得，则可以由残量关系式和(4.8)约化为更低阶。为此只需考虑时的情形。此时矩阵函数都是有意义的。利用形如和矩阵的可交换性质，

22、可以得出设所对应的Jordan块有t个，即根据Jordan块的排列规则有则 (4.9)定理1将建立GMRES残量和一个对照GMRES过程的残量之间的联系，这一对照过程定义如下: 定义1. 记，其中.相应地，记.对照GMRES过程是一个以作为初始残量的GMRES迭代过程(也就是在实际的第k步残量中去掉所对应的分量)，它的第步残量及残量多项式分别记为.如果，根据残量关系式，相应的特征值对GMRES迭代过程不会产生影响，我们就说在这里不参与迭代。定理1. 设，对任意的有 其中的证明. 根据定义1有 (4.10)根据(4.9)，显然 (4.11)由于，所以 (4.12)综合(4.10，4.11，4.1

23、2)和GMRES残量最小化性质(4.4)得到 类似于3，我们可以给出量的一个上界: (4.13)当是可对角化的，即对所有的时，(4.13)可以简化为 (4.14)在对照GMRES过程中，由于不再参与迭代，系数矩阵的谱分布将因此得到改善，从而这一过程会比实际GMRES过程有更快的收敛率。根据定理1和(4.13)，若在GMRES的第步迭代时，有个值能够很好地逼近，使得充分接近于1，则此后其迭代残量将和对照GMRES过程的残量是同阶收敛的，另外，根据求解矩阵特征值问题斜投影方法的一般理论，随着迭代步数的增加，由(4.5)所得的值总是收敛的，这样就解释了GMRES方法的加速收敛行为。由于计算量和存储量

24、的限制，实际使用的GMRES方法一般为其重新开始格式，所以应该选取重新开始GMRES方法的迭代步长够恰当，以使得在每次迭代循环中值都能够充分逼近一些对迭代过程有关键影响的特征值，尤指极端特征值，否则残量的收敛率可能非常慢。结合1对GMRES方法做出的收敛分析，应用定理1也可以建立残量基于的一个上界，这对分析GMRES方法具体的迭代过程是非常有用的。但得出的表达式含有因子，为了避免这一缺陷，定理2给出的表达式中仅含有因子.这里只考虑是可对角化的，即对任意的。所得出的结论可以直接推广到Jordan块时的情形。定理2. 设是可对角化的且定义 (4.15)则 其中如(4.14)中定义。证明. 根据7,

25、P.115,满足(4.15)的多项式是唯一存在，记为.注意到，由GMRES残量最小化性质(4.4)得到 定理2考虑了当被值逼近时GMRES方法的收敛情况，当有多个特征值同时被值逼近时，结论可作如下推广。定理3. 设是可对角化的且.记为一个由的个不同特征值所组成的集合，并记. 定义则 证明. 具体的证明过程和定理2的证明是一致的，只需将其中多项式和有理分式分别选择为 满足8在实际应用中，我们所使用的GMRES算法一般为其重新开始版本，所以我们有必要重点对重新开始GMRES算法的加速收敛特性进行讨论。其重新开始版本也就是使它其中的m固定不变，我们也称其为GMRRS(m)算法。算法: GMRES(m

26、)(1)开始:选取，计算(2)迭代:对做 (3)计算近似解: 其中极小化(4)重新开始:计算，如果满足精度要求则停止；否则令返回(2)重新开始GMRES算法的收敛性在理论上得不到保证，至今还不能像直接法那样形成一种通用的数学软件，在具体计算时需要密切地观察其收敛状态以做出相应的调整。因此对算法迭代过程的理论研究尤为重要。以下我们将描述重新开始GMRES算法的一种加速收敛特性。首先看一个算法加速收敛的数值例子，选取矩阵具有最简单的非对称形式，也即程序中我们所说的二对角矩阵 初始残量。选取，误差条件=1.0e-6, 进行GMRES(2)迭代，具体数值表明，在算法的前25次迭代循环中，残量收敛速度稳

27、定在；而从第26次循环迭代开始，收敛速度变为零。在算法收敛加速的过程中，我们来看一下它的收敛过程,如下表4-1:表4-1 收敛过程数据表Table4-1 Table of convergence process data迭代次数绝对残差残量收敛速度10.7992770.55500920.4436060.64952230.2881320.66342040.1911530.66603050.1273130.66654160.0848600.66664270.0565710.66666280.0377140.66666690.0251420.666666100.0167620.666667110.0

28、111740.666667120.0074500.666667130.0049660.666667140.0033110.666667150.0022070.666667160.0014720.666667170.0009810.666667180.0006540.666667190.0004360.666667200.0002910.666667210.0001940.666667220.0001290.666667230.0000860.666667240.0000570.666667250.0000380.666667260.0000260.000000表4-1的数据显示，在算法收敛加速

29、前的很长一段时间内，迭代值的变化是非常稳定的。由此可以设想，重新开始GMRES算法的这种加速收敛行为是特征性的，它仅依赖于循环迭代这个过程。为了验证这一设想，我们构造如下的一个混合多项式算法。HGMRES(m):(1) 执行GMRES(m)至迭代第s次迭代循环，形成残量多项式(2) 以作迭代多项式，执行Richardson迭代直至收敛 当选取s=1时，算法HGMRES(m)即9混合GMRES算法。对应于s的其它取值，Richardson迭代的迭代多项式可以是GMRES(m)在任一迭代循环中产生的残量多项式。根据表4-1的数据，我们取s=15，即HGMRES(2)的前14次迭代循环和GMRES(

30、2)完全相同，而从第15次循环开始，HGMRES(2)的迭代多项式固定为执行结果表明，新算法以完全的循环迭代过程很好地模拟了重新开始算法的加速收敛行为。根据以上讨论，下面借助HGMRES(m)算法，对这种加速特性做出分析。设存在N阶非奇异矩阵约化为约当标准形，即 其中阶子矩阵具有形式 在HGMRES(m)算法，则有 或其中 其子矩阵 很明显，只要便有此时，等价地，.但在算法收敛的过程中，考虑到不同的初始值，若在迭代开始的某一段时间内，矩阵的某些元素不是单调下降的，则残量的收敛速度就有可能受到很大的影响。而容易看出，这个问题可随着迭代循环次数k的增加而得到改善。对应地，残量收敛将呈加速趋势。下面

31、大致地给出算法收敛加速所需的循环次数。记，且设。考虑第个若当块，有(1) 主对角线元素: 当时，即单调收敛；(3) 次对角线元素:由可得，当时单调收敛；(4) 第三条对角线元素:当时，右端第一项单调收敛；由(2)知，此时第二项也不单调收敛；右上角元素: 当时，右端第一项单调收敛；根据(1)，(2)，以后各项都单调收敛。由以上讨论，一般当算法的迭代循环次数后可以期望残量的收敛加速。对于GMRES(m)，若从某个时刻开始各次迭代循环所产生的残量多项式都满足 其中为正常数，则可相应地讨论其加速收敛特性。学院理学学士论文 第五章 数值示例与算法实现第五章 数值示例与算法实现5.1 数值实验首先选取一个

32、三阶线性方程组来验证GMRES算法的正确性，选如下一个数值例子取 x 0=0, 0, 0，迭代步长 m=3, 误差=1.0e-6, 计算结果为x=0.518519, 0.222222, 0.666667, 把计算结果代入以上三阶方程组中，方程组成立。然后我们选取一个数值例子来看一下GMRES算法的加速收敛现象，选取如下方程组:取x 0=0, 0, 0, 0, 0, 0，迭代步长m=2，误差=0.0001，计算结果如表5-1所示。表5-1 计算结果表Table5-1 Table of calculating result迭代次数绝对残差残量收敛速度10.3123480.45365420.1416

33、980.61216530.0867420.62848140.0545160.64025350.0349040.64118760.0223700.64120270.0143430.64120880.0091970.64120990.0058970.641209100.0037810.641209110.0024250.641209120.0015550.641209130.0009970.000000分析表5-1可知，算法在迭代过程中残量收敛速度加快，并在第13次迭代后出精确解。而对同一方程组在其它条件不变的情况下，设m=3时，其结果在迭代第11次时，就已经收敛。设m=4时，其结果在迭代第26次

34、时，才收敛。所以我们在讨论GMRES(m)算法时，必须考虑一个重要问题是合理选取GMRES迭代步长m，以使得在其每次迭代循环中Ritz值都能够充分逼近一些对迭代过程有关键影响的特征值，否则残量的收敛速度可能非常慢。以下选择了一个具备明显特征值的数值例子来演绎收敛定理，为了进一步研究收敛特征，我们在计算过程中将m设定为20,误差条件设为10e-10。其中系数矩阵 形如： 并取，。迭代过程中Ritz值根据（4.7）求得。已知A是可对角化的，其特征值 1，2，3 = 1；i = i(i=4，5，.，20)全部落在正半实轴上。根据文献1 的定理5，容易证明在这种情况下GMRES方法有如下收敛估计： ，

35、其中max和min分别表示参与迭代过程的最大特征值和最小特征值。显然，A的谱中极小特征值是否参与迭代过程将对GMRES收敛率变化起关键作用。为此，考察迭代过程中极小Ritz值1(k)的收敛情况（见表5-2）。当k=8时，1(k)对1已有很好的逼近，使得成为一个适当的值。根据定理1，次后的GMRES过程将和1，2，3不参与迭代的对照GMRES过程有着同阶的收敛率。则 1k充分逼近1的前后残量的收敛因子分别是： 对应的，如表（5-3）所示GMRES在第八次迭代后表现出加速收敛行为，并在第十七次迭代收敛。其中第k步残量收敛的快慢由量的大小表示。表5-2 Ritz值收敛情况表Table5-2 Tabl

36、e of Ritz convergence instancek567891(k)5.503.312.121.351.1010.004.341.851.301.05表5-3 残量速度收敛表Table5-3 Table of residual variable convergence迭代次数绝对残差残量收敛速度12.238880.7115821.593160.9003431.434390.9942341.426111.0056151.434121.0002761.438101.0009571.439471.0002781.439511.0000691.439520.24699100.355540.

37、46473110.165230.51726120.089530.51726130.046310.53833140.024930.54231150.013520.63683160.008610.65621170.005650.000005.2 算法改进与实现从上面的数值示例可以看出，一般情况下，GMRES算法收敛速度较慢，为了提高GMRES算法的加速收敛速度，我们提出一种预处理技术以加快算法收敛。5.2.1 预处理技术我们来讨论一下怎样选择预处理子P来加速迭代过程。理想的情况是，应该能够找到一个P，使得的条件数小于A的，且方程很容易求解。一些简单的预处理子如下10：1. 对角化预处理子。令为一对

38、角矩阵，并且的主对角线与的主对角线相同；就是说，。如果的主对角线上有零元素，令其中。2. SSOR(对称连续超松弛)预处理子。令，其中，为松弛系数，为一个三对角矩阵，对角线元素为零，为一对角化矩阵，F为一上三角矩阵，对角线元素为零，这三个矩阵满足。注意，如果是共轭的，则SSOR预处理子也是共轭的。3. LU (LU分解) 预处理子。假设的LU分解形式A=LU，其中L为一单位下三角矩阵，U为一上三角矩阵。那么，令P=LU。不完整的LU预测器为，其中和是将A中的零元素对应的L和U中元素置为零而得到的。因此，A的稀疏结构可以用来保存和。下面举个数值例子来证明预处理技术加速GMRES算法收敛： =b对

39、于以上6阶方程组，选取迭代步长m=2，误差=1.0e-6, 在没有选择任何预处理子时，它迭代19次后，算法收敛。迭代过程如下表5-4： 表5-4 无预处理子残量收敛表Table5-4 No disposal residual variable convergence迭代次数绝对残差残量收敛速度10.7726670.56665420.4378350.64887630.2841010.66314340.1884000.66590750.1254570.66645760.0836110.40795570.0341100.15679980.0053480.63672290.0034050.180858

40、100.0006120.344478110.0002120.590836120.0001250.016349130.0000021.000000140.0000021.000000150.0000021.000000160.0000021.000000170.0000021.000000180.0000021.000000190.0000020.000000而在保持其它条件不变的情况下，选择SSOR预处理子时，GMRES算法迭代11次后收敛。迭代过程见表5-5：表5-5 SSOR预处理残量收敛表Table5-5 SSOR disposal residual variable convergen

41、ce迭代次数绝对残差残量收敛速度10.5945230.09076320.0539610.13058530.0070460.06807040.0004800.11489850.0000550.17465760.0000100.12259070.0000011.00000080.0000011.00000090.0000011.000000100.0000011.000000110.0000010.0000005.2.2 算法实现为了更好的实现算法，我选择c+作为编程语言（其编程环境是Windows XP + VC+ 6.0）。C+语言作为C语言的继承者，它首先继承了C语言的以下特点：丰富的运算符

42、和数据类型、结构化的程序设计方法、高效的机器代码、良好的可移植性等。在程序设计方法方面，C+是一种混合型程序设计语言，它既支持传统的面向过程程序设计（procedure oriented programming, POP），也支持当代的面向对象程序设计（object oriented programming, OOP）11。首先定义一个矩阵类SparseMtx, 其定义如下:Class SparseMtx int nrows;/矩阵的行数private: /压缩稀疏行格式 int lenth;/原矩阵中非零元个数 double* sra;/存储非零元的数组 int* clm;/存储在sra中的

43、元素的列下标 int* fnz; /sra的每一行中的第一个非零元素 Vector preconding(const Vector&,int i=0); /预处理函数public: double epsR; /输出时的绝对残差SparseMtx(int n,int m,double* t,int* c,int* f);/n: 原矩阵的行数(列数)/m: 非零元数组sra的长度/t: 原矩阵中的非零元/c: sra中元素的列下标/f: sra中每一行第一个非零元的下标 SparseMtx(int n,int m);/将所有元素初始化为0/n: 行数(列数)/m: 非零元的数目SparseMtx(

44、const SparseMtx&);/拷贝构造函数SparseMtx() /析构函数delete sra;delete fnz;delete clm;SparseMtx& operator=(const SparseMtx&); /重载=Vector operator*(const Vector&) const; /矩阵向量积double& operator(int i) const /下标取值return srai;int& getfnz(int i) const /每行中的第一个非零元return fnzi;int& getclm(int i) const /sra中元素的列下标retur

45、n clmi; int GMRES(Vector&,const Vector&,double&,int&,int);采用以上定义的矩阵是为了减少计算存储量，而且我们采用以上适用于任意编号节点生成的大型稀疏矩阵的按行存储非零元存储格式（压缩稀疏行格式）。在这种格式中，n乘n的矩阵A存储在三个一维数组中：数组sra逐行存储A非零元。假设A有nz个非零元，sra的元素个数就是nz。数组clm存储sra元素的列下标。即srai在矩阵A中处于列clmi。clm的元素个数为nz。数组fnz有n+1个元素，fnzn=nz。fnzi 存储了A的第i行的第一个非零元在sra中的位置，fnzn存储了A的最后一个非

46、零元素后面一个元素在sra中的位置。也就是说，srafnzi给出A的i行第一个非零元素，0in，而srafnzn-1代表了A的最后一个非零元素。这里以一个6 6 的矩阵A 为例。可用如下压缩稀疏行格式存储：sra=1，2，3，4，2，5，7，3，8，4，1，3, 5；clm=0，1，2，0，1，3，1，2，3，3, 4，4, 5；fnz=0，3，6，9，11，13。对于用压缩稀疏行格式存储的n乘n的矩阵和向量v来说，矩阵-向量乘法可以这样实现（假设结果存储在向量w中）：for(int i=0;in;i+)/矩阵向量乘法for(int j=fnzi;jfnzi+1;j+)wi+=sraj*vcl

47、mj;/wi存储积的第i行从上面矩阵定义的代码可以清楚地看出用稀疏存储是相当方便的。GMRES算法编程：int SparseMtx:GMRES(Vector& x , const Vector& b, double& eps, int& iter, int pcn, int m)/带重启的GMRES算法，求解Ax=b，称为GMRES(m) /算法如果成功返回0，如果失败返回1/算法输出近似解，残差和迭代次数/A:任何非奇异矩阵，可能不对称/x:输入时，初始向量。返回时，近似解/b:右边向量int j = 0;j nrows)error(In GMRES(m), m is bigger than

48、 number of rows);if(beta=stp) /如果初始条件满足解epr=beta; /返回iter=1;return 0;/Krylov子空间的正交基/vi指向它的第i个基向量Vector* v=new Vector* m+1;for (int i=0;i=m;i+)vi=new Vector(nrows); /第i个基向量/Hessenburg矩阵h,(m+1)*m;/hi存储h的第i列，有i+2项/只有h的非零部分被存储double* h=new double*m;for( i=0;im; i+) hi=new double i+2;iter=1;while (iter=m

49、axiter) /gmres(m)的迭代 *v0=r/beta;Vector g(m+1); /最小二乘问题的向量 g0=beta;Vector cs(m),sn(m); /Givens旋转for (k=0;km & iter=maxiter;k+,iter+) /正交化 Vector w=preconding(*this)*(*vk),pcn);double nmw=w.twonorm();for( i=0;i=k;i+)hki=dot(w,*vi);w -= hki*(*vi);hkk+1=w.twonorm(); /如果hkk+1很小，重新正交化if (nmw+1.0e-4*hkk+1=

50、nmw)for( i=0;i=k;i+)double hri=dot(w,*vi); /重新正交化hki+=hri;w -= hri*(*vi);hkk+1=w.twonorm();if (hkk+1=0)error(unexpected zero-divisor in GMRES();*vk+1=w/hkk+1;/对h的第k列作用Givens G0,G1,.G(k-1)for( i=0;ik;i+)double tmp=csi*hki+sni*hki+1;hkk+1=-sni*hki+csi*hki+1;hki=tmp;/从h的第k列产生Givens旋转Gkif(hkk+1=0)csk=1;

51、snk=0;elsedouble tpm=sqrt(hkk*hkk+hkk+1*hkk+1);csk=hkk/tpm;snk=hkk+1/tpm;/对h和g的第k列作用Givens旋转Gk hkk=csk*hkk+snk*hkk+1;double tmp=csk*gk+snk*gk+1;gk+1=-snk*gk+ csk*gk+1;gk=tmp;if(fabs(gk+1)=0;i-)for(int j=i+1;jk;j+) gi-=hji*gj;gi/=hii; /更新解x=x0+sum vi*yi for (i=0;ik;i+) x +=gi*(*vi);r =preconding(b-(*

52、this)*x, pcn); /计算残差，检查是否收敛beta=r.twonorm();if(beta=stp) /如果残差很小，终止conv=true;break;/end while epR= epr=(b-(*this)*x).twonorm(); /得到最后的残差/为v和h重新分配空间for (i=0; i=m;i+)delete vi;delete v;for(i=0;im;i+)delete hi;delete h;if(conv)return 0;else return 1; /gmres(m)结束以下是一个GMRES算法运算界面，做界面的目的是为了更好地分析GMRES算法。我做的是一个基于对话框的应用程序。运算主界面如图5-1:图5-1 运算界面图Fig5-1 Chart of calculating图5-1这个GMRES运算对话框功能主要分成三个部分，它们分别是输入、检验、输出。输入部分包括输入所要计算的方程组、输入参与GMRES的参数。检验部分主要是对稀疏压缩矩阵进行长度匹配检验和对重启GMRES迭代中的m进行检验。输出部分主要由列表输出结果和图形输出结果两个功能。在输入所要计算方程组中，有三种矩阵输入:一般矩阵输入、特殊二对角矩阵输入、特殊三对角矩阵输入。界面分别如下图5-2、5-3、5-4: