function x = csne(A,b)
%CSNE	Solution of sparse linear least squares problems.
%       @(#)csne.m Version 1.4 6/23/93
%       Pontus Matstoms, Linkoping University.
%       e-mail: pomat@math.liu.se
%
%       x=csne(A,b) solves the sparse linear least squares problem
%                    min ||Ax-b||
%                     x          2
%       using the semi-normal equation R'Rx=A'b. R is here
%       the upper triangular matrix R in the QR factorization of A.
%       The accuracy of the computed solution is improved by a few 
%       steps of iterative refinement.

% Minimum-degree ordering of A.

q=colmmd(A);
A=A(:,q);

% Compute the QR factorization of A.

[R,p]=sqr(A);

% Resulting column permutation.

A=A(:,p);
Pc=q(p);

% Compute the SNE-solution.

ATb=A'*b;
y=R'\ATb;
x=R\y;

% One step of iterative refinement.

r=b-A*x;
ATb=A'*r;
y=R'\ATb;
dx=R\y;
x=x+dx;

% Permute the computed solution.

x(Pc)=x;