%  Pivot the tableau (file is pivot.m)
%  Performs one pivot operation on given data
%
%  (LP)  min cx s.t. ax=b, x>=0
%
%
%  Data input by user or calling routine: a,b,c,bas where
%    a   - m by (m+n) matrix (containing identity)
%    b   - m-column vector assumed >=0
%    c   - m+n-row vector of costs which is priced out below using bas
%    bas - m-row vector with column indices corresp. to the
%          identity matrix in a, i.e. a(bas,:)=identity.
%    i    - pivot row
%    t    - pivot column
%
%
%  The matrix a is assumed to contain an m by m identity matrix
%      corresponding to the basic columns.
%  A relative accuracy of approximately 15 significant decimal
%      digits is assumed. This affects the values of 4 accuracy
%      parameters, epsi, i=0,1,2,3.
%  The pivot step is done using 'Gauss-Jordan' elimination.
%      No special factorizations to ensure stability are used.
%  Problem (LP) is assumed to be in canonical form, i.e. slacks
%    have been added and/or phase 1 has finished. However, we still
%    price out the cost vector c.
%  The original data is saved in cO,a0,b0,bas0.
%
% Example: Enter the data:
%   a=[1 -2 1 1 0 3
%      2 1 0 -1 1 1]
%   b=[2 3]'
%   c=[1 2 2 -1 3 3]
%   bas=[3 5]
%   t=4
%   i=1
% Now pivot the LP by executing the file pivot.m:
%   pivot
%
%   Program by Henry Wolkowicz, University of Waterloo, June. 1989
%
% The Program:
%
%
rnderr=0;
stop=1; % use to overcome the bug in the return statement
eps0=10^(-10);  %  numerical zero
eps1=10^(-5);   %  accuracy parameter for optimality check
eps2=10^(-8);   %  accuracy parameter pivot element (threshold test)
eps3=10^(-6);   %  accuracy parameter for final roundoff error check
a0=a;           %  save the matrix a for the final roundoff error test
b0=b;           %  save the vector b for the final roundoff error test
c0=c;bas0=bas;
[m,mn]=size(a); %  row and column size of a
% Ensure that the column entering the basis is not basic
if ~(all(t~=bas)),
     clc
     home,disp(['ERROR*****variable entering basis is already basic'])
     break
end
z=-c(bas)*b;    %  initial value for z
	      clc	
       	      home,disp([blanks(30)]),disp(['Initial tableau ' blanks(10)])
	      [a b
	       c z]	
            pause(2)
%  price out the cost vector
z=-c(bas)*b;
clc
for ii=1:m,
        c=c-c(bas(ii))*a(ii,:);
       	      home,disp([blanks(30)]),disp(['Price out the cost vector ' blanks(10) ])
	      [a b
	       c z]	
end
              pause(2)
		clc
       	      home,disp([blanks(30)]),disp(['Price out the cost vector ' blanks(10) ])
	      [a b
	       c z]	
              pause(2)
n=mn-m;         %  number of nonbasic variables
%  nbas - indices of the nonbasic variables
nbas=[];
for j=1:mn,
       if all(j~=bas),
                nbas=[nbas j];
        end
end
%  Update the basic and nonbasic vectors.
                nbas(nbas==t)=bas(i);
                bas(i)=t;
                alpha=a(i,t);   % pivot element
%  Store the data in ap,bp
                ap=a;
                bp=b;
%  Now pivot by row
		iter=iter+1;
	      clc	
             home,disp([blanks(30)]),disp(['pivot= a(' int2str(i) ',' int2str(t) ')' blanks(10)])
              disp(['tableau before pivot'])
	      [a b
	       c z]	
              pause(2)
%
		clc
                for k=1:m,
                      ratio=ap(k,t)/ap(i,t);
                      a(k,:)=ap(k,:)-ap(i,:)*ratio;
                      b(k)=bp(k)-bp(i)*ratio;
       	      home,disp([blanks(30)]),disp(['pivot= a(' int2str(i) ',' int2str(t) ')' blanks(10)])
              disp(['pivoting                 '])
	      [a b
	       c z]	
                end
%  Now for the objective row update
                ratio=c(t)/ap(i,t);
                c=c-ap(i,:)*ratio;
                z=z-bp(i)*ratio;
       	      home,disp([blanks(30)]),disp(['pivot= a(' int2str(i) ',' int2str(t) ')' blanks(10)])
              disp(['pivoting                 '])
	      [a b
	       c z]	
%  Now for the pivot row update
                a(i,:)=ap(i,:)/ap(i,t);
                b(i)=bp(i)/ap(i,t);
       	      home,disp([blanks(30)]),disp(['pivot= a(' int2str(i) ',' int2str(t) ')' blanks(10)])
              disp(['pivoting                 '])
	      [a b
	       c z]	
%
end
disp(['new tableau in:  [a b ; c z] '])
%   End of Program - June 1989
