%  Simplex Method  (Phase 1 and 2)
%  Solves the Linear Programming Problem
%  (LP)  min cx s.t. ax<=and=b, x>=0
%  The file reg.m is called to do phase 1 and phase 2.
%  See the file   templat1.m for an example.
%
%  Data input by user or calling routine: a,b,c,m1,m2 where
%    a   - m by n1 matrix containing inequality constraints in the
%               first m1 rows and equality constraints in the next m2 rows.
%    b   - m-column vector assumed >=0
%    c   - n1-row vector
%    m1,m2 - 2 nonnegative integers with the number of ineqality and equality
%               constraints respectively.
%    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.
%
%  (Program by Henry Wolkowicz - Dept of Comb & Opt, Univ of Waterloo- Oct/86)
%                                             (updated Mar, 1989)
%
nbas=[];         %  initialize to avoid compiler definition error
eps4=.00001;    %  accuracy parameter
eps0=10^(-10);  %  numerical zero
eps1=10^(-5);  %  accuracy parameter for optimality check
z=0;            %  initial objective value
if any(b<-eps0),
        disp('ERROR -the rhs vector b is *not* nonnegative!')
        return
end
[m,n1]=size(a);      %  number of rows and columns of a
clc
home,disp([blanks(20)]);
disp(['You have entered the following linear programming data'])
disp(['Constraint matrix a is ' int2str(m) ' by ' int2str(n1) ])
disp(['There are ' int2str(m1) ' inequality (<=) and ' int2str(m2) ' equality constraints'])
if pt <= 20,
	disp(['The pause time between iterations is ' int2str(pt) ' seconds'])
else
	disp(['The program stops after each iteration'])
end
disp(['Hit the return key to begin'])
%pause
if m1~=0,       % inequality constraints exist
        a=[a [eye(m1) zeros(m1,m2)]'];
        c=[c zeros(1,m1)];
	disp(['Slack variables are added'])
end
if m2==0,       % if no equality constraints
        bas=[n1+1:n1+m1];      % slack variables are basic
        disp(['start phase2 '])
        reg;
	clc
	home,disp([blanks(20)]),disp(['Final tableau'])
	[a b
         c z]
	pause(3)
        return
else
        disp(['solve phase 1 first '])
        corig=c;
        c=[zeros(1,n1+m1) ones(1,m2)];
        a=[a [zeros(m2,m1) eye(m2)]'];
        bas=[n1+1:n1+m1+m2];
	disp(['Artificial variables are added.'])
        reg;        % solve phase 1 using the reg.m file
        if z<-eps4,
                disp(['optimal value from phase 1 is: ' num2str(z) ])
                disp(['the above shows that the problem is infeasible'])
		pause(3)
		disp(['Final tableau'])
		[a b
         	c z]
                return
        else
                disp(['phase 1 is solved now solve phase 2'])
		pause(3)
                a=a(:,1:n1+m1);
                c=c(1:n1+m1);
%  Ensure that no artificial variable remains in the basis.
                while ~all(bas<n1+m1+1),
                disp(['an artificial variable remains in the basis after phase1'])
                disp(['pivot to remove the remaining artificial variables'])
                        mto1=[1:m];
                        i=mto1(bas>n1+m1); % pivot row
                        i=i(1)
                        n1m1to1=[1:n1+m1];
                        t=n1m1to1(abs(a(i,:))>eps1);    % pivot column
                        t=t(1)
                        v=nbas(nbas==t);  % variable entering the basis
                        nbas(v)=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
%
		clc,home,disp([blanks(20)])
                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(20)]),disp(['pivot= a(' int2str(i) ',' int2str(t) ')' blanks(10)])
		      [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(20)]),disp(['pivot= a(' int2str(i) ',' int2str(t) ')' blanks(10)])
	      [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(20)]),disp(['pivot= a(' int2str(i) ',' int2str(t) ')' blanks(10)])
	      [a b
	       c z]	
		if pt > 20,
			disp(['Hit return key to continue'])
			pause
		else
			pause(pt)
		end
              end
		disp(['The artificial variables have been removed'])
		pause(3)
                c=corig;
                reg  % solve the problem using file reg.m
		disp(['Final tableau'])
		[a b
         	c z]
        end
end
%   End of Program - Oct 16, 1986.
