cgls_mod

function [X,rho,eta,F] = cgls(A,b,k,reorth,s)

0001 %function [X,rho,eta,F] = cgls(A,b,k,reorth,s)
0002 function [x,rho,eta,F] = cgls_mod(A,b,k,reorth,s)
0003 
0004 %CGLS Conjugate gradient algorithm applied implicitly to the normal equations.
0005 %
0006 % [X,rho,eta,F] = cgls(A,b,k,reorth,s)
0007 %
0008 % Performs k steps of the conjugate gradient algorithm applied
0009 % implicitly to the normal equations A'*A*x = A'*b.
0010 %
0011 % The routine returns all k solutions, stored as columns of
0012 % the matrix X.  The corresponding solution and residual norms
0013 % are returned in the vectors eta and rho, respectively.
0014 %
0015 % If the singular values s are also provided, cgls computes the
0016 % filter factors associated with each step and stores them
0017 % columnwise in the matrix F.
0018 %
0019 % Reorthogonalization of the normal equation residual vectors
0020 % A'*(A*X(:,i)-b) is controlled by means of reorth:
0021 %    reorth = 0 : no reorthogonalization (default),
0022 %    reorth = 1 : reorthogonalization by means of MGS.
0023 
0024 % References: A. Bjorck, "Numerical Methods for Least Squares Problems",
0025 % SIAM, Philadelphia, 1996.
0026 % C. R. Vogel, "Solving ill-conditioned linear systems using the
0027 % conjugate gradient method", Report, Dept. of Mathematical
0028 % Sciences, Montana State University, 1987.
0029 
0030 % Per Christian Hansen, IMM, July 23, 2007.
0031 
0032 % modified by Wolfgang Förstner: eliminated the storage of X.
0033 
0034 % The fudge threshold is used to prevent filter factors from exploding.
0035 fudge_thr = 1e-4;
0036 
0037 % Initialization.
0038 if (k < 1), error('Number of steps k must be positive'), end
0039 if (nargin==3), reorth = 0; end
0040 %if (nargout==4 & nargin<5), error('Too few input arguments'), end
0041 if (nargout==4 & nargin<4), error('Too few input arguments'), end
0042 if (reorth<0 | reorth>1), error('Illegal reorth'), end
0043 [m,n] = size(A); 
0044 %X = zeros(n,k);
0045 if (reorth==1), ATr = zeros(n,k+1); end
0046 if (nargout > 1)
0047   eta = zeros(k,1); rho = eta;
0048 end
0049 F=0;
0050 if (nargin==5)
0051   F = zeros(n,k); Fd = zeros(n,1); s2 = s.^2;
0052 end
0053 
0054 % Prepare for CG iteration.
0055 x = zeros(n,1);
0056 d = A'*b;
0057 r = b;
0058 normr2 = d'*d;
0059 if (reorth==1), ATr(:,1) = d/norm(d); end
0060 
0061 % Iterate.
0062 for j=1:k
0063 
0064   % Update x and r vectors.
0065   Ad = A*d; alpha = normr2/(Ad'*Ad);
0066   x  = x + alpha*d;
0067   r  = r - alpha*Ad;
0068   s  = A'*r;
0069 
0070   % Reorthogonalize s to previous s-vectors, if required.
0071   if (reorth==1)
0072     for i=1:j, s = s - (ATr(:,i)'*s)*ATr(:,i); end
0073     ATr(:,j+1) = s/norm(s);
0074   end
0075 
0076   % Update d vector.
0077   normr2_new = s'*s;
0078   beta = normr2_new/normr2;
0079   normr2 = normr2_new;
0080   d = s + beta*d;
0081   %X(:,j) = x;
0082 
0083   % Compute norms, if required.
0084   if (nargout>1), rho(j) = norm(r); end
0085   if (nargout>2), eta(j) = norm(x); end
0086 
0087   % Compute filter factors, if required.
0088   if (nargin==5)
0089     if (j==1)
0090       F(:,1) = alpha*s2;
0091       Fd = s2 - s2.*F(:,1) + beta*s2;
0092     else
0093       F(:,j) = F(:,j-1) + alpha*Fd;
0094       Fd = s2 - s2.*F(:,j) + beta*Fd;
0095     end
0096     if (j > 2)
0097       f = find(abs(F(:,j-1)-1) < fudge_thr & abs(F(:,j-2)-1) < fudge_thr);
0098       if ~isempty(f), F(f,j) = ones(length(f),1); end
0099     end
0100   end
0101 
0102 end