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; endif 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;endx=rx((m+1):(m+n));