function [X,rho,eta,F] = pcgls(A,L,W,b,k,sm)
%PCGLS "Preconditioned" conjugate gradients appl. implicitly to normal equations.
% [X,rho,eta,F] = pcgls(A,L,W,b,k,sm)
%
% Performs k steps of the `preconditioned' conjugate gradient
% algorithm applied implicitly to the normal equations
%    (A*L_p)'*(A*L_p)*x = (A*L_p)'*b ,
% where L_p is the A-weighted generalized inverse of L.  Notice
% that the matrix W holding a basis for the null space of L must
% also be specified.
%
% The routine returns all k solutions, stored as columns of the
% matrix X.  The solution seminorm and residual norm are returned in eta
% and rho, respectively.
%
% If the generalized singular values sm of (A,L) are also provided,
% pcgls computes the filter factors associated with each step and
% stores them columnwise in the matrix F.

% References: A. Bjorck, "Least Squares Methods", in P. G.
% Ciarlet & J. L Lions (Eds.), "Handbook of  Numerical Analysis,
% Vol. I", Elsevier, Amsterdam, 1990; p. 560.
% M. Hanke, "Regularization with differential operators.  An itera-
% tive approach", J. Numer. Funct. Anal. Optim. 13 (1992), 523-540.
% 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, UNI-C and Martin Hanke, Institut fuer
% Praktische Mathematik, Universitaet Karlsruhe, 11/05/92.

% 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
[m,n] = size(A); [p,n1] = size(L); X = zeros(n,k);
if (nargout > 1)
  eta = zeros(k,1); rho = eta;
end
if (nargout>3 & nargin==5), error('Too few imput arguments'), end
if (nargin==6)
  F = zeros(p,k); Fd = zeros(p,1); s = (sm(:,1)./sm(:,2)).^2;
end

% Prepare for computations with L_p.
[NAA,x_0] = pinit(W,A,b);

% Prepare for CG iteartion.
x  = x_0;
r  = b - A*x_0; s = A'*r;
q1 = ltsolve(L,s);
q  = lsolve(L,q1,W,NAA);
z  = q;
dq = s'*q;
if (nargout>2), z1 = q1; x1 = zeros(p,1); end

% Iterate.
for j=1:k

  Az  = A*z; alpha = dq/(Az'*Az);
  x   = x + alpha*z;
  r   = r - alpha*Az; s = A'*r;
  q1  = ltsolve(L,s);
  q   = lsolve(L,q1,W,NAA);
  dq2 = s'*q; beta = dq2/dq;
  dq  = dq2;
  z   = q + beta*z;
  X(:,j) = x;
  if (nargout>1), rho(j) = norm(r); end
  if (nargout>2)
    x1 = x1 + alpha*z1; z1 = q1 + beta*z1; eta(j) = norm(x1);
  end

  if (nargin==6)
    if (j==1)
      F(:,1) = alpha*s;
      Fd = s - s.*F(:,1) + beta*s;
    else
      F(:,j) = F(:,j-1) + alpha*Fd;
      Fd = s - s.*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
