function [x,y,z,A,b,c,n0]=ProgLineal_InP_MPS_sc(file)
% Resuelve min c'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
%
 eval(['!copy ' file ' ' 'D:\MATLAB2013a\work\tmp\in.mps']);   % Lee file(MPS)
 LIPSO = 'D:\MATLAB2013a\work'; ar=mps2mat(LIPSO); fprintf('\n');
 load D:\MATLAB2013a\work\tmp\default;
 [m,n]=size(A); x=sqrt(n)*ones(n,1); y=zeros(m,1);
 z = x; bc = 1+max(norm(b),norm(c)); t0=cputime;
 p = symamd(A*A');                      % Reordenación grado mínimo
 fprintf('\n              No fac.    No fac.    Brecha     Error');
 fprintf('\n              primal      dual       dual      relati.\n');
 fprintf('  Iter        A*x-b     A''*y+z-c     x''*z      total\n');
 fprintf('  ---------------------------------------------------\n');
%
 for iter = 1:100
    Rd=A'*y+z-c; rrd=norm(Rd);       % No factibilidad dual
    Rp=A*x-b;    rrp=norm(Rp);       % No factibilidad priaml
    Rc=x.*z;                         % No factibilidad complementariedad
    residuo=norm([Rd;Rp;Rc])/bc;
    gap=mean(Rc);                    % Brecha dual
    fprintf('%5i %15.2e %10.2e %11.2e %10.2e\n',iter,rrd,rrp,gap,residuo);
    if residuo<5.e-10, break, end
    Rc=Rc-min(.1,100*gap)*gap;
    d=min(5.e+15,x./z);
    B=A*sparse(1:n,1:n,d)*A'; 
    R=chol(sparse(B(p,p)));          % Cholesky 
    t1=x.*Rd - Rc; 
    t2=-(Rp+A*(t1./z)); 
    dy=zeros(m,1); 
    dy(p)=R\(R'\t2(p));              % Resuelve sistema 
    dx=(x.*(A'*dy)+t1)./z;
    dz=-(z.*dx+Rc)./x;
    tau=max(.9995,1-gap);            
    ap=-1/min(min(dx./x),-1);        % Amplitud de paso primal
    ad=-1/min(min(dz./z),-1);        % Amplitud de paso dual
    ap=tau*ap;
    ad=tau*ad;
    x=x+ap*dx; z=z+ad*dz; y=y+ad*dy;
 end
 fprintf('Termina bien!\t[m n] = [%g %g]\tCPU = %g\n',m,n,cputime-t0);
 f=c'*x;
 fprintf('Fun. Obje.= %18.14f\n',f);
end