function x = lsqr(a,b)
% Linear least squares solver
% Routine LSQR from Paige and Saunders to solve linear
% least-squares problems.
% ref: ACM Trans. Math. Software, v 8, pp 43 - 71, 1982.
%
% By D. Holmgren, Mar 13 94. holmgren@phobos.astro.uwo.ca
%
% works - tested on sample problem from Forsythe et al.
% minor error corrected 23/3/94
% need a tolerance for breaking out of iteration loop...
tol = 5.e-05;
nl = length(b);
ni = input('Number of iterations ');
nu = input('Number of unknowns ');
% initialization...
% most of these are work variables
x = zeros(nu,1);
u = zeros(nl,1);
v = zeros(nu,1);
u = b;
beta = norm(u);
u = u ./ beta;
v = a' * u;
alpha = norm(v);
v = v ./ alpha;
w = v;
phib = beta;
rhob = alpha;
% start bidiagonalization...
for i = 1:ni
   u1 = a * v - alpha .* u;
   beta = norm(u1);
   u = u1 ./ beta;
   v1 = a' * u - beta .* v;
   alpha = norm(v1);
   v = v1 ./ alpha;
% construct and apply next orthogonal transformation...
   rho = sqrt(rhob^2 + beta^2);
   c = rhob/rho;
   s = beta/rho;
   theta = s * alpha;
   rhob = -c * alpha;
   phi = c * phib;
   phib = s * phib;
% update x, w...
   x = x + (phi/rho) .* w;
   w = v - (theta/rho) .* w;
% compute sum of squares...   
   r = a * x - b;
   sumsq = norm(r)^2;   
   disp('Sum of squares of residuals:')
   disp( sumsq )
% break out if r is small enough.
% continuation beyond this can cause problems.
   if sumsq <= tol, break, end
end