function x = cne(A,b)
%SSNE	Solution of sparse linear least squares problems.
%       @(#)cne.m Version 1.3 5/29/92
%       Pontus Matstoms, Linkoping University.
%       e-mail: pomat@math.liu.se
%
%       x=cne(A,b) solves the sparse linear least squares problem
%                    min ||Ax-b||
%                     x          2
%       using the normal equation A'Ax=A'b by Cholesky factorization
%       of A'A. The accuracy of the computed solution is improved by
%       a few steps of iterative refinement.

% Compute a minimum degree ordering and the Cholesky factor.

p=colmmd(A);
R=chol(A(:,p)'*A(:,p));

% Compute the ne-solution and do iterative refinement.

ATb=A(:,p)'*b;
y=R'\ATb;
x=R\y;

for i=1:1,
   r=b-A(:,p)*x;
   ATb=A(:,p)'*r;
   y=R'\ATb;
   dx=R\y;
   x=x+dx;
end

% Permute the computed solution

n=size(A,2);
pinv(p)=1:n;
x=x(pinv);
