gmres源程序matlab.docx

gmres源程序matlab.docx

《gmres源程序matlab.docx》由会员分享，可在线阅读，更多相关《gmres源程序matlab.docx（23页珍藏版）》请在冰豆网上搜索。

gmres源程序matlab

function[x,flag,relres,iter,resvec]=gmres（A,b,restart,tol,maxit,M1,M2,x,varargin）

%GMRESGeneralizedMinimumResidualMethod.

%X=GMRES（A,B）attemptstosolvethesystemoflinearequationsA*X=B

%forX.TheN-by-NcoefficientmatrixAmustbesquareandtheright

%handsidecolumnvectorBmusthavelengthN.Thisusestheunrestarted

%methodwithMIN（N,10）totaliterations.

%

%X=GMRES（AFUN,B）acceptsafunctionhandleAFUNinsteadofthematrix

%A.AFUN（X）acceptsavectorinputXandreturnsthematrix-vector

%productA*X.Inallofthefollowingsyntaxes,youcanreplaceAby

%AFUN.

%

%X=GMRES（A,B,RESTART）restartsthemethodeveryRESTARTiterations.

%IfRESTARTisNor[]thenGMRESusestheunrestartedmethodasabove.

%

%X=GMRES（A,B,RESTART,TOL）specifiesthetoleranceofthemethod.If

%TOLis[]thenGMRESusesthedefault,1e-6.

%

%X=GMRES（A,B,RESTART,TOL,MAXIT）specifiesthemaximumnumberofouter

%iterations.Note:

thetotalnumberofiterationsisRESTART*MAXIT.If

%MAXITis[]thenGMRESusesthedefault,MIN（N/RESTART,10）.IfRESTART

%isNor[]thenthetotalnumberofiterationsisMAXIT.

%

%X=GMRES（A,B,RESTART,TOL,MAXIT,M）and

%X=GMRES（A,B,RESTART,TOL,MAXIT,M1,M2）usepreconditionerMorM=M1*M2

%andeffectivelysolvethesysteminv（M）*A*X=inv（M）*BforX.IfMis

%[]thenapreconditionerisnotapplied.Mmaybeafunctionhandle

%returningM\X.

%

%X=GMRES（A,B,RESTART,TOL,MAXIT,M1,M2,X0）specifiesthefirstinitial

%guess.IfX0is[]thenGMRESusesthedefault,anallzerovector.

%

%[X,FLAG]=GMRES（A,B,...）alsoreturnsaconvergenceFLAG:

%0GMRESconvergedtothedesiredtoleranceTOLwithinMAXITiterations.

%1GMRESiteratedMAXITtimesbutdidnotconverge.

%2preconditionerMwasill-conditioned.

%3GMRESstagnated（twoconsecutiveiterateswerethesame）.

%

%[X,FLAG,RELRES]=GMRES（A,B,...）alsoreturnstherelativeresidual

%NORM（B-A*X）/NORM（B）.IfFLAGis0,thenRELRES<=TOL.Notewith

%preconditionersM1,M2,theresidualisNORM（M2\（M1\（B-A*X）））.

%

%[X,FLAG,RELRES,ITER]=GMRES（A,B,...）alsoreturnsboththeouterand

%inneriterationnumbersatwhichXwascomputed:

0<=ITER

（1）<=MAXIT

%and0<=ITER

（2）<=RESTART.

%

%[X,FLAG,RELRES,ITER,RESVEC]=GMRES（A,B,...）alsoreturnsavectorof

%theresidualnormsateachinneriteration,includingNORM（B-A*X0）.

%NotewithpreconditionersM1,M2,theresidualisNORM（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:

111:

10]）;

%x=gmres（A,b,10,tol,maxit,M）;

%Or,usethismatrix-vectorproductfunction

%%-----------------------------------------------------------------%

%functiony=afun（x,n）

%y=[0;x（1:

n-1）]+[（（n-1）/2:

-1:

0）';（1:

（n-1）/2）'].*x+[x（2:

n）;0];

%%-----------------------------------------------------------------%

%andthispreconditionerbacksolvefunction

%%------------------------------------------%

%functiony=mfun（r,n）

