function x=SimplexLP(c,A,b)

% Solves the LP
%
%  min c'*x s.t. A*x=b, x>=0
%
% by using the Simplex method
%
% S. Ulbrich, December 2006

n=length(c);
m=length(b);
msk=(b<0);
if any(msk)
 b(msk)=-b(msk);
 A(msk,:)=-A(msk,:);
end
Ah=[A,speye(m,m)];
ch=[0*c;ones(size(b))];
xh0=[0*c;b];
B0=[n+1:n+m];
[xhopt,Bopt]=simplex(Ah,b,ch,B0,xh0);
if ch'*xhopt>1e-12
 disp('Problem infeasible!');
 stop
end
if max(Bopt)<=n
    x0=xhopt(1:n);
    B0=Bopt;
else
    Bopt=[n+1:n+m];
    tmp=find(Bopt>n);
    N = setdiff(1:n,Bopt);
    jN=length(N)+1;
    for j=1:length(tmp)
     jN=jN-1;
     for k=1:jN
      Bopt(tmp(j))=N(k);
      if (rank(full(Ah(:,Bopt)))==m)
       N(k)=N(jN);
       break;
      end
     end
    end
end
[x,Bopt]=simplex(A,b,c,B0,x0);