clear()
////////////////////////////////////////////////
///
/// simplex programming: benefit maximization
///
////////////////////////////////////////////////
//
//
m=3 // number of ingredients
p=2 // number of sold products
B=[5,4;1,2;3,2]; // amount needed for each ingredient and each product
c1=[40;50];// unitary margin of each product
b=[80;24;36]; // total available amount of each ingredient
//
//
//////////the case of 2012 exam /////////////////
B=[1,-1,1;2,2,1;4,1,5];
b=[1;10;12];
c1=[5;-2;1];
m=3;p=3;
//////////////////////////////////////////////////
//
// canonical form
//
n=m+p
A=[B,eye(m,m)];
c=[c1;zeros(m,1)];
//
//
disp('successive extreme point and associated optimal value:');
Jout=1:p;J=(p+1):n;
X=zeros(n,1);X(J)=inv(A(:,J))*b;
disp(X);disp(0);
//
d=c(Jout)'-(c(J)')*inv(A(:,J))*A(:,Jout)
//
while max(d)>0
  [dmax,iin]=max(d);A1=inv(A(:,J))*A(:,Jout);
  btilde=A1(:,iin);
  bprime=inv(A(:,J))*b;
  bhat=bprime./(btilde+1E-12);      // to avoid division by 0
  [d,ii]=find(bhat<0);bhat(ii)=1E10;// to avoid negative values
  [alpha,iout]=min(bhat);
  //J(iout) goes out,  Jout(iin) goes in
  buffer=J(iout);J(iout)=Jout(iin);Jout(iin)=buffer;
  d=c(Jout)'-(c(J)')*inv(A(:,J))*A(:,Jout);
  X=zeros(n,1);X(J)=inv(A(:,J))*b;
  disp(X);disp(c(J)'*X(J)); 
end  