function [fobj,niters,x] = Rsimplex_1(file,A,b,c,lb,ub)
% Método Simplex revisado en dos FASES
%    Resuelve min c'x st Ax=b, lb<=x<=ub.
%
 if nargin>2
    if nargin<4, error('Si no se lee de fichero, se requieren al menos',...
                      ' cuatro argumentos en entrada'); end
    if nargin<6, ub = Inf*ones(size(c)); end
    if nargin<5, lb = zeros(size(c)); end
    A = sparse(A); c=c(:); b=b(:); lb=lb(:); ub=ub(:);
 else                          % Lectura de los datos desde fichero 'file' en MPS
    eval(['!copy ' file ' d:\Matlab2013a\work\tmp\in.mps']); clear all
    LIPSO = 'd:\Matlab2013a\work'; fin=mps2mat(LIPSO); fprintf('\n');
    load d:\Matlab2013a\work\tmp\default;   % En default.mat los datos 
    lb=lbounds; ub=ubounds;
 end
 [m,n] = size(A); zer_tol = 1.0e-5; piv_tol = 1.0e-8; niters=0; t=cputime;
 f=fopen('Salida_Rsimplex','w');
 fprintf(f,' It.  Fase    Sinf/Fobj      Nopt\n');
 fprintf(f,'----------------------------------\n');
% FASE I
 v    = zeros(n,1);
 lbdd = find(lb>-Inf);      v(lbdd) = lb(lbdd); free = setdiff(1:n,lbdd);
 ubdd = free(ub(free)<Inf); v(ubdd) = ub(ubdd);
 if isempty(ubdd), N = -lbdd';      % Var no en la base: en cota inf. No hay en u
     else N = [-lbdd' ubdd]; free = setdiff(free,ubdd);
 end
 f = length(free);
 if f>0                              % ¡OJO! Hay variables libres
    [L,U,P,Q] = lu(A(:,free));       % Q es matriz de permutación
    free = free*Q';
    i = length(find(abs(diag(U))>piv_tol)); j = free(i+1:f);
    if ~isempty(j)
       N = [N j]; free = setdiff(free,j); f = length(free);
       lb(j) = zeros(size(j)); ub(j) = zeros(size(j));
    end;
    [k,~] = find(P(1:f,:));          % relies on diag(U) having zeros at end
    v(free) = U(1:f,1:f)\(L(1:f,1:f)\(P(1:f,:)*(b - A*v)));
    k = setdiff(1:m,k); m = m-f; B = [free n+1:n+m];
    A = [A sparse(k,1:m,sign(b(k)-A(k,:)*v+eps),m+f,m)];
 else                                % No hay variables libres 
    B = n+1:n+m; j = []; A = [A sparse(1:m,1:m,sign(b-A*v+eps*ones(size(b))))];
 end
 lb = [lb; zeros(m,1)]; ub = [ub; Inf*ones(m,1)]; 
 w  = [zeros(n,1); ones(m,1)];       % Variables + artificiales para FASE I
 [x,niters,B,N,L,U] = rsmbdd(A,b,w,lb,ub,B,N,1,niters);            % FASE I
 if w'*x>zer_tol, error('El problema no es factible'); end
% Adaptar para FASE II
 ub(n+1:n+m) = zeros(m,1); c = [c; zeros(m,1)];
 if ~isempty(j) 
    c_r = c(j)'-c(B)'*(U\(L\A(:,j)));  
    if norm(c_r,inf)>zer_tol, error('Problema no acotado'); end
 end
 [x,niters] = rsmbdd(A,b,c,lb,ub,B,N,2,niters);          % FASE II
 x = x(1:n); fobj=c(1:n)'*x;
 fclose('all');
 fprintf('[m n]=[%g %g]. Tiempo CPU:%10.4f; fun. obj.=%17.4e\n',m,n,cputime-t,fobj);
end

function [x,niters,B,N,L,U] = rsmbdd(A,b,c,lb,ub,B,N,fase,niters)
% Rutina de SIMPLEX revisado para min c'x s.a Ax=b, lb<=x<=ub.
% B vector 1xm con los índices de las variables (columnas) en la base.
% N vector 1x(l-m) con los índices de las variables no básicas,
%      es decir, en sus límites (+/-1).
%
 f=fopen('Salida_Rsimplex','a');
 zer_tol = 1.0e-5; piv_tol = 1.0e-8;
 [L,U] = lu(A(:,B));
 if any(abs(diag(U))<piv_tol), error('La base inicial no es invertible'); end
 x = lb;  upper = N(N>0); x(upper) = ub(upper);
 x(B) = U\(L\(b-A(:,abs(N))*x(abs(N)))); 
 while 1
    if any(x(B)<lb(B)-zer_tol | x(B)>ub(B)+zer_tol), error('Punto no factible'); end
    u = L'\(U'\c(B)); 
    c_r = c(abs(N))'-u'*A(:,abs(N));                % Costes reducidos no básicas
    if ~any(c_r.*N>zer_tol), return, end            % SE HA LLEGADO A UNA SOLUCIÓN
    Nopt=sum(c_r.*N>zer_tol);
    [~,s] = max(c_r.*sign(N));               % Entra var. s no básica con>coste reducido
    blk = abs(N(s)); 
    min_ratio = ub(blk) - lb(blk); order = -sign(N(s)); % Puede s ir de u a l o l a u
    d = order*U\(L\A(:,blk));
    block_lb = find(d>=piv_tol & lb(B)>-1e32); % Comprobar básica que llega a lb
    if ~isempty(block_lb)
      [min_lb,index_r] = min((x(B(block_lb))-lb(B(block_lb)))./d(block_lb));
      if min_lb<min_ratio
         r=block_lb(index_r);
         min_ratio = min_lb; 
         blk = B(r); 
      end
    end
    block_ub = find(d<=-piv_tol & ub(B)<1e32); % Comprobar básica que llega a ub
    if ~isempty(block_ub)
      [min_ub,index_r] = min((x(B(block_ub))-ub(B(block_ub)))./d(block_ub));
      if min_ub<min_ratio
         r=block_ub(index_r); 
         min_ratio = min_ub; 
         blk = B(r); 
      end
    end
    if min_ratio>1e32, error('Problema no acotado'); end
  % Readaptar solución y datos
    x(B) = x(B)-min_ratio*d;
    x(abs(N(s))) = x(abs(N(s))) + order*min_ratio;
    if blk == abs(N(s)), N(s) = -N(s);     % Se bloquea la no básica
       else swap = B(r); B(r) = abs(N(s)); N(s) = -sign(d(r))*swap;
          [L,U] = lu(A(:,B)); 
    end
    niters=niters+1;
    fprintf(f,'%5d%4d%18.9e%6d\n',niters,fase,c'*x,Nopt);
 end
end