%  Regular Phase 2 of the Simplex Method  (file is reg.m)
%  Solves 'small' Linear Programming Problems (in canonical form)
%
%  (LP)  min cx s.t. ax=b, x>=0
%
%
%  Data input by user or calling routine: a,b,c,bas,pt 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.
%    pt  - pause time after each iteration, 0 denotes no pause, while
%                   any number > 20 denotes infinite time, i.e. you
%                   are prompted to hit return to continue.
%  The user can change the upper bound of 100 iterations - see iterm below.
%
%
%  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.
%      We do **not** use the revised simplex method.
%      A final check on roundoff error is made.
%      We use a threshold value when finding the pivot element.
%  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.
%
% 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]
% Now solve the LP by executing the file reg.m:
%   reg
%
%   Program by Henry Wolkowicz, University of Waterloo, Oct. 1986
%
% The Program:
%
%
rnderr=0;
iterm=100;
stop=1; % use to overcome the bug in the return statement
% Error tolerances (from 'Advanced Linear Progr.' by B.A.Murtaugh, pg 34.)
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
z=-c(bas)*b;    %  initial value for z
	      clc	
       	      home,disp([blanks(30)]),disp(['Initial tableau ' blanks(10)])
	      [a b
	       c z]	
		if pt > 20,
			disp(['Hit return to continue'])
			pause
		else
			pause(pt)
		end
%  price out the cost vector
z=-c(bas)*b;
clc
for i=1:m,
        c=c-c(bas(i))*a(i,:);
       	      home,disp([blanks(30)]),disp(['Price out the cost vector ' blanks(10) ])
	      [a b
	       c z]	
end
		clc
       	      home,disp([blanks(30)]),disp(['Price out the cost vector ' blanks(10) ])
	      [a b
	       c z]	
		if pt > 20,
			disp(['Hit return to continue'])
			pause
		else
			pause(pt)
		end
iter=0;         %  initialize the iteration count
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
%  Perform simplex iterations as long as there is a neg cost
while iter<iterm,
%  Find a negative reduced cost.
        ctemp=c;        % temporary work vector
        neg=[];
        for j=1:n,
                 if ctemp(nbas(j))<-eps1,
                        neg=[neg nbas(j)];
                 end
        end
        ct=-1;
   if length(neg)==0,
                          disp(['This phase is completed   - current basis is: '])
                           bas=bas
                          disp(['The current basic variable values are : '])
                          b
                          disp(['The current objective value is:'])
                          c0(bas)*b
                           disp(['The number of iterations is ' int2str(iter) ])
                           if norm(a0(:,bas)*b-b0,inf)>eps3, % check solution
                           disp(['**WARNING** roundoff error is significant'])
                           end
                           if any(b<-eps0), % check positive final solution
                             disp(['**WARNING** final b not nonnegative'])
                           end
			   pause(3)
              		   clc
       	      		   home,disp([blanks(30)]),disp(['Final tableau in this phase' blanks(10)])
              		   disp(['           ' blanks(10)])
			   [a b
	       		   c z]	
                           stop=0;
                           return
   else
              while ct<-eps1,   % continue till we find a suitable pivot
                [ct,i]=min(ctemp(neg));
                if ct>=-eps1,   % no suitable pivot columns are left
                        disp(['a suitable pivot element cannot be found'])
                        disp(['probable cause: roundoff error or ill-cond prob'])
                        disp(['equilibrate problem before solving'])
                        stop=0;
                        return
                end
                t=neg(i);       % index of the most neg reduced cost
%  Now, let x sub t enter the basis
%
%  First, we need to find the variable which leaves the basis
                pos=[];
                ind=[];
                for i=1:m,
                        if a(i,t)>eps0,
                                ind=[ind i]; % suitable rows
                        end
                end
                if length(ind)==0,
                        disp(['The problem is unbounded '])
                        stop=0;
                        return
                end
                [alpha,i]=min(b(ind)./a(ind,t));
                i=ind(i);       % pivot row
                if a(i,t)>eps2, % a suitable pivot element is found
                        ct=0;
                else
                        ctemp(t)=0;     % column t is unsuitable pivot col.
                end
              end
              if stop==0,
                return  % Ensure that we return
              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]	
		if pt > 20,
			disp(['Hit return to continue'])
			pause
		else
			pause(pt)
		end
%
		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]	
		if pt <= 20,
		pause(pt)
		else
		disp(['Hit return to continue'])
		pause
		end
    end
%
end
if iter>=100,
        text='Iteration bound has been exceeded ****** '
end
%   End of Program - Oct 16, 1986.
