gmres源程序matlab.doc

《gmres源程序matlab.doc》由会员分享，可在线阅读，更多相关《gmres源程序matlab.doc（13页珍藏版）》请在课桌文档上搜索。

1、-function *,flag,relres,iter,resvec = gmres(A,b,restart,tol,ma*it,M1,M2,*,varargin)%GMRES Generalized Minimum Residual Method.% * = GMRES(A,B) attempts to solve the system of linear equations A* = B% for *. The N-by-N coefficient matri* A must be square and the right% hand side column vector B must

2、have length N. This uses the unrestarted% method with MIN(N,10) total iterations.% * = GMRES(AFUN,B) accepts a function handle AFUN instead of the matri*% A. AFUN(*) accepts a vector input * and returns the matri*-vector% product A*. In all of the following synta*es, you can replace A by% AFUN.% * =

3、 GMRES(A,B,RESTART) restarts the method every RESTART iterations.% If RESTART is N or then GMRES uses the unrestarted method as above.% * = 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,MA*IT) specifies the ma*i

4、mum number of outer% iterations. Note: the total number of iterations is RESTART*MA*IT. If% MA*IT is then GMRES uses the default, MIN(N/RESTART,10). If RESTART% is N or then the total number of iterations is MA*IT.% * = GMRES(A,B,RESTART,TOL,MA*IT,M) and% * = GMRES(A,B,RESTART,TOL,MA*IT,M1,M2) use p

5、reconditioner M or M=M1*M2% and effectively solve the system inv(M)*A* = inv(M)*B for *. If M is% then a preconditioner is not applied. M may be a function handle% returning M*.% * = GMRES(A,B,RESTART,TOL,MA*IT,M1,M2,*0) specifies the first initial% guess. If *0 is then GMRES uses the default, an al

6、l zero vector.% *,FLAG = GMRES(A,B,.) also returns a convergence FLAG:% 0 GMRES converged to the desired tolerance TOL within MA*IT iterations.% 1 GMRES iterated MA*IT times but did not converge.% 2 preconditioner M was ill-conditioned.% 3 GMRES stagnated (two consecutive iterates were the same).% *

