function [x,eta,d,B,N] = simplex(d)
    % This function solves the LP problem
    % maximize eta = c' * x, subject to Ax <= b, x >= 0,
    % using the simplex algorithm.
    % Outputs:
    % x: solution vector, x = [x_1,...,x_n]'.
    % eta: optimal value of objective function
    % d: Optimal dictionary
    % B: Indices of basic variables in d
    % N: Indices of nonbasic variables in d
    % Test from section 5.6:  [x,eta,d,B,N]=simplex([-2 -1;-2 4;-1 3],[4; -8; -7],[-1;-1]);
  
    m = size(d,1)-1; n = size(d,2)-1;
    c = d(1,2:end)';
    if sum(d(2:end,1)<0) == 0
        % We can make an initial feasible dictionary.
        B=[n+1:n+m]';
        N=[1:n]';
    else
        % Replace the objective with -x_1-...-x_n, to obtain a dual feasible dictionary
        dualdict = -d';
        dualdict(2:end,1) = 1;
        [x_0,eta_0,d_0,B_0,N_0]=simplex(dualdict);
        assert(eta_0 ~= inf, 'The problem is infeasible') % Unbounded dual implies infeasible primal
      
        % Convert dual dictionary to primal
        d = -d_0'; % The negative transpose
        I=[(n+1:n+m) (1:n)]'; % Used to map dual variables to their complementary primal variables. 
        B=I(N_0); N=I(B_0);
        % Replace the objective function by the original one, see the notes
        c_ext = [c;zeros(m,1)];
        C_B = c_ext(B);
        C_N = c_ext(N);
        d(1,1) = C_B'*d(2:(m+1),1);
        d(1,2:n+1) = C_B'*d(2:(m+1),2:(n+1)) + C_N';
    end

    MaxIt=10000; it = 1;
    tol = 1.0e-12; % tolerance used for sign tests.
    while it < MaxIt
        inds = find(d(1,2:end)>tol); % Choice set for entering variable
        if isempty(inds) % Dictionary is optimal.
            break;
        else
            [t,q] = min(N(inds)); % Blands rule for choosing entering variable
            r = inds(q); % r is entering
            
            ratios = -d(2:end,r+1)./d(2:end,1);
            % Apply the covention 0/0=0. Matlab computes this to NaN otherwise
            inds = find(( abs(d(2:end,r+1)) < tol) & ( abs(d(2:end,1)) < tol));
            ratios(inds)=0;
            
            val= max(ratios);
            if val<=tol % The problem is unbounded
                eta = inf; break;
            else
                inds = find(ratios==val); % Choice set for leaving variable
                [t,q] = min(B(inds)); % Blands rule for choosing leaving variable
                s = inds(q); % s is leaving
            
                d=pivot(d,r,s);
                temp=B(s);
                B(s)=N(r);
                N(r)=temp;
                it = it+1;
            end
        end
    end

    assert(it < MaxIt, 'Maximum number of iterations exceeded.')

    x_ext = zeros(n+m,1);
    x_ext(B) = d(2:end,1);
    x = x_ext(1:n);
    eta = d(1,1);
end