gmres源程序matlabWord文档格式.docx
《gmres源程序matlabWord文档格式.docx》由会员分享，可在线阅读，更多相关《gmres源程序matlabWord文档格式.docx（22页珍藏版）》请在冰点文库上搜索。
%[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:
0）'
;
（1:
（n-1）/2）'
].*x+[x（2:
n）;
0];
%andthispreconditionerbacksolvefunction
%%------------------------------------------%
%functiony=mfun（r,n）
%y=r./[（（n-1）/2:
1）'
1;
];
%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.
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）
SquareMatrix'
Matrixmustbesquare.'
end
if~isequal（size（b）,[m,1]）
VectorSize'
'
%s%d%s'
...
'
Righthandsidemustbeacolumnvectoroflength'
m,...
tomatchthecoefficientmatrix.'
else
m=size（b,1）;
n=m;
if~iscolumn（b）
Vector'
Righthandsidemustbeacolumnvector.'
%Assigndefaultvaluestounspecifiedparameters
3）||isempty（restart）||（restart==n）
restarted=false;
restarted=true;
4）||isempty（tol）
tol=1e-6;
warned=0;
iftol<
eps
warning（'
tooSmallTolerance'
strcat（'
Inputtolissmallerthanepsandmaynotbeachieved'
...
byGMRES\n'
Trytouseabiggertolerance'
））;
warned=1;
tol=eps;
elseiftol>
=1
tooBigTolerance'
Inputtolisbiggerthan1\n'
Trytouseasmallertolerance'
tol=1-eps;
5）||isempty（maxit）
ifrestarted
maxit=min（ceil（n/restart）,10）;
else
maxit=min（n,10）;
ifrestarted
outer=maxit;
ifrestart>
n
tooManyInnerIts'
RESTARTis%d%s\n%s%d.'
restart,'
butitshouldbeboundedbySIZE（A,1）.'
SettingRESTARTto'
n）;
restart=n;
inner=restart;
outer=1;
ifmaxit>
MAXITis%d%s\n%s%d.'
maxit,'
SettingMAXITto'
maxit=n;
inner=maxit;
%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<
itermsg（'
gmres'
tol,maxit,0,flag,iter,NaN）;
return
if（（nargin>
=6）&
&
~isempty（M1））
existM1=1;
[m1type,m1fun,m1fcnstr]=iterchk（M1）;
ifstrcmp（m1type,'
if~isequal（size（M1）,[m,m]）
PreConditioner1Size'
Preconditionermustbeasquarematrixofsize'
m,...
tomatchtheproblemsize.'
existM1=0;
m1type='
=7）&
~isempty（M2））
existM2=1;
[m2type,m2fun,m2fcnstr]=iterchk（M2）;
ifstrcmp（m2type,'
if~isequal（size（M2）,[m,m]）
PreConditioner2Size'
existM2=0;
m2type='
=8）&
~isempty（x））
if~isequal（size（x）,[n,1]）
XoSize'
Initialguessmustbeacolumnvectoroflength'
n,...
8）&
strcmp（atype,'
）&
...
strcmp（m1type,'
strcmp（m2type,'
））
TooManyInputs'
Toomanyinputarguments.'
%Setupforthemethod
flag=1;
xmin=x;
%Iteratewhichhasminimalresidualsofar
imin=0;
%"
Outer"
iterationatwhichxminwascomputed
jmin=0;
Inner"
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
relres=normr/n2b;
resvec=normr;
tol,maxit,[00],flag,iter,relres）;
minv_b=b;
ifexistM1
r=iterapp（'
mldivide'
m1fun,m1type,m1fcnstr,r,varargin{:
if~x0iszero
minv_b=iterapp（'
m1fun,m1type,m1fcnstr,b,varargin{:
minv_b=r;
if~all（isfinite（r））||~all（isfinite（minv_b））
flag=2;
x=xmin;
ifexistM2
m2fun,m2type,m2fcnstr,r,varargin{:
m2fun,m2type,m2fcnstr,minv_b,varargin{:
%normofthepreconditionedresidual
n2minv_b=norm（minv_b）;
%normofthepreconditionedrhs
clearminv_b;
tolb=tol*n2minv_b;
relres=normr/n2minv_b;
resvec=n2minv_b;
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
v=v-U（:
k）*（2*（U（:
k）'
*v））;
%Explicitlynormalizevtoreducetheeffectsofround-off.
v=v/norm（v）;
%ApplyAtov.
v=iterapp（'
afun,atype,afcnstr,v,varargin{:
%
升级会员 