function x = simplex(A,b,c,I)

% function x = simplex(A,b,c,I)
% solvex the linear program 
% c'x -> min, Ax=b, x>=0
% starting with an admissible basic solution 
% with indices I

[m,n] = size(A); 
E = eye(m,m);
T = inv(A(:,I)); 
T = [T, T*b];

for i = 1:nchoosek(n,m);
   K = setdiff(1:n,I);
   d = c(K) - A(:,K)'*T(:,1:m)'*c(I); 
   [dmin,k] = min(d);
   if d(k) >= 0;
      x = zeros(n,1); x(I) = T(:,m+1);
      return;
   end;
   k = K(k);
   z = T(:,1:m)*A(:,k);
   J = find(z>0);
   if length(J) == 0;
      display('no solution');
      x = [];
      return;
   end;
   [zxmin,j] = min(T(J,m+1)./z(J));
   j = J(j);
   T = (E + (E(:,j)-z)*E(j,:)/z(j)) * T;
   I(j) = k;  
end;

display('no solution after maximal number of steps');

