function x = augm(A,b,alpha,N)
% AUGM  Solution of sparse linear least squares problems.
%
%       x=augm(A,b,alpha,N) solves the sparse linear least squares problem
%                    min ||Ax-b||
%                     x          2
%       using the augmented system approach with N steps of iterative
%       refinement. alpha is the scaling parameter used.

if nargin < 3, alpha=max(max(abs(A)))/1000; end
if nargin < 4, N=0; end

[m,n]=size(A);

C=[alpha*speye(m) A ; A' zeros(n)];
fg=[b; zeros(n,1)];

p=symmmd(C);
[L,U]=lu(C(p,p));
y=L\fg(p);
rx(p,1)=U\y;

for i=1:N,
  r=alpha*rx(1:m); x=rx((m+1):(m+n));
  fg=[b-A*x-r ; -A'*r/alpha];
  y=L\fg(p);
  drdx(p,1)=U\y;
  rx=rx+drdx;
end

x=rx((m+1):(m+n));