gmres源程序matlab

合集下载

GM(1,1)模型中的MATLAB程序

GM （1,1）模型中的MATLAB 程序一、GM （1,1）模型的建立：（1）、一次累加生成序列的MA TLAB 命令：>> X0=[142 340 200 500 900 800 490 980 463 1100];>> X1(1)=X0(1)X1 =142>> for k=2:10X1(k)=X1(k-1)+X0(k)endX1 =142 482 682 1182 20822882 3372 4352 4815 5915（2）、由一次累加生成序列紧邻均值生成的)1(Z 的MA TLAB 命令：>> X0=[142 340 200 500 900 800 490 980 463 1100];>> X1(1)=X0(1);>> for k=2:10X1(k)=X1(k-1)+X0(k)Z(k)=(1/2)*(X1(k)+X1(k-1))EndX1 =142 482 682 1182 20822882 3372 4352 4815 5915Z =1.0e+003 *0 0.3120 0.5820 0.9320 1.6320 2.4820 3.12703.86204.58355.3650（3）、GM （1，1）的灰微分方程模型为：b k aZ k X =+)()()1()0(。

设∧α为待估计参数向量，⎥⎦⎤⎢⎣⎡=∧b a α。

利用最小二乘法得到Y B B B ')'(1-∧=α，MA TLAB 程序如下：>> B=[-Z(2:10)',ones(9,1)];>> Y=(X0(2:10))';>> alfa=inv(B'*B)*B'*Yalfa =-0.1062371.6018（4）、GM （1，1）的灰微分方程模型 b k aZ k X =+)()()1()0(的时间相应序列为：ab e a b X k X ak +⋅-=+-∧))1(()1()0()1( 由.6018.371,1062.0=-=b a 令.)1(,)0(u X v a b u -== 计算得到1.3499-=u ， 1.3641=v 。

gmres源程序 matlab

function [x，flag,relres,iter,resvec] = gmres(A,b,restart,tol,maxit,M1，M2,x,varargin)%GMRES Generalized Minimum Residual Method。

％X = GMRES(A,B) attempts to solve the system of linear equations A*X = B% for X。

The N-by—N coefficient matrix A must be square and the right％hand side column vector B must have length N。

This uses the unrestarted% method with MIN(N,10) total iterations.%％X = GMRES（AFUN,B）accepts a function handle AFUN instead of the matrix% A。

AFUN(X）accepts a vector input X and returns the matrix-vector% product A*X. In all of the following syntaxes, you can replace A by% AFUN.％％X = GMRES(A，B,RESTART）restarts the method every RESTART iterations.％If RESTART is N or [］then GMRES uses the unrestarted method as above.％% X = GMRES(A，B,RESTART,TOL）specifies the tolerance of the method。

If% TOL is ［] then GMRES uses the default, 1e-6.％% X = GMRES（A，B,RESTART，TOL，MAXIT）specifies the maximum number of outer % iterations。

matlab经典源程序带有注释(详细经典)

2.1set 与get 函数 (1)2.2callback函数 (2)2.3元胞数组 (4)2.4结构数组 (6)2.5矩阵操作 (9)2.6字符串操作 (13)2.7判断函数使用大全 (16)2.11打开外部程序 (21)2.11程序运行时间 (22)2.14动画 (23)2.12动画 (24)2.23显示多行内容 (26)2.24 uitable 使用 (26)2.27鼠标操作 (27)2.28键盘操作 (27)2.32粘贴板 (28)2.1set 与get 函数set(edit_handle,'String','my value!'); %String为Edit控件的属性%%%2.1-1%创建figure对象hfig=figure(1);%创建坐标轴对象，指定其父对象为figure 1haxes1=axes('parent',hfig);prop.Color='b';prop.FontSize=12;set(haxes1,prop);%%%2.1-2hfig=figure(1);%查询其Units属性值get(hfig,'units')%其Units属性值为pixels（像素）% ans=% pixels%%%2.1-3%figure的Pointer属性标识了鼠标指针的形状set(gcf,'pointer');% 返回值为：[ crosshair | fullcrosshair | {arrow} | ibeam | watch | topl | topr | botl | botr | left | top | right | bottom | circle | cross | fleur | custom | hand ]%%%2.1-4%首先取得标识电脑屏幕大小的度量单位get(0,'units')% ans =% pixels%取得屏幕的尺寸get(0,'screensize')% ans =% 1 1 1280 8002.2callback函数%定义M文件的主函数名称为DefineCallback，不带输入和输出参数function DefineCallbackhFig= figure('units','normalize',...'position',[0.4 0.4 0.3 0.2]);%在窗口中创建按钮控件，并定义其Callback属性uicontrol('parent',hFig,...'style','pushbutton',...'String','Execute Callback',...'units','normalize',...'position',[0.4 0.4 0.3 0.2],...'callback',['figure;',...'x = 0:pi/20:2*pi;',...'y = sin(x);',...'plot(x,y);']);%定义M文件的主函数名称为DefineCallback，不带输入和输出参数function DefineCallback%创建界面窗口hFig= figure('units','normalize',...'position',[0.4 0.4 0.3 0.2]);%在窗口中创建按钮控件hpush=uicontrol('parent',hFig,...'style','pushbutton',...'String','Execute Callback',...'units','normalize',...'position',[0.4 0.4 0.3 0.2]);%设置按钮的Callback属性set(hpush,'callback',@mycallback);%定义回调函数为子函数function mycallback(hobj,event)figure;x = 0:pi/20:2*pi;y = sin(x);plot(x,y);2.3元胞数组a={'hello' [1 2 3;4 5 6];1 {'1''2'}}a ='hello' [2x3 double][ 1] {1x2 cell }%示例2：将元胞数组a中的元胞逐一赋值>> a{1,1}='hello';a{1,2}=[1 2 3;4 5 6];a{2,1}=1;a{2,2}={'1' '2'};>> aa ='hello' [2x3 double][ 1] {1x2 cell }%示例3：使用cell函数来创建元胞数组%生成2x3的元素为空数组的元胞数组>> a=cell(2,3)a =[] [] [][] [] []%示例4：判断数组A是否为元胞数组%定义一个元胞数组A>> A={1 2 3};%判断A是否为元胞数组，如果为元胞数组，则函数>> tf = iscell(A)tf =1%示例5：显示元胞数组C中的内容>> clear>> C={'Smith' [1 2;3 4] [12]};%直接显示元胞数组C中的内容>> celldisp(C)C{1} =SmithC{2} =1 23 4C{3} =12%显示元胞数组C中的内容，数组的名称用cellcontent代替>> celldisp(C,'cellcontent')cellcontent{1} =Smithcellcontent{2} =1 23 4cellcontent{3} =12%示例6：将字符数组转换为元胞数组>> S = ['abc '; 'defg'; 'hi m'];>> cellstr(S)ans ='abc'%原先abc后面的空格被清除'defg''hi m'%i和m之间的空格仍然保留%示例7：显示元胞数组S中的内容（包括空格和字符）>> S = {'abc ', 'defg','hi m'};>> cellplot(S)%示例8：将数字数组A按行或按列转换为元胞数组%A是4x3的数组>> A=[1 2 3;4 5 6;7 8 9;10 11 12];%把A的每一列转换为一个元胞，得到的C是1×3的元胞数组>> C=num2cell(A,1)C =[4x1 double] [4x1 double] [4x1 double] %把A的每一行转换为一个元胞，得到的C是4×1的元胞数组>> C=num2cell(A,2)C =[1x3 double][1x3 double][1x3 double][1x3 double]2.4结构数组%示例1：使用直接法来创建结构数组>> A(1).name = 'Pat';A(1).number = 176554;A(2).name = 'Tony';A(2).number = 901325;>> AA =1x2 struct array with fields:namenumber%示例2：利用struct函数来创建结构数组>> A(1)=struct('name','Pat','number',176554);A(2)=struct('name','Tony','number',901325);>> AA =1x2 struct array with fields:namenumber%示例3：使用deal函数来得到结构体中各结构域的值%定义结构数组A>> = 'Pat'; A.number = 176554;A(2).name = 'Tony';A(2).number = 901325;%得到结构数组中所有name结构域的数据>>[name1,name2] = deal(A(:).name)name1 =Patname2 =Tony%示例4：使用getfield函数来取得结构体中结构域的值%定义mystr结构数组>> mystr(1,1).name = 'alice';mystr(1,1).ID = 0;mystr(2,1).name = 'gertrude';mystr(2,1).ID = 1;%取得mystr(2,1)的结构域name的值>> f = getfield(mystr, {2,1}, 'name')f =gertrude%示例5：删除结构数组中的指定结构域%定义结构数组s>> s.field1=[1 2 3];s.field2='string';s.field3={1 2 3;4 5 6}; %删除结构域field1>> s=rmfield(s,'field1')s =field2: 'string'field3: {2x3 cell}%删除结构域'field2','field3'>> s=rmfield(s,{'field2','field3'})s =field1: [1 2 3]%示例6：%定义结构数组s>> s.field1=[1 2 3];s.field2='string';s.field3={1 2 3;4 5 6}; >> ss =field1: [1 2 3]field2: 'string'field3: {2x3 cell}%将结构数组转换为元胞数组>> c=struct2cell(s)c =[1x3 double]'string'{2x3 cell }%示例7：>>c = {'birch', 'betula', 65; 'maple', 'acer', 50}c ='birch''betula' [65]'maple''acer' [50]>>fields = {'name', 'genus', 'height'}; %fields包含struct中的结构域名>>s = cell2struct(c, fields, 2); %dim=2表示把c中的各行转换为struct数组s =2x1 struct array with fields:namegenusheight>> s(1)ans =name: 'birch'genus: 'betula'height: 65>> s(2)ans =name: 'maple'genus: 'acer'height: 50>> fields = {'field1', 'field2'};>> s = cell2struct(c, fields, 1); %dim=1表示把c中的各列转换为struct数组>> s(1)ans =field1: 'birch'field2: 'maple'>> s(2)ans =field1: 'betula'field2: 'acer'>> s(3)ans =field1:65field2:502.5矩阵操作%示例1： find函数的使用方法。

MATLAB源程序

1、实验目的（1）学习MATLAB的各种二维绘图；（2）学习MATLAB的三维绘图；（3）M ATLAB的绘图修饰（多种、图形注释、绘图颜色、色图矩阵）。

2、实验内容（1）基本二维绘图1）向量绘图x=0:2*pi/100:2*pi;y1=sin(2*x);y2=cos(2*x);plot(x,y1)plot(x,y2)plot(x,y1,x,y2)plot(x,y1);hold on;plot(x,y2);hold off;plot(x',[y1',y2'])设定颜色与线型plot(x,y1,'c:',x,y2,'wo')多窗口绘图figure(1);plot(x,y1);figure(2);plot(x,y2);子图绘图subplot(221);plot(x,y1);subplot(222);plot(x,y2);subplot(223);plot(x,y1,x,y1+y2);subplot(224);plot(x,y2,x,y1-y2);复变函数绘图w=0.01:0.01:10;G=1./(1+2*w*i);subplot(121);plot(G);subplot(122);plot(real(G),imag(G));插值绘图x=0:2*pi/8:2*pi;y=sin(x); plot(x,y,'o');hold on;xi=0:2*pi/100:2*pi;yi=spline(x,y,xi);plot(xi,yi,'m');反白绘图与绘图背景颜色设定whitebgwhitebg('b')whitebg('k')2）函数绘图fplot('sin',[0 4*pi])f='sin(x)';fplot(f,[0 4*pi])fplot('sin(1/x)',[0.01 0.1],1e-3)fplot('[tan(x),sin(x),cos(x)]',[-2*pi,2*pi,-2*pi,2*pi])3)符号函数快捷绘图syms xf=exp(-0.5*x)*sin(x);ezplot(f,[0,10])（2）多种二维绘图1）半对数绘图（频率特性绘图）w=logspace(-1,1);g=20*log10(1./(1+2*w*i));p=angle(1./(1+2*w*i))*180/pi;subplot(211);semilogx(w,g);grid;subplot(212);semilogx(w,p);grid;2）极坐标绘图t=0:2*pi/180:2*pi;mo=cos(2*t);polar(t,mo);t=0:2*pi/8:2*pi; y=sin(t);bar(t,y);t=0:2*pi/8:2*pi; y=sin(t);stem(t,y);t=0:2*pi/8:2*pi; y=sin(t);stairs(t,y);t=-pi:pi/200:pi;comet(t,tan(sin(t))-sin(tan(t)));（3）图形注释y1=dsolve('D2u+2*Du+10*u=0','Du(0)=1,u(0)=0','x');y2=dsolve('D2u+2*Du+10*u=1','Du(0)=1,u(0)=0','x');ezplot(y1,[0,5]);hold on;ezplot(y2,[0,5]);axis([0,5,-0.1,0.2]);title('二阶系统时间响应');axis([0,5,-0.1,0.2]);title('二阶系统时间响应');xlabel('时间t');xlabel('时间t');title('二阶系统时间响应');xlabel('时间t');ylabel('响应幅值y');gtext('零输入响应');gtext('零状态响应');grid;hold off;（4）三维绘图1）三维线图t=0:pi/50:10*pi;plot3(sin(t),cos(t),t);comet3(sin(t),cos(t),t);Z2=[1 1;1 -1];Z4=[Z2 Z2;Z2 -Z2]; Z8=[Z4 Z4;Z4 -Z4]; mesh(Z8);x=-4:0.5:4;y=x; [X,Y]=meshgrid(x,y); Z=X.^2-Y.^2;mesh(X,Y,Z);4)圆锥面网线图t1=0:0.1:0.9;t2=1:0.1:2;r=[t1 -t2+2];[X,Y,Z]=cylinder(r,40);mesh(X,Y,Z);5）视角修饰t1=0:0.1:0.9;t2=1:0.1:2;r=[t1 -t2+2];[X,Y,Z]=cylinder(r,30);mesh(X,Y,Z);subplot(2,2,1);mesh(X,Y,Z);view(0,0);subplot(2,2,2);mesh(X,Y,Z);view(-20,20);subplot(2,2,3);mesh(X,Y,Z);view(-30,30);subplot(2,2,4);mesh(X,Y,Z);view(-40,40);6）MATLAB暖色（hot）函数色图peaks(20);z = 3*(1-x).^2.*exp(-(x.^2) - (y+1).^2) ...- 10*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2) ...- 1/3*exp(-(x+1).^2 - y.^2)axis('off')colormap(hot);colorbar('horiz');7）光照修饰surfl(peaks(20));colormap(hot);shading interp;8）表面修饰subplot(131);peaks(20);shading flat;z = 3*(1-x).^2.*exp(-(x.^2) - (y+1).^2) ...- 10*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2) ...- 1/3*exp(-(x+1).^2 - y.^2)subplot(132);peaks(20);shading interp;z = 3*(1-x).^2.*exp(-(x.^2) - (y+1).^2) ... - 10*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2) ... - 1/3*exp(-(x+1).^2 - y.^2)subplot(133);peaks(20);shading faceted;z = 3*(1-x).^2.*exp(-(x.^2) - (y+1).^2) ... - 10*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2) ... - 1/3*exp(-(x+1).^2 - y.^2)9）透视与消隐P=peaks(30);subplot(121);mesh(P);hidden off;subplot(122);mesh(P);hidden on;10）裁剪修饰P=peaks(30);P(20:23,9:15)=NaN*ones(4,7);subplot(121);meshz(P);P=peaks(30);subplot(121);mesh(P);hidden off;subplot(122);mesh(P);hidden on;P=peaks(30);P(20:23,9:15)=NaN*ones(4,7);subplot(121);meshz(P);subplot(122);meshc(P);11）水线修饰waterfall(peaks(20))colormap([1 0 0])12)等高线修饰contour(peaks(20),6)contour3(peaks(20),10)clabel( contour(peaks(20),4))clabel( contour3(peaks(20),3))。

以下是matlab源程序

以下是matlab源程序=========================================== %%%%%%%%%%%%%%人造地震波程序%%%%%%%%%%%%%% 本程序利用反应谱转换为功率谱方法模拟人造地震波% 变量说明：% s0: 谱强度因子dt: 地震波采样时间间隔% kao: 地震波总时间长度% Ts: 平稳地震动的持续时间% t1,t2,C: 控制主震段首末时间和衰减快慢的参数% keceg,omigag: 地表覆盖层的阻尼比和卓越频率% omigah: 反映基岩特性的谱参数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%xg(1:1500,1)=zeros(1500,1);dt=0.02;kao=21.9;t1=1.5;t2=15;c=0.18;Domiga=2*pi*(25-1/6)/100;Tg=0.25;kece=0.05;alfa=0.9;for i=1:100omiga=(2*pi/6)+Domiga*(i-0.5);T=2*pi/omiga;if T<0.1r=5.5*alfa*T+0.45*alfa;elseif (T>=0.1&T<=Tg)r=alfa;elser=alfa*(Tg/T)^0.9;endgama=omiga/pi;M1=(t2-0.3*t1)/kao;M2=(exp(-2*kece*omiga*t2)-exp(-2*kece*omiga*t1))/(2*kece*omiga*kao);M3=(1-exp(-c*(kao+0.5*t1-t2)))/(c*kao);M=M1+M2+M3;Gomiga=4*kece*r^2/(pi*omiga*M*(sqrt(2*log(gama*kao))+0.5772/sqrt(2*log(gama*kao)))^2); ak=2*sqrt(Gomiga*Domiga);qrand=rand(1)*2*pi;for j=1:1500xg(j)=xg(j)+ak*sin(omiga*j*dt+qrand);endendfor i=1:1500if i*dt<t1fai=(i*dt/t1)^2;elseif (i*dt>=t1&i*dt<t2)fai=1;elsefai=exp(-c*(i*dt-t2));endxg(i)=fai*xg(i);endsubplot(2,1,1);plot(xg);%xg=xg./(max(abs(xg)));save ren xg -ascii;fans=fft(xg);fans=fans.*conj(fans);subplot(2,1,2);plot(fans);。

Matlab源程序代码

正弦波的源程序：(一)，用到的函数1,f2t函数function x=f2t(X)global dt df t f T N%x=f2t(X)%x为时域的取样值矢量%X为x的傅氏变换%X与x长度相同并为2的整幂%本函数需要一个全局变量dt(时域取样间隔) X=[X(N/2+1:N),X(1:N/2)];x=ifft(X)/dt;end2,t2f函数。

function X=t2f(x)global dt df N t f T%X=t2f(x)%x为时域的取样值矢量%X为x的傅氏变换%X与x长度相同,并为2的整幂。

%本函数需要一个全局变量dt(时域取样间隔) H=fft(x);X=[H(N/2+1:N),H(1:N/2)]*dt;end(二),主程序。

1，%（1）绘出正弦信号波形及频谱global dt df t f Nclose allk=input('取样点数＝2^k, k取10左右');if isempty(k), k=10; endf0=input('f0=取1(kz)左右');if isempty(f0), f0=1; endN=2^k;dt=0.01; %msdf=1/(N*dt); %KHzT=N*dt; %截短时间Bs=N*df/2; %系统带宽f=[-Bs+df/2:df:Bs]; %频域横坐标t=[-T/2+dt/2:dt:T/2]; %时域横坐标s=sin(2*pi*f0*t); %输入的正弦信号S=t2f(s); %S是s的傅氏变换a=f2t(S); %a是S的傅氏反变换a=real(a);as=abs(S);subplot(2,1,1) %输出的频谱plot(f,as,'b');gridaxis([-2*f0,+2*f0,min(as),max(as)])xlabel('f (KHz)')ylabel('|S(f)| (V/KHz)') %figure(2)subplot(2,1,2)plot(t,a,'black') %输出信号波形画图gridaxis([-2/f0,+2/f0,-1.5,1.5])xlabel('t(ms)')ylabel('a(t)(V)')gtext('频谱图')最佳基带系统的源程序：(一)，用到的函数f2t函数和t2f函数。

高等数值分析第二章答案

第二章习题参考答案1.解：由于20Ax b−≥，极小化2b Ax −与极小化22Ax b −是等价的。

令22()(,)(,)2(,)x Ax b Ax Ax b b Ax b ϕ=−=+−，对于任意的n R y x ∈,和实数α，)()(),()()(,*222*2****x Ay a x Ay Ay a x ay x b Ax x ϕϕϕϕ≥+=+=+=则有满足若这表示处达到极小值。

在*)(x x ϕ反之，若必有处达到极小，则对任意在nR y x ay x ∈+*)(ϕ0),(2),(2),(20)(**0*=−=+−=+=Ay b Ax Ay Ay a Ay b Ax daay x d a 即ϕ故有 b Ax =*成立。

以上证明了求解,22b Ax b Ax −=等价于极小化即。

等价于极小化2b Ax b Ax −= 推导最速下降法过程如下：),/(),(0),(),(,0),,2)(222)()(11k T k T k T k k T k T k T k k T k k k T k k kT k T k T T x x k r AA r AA r AA r a r AA r AA a r AA r r aA x da dx a r aA x x r A Ax b A Ax A b A x grad x x k==+−=++==−=−=−++=最终得到得出（由取得极小值。

使求出取的负梯度方向，且下降最快的方向是该点在ϕϕϕ给出的算法如下：1）)(000Ax b A r A R x T T n −=∈，计算给定； 2）L ,2,1,0=k 对于）转到否则数。

为一事先给定的停机常则停止；其中若2),/(),(10,11kT k k k k T k k k k k k k k k r A p Ax b r r A a x x Ap Ap p p a k k r =−=+==+=>≤−−εε2.证明 1）正定性由对称正定矩阵的性质，(),0x Ax ≥（当且仅当x ＝0时取等号），所以 ()12,0Axx Ax =≥（当且仅当x ＝0时取等号）2）齐次性()()()121122,(),,AA xx A x x Ax x Ax x αααααα⎡⎤====⎣⎦3）o1方法（一）A 是对称正定矩阵，得到(,())0x y A x y λλ++≥，把它展开如下2(,)(,)(,)(,)0y Ay x Ay y Ax x Ax λλλ+++≥考虑到(,)(,)(,)x Ay Ax y y Ax ==，把上式看成关于λ的一元二次方程，则式子等价于24(,)4(,)(,)0x Ay x Ax y Ay ∆=−≤因此1/21/2(,)(,)(,)x Ay x Ax y Ay ≤所以1/21/221/21/2((,)(,))(,)(,)2(,)(,)(,)(,)2(,)(,)(,)(,)(,)((),())x Ax y Ay x Ax y Ay x Ax y Ay x Ax y Ay x Ay x Ax y Ay x Ay y Ax x y A x y +=++≥++=+++=++两边开平方即可得到AA A x yx y +≤+因此，1/2(,)A x Ax x =是一种向量范数。

GMRES_matlab

gmresGeneralized minimum residual method (with restarts)expand all in page Syntaxx = gmres(A,b)gmres(A,b,restart)gmres(A,b,restart,tol)gmres(A,b,restart,tol,maxit)gmres(A,b,restart,tol,maxit,M)gmres(A,b,restart,tol,maxit,M1,M2)gmres(A,b,restart,tol,maxit,M1,M2,x0)[x,flag] = gmres(A,b,...)[x,flag,relres] = gmres(A,b,...)[x,flag,relres,iter] = gmres(A,b,...)[x,flag,relres,iter,resvec] = gmres(A,b,...)Descriptionx = gmres(A,b) attempts to solve the system of linear equations A*x = b for x. The n-by-n coefficient matrix A must be square and should be large and sparse. The column vector b must have length n. A can be a function handle, afun, such that afun(x) returns A*x. For this syntax, gmres does not restart; the maximum number of iterations is min(n,10).Parameterizing Functions explains how to provide additional parameters to the function afun, as well as the preconditioner function mfun described below, if necessary.If gmres converges, a message to that effect is displayed. If gmres fails to converge after the maximum number of iterations or halts for any reason, a warning message is printed displaying the relative residual norm(b-A*x)/norm(b) and the iteration number at which the method stopped or failed. gmres(A,b,restart) restarts the method every restart inner iterations. The maximum number of outer iterations is min(n/restart,10). The maximum number of total iterationsis restart*min(n/restart,10). If restart is n or [], then gmres does not restart and the maximum number of total iterations is min(n,10).gmres(A,b,restart,tol) specifies the tolerance of the method. If tol is [], then gmres uses the default, 1e-6.gmres(A,b,restart,tol,maxit) specifies the maximum number of outer iterations, i.e., the total number of iterations does not exceed restart*maxit. If maxit is [] then gmres uses the default, min(n/restart,10). If restart is n or [], then the maximum number of total iterationsis maxit (instead of restart*maxit).gmres(A,b,restart,tol,maxit,M) and gmres(A,b,restart,tol,maxit,M1,M2) use preconditioner M or M = M1*M2and effectively solve the system inv(M)*A*x = inv(M)*b for x.If M is [] then gmres applies no preconditioner. M can be a function handle mfun suchthat mfun(x) returns M\x.gmres(A,b,restart,tol,maxit,M1,M2,x0) specifies the first initial guess. If x0 is [],then gmres uses the default, an all-zero vector.[x,flag] = gmres(A,b,...) also returns a convergence flag:flag = 0gmres converged to the desired tolerance tol within maxit outer iterations.flag = 1gmres iterated maxit times but did not converge.flag = 2Preconditioner M was ill-conditioned.flag = 3gmres stagnated. (Two consecutive iterates were the same.)Whenever flag is not 0, the solution x returned is that with minimal norm residual computed over all the iterations. No messages are displayed if the flag output is specified.[x,flag,relres] = gmres(A,b,...) also returns the relative residual norm(b-A*x)/norm(b).If flag is 0, relres <= tol. The third output, relres, is the relative residual of the preconditioned system.[x,flag,relres,iter] = gmres(A,b,...) also returns both the outer and inner iterationnumbers at which x was computed, where 0 <= iter(1) <= maxit and 0 <= iter(2) <= restart.[x,flag,relres,iter,resvec] = gmres(A,b,...) also returns a vector of the residual norms at each inner iteration. These are the residual norms for the preconditioned system.ExamplesUsing gmres with a Matrix InputA = gallery('wilk',21);b = sum(A,2);tol = 1e-12;maxit = 15;M1 = diag([10:-1:1 1 1:10]);x = gmres(A,b,10,tol,maxit,M1);displays the following message:gmres(10) converged at outer iteration 2 (inner iteration 9) toa solution with relative residual 3.3e-013Using gmres with a Function HandleThis example replaces the matrix A in the previous example with a handle to a matrix-vector productfunction afun, and the preconditioner M1 with a handle to a backsolve function mfun. The example iscontained in a function run_gmres that∙Calls gmres with the function handle @afun as its first argument.∙Contains afun and mfun as nested functions, so that all variables in run_gmres are available to afun and mfun.The following shows the code for run_gmres:function x1 = run_gmresn = 21;b = afun(ones(n,1));tol = 1e-12; maxit = 15;x1 = gmres(@afun,b,10,tol,maxit,@mfun);function y = afun(x)y = [0; x(1:n-1)] + ...[((n-1)/2:-1:0)'; (1:(n-1)/2)'].*x + ...[x(2:n); 0];endfunction y = mfun(r)y = r ./ [((n-1)/2:-1:1)'; 1; (1:(n-1)/2)'];endendWhen you enterx1 = run_gmres;MATLAB® software displays the messagegmres(10) converged at outer iteration 2 (inner iteration 10) to a solution with relative residual 1.1e-013.Using a Preconditioner without RestartThis example demonstrates the use of a preconditioner without restarting gmres.Load west0479, a real 479-by-479 nonsymmetric sparse matrix.load west0479;A = west0479;Set the tolerance and maximum number of iterations.tol = 1e-12; maxit = 20;Define b so that the true solution is a vector of all ones.b = full(sum(A,2));[x0,fl0,rr0,it0,rv0] = gmres(A,b,[],tol,maxit);fl0 is 1 because gmres does not converge to the requested tolerance 1e-12 within the requested 20 iterations. The best approximate solution that gmres returns is the last one (as indicated by it0(2) = 20). MATLAB stores the residual history in rv0.Plot the behavior of gmres.semilogy(0:maxit,rv0/norm(b),'-o');xlabel('Iteration number');ylabel('Relative residual');The plot shows that the solution converges slowly. A preconditioner may improve the outcome.Use ilu to form the preconditioner, since A is nonsymmetric.[L,U] = ilu(A,struct('type','ilutp','droptol',1e-5));Error using iluThere is a pivot equal to zero. Consider decreasingthe drop tolerance or consider using the 'udiag' option.Note MATLAB cannot construct the incomplete LU as it would result in a singular factor, which is useless as a preconditioner.As indicated by the error message, try again with a reduced drop tolerance.[L,U] = ilu(A,struct('type','ilutp','droptol',1e-6));[x1,fl1,rr1,it1,rv1] = gmres(A,b,[],tol,maxit,L,U);fl1 is 0 because gmres drives the relative residual to 9.5436e-14 (the value of rr1). The relative residual is less than the prescribed tolerance of 1e-12 at the sixth iteration (the value of it1(2)) when preconditioned by the incomplete LU factorization with a drop tolerance of 1e-6. The output, rv1(1), is norm(M\b), where M = L*U. The output, rv1(7), is norm(U\(L\(b-A*x1))).Follow the progress of gmres by plotting the relative residuals at each iteration starting from the initial estimate (iterate number 0).semilogy(0:it1(2),rv1/norm(b),'-o');xlabel('Iteration number');ylabel('Relative residual');Using a Preconditioner with RestartThis example demonstrates the use of a preconditioner with restarted gmres.Load west0479, a real 479-by-479 nonsymmetric sparse matrix.load west0479;A = west0479;Define b so that the true solution is a vector of all ones.b = full(sum(A,2));Construct an incomplete LU preconditioner as in the previous example.[L,U] = ilu(A,struct('type','ilutp','droptol',1e-6));The benefit to using restarted gmres is to limit the amount of memory required to execute the method. Without restart, gmres requires maxit vectors of storage to keep the basis of the Krylov subspace. Also,gmres must orthogonalize against all of the previous vectors at each step. Restarting limits the amount of workspace used and the amount of work done per outer iteration. Note that even though preconditioned gmres converged in six iterations above, the algorithm allowed for as many as twenty basis vectors and therefore, allocated all of that space up front.Execute gmres(3), gmres(4), and 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, and fl5 are all 0 because in each case restarted gmres drives the relative residual to less than the prescribed tolerance of 1e-12.The following plots show the convergence histories of eachrestarted gmres method. gmres(3) converges at outer iteration 5, inner iteration 3 (it3 = [5, 3]) which would be the same as outer iteration 6, inner iteration 0, hence the marking of 6 on the final tick mark.figuresemilogy(1:1/3:6,rv3/norm(b),'-o');set(gca,'XTick',[1:1/3:6],...'XTickLabel',...['1';' ';' ';'2';' ';' ';'3';' ';' ';'4';' ';' ';'5';' ';' ';'6';])title('gmres(3)')xlabel('Iteration number');ylabel('Relative residual');figuresemilogy(1:1/4:3,rv4/norm(b),'-o');set(gca,'XTick',[1:1/4:3],...'XTickLabel',['1';' ';' ';' ';'2';' ';' ';' ';'3'])title('gmres(4)')xlabel('Iteration number');ylabel('Relative residual');figuresemilogy(1:1/5:2.8,rv5/norm(b),'-o');set(gca,'XTick',[1:1/5:2.8],...'XTickLabel',['1';' ';' ';' ';' ';'2';' ';' ';' ';' ']) title('gmres(5)')xlabel('Iteration number');ylabel('Relative residual');ReferencesBarrett, R., M. Berry, T. F. Chan, et al., Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, SIAM, Philadelphia, 1994.Saad, Youcef and Martin H. Schultz, "GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems," SIAM J. Sci. Stat. Comput., July 1986, Vol. 7, No. 3, pp. 856-869.。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

f u n c t i o n[x,f l a g,r e l r e s,i t e r,r e s v e c]=g m r e s(A,b,r e s t a r t,t o l,m a x i t,M1,M2,x,v a r a r g i n)%GMRES Generalized Minimum Residual Method.% X = GMRES(A,B) attempts to solve the system of linear equations A*X = B% for X. The N-by-N coefficient matrix A must be square and the right% hand side column vector B must have length N. This uses the unrestarted% method with MIN(N,10) total iterations.%% X = GMRES(AFUN,B) accepts a function handle AFUN instead of the matrix% A. AFUN(X) accepts a vector input X and returns the matrix-vector% product A*X. In all of the following syntaxes, you can replace A by% AFUN.%% X = GMRES(A,B,RESTART) restarts the method every RESTART iterations.% If RESTART is N or [] then GMRES uses the unrestarted method as above.%% X = GMRES(A,B,RESTART,TOL) specifies the tolerance of the method. If% TOL is [] then GMRES uses the default, 1e-6.%% X = GMRES(A,B,RESTART,TOL,MAXIT) specifies the maximum number of outer% iterations. Note: the total number of iterations is RESTART*MAXIT. If% MAXIT is [] then GMRES uses the default, MIN(N/RESTART,10). If RESTART% is N or [] then the total number of iterations is MAXIT.%% X = GMRES(A,B,RESTART,TOL,MAXIT,M) and% X = GMRES(A,B,RESTART,TOL,MAXIT,M1,M2) use preconditioner M or M=M1*M2% and effectively solve the system inv(M)*A*X = inv(M)*B for X. If M is% [] then a preconditioner is not applied. M may be a function handle% returning M\X.%% X = GMRES(A,B,RESTART,TOL,MAXIT,M1,M2,X0) specifies the first initial% guess. If X0 is [] then GMRES uses the default, an all zero vector.%% [X,FLAG] = GMRES(A,B,...) also returns a convergence FLAG:% 0 GMRES converged to the desired tolerance TOL within MAXIT iterations.% 1 GMRES iterated MAXIT times but did not converge.% 2 preconditioner M was ill-conditioned.% 3 GMRES stagnated (two consecutive iterates were the same).%% [X,FLAG,RELRES] = GMRES(A,B,...) also returns the relative residual% NORM(B-A*X)/NORM(B). If FLAG is 0, then RELRES <= TOL. Note with% preconditioners M1,M2, the residual is NORM(M2\(M1\(B-A*X))).%% [X,FLAG,RELRES,ITER] = GMRES(A,B,...) also returns both the outer and% inner iteration numbers at which X was computed: 0 <= ITER(1) <= MAXIT% and 0 <= ITER(2) <= RESTART.%% [X,FLAG,RELRES,ITER,RESVEC] = GMRES(A,B,...) also returns a vector of% the residual norms at each inner iteration, including NORM(B-A*X0).% Note with preconditioners M1,M2, the residual is NORM(M2\(M1\(B-A*X))).%% Example:% n = 21; A = gallery('wilk',n); b = sum(A,2);% tol = 1e-12; maxit = 15; M = diag([10:-1:1 1 1:10]);% x = gmres(A,b,10,tol,maxit,M);% Or, use this matrix-vector product function% %-----------------------------------------------------------------%% function y = afun(x,n)% y = [0; x(1:n-1)] + [((n-1)/2:-1:0)'; (1:(n-1)/2)'].*x+[x(2:n); 0];% %-----------------------------------------------------------------%% and this preconditioner backsolve function% %------------------------------------------%% function y = mfun(r,n)% y = r ./ [((n-1)/2:-1:1)'; 1; (1:(n-1)/2)'];% %------------------------------------------%% as inputs to GMRES:% x1 = gmres(@(x)afun(x,n),b,10,tol,maxit,@(x)mfun(x,n));%% Class support for inputs A,B,M1,M2,X0 and the output of AFUN:% float: double%% See also BICG, BICGSTAB, BICGSTABL, CGS, LSQR, MINRES, PCG, QMR, SYMMLQ, % TFQMR, ILU, FUNCTION_HANDLE.% References% H.F. Walker, "Implementation of the GMRES Method Using Householder% Transformations", SIAM J. Sci. Comp. Vol 9. No 1. January 1988.% Copyright 1984-2010 The MathWorks, Inc.if (nargin < 2)error('MATLAB:gmres:NumInputs','Not enough input arguments.');end% Determine whether A is a matrix or a function.[atype,afun,afcnstr] = iterchk(A);if strcmp(atype,'matrix')% Check matrix and right hand side vector inputs have appropriate sizes[m,n] = size(A);if (m ~= n)error('MATLAB:gmres:SquareMatrix','Matrix must be square.');endif ~isequal(size(b),[m,1])error('MATLAB:gmres:VectorSize', '%s %d %s', ...'Right hand side must be a column vector of length', m,...'to match the coefficient matrix.');endn = m;if ~iscolumn(b)error('MATLAB:gmres:Vector','Right hand side must be a column vector.');endend% Assign default values to unspecified parametersif (nargin < 3) || isempty(restart) || (restart == n)restarted = false;elserestarted = true;endif (nargin < 4) || isempty(tol)tol = 1e-6;endwarned = 0;if tol < epswarning('MATLAB:gmres:tooSmallTolerance', ...strcat('Input tol is smaller than eps and may not be achieved',...' by GMRES\n',' Try to use a bigger tolerance'));warned = 1;tol = eps;elseif tol >= 1warning('MATLAB:gmres:tooBigTolerance', ...strcat('Input tol is bigger than 1 \n',...' Try to use a smaller tolerance'));warned = 1;tol = 1-eps;endif (nargin < 5) || isempty(maxit)if restartedmaxit = min(ceil(n/restart),10);elsemaxit = min(n,10);endendif restartedouter = maxit;if restart > nwarning('MATLAB:gmres:tooManyInnerIts', 'RESTART is %d %s\n%s %d.', ...restart, 'but it should be bounded by SIZE(A,1).', ...' Setting RESTART to', n);restart = n;endinner = restart;if maxit > nwarning('MATLAB:gmres:tooManyInnerIts', 'MAXIT is %d %s\n%s %d.', ...maxit, 'but it should be bounded by SIZE(A,1).', ...' Setting MAXIT to', n);maxit = n;endinner = maxit;end% Check for all zero right hand side vector => all zero solutionn2b = norm(b); % Norm of rhs vector, bif (n2b == 0) % if rhs vector is all zerosx = zeros(n,1); % then solution is all zerosflag = 0; % a valid solution has been obtainedrelres = 0; % the relative residual is actually 0/0iter = [0 0]; % no iterations need be performedresvec = 0; % resvec(1) = norm(b-A*x) = norm(0)if (nargout < 2)itermsg('gmres',tol,maxit,0,flag,iter,NaN);endreturnendif ((nargin >= 6) && ~isempty(M1))existM1 = 1;[m1type,m1fun,m1fcnstr] = iterchk(M1);if strcmp(m1type,'matrix')if ~isequal(size(M1),[m,m])error('MATLAB:gmres:PreConditioner1Size', '%s %d %s',...'Preconditioner must be a square matrix of size', m, ...'to match the problem size.');endendelseexistM1 = 0;m1type = 'matrix';endif ((nargin >= 7) && ~isempty(M2))existM2 = 1;[m2type,m2fun,m2fcnstr] = iterchk(M2);if strcmp(m2type,'matrix')if ~isequal(size(M2),[m,m])error('MATLAB:gmres:PreConditioner2Size', '%s %d %s', ...'Preconditioner must be a square matrix of size', m, ...'to match the problem size.');endendelseexistM2 = 0;m2type = 'matrix';endif ((nargin >= 8) && ~isempty(x))if ~isequal(size(x),[n,1])error('MATLAB:gmres:XoSize', '%s %d %s', ...'Initial guess must be a column vector of length', n, ...'to match the problem size.');endelsex = zeros(n,1);endif ((nargin > 8) && strcmp(atype,'matrix') && ...strcmp(m1type,'matrix') && strcmp(m2type,'matrix'))error('MATLAB:gmres:TooManyInputs','Too many input arguments.');end% Set up for the methodflag = 1;xmin = x; % Iterate which has minimal residual so farimin = 0; % "Outer" iteration at which xmin was computed jmin = 0; % "Inner" iteration at which xmin was computed tolb = tol * n2b; % Relative toleranceevalxm = 0;stag = 0;moresteps = 0;maxmsteps = min([floor(n/50),5,n-maxit]);maxstagsteps = 3;minupdated = 0;x0iszero = (norm(x) == 0);r = b - iterapp('mtimes',afun,atype,afcnstr,x,varargin{:});normr = norm(r); % Norm of initial residualif (normr <= tolb) % Initial guess is a good enough solution flag = 0;relres = normr / n2b;iter = [0 0];resvec = normr;if (nargout < 2)itermsg('gmres',tol,maxit,[0 0],flag,iter,relres);endreturnendminv_b = b;if existM1r = iterapp('mldivide',m1fun,m1type,m1fcnstr,r,varargin{:});if ~x0iszerominv_b = iterapp('mldivide',m1fun,m1type,m1fcnstr,b,varargin{:});elseminv_b = r;endif ~all(isfinite(r)) || ~all(isfinite(minv_b))flag = 2;x = xmin;relres = normr / n2b;iter = [0 0];resvec = normr;returnendendif existM2r = iterapp('mldivide',m2fun,m2type,m2fcnstr,r,varargin{:});if ~x0iszerominv_b = iterapp('mldivide',m2fun,m2type,m2fcnstr,minv_b,varargin{:});elseminv_b = r;endif ~all(isfinite(r)) || ~all(isfinite(minv_b))flag = 2;x = xmin;relres = normr / n2b;iter = [0 0];resvec = normr;returnendendnormr = norm(r); % norm of the preconditioned residualn2minv_b = norm(minv_b); % norm of the preconditioned rhsclear minv_b;tolb = tol * n2minv_b;if (normr <= tolb) % Initial guess is a good enough solution flag = 0;relres = normr / n2minv_b;iter = [0 0];resvec = n2minv_b;if (nargout < 2)itermsg('gmres',tol,maxit,[0 0],flag,iter,relres);endreturnendresvec = zeros(inner*outer+1,1); % Preallocate vector for norm of residuals resvec(1) = normr; % resvec(1) = norm(b-A*x0)normrmin = normr; % Norm of residual from xmin% Preallocate J to hold the Given's rotation constants.J = zeros(2,inner);U = zeros(n,inner);R = zeros(inner,inner);w = zeros(inner+1,1);for outiter = 1 : outer% Construct u for Householder reflector.% u = r + sign(r(1))*||r||*e1u = r;normr = norm(r);beta = scalarsign(r(1))*normr;u(1) = u(1) + beta;u = u / norm(u);U(:,1) = u;% Apply Householder projection to r.% w = r - 2*u*u'*r;w(1) = -beta;for initer = 1 : inner% Form P1*P2*P3...Pj*ej.% v = Pj*ej = ej - 2*u*u'*ejv = -2*(u(initer)')*u;v(initer) = v(initer) + 1;% v = P1*P2*...Pjm1*(Pj*ej)for k = (initer-1):-1:1v = v - U(:,k)*(2*(U(:,k)'*v));end% Explicitly normalize v to reduce the effects of round-off.v = v/norm(v);% Apply A to v.v = iterapp('mtimes',afun,atype,afcnstr,v,varargin{:});% Apply Preconditioner.if existM1v = iterapp('mldivide',m1fun,m1type,m1fcnstr,v,varargin{:});if ~all(isfinite(v))flag = 2;breakendendif existM2v = iterapp('mldivide',m2fun,m2type,m2fcnstr,v,varargin{:});if ~all(isfinite(v))flag = 2;breakendend% Form Pj*Pj-1*...P1*Av.for k = 1:initerv = v - U(:,k)*(2*(U(:,k)'*v));end% Determine Pj+1.if (initer ~= length(v))% Construct u for Householder reflector Pj+1.u = [zeros(initer,1); v(initer+1:end)];alpha = norm(u);if (alpha ~= 0)alpha = scalarsign(v(initer+1))*alpha;% u = v(initer+1:end) +% sign(v(initer+1))*||v(initer+1:end)||*e_{initer+1)u(initer+1) = u(initer+1) + alpha;u = u / norm(u);U(:,initer+1) = u;% Apply Pj+1 to v.% v = v - 2*u*(u'*v);v(initer+2:end) = 0;v(initer+1) = -alpha;endend% Apply Given's rotations to the newly formed v.for colJ = 1:initer-1tmpv = v(colJ);v(colJ) = conj(J(1,colJ))*v(colJ) + conj(J(2,colJ))*v(colJ+1);v(colJ+1) = -J(2,colJ)*tmpv + J(1,colJ)*v(colJ+1);end% Compute Given's rotation Jm.if ~(initer==length(v))rho = norm(v(initer:initer+1));J(:,initer) = v(initer:initer+1)./rho;w(initer+1) = -J(2,initer).*w(initer);w(initer) = conj(J(1,initer)).*w(initer);v(initer) = rho;v(initer+1) = 0;endR(:,initer) = v(1:inner);normr = abs(w(initer+1));resvec((outiter-1)*inner+initer+1) = normr;normr_act = normr;if (normr <= tolb || stag >= maxstagsteps || moresteps)if evalxm == 0ytmp = R(1:initer,1:initer) \ w(1:initer);additive = U(:,initer)*(-2*ytmp(initer)*conj(U(initer,initer)));additive(initer) = additive(initer) + ytmp(initer);for k = initer-1 : -1 : 1additive(k) = additive(k) + ytmp(k);additive = additive - U(:,k)*(2*(U(:,k)'*additive));endif norm(additive) < eps*norm(x)stag = stag + 1;elsestag = 0;endxm = x + additive;evalxm = 1;elseif evalxm == 1addvc = [-(R(1:initer-1,1:initer-1)\R(1:initer-1,initer))*...(w(initer)/R(initer,initer)); w(initer)/R(initer,initer)];if norm(addvc) < eps*norm(xm)stag = stag + 1;elsestag = 0;endadditive = U(:,initer)*(-2*addvc(initer)*conj(U(initer,initer)));additive(initer) = additive(initer) + addvc(initer);for k = initer-1 : -1 : 1additive(k) = additive(k) + addvc(k);additive = additive - U(:,k)*(2*(U(:,k)'*additive));endxm = xm + additive;endr = b - iterapp('mtimes',afun,atype,afcnstr,xm,varargin{:});if norm(r) <= tol*n2bx = xm;flag = 0;iter = [outiter, initer];breakendminv_r = r;if existM1minv_r = iterapp('mldivide',m1fun,m1type,m1fcnstr,r,varargin{:});if ~all(isfinite(minv_r))flag = 2;breakendendif existM2minv_r = iterapp('mldivide',m2fun,m2type,m2fcnstr,minv_r,varargin{:});if ~all(isfinite(minv_r))flag = 2;breakendendnormr_act = norm(minv_r);resvec((outiter-1)*inner+initer+1) = normr_act;if normr_act <= normrminnormrmin = normr_act;imin = outiter;jmin = initer;xmin = xm;minupdated = 1;endif normr_act <= tolbx = xm;flag = 0;iter = [outiter, initer];breakelseif stag >= maxstagsteps && moresteps == 0stag = 0;endmoresteps = moresteps + 1;if moresteps >= maxmstepsif ~warnedwarning('MATLAB:gmres:tooSmallTolerance', ...strcat('Input tol might be obviously smaller than',...' eps*cond(A) and may not be achieved by GMRES\n',...' Try to use a bigger tolerance'));endflag = 3;iter = [outiter, initer];break;endendendif normr_act <= normrminnormrmin = normr_act;imin = outiter;jmin = initer;minupdated = 1;endif stag >= maxstagstepsflag = 3;break;endend % ends inner loopevalxm = 0;if flag ~= 0if minupdatedidx = jmin;elseidx = initer;endy = R(1:idx,1:idx) \ w(1:idx);additive = U(:,idx)*(-2*y(idx)*conj(U(idx,idx)));additive(idx) = additive(idx) + y(idx);for k = idx-1 : -1 : 1additive(k) = additive(k) + y(k);additive = additive - U(:,k)*(2*(U(:,k)'*additive));endx = x + additive;xmin = x;r = b - iterapp('mtimes',afun,atype,afcnstr,x,varargin{:});minv_r = r;if existM1minv_r = iterapp('mldivide',m1fun,m1type,m1fcnstr,r,varargin{:});if ~all(isfinite(minv_r))flag = 2;breakendendif existM2minv_r = iterapp('mldivide',m2fun,m2type,m2fcnstr,minv_r,varargin{:});if ~all(isfinite(minv_r))flag = 2;breakendendnormr_act = norm(minv_r);r = minv_r;endif normr_act <= normrminxmin = x;normrmin = normr_act;imin = outiter;jmin = initer;endif flag == 3break;endif normr_act <= tolbflag = 0;iter = [outiter, initer];break;endminupdated = 0;end % ends outer loop% returned solution is that with minimum residualif flag == 0relres = normr_act / n2minv_b;elsex = xmin;iter = [imin jmin];relres = normr_act / n2minv_b;end% truncate the zeros from resvecif flag <= 1 || flag == 3resvec = resvec(1:(outiter-1)*inner+initer+1);indices = resvec==0;resvec = resvec(~indices);elseif initer == 0resvec = resvec(1:(outiter-1)*inner+1);elseresvec = resvec(1:(outiter-1)*inner+initer);endend% only display a message if the output flag is not usedif nargout < 2if restarteditermsg(sprintf('gmres(%d)',restart),tol,maxit,[outiter initer],flag,iter,relres);elseitermsg(sprintf('gmres'),tol,maxit,initer,flag,iter(2),relres);endendfunction sgn = scalarsign(d)sgn = sign(d);if (sgn == 0)sgn = 1;end。