%lpnewt1: Solve primal LP: min c'x s.t. Ax=<b 
%Via Newton for least 2-norm dual LP: max -b'v s.t. -A'v=c, v>=0
%Input: c,A,b,epsi,delta,tol,itmax;Output: v (l2norm dual sol), z primal sol
epsi=1e-3;tol=1e-12;delta=1e-4;itmax=100;%default inputs
pl=inline('(abs(x)+x)/2');%pl(us) function
tic;i=0;z=0;y=((A(1:n,:))'*A(1:n,:)+epsi*eye(n))\(A(1:n,:))'*b(1:n);%initial y
while (i<itmax & norm(y-z,inf)>tol & toc<1800)
    df=A'*pl((A*y-b))+epsi*c;
    d2f=A'*spdiags(sign(pl(A*y-b)),0,m,m)*A+delta*speye(n);
    z=y;y=y-d2f\df;
    i=i+1;
end
toc1=toc;v=pl(A*y-b)/epsi;t=find(v);z=A(t,:)\b(t);toc2=toc;
format short e;[epsi delta  tol i-1 toc1 toc2 norm(x-y,inf) norm(x-z,inf)]