52.gmres函数

在MATLAB中,提供了gmres函数实现广义最小残差法。函数的语法格式为:

x=gmres(A,b):针对x对线性方程组A∗x=b求解。n×n系数矩阵A必须是方阵,并且应为大型稀疏矩阵。列向量b的长度必须为n。A可以是函数句柄afun,这样,afun(x)返回A∗x。对于此语法,gmres不会重新启动;最大迭代次数为min(n,10)。

如果gmres收敛,则会显示一条有关该结果的消息。如果gmres无法在达到最大迭代次数后收敛或出于任何原因暂停,则会输出一条显示相对残差norm(b-A∗x)/norm(b)以及该方法停止或失败时所达到的迭代数的警告消息。

gmres(A,b,restart):每restart次内迭代重新启动该方法一次。最大外迭代次数为min(n/restart,10)。最大总迭代次数为restart∗min(n/restart,10)。如果restart为n或[],则gmres不重新启动,最大总迭代次数为min(n,10)。

gmres(A,b,restart,tol):指定该方法的容差。如果tol为[],gmres使用默认值1e-6。

gmres(A,b,restart,tol,maxit):指定最大外迭代次数,即,总迭代次数不超过restart∗maxit。如果maxit为[],gmres使用默认值min(n/restart,10)。如果restart为n或[],则最大总迭代次数为maxit(而非restart∗maxit)。

gmres(A,b,restart,tol,maxit,M)和gmres(A,b,restart,tol,maxit,M1,M2):使用预设子条件M或M=M1∗M2,并高效求解关于x的方程组inv(M)∗A∗x=inv(M)∗b。如果M为[],gmres不会应用预设子条件。M可以是函数句柄mfun,这样,mfun(x)返回M\x。

gmres(A,b,restart,tol,maxit,M1,M2,x0):指定第一个初始估计值。如果x0为[],gmres使用默认值(即全部为零的向量)。

[x,flag]=gmres(A,b,…):也返回一个收敛标志,其取值及含义见表1-5。

表1-5 flag取值

如果flag不为0,返回的解x具有在所有迭代中最小的范数残差。如果指定flag输出,则不会显示消息。

[x,flag,relres]=gmres(A,b,…):还返回相对残差norm(b-A∗x)/norm(b)。如果flag=0,则relres≤tol。第三个输出relres是预调节方程组的相对残差。

[x,flag,relres,iter]=gmres(A,b,…):还返回计算x时所达到的外和内迭代数,其中0≤iter(1)≤maxit且0≤iter(2)≤restart。

[x,flag,relres,iter,resvec]=gmres(A,b,…):还返回每次内迭代中的残差范数的向量。这些是预调节方程组的残差范数。

【例1-53】此示例演示如何使用预设子条件并重新启动gmres。

     %加载west0479,它是一个非对称的479×479实稀疏矩阵
     load west0479;
     A = west0479;
     %定义b以使实际解是全为1的向量
     b = full(sum(A,2));
     %构造一个不完全LU预设子条件
     [L,U] = ilu(A,struct('type','ilutp','droptol',1e-6));

使用重新启动的gmres的好处是限制执行该方法所需的内存量。如果不重新启动,gmres需要存储maxit向量来保存基本Krylov子空间。此外,gmres必须在每一步与以前的所有向量正交。重新启动可限制使用的工作区间以及每次外迭代执行的工作量。请注意,即使预调节的gmres在以上6次迭代中收敛,该算法也允许多达20个基向量,因此应事先分配所有这些空间。

     %请执行gmres(3)、gmres(4)和gmres(5)
     tol = 1e-12;
     maxit = 20;
     re3 = 3;
     [x3,fl3,rr3,it3,rv3] = gmres(A,b,re3,tol,maxit,L,U);
     re4 = 4;
     [x4,fl4,rr4,it4,rv4] = gmres(A,b,re4,tol,maxit,L,U);
     re5 = 5;
     [x5,fl5,rr5,it5,rv5] = gmres(A,b,re5,tol,maxit,L,U);

fl3、fl4和fl5均为0,因为在每种情况下重新启动的gmres使相对残差趋向小于1e-12的规定公差。

     %下列绘图显示了每个重新启动的gmres方法的收敛历史记录。gmres(3)收敛于外迭代5、
     %内迭代3(it3 = [5,3]),这将与外迭代6、内迭代0相同,因此最终刻度线上的标记是6
     subplot(3,1,1);semilogy(1:1/3:6,rv3/norm(b),'-o');
     h1 = gca;
     h1.XTick = [1:1/3:6];
     h1.XTickLabel = ['1';'';'';'2';'';'';'3';'';'';'4';'';'';'5';'';'';'6';];
     title('gmres(3)')
     xlabel('迭代次数');
     ylabel('相对残差');
     subplot(3,1,2);semilogy(1:1/4:3,rv4/norm(b),'-o');
     h2 = gca;
     h2.XTick = [1:1/4:3];
     h2.XTickLabel = ['1';'';'';'';'2';'';'';'';'3'];
     title('gmres(4)')
     xlabel('迭代次数');
     ylabel('相对残差');
     subplot(3,1,3);semilogy(1:1/5:2.8,rv5/norm(b),'-o');
     h3 = gca;
     h3.XTick = [1:1/5:2.8];
     h3.XTickLabel = ['1';'';'';'';'';'2';'';'';'';''];
     title('gmres(5)')
     xlabel('迭代次数');
     ylabel('相对残差');

运行程序,效果如图1-21所示。

图1-21 每个重新启动的gmres方法的收敛曲线