>> 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]

b =

     3
     2
     5
     1

>> c=[0;0;0;0;1;1;1;1]

c =

     0
     0
     0
     0
     1
     1
     1
     1

>> c=-c

c =

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

>> x = zeros(8,1)

x =

     0
     0
     0
     0
     0
     0
     0
     0

>> 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(sB) = B^(-1)*b

x =

     0
     0
     0
     0
     3
     2
     5
     1

>> b

b =

     3
     2
     5
     1

>> d = zeros(8,1);
>> cbar = c' - c(sB)'* B^(-1)*A

cbar =

     0     8    21     1     0     0     0     0

>> p=3

;

p =

     3

>> d(p)=1

d =

     0
     0
     1
     0
     0
     0
     0
     0

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

d =

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

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

ans =

    1.0000
    0.3333
    0.5556
    0.3333

>> lambda = ans(2)

lambda =

    0.3333

>> l = 2

l =

     2

>> q = sB(2)

q =

     6

>> x = x + lambda*d

x =

         0
         0
    0.3333
         0
    2.0000
         0
    2.0000
         0

>> sB(l) = p

sB =

     5     3     7     8

>> x = zeros

x =

     0

>> 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(sB) = B^(-1)*b

x =

  Columns 1 through 7

         0         0    0.3333         0    2.0000         0    2.0000

  Column 8

         0

>> x

x =

  Columns 1 through 7

         0         0    0.3333         0    2.0000         0    2.0000

  Column 8

         0

>> x=x'

x =

         0
         0
    0.3333
         0
    2.0000
         0
    2.0000
         0

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

cbar =

  Columns 1 through 7

    3.5000    1.0000         0    1.0000         0   -3.5000         0

  Column 8

         0

>> p = 1

p =

     1

>> d = zeros(8,1);
>> d(p)=1

d =

     1
     0
     0
     0
     0
     0
     0
     0

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

d =

    1.0000
         0
    0.1667
         0
   -1.5000
         0
   -1.5000
   -0.5000

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

ans =

    1.3333
   -2.0000
    1.3333
         0

>> x(3) 

ans =

    0.3333

>> x(3)+d(3)

ans =

    0.5000

>> d

d =

    1.0000
         0
    0.1667
         0
   -1.5000
         0
   -1.5000
   -0.5000

>> x(3)+1000*d(3)

ans =

   167

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

ans =

    1.3333
   -2.0000
    1.3333
         0

>> lambda = ans(4)

lambda =

     0

>> l = 4

l =

     4

>> q = sB(4)

q =

     8

>> x

x =

         0
         0
    0.3333
         0
    2.0000
         0
    2.0000
         0

>>