截断共轭梯度法MATLAB代码示例

以下是截断共轭梯度法的MATLAB代码示例：

function [x,flag,relres,iter,resvec] = tCG(A,b,tol,maxit,M,x0)
% tCG truncated Conjugate Gradient Method
%   x = tCG(A,b) attempts to solve the system of linear equations A*x = b
%   for x using the truncated Conjugate Gradient method. The matrix A 
%   should be symmetric and positive definite. If A is not symmetric, then 
%   the method automatically solves the symmetric part (A+A')/2.
%
%   x = tCG(A,b,tol) specifies the tolerance of the method. If tol is []
%   then tCG uses the default, 1e-6.
%
%   x = tCG(A,b,tol,maxit) specifies the maximum number of iterations.
%   If maxit is [] then tCG uses the default, min(n,20).
%
%   x = tCG(A,b,tol,maxit,M) and x = tCG(A,b,tol,maxit,M,x0) use the 
%   preconditioner M. If M is [] then tCG uses no preconditioner. If M is a
%   function handle then M(x) returns the product of the matrix inverse of
%   the preconditioner times x. If M is a matrix, then tCG applies the
%   preconditioner M\y to each column of b.
%
%   [x,flag] = tCG(A,b,...) also returns a convergence flag:
%    0 tCG converged to the desired tolerance tol within maxit iterations.
%    1 tCG iterated maxit times but did not converge.
%    2 preconditioner M was ill-conditioned.
%
%   [x,flag,relres] = tCG(A,b,...) also returns the relative residual
%   norm(b-A*x)/norm(b).
%
%   [x,flag,relres,iter] = tCG(A,b,...) also returns the iteration number
%   at which x was computed: 0 <= iter <= maxit.
%
%   [x,flag,relres,iter,resvec] = tCG(A,b,...) also returns a vector of
%   the residual norms at each iteration: resvec(iter+1) = norm(b-A*x)/norm(b).
%
%   Example:
%      n = 500; A = delsq(numgrid('S',n+2)); b = sum(A,2);
%      tol = 1e-10; maxit = 500; M = spdiags(spdiags(A,0),0,n,n);
%      [x,flag,relres,iter,resvec] = tCG(A,b,tol,maxit,M);
%      semilogy(0:iter,resvec,'.-'), xlabel('iteration'), ylabel('relative residual')
%
%   See also CG, PCG, SYMMLQ, MINRES.

n = size(A,1); % length of the system
if isempty(tol), tol = 1e-6; end % default tolerance
if isempty(maxit), maxit = min(n,20); end % default maximum iterations
if isempty(x0), x0 = zeros(n,1); end % default initial guess
if isempty(M), M = speye(n,n); end % default preconditioner

flag = 1; % initialize flag to failed
x = x0; % initialize solution
r = b - A*x; % initialize residual
normr = norm(r); % initialize norm of residual
resvec = zeros(maxit+1,1); % initialize residual vector
resvec(1) = normr; % save initial residual
if (normr/norm(b)) <= tol, flag = 0; return, end % bail out early if already converged

p = M\r; % initialize direction
if norm(p) == 0, flag = 2; return, end % bail out if preconditioner is ill-conditioned
rr = r'*p; % initialize square of residual norm
for iter = 1:maxit
    Ap = A*p;
    alpha = rr/(p'*Ap);
    x = x + alpha*p;
    r = r - alpha*Ap;
    normr = norm(r);
    resvec(iter+1) = normr;
    if (normr/norm(b)) <= tol, flag = 0; break, end % success
    if mod(iter,50) == 0, % restart
        r = b - A*x; % reset residual
        normr = norm(r); % reset norm of residual
        resvec(iter+1) = normr; % save new residual
        if (normr/norm(b)) <= tol, flag = 0; break, end % success
    else
        beta = (normr/normr_old)^2;
        p = M\r + beta*p;
        normr_old = normr;
        rr = r'*p;
    end
end
resvec(iter+2:end) = []; % truncate residual vector
relres = resvec(end)/resvec(1); % relative residual