cgls.m

这是在网上下的一个东东

function [X,rho,eta,F] = cgls(A,b,k,reorth,s) 
%CGLS Conjugate gradient algorithm applied implicitly to the normal equations. 
% 
% [X,rho,eta,F] = cgls(A,b,k,reorth,s) 
% 
% Performs k steps of the conjugate gradient algorithm applied 
% implicitly to the normal equations A'*A*x = A'*b. 
% 
% The routine returns all k solutions, stored as columns of 
% the matrix X.  The corresponding solution and residual norms 
% are returned in the vectors eta and rho, respectively. 
% 
% If the singular values s are also provided, cgls computes the 
% filter factors associated with each step and stores them 
% columnwise in the matrix F. 
% 
% Reorthogonalization of the normal equation residual vectors 
% A'*(A*X(:,i)-b) is controlled by means of reorth: 
%    reorth = 0 : no reorthogonalization (default), 
%    reorth = 1 : reorthogonalization by means of MGS. 
 
% References: A. Bjorck, "Numerical Methods for Least Squares Problems", 
% SIAM, Philadelphia, 1996. 
% C. R. Vogel, "Solving ill-conditioned linear systems using the 
% conjugate gradient method", Report, Dept. of Mathematical 
% Sciences, Montana State University, 1987. 
  
% Per Christian Hansen, IMM, April 8, 2001. 
 
% The fudge threshold is used to prevent filter factors from exploding. 
fudge_thr = 1e-4; 
  
% Initialization. 
if (k < 1), error('Number of steps k must be positive'), end 
if (nargin==3), reorth = 0; end 
if (nargout==4 & nargin<5), error('Too few input arguments'), end 
if (reorth<0 | reorth>1), error('Illegal reorth'), end 
[m,n] = size(A); X = zeros(n,k); 
if (reorth==1), ATr = zeros(n,k); end 
if (nargout > 1) 
  eta = zeros(k,1); rho = eta; 
end 
if (nargin==5) 
  F = zeros(n,k); Fd = zeros(n,1); s2 = s.^2; 
end 
 
% Prepare for CG iteration. 
x = zeros(n,1); 
d = (b'*A)';   % A'*b; 
r = b; 
normr2 = d'*d; 
if (reorth==1), ATr(:,1) = d/norm(d); end 
 
% Iterate. 
for j=1:k 
 
  % Update x and r vectors. 
  Ad = A*d; alpha = normr2/(Ad'*Ad); 
  x  = x + alpha*d; 
  r  = r - alpha*Ad; 
  s  = (r'*A)';   % A'*r; 
 
  % Reorthogonalize s to previous s-vectors, if required. 
  if (reorth==1) 
    for i=1:j-1, s = s - (ATr(:,i)'*s)*ATr(:,i); end 
    ATr(:,j) = s/norm(s); 
  end 
 
  % Update d vector. 
  normr2_new = s'*s; 
  beta = normr2_new/normr2; 
  normr2 = normr2_new; 
  d = s + beta*d; 
  X(:,j) = x; 
   
  % Compute norms, if required. 
  if (nargout>1), rho(j) = norm(r); end 
  if (nargout>2), eta(j) = norm(x); end 
 
  % Compute filter factors, if required. 
  if (nargin==5) 
    if (j==1) 
      F(:,1) = alpha*s2; 
      Fd = s2 - s2.*F(:,1) + beta*s2; 
    else 
      F(:,j) = F(:,j-1) + alpha*Fd; 
      Fd = s2 - s2.*F(:,j) + beta*Fd; 
    end 
    if (j > 2) 
      f = find(abs(F(:,j-1)-1) < fudge_thr & abs(F(:,j-2)-1) < fudge_thr); 
      if (length(f) > 0), F(f,j) = ones(length(f),1); end 
    end 
  end 
 
end