gmres源程序 matlab(13页).doc

《gmres源程序 matlab(13页).doc》由会员分享，可在线阅读，更多相关《gmres源程序 matlab(13页).doc（13页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、-gmres源程序 matlab-第 13 页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 sid

2、e 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

3、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,M

4、AXIT) 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,RESTAR

5、T,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 MX.% X = GMRES(A,B,RESTART,TOL,MAXIT,M1,M2,X0) specifies the first initial% guess. If X0 is then GMRE

6、S 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 it

7、erates 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 numb

8、ers 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;

9、 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% functi

10、on 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, FUN

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

12、ugh 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.); end if isequal(s

13、ize(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.); 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 unspecif

14、ied parametersif (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 =

15、1-eps;endif (nargin 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 maxit n warning(MATLAB:gmres:tooManyInnerIts, MAXIT is %d %sn%s %d., . maxit, but it s

16、hould be bounded by SIZE(A,1)., . Setting MAXIT to, n); maxit = n; end inner = 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 zeros x = zeros(n,1); % then solution is all zeros flag = 0; % a valid solu

17、tion has been 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*x) = norm(0) if (nargout = 6) & isempty(M1) existM1 = 1; m1type,m1fun,m1fcnstr = iterchk(M1); if strcmp(m1type,matrix) if isequal(size(M1),m,m)

18、error(MATLAB:gmres:PreConditioner1Size, %s %d %s,. Preconditioner must be a square matrix of size, m, . to match the problem size.); end endelse existM1 = 0; m1type = matrix;endif (nargin = 7) & isempty(M2) existM2 = 1; m2type,m2fun,m2fcnstr = iterchk(M2); if strcmp(m2type,matrix) if isequal(size(M2

19、),m,m) error(MATLAB:gmres:PreConditioner2Size, %s %d %s, . Preconditioner must be a square matrix of size, m, . to match the problem size.); end endelse existM2 = 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

20、 column vector of length, n, . to match the problem size.); endelse x = 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 mi

21、nimal residual so farimin = 0; % Outer iteration at which xmin was computedjmin = 0; % Inner iteration at which xmin was computedtolb = 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 =

22、 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); end returnendminv_b = b

23、;if existM1 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 = 0 0; resvec = normr; return endendif existM2

24、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 = 0 0; resvec = normr; return endendnormr = norm(r); %

25、 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,f

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

27、nner,inner);w = 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; fo

28、r initer = 1 : inner % 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 % Explicitly normalize v to reduce the effects of round-off. v = v/norm(v); % Apply A to v. v =

29、iterapp(mtimes,afun,atype,afcnstr,v,varargin:); % Apply Preconditioner. if existM1 v = iterapp(mldivide,m1fun,m1type,m1fcnstr,v,varargin:); if all(isfinite(v) flag = 2; break end end if existM2 v = iterapp(mldivide,m2fun,m2type,m2fcnstr,v,varargin:); if all(isfinite(v) flag = 2; break end end % Form

30、 Pj*Pj-1*.P1*Av. 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(init

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

32、(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 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);

33、 v(initer) = rho; 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 = maxstagsteps | moresteps) if evalxm = 0 ytmp = R(1:initer,1:initer) w(1:initer); additive = U(:,initer)*(-2*ytmp(initer)*conj(U(initer,i

34、niter); additive(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(x) stag = stag + 1; else stag = 0; end xm = x + additive; evalxm = 1; elseif evalxm = 1 addvc =

35、 -(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; 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 xm = xm + additive; end r = b - iterapp(mtimes,afun,atype,afcnstr,xm,varargin:); if norm(r) = tol*n2b x = xm; flag = 0; iter = outiter, initer; break end minv_r = r; if exist