function IntPointLP_2(A,b,c)
% Resuelve min p'x s.t. Ax=b, x>=0,
%    cuyo dual es max b'y  s.t. A'y+z=c,  z>=0. Método Primal-Dual IPt
%
 [m,n]=size(A); x=sqrt(n)*ones(n,1); y=zeros(m,1); % Solución inicial factible
 z = x; bc = 1+max(norm(b),norm(c)); t0=cputime;
 for iter = 1:100
    Rd=A'*y+z-c;
    Rp=A*x-b;
    Rc=x.*z;
    residual=norm([Rd;Rp;Rc])/bc;
    gap=mean(Rc);
    fprintf('iter %2i: medidor duali.=%9.2e, error tot. rel.=%9.2e',iter,gap,residual);
    fprintf('\tobj=%14.6e\n',c'*x);
    if residual<5.e-8, break, end
    Rc=Rc-min(.1,100*gap)*gap;
    d=min(5.e+15,x./z);
    B=A*diag(d)*A';
    R=chol(B);
    t1=x.*Rd-Rc;
    t2=-(Rp+A*(t1./z));
    dy=R\(R'\t2);
    dx=(x.*(A'*dy)+t1)./z;
    dz=-(z.*dx+Rc)./x;
    tau=max(.9995,1-gap);
    ap=-1/min(min(dx./x),-1);
    ad=-1/min(min(dz./z),-1);
    ap=tau*ap;
    ad=tau*ad;
    x=x+ap*dx; z=z+ad*dz; y=y+ad*dy;
 end
 fprintf('Termina!\t[m n] = [%g %g]\tCPU = %g\n',m,n,cputime-t0);
 x,y,z
 fprintf('Fun, Obj= %14.6f\n',c'*x);
end
