>> A = [1 2 3 0 1 0 0 0;-1 2 6 0 0 1 0 0;0 4 9 0 0 0 1 0;0 0 3 1 0 0 0 1]

A =

     1     2     3     0     1     0     0     0
    -1     2     6     0     0     1     0     0
     0     4     9     0     0     0     1     0
     0     0     3     1     0     0     0     1

>> b=[3;2;5;1];
>> c=[0;0;0;0;1;1;1;1];
>> c=-c

c =

     0
     0
     0
     0
    -1
    -1
    -1
    -1

>> sB=[5 6 7 8]

sB =

     5     6     7     8

>> B = A(:,sB)

B =

     1     0     0     0
     0     1     0     0
     0     0     1     0
     0     0     0     1

>> x=zeros(8,1);x(sB)=B^-1*b

x =

     0
     0
     0
     0
     3
     2
     5
     1

>> cbar = c'-c(sB)'*B^-1*A

cbar =

     0     8    21     1     0     0     0     0

>> p = 3

p =

     3

>> d = zeros(8,1); d(sB)=-B^-1*A(:,p)

d =

     0
     0
     0
     0
    -3
    -6
    -9
    -3

>> x(sB) ./ -d(sB)

ans =

    1.0000
    0.3333
    0.5556
    0.3333

>> l=2;q=sB(l)

q =

     6

>> sB(l)=p

sB =

     5     3     7     8

>> x=zeros(8,1);x(sB)=B^-1*b

x =

     0
     0
     2
     0
     3
     0
     5
     1

>> B = A(:,sB)

B =

     1     3     0     0
     0     6     0     0
     0     9     1     0
     0     3     0     1

>> x=zeros(8,1);x(sB)=B^-1*b

x =

         0
         0
    0.3333
         0
    2.0000
         0
    2.0000
         0

>> cbar = c'-c(sB)'*B^-1*A

cbar =

    3.5000    1.0000         0    1.0000         0   -3.5000         0         0

>> p=1

p =

     1

>> d = zeros(8,1); d(sB)=-B^-1*A(:,p)

d =

         0
         0
    0.1667
         0
   -1.5000
         0
   -1.5000
   -0.5000

>> x(sB) ./ -d(sB)

ans =

    1.3333
   -2.0000
    1.3333
         0

>> l = 4;q=sB(l)

q =

     8

>> sB(l)=p

sB =

     5     3     7     1

>> B = A(:,sB)

B =

     1     3     0     1
     0     6     0    -1
     0     9     1     0
     0     3     0     0

>> x=zeros(8,1);x(sB)=B^-1*b

x =

         0
         0
    0.3333
         0
    2.0000
         0
    2.0000
         0

>> B = A(:,sB)

B =

     1     3     0     1
     0     6     0    -1
     0     9     1     0
     0     3     0     0

>> cbar = c'-c(sB)'*B^-1*A

cbar =

     0     8     0    -6     0     0     0    -7

>> p=2

p =

     2

>> d = zeros(8,1); d(sB)=-B^-1*A(:,p)

d =

     2
     0
     0
     0
    -4
     0
    -4
     0

>> x(sB) ./ -d(sB)

ans =

    0.5000
      -Inf
    0.5000
         0

>> l=1;q=sB(l)

q =

     5

>> sB(l)=p

sB =

     2     3     7     1

>> B = A(:,sB)

B =

     2     3     0     1
     2     6     0    -1
     4     9     1     0
     0     3     0     0

>> x=zeros(8,1);x(sB)=B^-1*b

x =

    1.0000
    0.5000
    0.3333
         0
         0
         0
         0
         0

>> cbar = c'-c(sB)'*B^-1*A

cbar =

     0     0     0     0    -2    -2     0    -1

>>