gmres源程序matlab.doc

-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 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.% * = 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*imum 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 preconditioner 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 all 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).% *,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: 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(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)% 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.% 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 arguments.');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) error('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 parametersif (nargin < 3) | isempty(restart) | (restart = n) restarted = false;else restarted = true;endif (nargin < 4) | isempty(tol) tol = 1e-6;endwarned = 0;if tol < eps warning('MATLAB:gmres:tooSmallTolerance', . strcat('Input tol is smaller than eps and may not be achieved',. ' by GMRESn',' Try to use a bigger tolerance'); warned = 1; tol = eps;elseif 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 (nargin < 5) | isempty(ma*it) if restarted ma*it = min(ceil(n/restart),10); else ma*it = min(n,10); endendif restarted outer = ma*it; if restart > 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 bounded 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 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 < 2) itermsg('gmres',tol,ma*it,0,flag,iter,NaN); end returnendif (nargin >= 6) && isempty(M1) e*istM1 = 1; m1type,m1fun,m1fcnstr = iterchk(M1); if strcmp(m1type,'matri*') if isequal(size(M1),m,m) error('MATLAB:gmres: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('MATLAB: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 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 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('mtimes',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 = 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('mldivide',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 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,ma*it,0 0,flag,iter,relres); 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 Given's 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(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 : 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 % E*plicitly 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 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. 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+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 Given's 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+1); end % Compute Given's 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); resvec(outiter-1)*inner+initer+1) = normr; normr_act = normr; if (normr <= tolb | stag >= 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(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,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)