% doing revised simplex with bounds "manually".
% problem is: min 3x_1 + 4x_2 - 5x_3 + 7x_4
% subject to x_1 - x_2 + 2x_3 = 6,  -x_1 + 2x_2 + x_4 = 5, 
% 1 <= x_1 <= 2,
% 2 <= x_2 <= 6,
% -1 <= x_3 <= 4,
% 1 <= x_4
%
p=[3 4 -5 7]'; A=[1 -1 2 0; -1 2 0 1]; b=[6 5]';
% choose initial nonbasic set to be {1,2}, 
% with x_1 at upper bound  and x_2 at lower bound
N=[1 -2]; x(1)=2; x(2)=2;
% The nonbasic set is thus {3,4}, and we can calculate the
% corresponding values of x(3) and x(4) from the equality constraints
[L,U]=lu(A(:,B));
B=[3 4];
x(B)=U\(L\(b-A(:,abs(N))*x(abs(N))))

x =

     2
     2
     3
     3

% now compute reduced cost vector
[L,U]=lu(A(:,B));
u=L'\(U'\p(B)); c=p(abs(N))-A(:,abs(N))'*u

c =

   12.5000
  -12.5000


% BOTH nonbasic elements are eligible pivots. Since c_1>0, objective will
% decrease as x_1 is decreased below its upper bound, and since c_2<0, the
% objective will also decrease when x_2 is increased away from its lower bound.
% We set s=1 to try decreasing x_1.
s=1;
% calculate d, the negative of column 1 of the tableau.
d=U\(L\A(:,abs(N(s))))

d =

    0.5000
   -1.0000

% if x1=2-lambda, x3=3+.5*lambda and x4=3-lambda.  the first of these
% three variables to hit a bound is actually x_1, which hits its lower
% bound when lambda=1.  Let's take this step in x_1 and recalculate
% the basic variables, and make the relevant changes to N

N(1)=-1;
lambda=1; x(1)=2-lambda;
x(B)=U\(L\(b-A(:,abs(N))*x(abs(N))))

x =

    1.0000
    2.0000
    3.5000
    2.0000

% Now to the next step. Since B and abs(N) have not changed, c and u are the same.
c

c =

   12.5000
  -12.5000

% now the only eligible pivot column is s=2, since we can decrease the
% objective by allowing x_2 to increase from its lower bound, since
% c_N(2)<0.

s=2;
% calculate the negative of this tableau column
d=U\(L\(A(:,abs(N(s)))));
d

d =

   -0.5000
    2.0000

% when x_2=2+lambda, basic variables become x_3=3.5+.5*lambda and
% x4=2-2*lambda.  The first variable to hit a bound is x_2, when
% lambda=.5. So set r=2, since B(2)=4;

r=2;
lambda=0.5;
x(2)=x(2)+lambda;
x(B)=x(B)-lambda*d;
x

x =

    1.0000
    2.5000
    3.7500
    1.0000

% update B and N
B(2)=2; N(2)=-4;
B

B =

     3     2

N

N =

    -1    -4

% start another iteration. 
[L,U]=lu(A(:,B));
u=L'\(U'\p(B));
c=p(abs(N))-A(:,abs(N))'*u

c =

    6.2500
    6.2500

% Current x is optimal! Both nonbasic variables are at their lower
% bound, and in both cases the objective will increase if we increase
% them away from their lower bound.  calculate optimal objective

z=p'*x

z =

    1.2500

diary off