gmres源程序matlabWord版

1、傳播優(yōu)秀Word版文檔 ，希望對(duì)您有幫助，可雙擊去除！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

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

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

4、OL,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,RE

5、START,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

6、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 consecutiv

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

8、numbers at which X was computed: 0 = ITER(1) = MAXIT傳播優(yōu)秀Word版文檔 ，希望對(duì)您有幫助，可雙擊去除！% 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(M

9、2(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 pr

10、econditioner 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, L

11、SQR, 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.% $Revision: 1.21.4.13 $ $Date: 2010/02/25 08:11:48 $if

12、(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傳播優(yōu)秀Word版文檔 ，希望對(duì)您有幫助，可雙擊去除！ m,n = size(A); if (m = n)

13、 error(MATLAB:gmres:SquareMatrix,Matrix must be square.); end if 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.); endelse m = size(b,1); n = m; if iscolumn(b) error(MATLAB:gmres:Vector,Right hand

14、 side must be a column vector.); endend% Assign default values to unspecified 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 to

15、l is bigger than 1 n,. Try to use a smaller tolerance); warned = 1; tol = 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 warn

16、ing(MATLAB:gmres:tooManyInnerIts, MAXIT is %d %sn%s %d., . maxit, but it should 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 z

17、eros x = zeros(n,1); % then solution is all zeros flag = 0; % a valid solution 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,

18、m1fcnstr = iterchk(M1); if strcmp(m1type,matrix) if isequal(size(M1),m,m) error(MATLAB:gmres:PreConditioner1Size, %s %d %s,.傳播優(yōu)秀Word版文檔 ，希望對(duì)您有幫助，可雙擊去除！ Preconditioner must be a square matrix of size, m, . to match the problem size.); end endelse existM1 = 0; m1type = matrix;endif (nargin = 7) & isem

19、pty(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.); end endelse existM2 = 0; m2type = matrix;endif (nargin = 8

20、) & 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.); endelse x = zeros(n,1);endif (nargin 8) & strcmp(atype,matrix) & . strcmp(m1type,matrix) & strcmp(m2type,matrix) error(MATLAB:gmres:TooMany

21、Inputs,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 computedjmin = 0; % Inner iteration at which xmin was computed傳播優(yōu)秀Word版文檔 ，希望對(duì)您有幫助，可雙擊去除！tolb = tol * n2b; % Relative toleranceevalx

22、m = 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

23、= normr / n2b; iter = 0 0; resvec = normr; if (nargout 2) itermsg(gmres,tol,maxit,0 0,flag,iter,relres); end returnendminv_b = b;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(i

24、sfinite(r) | all(isfinite(minv_b) flag = 2; x = xmin; relres = normr / n2b; iter = 0 0; resvec = normr; return endendif existM2 r = iterapp(mldivide,m2fun,m2type,m2fcnstr,r,varargin:); if x0iszero minv_b = iterapp(mldivide,m2fun,m2type,m2fcnstr,minv_b,varargin:);傳播優(yōu)秀Word版文檔 ，希望對(duì)您有幫助，可雙擊去除！ else minv

25、_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); % norm of the preconditioned residualn2minv_b = norm(minv_b); % norm of the preconditioned rhsclear minv_b;tolb = tol * n2minv_b;if (normr = tolb)

26、 % 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); end returnendresvec = zeros(inner*outer+1,1); % Preallocate vector for norm of residualsresvec(1) = normr; % resvec(1) = norm(b

27、-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(inner,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(

28、r);傳播優(yōu)秀Word版文檔 ，希望對(duì)您有幫助，可雙擊去除！ 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*ej v = -2*(u(initer)*u; v(initer) = v(initer) + 1

29、; % 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 = iterapp(mtimes,afun,atype,afcnstr,v,varargin:); % Apply Preconditioner. if existM1 v = iterapp(mldivide,m1fun,m1type,m

30、1fcnstr,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 Pj*Pj-1*.P1*Av. for k = 1:initer v = v - U(:,k)*(2*(U(:,k)*v); end傳播優(yōu)秀Word版文檔 ，希望對(duì)您有幫助，可雙擊去除！ % Determine Pj+1. if (i

31、niter = 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(:

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

33、+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; v(initer+1) = 0; end R(:,initer) = v(1:inner); normr = abs(w(initer+1); r

34、esvec(outiter-1)*inner+initer+1) = normr; normr_act = normr; 傳播優(yōu)秀Word版文檔 ，希望對(duì)您有幫助，可雙擊去除！ 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,initer); additive(initer) = additive(initer) + ytmp(initer); for

35、 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 = -(R(1:initer-1,1:initer-1)R(1:initer-1,initer)*. (w(initer)/R(

36、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 =

37、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 existM1 minv_r = iterapp(mldivide,m1fun,m1type,m1fcnstr,r,varargin:); if all(isfinite(m

38、inv_r) flag = 2;傳播優(yōu)秀Word版文檔 ，希望對(duì)您有幫助，可雙擊去除！ break end end if existM2 minv_r = iterapp(mldivide,m2fun,m2type,m2fcnstr,minv_r,varargin:); if all(isfinite(minv_r) flag = 2; break end end normr_act = norm(minv_r); resvec(outiter-1)*inner+initer+1) = normr_act; if normr_act = normrmin normrmin = normr_ac

39、t; imin = outiter; jmin = initer; xmin = xm; minupdated = 1; end if normr_act = maxstagsteps & moresteps = 0 stag = 0; end moresteps = moresteps + 1; if moresteps = maxmsteps if warned warning(MATLAB:gmres:tooSmallTolerance, . strcat(Input tol might be obviously smaller than,. eps*cond(A) and may no

40、t be achieved by GMRESn,. Try to use a bigger tolerance); end flag = 3; iter = outiter, initer; break; end end傳播優(yōu)秀Word版文檔 ，希望對(duì)您有幫助，可雙擊去除！ end if normr_act = maxstagsteps flag = 3; break; end end % ends inner loop evalxm = 0; if flag = 0 if minupdated idx = jmin; else idx = initer; end y = 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 :