function [x_opt,B] = simplex (A, b, c, B, x)

% Solves the LP
%
%  min c'*x s.t. A*x=b, x>=0
%
% for an admissible starting basis B
% and an admissible solution x
%
% by using the Simplex method
%
% Ch. Brandenburg, December 2005
% S. Ulbrich, December 2006


tol=1e-12;
N = setdiff(1:size(x,1),B);

% standard simplex algorithm
maxit = 1000000;
it = 1;
z = zeros (length(x),1);

while (it <= maxit)

    cTx=c'*x

    % (1) BTRAN
    % y' * A_B = c_B'

    y = A(:,B)'\c(B);


    % (2) Pricing
    % z_N = c_N - A_N'*y

    z(N) = c(N) - A(:,N)'*y;

    % check for optimality
    if (sum (z(N) < -tol) == 0)
        x_opt = x;
        break;
    end

    % choose entering variable
    % j in N with z_j < 0
    minz=min(z(N));
    tmp = find(z(N)==minz);
    pos_j = find (N == min(N(tmp)));
    j = N(pos_j);
    

    % (3) FTRAN
    % A_B*w = A_.j

    w = A(:,B)\full(A(:,j));


    % (4) Ratio-Test

    % check if LP is bounded
    % ie w <= 0
    msk=(abs(w)<tol);
    w(msk)=0;
    if (sum (w > 0) == 0)
        disp ('LP is unbounded');
        break;
    end

    % find gamma
    % C : minimal values
    % I : indeces

    % find gamma and i
    tmp = find (w>0); % indeces k of w with w(k) > 0
    [C, I] = min (x(B(tmp)) ./ w(tmp));
    I = tmp(I); % indeces i of w for which x_i/w_i is minimal
    gamma = C(1);
    
    % selection rule:
    i = find(B == min(B(I)));


    % (5) Update

    x(B) = x(B) - gamma*w;
    x(B(i))=0;
    N(pos_j) = B(i);
    B(i) = j;
    x(j) = gamma;
    msk=(x(B)<tol);
    x(B(msk))=0;

    it = it + 1
end