% Conjugate gradient least-squares routine.
% Applies CG method to a linear least-squares problem.
% Based on algorithm given in paper by Paige and Saunders.
% ref: ACM Trans. Math. Software, v 8 , pp 43 - 71, 1982.
%
% Syntax: x = cgls(a,b)
%
% By D. Holmgren, Mar 13 94. holmgren@phobos.astro.uwo.ca
% works - tested on sample problem from Forsythe et al.
function x = cgls(a,b)
ni = input('Number of iterations ');
nu = input('Number of unknowns ');
% set up starting data...
% x = solution vector
% r = residual vector
nl = length(b);
r = zeros(nl,1);
s = zeros(nu,1);
x = zeros(nu,1);
p = zeros(nu,1);
q = zeros(nu,1);
r = b;
s = a' * b;
p = s;
gamma = norm(s)^2;
% CG loop ...
for i = 1:ni
   q = a * p;
   alpha = gamma / norm(q)^2;
   gamold = gamma;
   x = x + alpha * p;
   r = r - alpha * q;
   s = a' * r;
   gamma = norm(s)^2;
   beta = gamma/gamold;
   p = s + beta .* p;
   disp(' Sum of squares of residuals:')
   disp( norm(r)^2 )
end