%y=r./[（（n-1）/2:

-1:

1）';1;（1:

（n-1）/2）'];

%%------------------------------------------%

%asinputstoGMRES:

%x1=gmres（@（x）afun（x,n）,b,10,tol,maxit,@（x）mfun（x,n））;

%

%ClasssupportforinputsA,B,M1,M2,X0andtheoutputofAFUN:

%float:

double

%

%SeealsoBICG,BICGSTAB,BICGSTABL,CGS,LSQR,MINRES,PCG,QMR,SYMMLQ,

%TFQMR,ILU,FUNCTION_HANDLE.

%References

%H.F.Walker,"ImplementationoftheGMRESMethodUsingHouseholder

%Transformations",SIAMJ.Sci.Comp.Vol9.No1.January1988.

%Copyright1984-2010TheMathWorks,Inc.

%$Revision:

1.21.4.13$$Date:

2010/02/2508:

11:

48$

if（nargin<2）

error（'MATLAB:

gmres:

NumInputs','Notenoughinputarguments.'）;

end

%DeterminewhetherAisamatrixorafunction.

[atype,afun,afcnstr]=iterchk（A）;

ifstrcmp（atype,'matrix'）

%Checkmatrixandrighthandsidevectorinputshaveappropriatesizes

[m,n]=size（A）;

if（m~=n）

error（'MATLAB:

gmres:

SquareMatrix','Matrixmustbesquare.'）;

end

if~isequal（size（b）,[m,1]）

error（'MATLAB:

gmres:

VectorSize','%s%d%s',...

'Righthandsidemustbeacolumnvectoroflength',m,...

'tomatchthecoefficientmatrix.'）;

end

else

m=size（b,1）;

n=m;

if~iscolumn（b）

error（'MATLAB:

gmres:

Vector','Righthandsidemustbeacolumnvector.'）;

end

end

%Assigndefaultvaluestounspecifiedparameters

if（nargin<3）||isempty（restart）||（restart==n）

restarted=false;

else

restarted=true;

end

if（nargin<4）||isempty（tol）

tol=1e-6;

end

warned=0;

iftol

warning（'MATLAB:

gmres:

tooSmallTolerance',...

strcat（'Inputtolissmallerthanepsandmaynotbeachieved',...

'byGMRES\n','Trytouseabiggertolerance'））;

warned=1;

tol=eps;

elseiftol>=1

warning（'MATLAB:

gmres:

tooBigTolerance',...

strcat（'Inputtolisbiggerthan1\n',...

'Trytouseasmallertolerance'））;

warned=1;

tol=1-eps;

end

if（nargin<5）||isempty（maxit）

ifrestarted

maxit=min（ceil（n/restart）,10）;

else

maxit=min（n,10）;

end

end

ifrestarted

outer=maxit;

ifrestart>n

warning（'MATLAB:

gmres:

tooManyInnerIts','RESTARTis%d%s\n%s%d.',...

restart,'butitshouldbeboundedbySIZE（A,1）.',...

'SettingRESTARTto',n）;

restart=n;

end

inner=restart;

else

outer=1;

ifmaxit>n

warning（'MATLAB:

gmres:

tooManyInnerIts','MAXITis%d%s\n%s%d.',...

maxit,'butitshouldbeboundedbySIZE（A,1）.',...

'SettingMAXITto',n）;

maxit=n;

end

inner=maxit;

end

%Checkforallzerorighthandsidevector=>allzerosolution

n2b=norm（b）;%Normofrhsvector,b

if（n2b==0）%ifrhsvectorisallzeros

x=zeros（n,1）;%thensolutionisallzeros

flag=0;%avalidsolutionhasbeenobtained

relres=0;%therelativeresidualisactually0/0

iter=[00];%noiterationsneedbeperformed

resvec=0;%resvec

（1）=norm（b-A*x）=norm（0）

if（nargout<2）

itermsg（'gmres',tol,maxit,0,flag,iter,NaN）;

end

return

end

if（（nargin>=6）&&~isempty（M1））

existM1=1;

[m1type,m1fun,m1fcnstr]=iterchk（M1）;

ifstrcmp（m1type,'matrix'）

if~isequal（size（M1）,[m,m]）

error（'MATLAB:

gmres:

PreConditioner1Size','%s%d%s',...

'Preconditionermustbeasquarematrixofsize',m,...

'tomatchtheproblemsize.'）;

end

end

else

existM1=0;

m1type='matrix';

end

if（（nargin>=7）&&~isempty（M2））

existM2=1;

[m2type,m2fun,m2fcnstr]=iterchk（M2）;

ifstrcmp（m2type,'matrix'）

if~isequal（size（M2）,[m,m]）

error（'MATLAB:

gmres:

PreConditioner2Size','%s%d%s',...

'Preconditionermustbeasquarematrixofsize',m,...

'tomatchtheproblemsize.'）;

end

end

else

existM2=0;

m2type='matrix';

end

if（（nargin>=8）&&~isempty（x））

if~isequal（size（x）,[n,1]）

error（'MATLAB:

gmres:

XoSize','%s%d%s',...

'Initialguessmustbeacolumnvectoroflength',n,...

'tomatchtheproblemsize.'）;

end

else

x=zeros（n,1）;

end

if（（nargin>8）&&strcmp（atype,'matrix'）&&...

strcmp（m1type,'matrix'）&&strcmp（m2type,'matrix'））

error（'MATLAB:

gmres:

TooManyInputs','Toomanyinputarguments.'）;

end

%Setupforthemethod

flag=1;

xmin=x;%Iteratewhichhasminimalresidualsofar

imin=0;%"Outer"iterationatwhichxminwascomputed

jmin=0;%"Inner"iterationatwhichxminwascomputed

tolb=tol*n2b;%Relativetolerance

evalxm=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）;%Normofinitialresidual

if（normr<=tolb）%Initialguessisagoodenoughsolution

flag=0;

relres=normr/n2b;

iter=[00];

resvec=normr;

if（nargout<2）

itermsg（'gmres',tol,maxit,[00],flag,iter,relres）;

end

return

end

minv_b=b;

ifexistM1

r=iterapp（'mldivide',m1fun,m1type,m1fcnstr,r,varargin{:

}）;

if~x0iszero

minv_b=iterapp（'mldivide',m1fun,m1type,m1fcnstr,b,varargin{:

}）;

else

minv_b=r;

end

if~all（isfinite（r））||~all（isfinite（minv_b））

flag=2;

x=xmin;

relres=normr/n2b;

iter=[00];

resvec=normr;

return

end

end

ifexistM2

r=iterapp（'mldivide',m2fun,m2type,m2fcnstr,r,varargin{:

}）;

if~x0iszero

minv_b=iterapp（'mldivide',m2fun,m2type,m2fcnstr,minv_b,varargin{:

}）;

else

minv_b=r;

end

if~all（isfinite（r））||~all（isfinite（minv_b））

flag=2;

x=xmin;

relres=normr/n2b;

iter=[00];

resvec=normr;

return

end

end

normr=norm（r）;%normofthepreconditionedresidual

n2minv_b=norm（minv_b）;%normofthepreconditionedrhs

clearminv_b;

tolb=tol*n2minv_b;

if（normr<=tolb）%Initialguessisagoodenoughsolution

flag=0;

relres=normr/n2minv_b;

iter=[00];

resvec=n2minv_b;

if（nargout<2）

itermsg（'gmres',tol,maxit,[00],flag,iter,relres）;

end

return

end

resvec=zeros（inner*outer+1,1）;%Preallocatevectorfornormofresiduals

resvec

（1）=normr;%resvec

（1）=norm（b-A*x0）

normrmin=normr;%Normofresidualfromxmin

%PreallocateJtoholdtheGiven'srotationconstants.

J=zeros（2,inner）;

U=zeros（n,inner）;

R=zeros（inner,inner）;

w=zeros（inner+1,1）;

foroutiter=1:

outer

%ConstructuforHouseholderreflector.

%u=r+sign（r

（1））*||r||*e1

u=r;

normr=norm（r）;

beta=scalarsign（r

（1））*normr;

u

（1）=u

（1）+beta;

u=u/norm（u）;

U（:

1）=u;

%ApplyHouseholderprojectiontor.

%w=r-2*u*u'*r;

w

（1）=-beta;

foriniter=1:

inner

%FormP1*P2*P3...Pj*ej.

%v=Pj*ej=ej-2*u*u'*ej

v=-2*（u（initer）'）*u;

v（initer）=v（initer）+1;

%v=P1*P2*...Pjm1*（Pj*ej）

fork=（initer-1）:

-1:

1

v=v-U（:

k）*（2*（U（:

k）'*v））;

end

%Explicitlynormalizevtoreducetheeffectsofround-off.

v=v/norm（v）;

%ApplyAt