7、,FLAG,RELRES = GMRES(A,B,.) also returns the relative residual% NORM(B-A*)/NORM(B). If FLAG is 0, then RELRES = TOL. Note with% preconditioners M1,M2, the residual is NORM(M2(M1(B-A*).% *,FLAG,RELRES,ITER = GMRES(A,B,.) also returns both the outer and% inner iteration numbers at which * was computed

8、: 0 = ITER(1) = MA*IT% and 0 = ITER(2) = RESTART.% *,FLAG,RELRES,ITER,RESVEC = GMRES(A,B,.) also returns a vector of% the residual norms at each inner iteration, including NORM(B-A*0).% Note with preconditioners M1,M2, the residual is NORM(M2(M1(B-A*).% E*ample:% n = 21; A = gallery(wilk,n); b = sum

9、(A,2);% tol = 1e-12; ma*it = 15; M = diag(10:-1:1 1 1:10);% * = gmres(A,b,10,tol,ma*it,M);% Or, use this matri*-vector product function% %-% function y = afun(*,n)% y = 0; *(1:n-1) + (n-1)/2:-1:0); (1:(n-1)/2).*+*(2:n); 0;% %-% and this preconditioner backsolve function% %-% function y = mfun(r,n)%

10、y = r ./ (n-1)/2:-1:1); 1; (1:(n-1)/2);% %-% as inputs to GMRES:% *1 = gmres(*)afun(*,n),b,10,tol,ma*it,(*)mfun(*,n);% Class support for inputs A,B,M1,M2,*0 and the output of AFUN:% float: double% See also BICG, BICGSTAB, BICGSTABL, CGS, LSQR, MINRES, PCG, QMR, SYMMLQ,% TFQMR, ILU, FUNCTION_HANDLE.%

11、 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.% $Revision: 1.21.4.13 $ $Date: 2010/02/25 08:11:48 $if (nargin 2) error(MATLAB:gmres:NumInputs,Not enough input argu

12、ments.);end% Determine whether A is a matri* or a function.atype,afun,afcnstr = iterchk(A);if strcmp(atype,matri*) % Check matri* and right hand side vector inputs have appropriate sizes m,n = size(A); if (m = n) error(MATLAB:gmres:SquareMatri*,Matri* must be square.); end if isequal(size(b),m,1) er

13、ror(MATLAB:gmres:VectorSize, %s %d %s, . Right hand side must be a column vector of length, m,. to match the coefficient matri*.); endelse m = size(b,1); n = m; if iscolumn(b) error(MATLAB:gmres:Vector,Right hand side must be a column vector.); endend% Assign default values to unspecified parameters

14、if (nargin 3) | isempty(restart) | (restart = n) restarted = false;else restarted = true;endif (nargin 4) | isempty(tol) tol = 1e-6;endwarned = 0;if tol = 1 warning(MATLAB:gmres:tooBigTolerance, . strcat(Input tol is bigger than 1 n,. Try to use a smaller tolerance); warned = 1; tol = 1-eps;endif (n

15、argin n warning(MATLAB:gmres:tooManyInnerIts, RESTART is %d %sn%s %d., . restart, but it should be bounded by SIZE(A,1)., . Setting RESTART to, n); restart = n; end inner = restart;else outer = 1; if ma*it n warning(MATLAB:gmres:tooManyInnerIts, MA*IT is %d %sn%s %d., . ma*it, but it should be bound

16、ed by SIZE(A,1)., . Setting MA*IT to, n); ma*it = n; end inner = ma*it;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 zeros * = zeros(n,1); % then solution is all zeros flag = 0; % a valid solution has been

17、obtained relres = 0; % the relative residual is actually 0/0 iter = 0 0; % no iterations need be performed resvec = 0; % resvec(1) = norm(b-A*) = norm(0) if (nargout = 6) & isempty(M1) e*istM1 = 1; m1type,m1fun,m1fcnstr = iterchk(M1); if strcmp(m1type,matri*) if isequal(size(M1),m,m) error(MATLAB:gm

18、res:PreConditioner1Size, %s %d %s,. Preconditioner must be a square matri* of size, m, . to match the problem size.); end endelse e*istM1 = 0; m1type = matri*;endif (nargin = 7) & isempty(M2) e*istM2 = 1; m2type,m2fun,m2fcnstr = iterchk(M2); if strcmp(m2type,matri*) if isequal(size(M2),m,m) error(MA

19、TLAB:gmres:PreConditioner2Size, %s %d %s, . Preconditioner must be a square matri* of size, m, . to match the problem size.); end endelse e*istM2 = 0; m2type = matri*;endif (nargin = 8) & isempty(*) if isequal(size(*),n,1) error(MATLAB:gmres:*oSize, %s %d %s, . Initial guess must be a column vector

20、of length, n, . to match the problem size.); endelse * = zeros(n,1);endif (nargin 8) & strcmp(atype,matri*) & . strcmp(m1type,matri*) & strcmp(m2type,matri*) error(MATLAB:gmres:TooManyInputs,Too many input arguments.);end% Set up for the methodflag = 1;*min = *; % Iterate which has minimal residual

21、so farimin = 0; % Outer iteration at which *min was computedjmin = 0; % Inner iteration at which *min was computedtolb = tol * n2b; % Relative toleranceeval*m = 0;stag = 0;moresteps = 0;ma*msteps = min(floor(n/50),5,n-ma*it);ma*stagsteps = 3;minupdated = 0;*0iszero = (norm(*) = 0);r = b - iterapp(mt

22、imes,afun,atype,afcnstr,*,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,ma*it,0 0,flag,iter,relres); end returnendminv_b = b;if e*istM1 r =

23、 iterapp(mldivide,m1fun,m1type,m1fcnstr,r,varargin:); if *0iszero minv_b = iterapp(mldivide,m1fun,m1type,m1fcnstr,b,varargin:); else minv_b = r; end if all(isfinite(r) | all(isfinite(minv_b) flag = 2; * = *min; relres = normr / n2b; iter = 0 0; resvec = normr; return endendif e*istM2 r = iterapp(mld

24、ivide,m2fun,m2type,m2fcnstr,r,varargin:); if *0iszero 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; * = *min; relres = normr / n2b; iter = 0 0; resvec = normr; return endendnormr = norm(r); % norm of the pr

25、econditioned 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,ma*it,0 0,flag,iter,relres

26、); end returnendresvec = zeros(inner*outer+1,1); % Preallocate vector for norm of residualsresvec(1) = normr; % resvec(1) = norm(b-A*0)normrmin = normr; % Norm of residual from *min% Preallocate J to hold the Givens rotation constants.J = zeros(2,inner);U = zeros(n,inner);R = zeros(inner,inner);w =

27、zeros(inner+1,1);for outiter = 1 : outer % Construct u for Householder reflector. % 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; % Apply Householder projection to r. % w = r - 2*u*u*r; w(1) = -beta; for initer = 1 : i

28、nner % Form P1*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) for k = (initer-1):-1:1 v = v - U(:,k)*(2*(U(:,k)*v); end % E*plicitly normalize v to reduce the effects of round-off. v = v/norm(v); % Apply A to v. v = iterapp(mtimes,a

29、fun,atype,afcnstr,v,varargin:); % Apply Preconditioner. if e*istM1 v = iterapp(mldivide,m1fun,m1type,m1fcnstr,v,varargin:); if all(isfinite(v) flag = 2; break end end if e*istM2 v = iterapp(mldivide,m2fun,m2type,m2fcnstr,v,varargin:); if all(isfinite(v) flag = 2; break end end % Form Pj*Pj-1*.P1*Av.

30、 for k = 1:initer v = 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+

31、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; end end % Apply Givens rotations to the newly formed v. for colJ = 1:initer-1 tmpv = v(colJ); v(colJ) = conj(J(1,colJ)*v(colJ)

32、 + conj(J(2,colJ)*v(colJ+1); v(colJ+1) = -J(2,colJ)*tmpv + J(1,colJ)*v(colJ+1); end % Compute Givens 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

33、; v(initer+1) = 0; end R(:,initer) = v(1:inner); normr = abs(w(initer+1); resvec(outiter-1)*inner+initer+1) = normr; normr_act = normr; if (normr = ma*stagsteps | moresteps) if eval*m = 0 ytmp = R(1:initer,1:initer) w(1:initer); additive = U(:,initer)*(-2*ytmp(initer)*conj(U(initer,initer); additive

34、(initer) = additive(initer) + ytmp(initer); for k = initer-1 : -1 : 1 additive(k) = additive(k) + ytmp(k); additive = additive - U(:,k)*(2*(U(:,k)*additive); end if norm(additive) eps*norm(*) stag = stag + 1; else stag = 0; end *m = * + additive; eval*m = 1; elseif eval*m = 1 addvc = -(R(1:initer-1,

35、1:initer-1)R(1:initer-1,initer)*. (w(initer)/R(initer,initer); w(initer)/R(initer,initer); if norm(addvc) eps*norm(*m) stag = stag + 1; else stag = 0; end additive = U(:,initer)*(-2*addvc(initer)*conj(U(initer,initer); additive(initer) = additive(initer) + addvc(initer); for k = initer-1 : -1 : 1 additive(k) = additive(k) + addvc(k); additive = additive - U(:,k)*(2*(U(:,k)*additive); end *m = *m + additive; end r = b - iterapp(mtimes,afun,atype,afcnstr,*m,varargin:); if norm(r)