% Title      : Linear Programming
%%  Originally, this code was programmed by Prof. Henry Wolkowicz 
%  (Univ. of Waterloo) and modified by authors.

% Modified by   Hojong Baik   (231-3363, hbaik@vt.edu)
% 				Dr. Antonio A. Trani (231-4418, vuela@vt.edu)
%				Dr. Donald R. Drew (drdrew@vt.edu)
% Date       : Nov. 29, 1998


%  The Simplex Method 
%  Solves 'small' Linear Programming Problems (in canonical form)
%

%
%  (LP)  max cx s.t. ax=b, x>=0
%
%
%  Data input by user or calling routine: a,b,c,bas,pt where
%    minmax=1, if minimizing problem
%                   0, if maximizing problem
%    a   - m by (m+n) matrix (containing identity)
%    b   - m-column vector assumed >=0
%    c   - m+n-row vector of costs which is placed left hand side(LHS)
%    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.
%
%
% datafile :
osaka2

disp(['If the input data is correct, hit return to continue'])
                 pause
if minmax==1,
     c=-c;
end

echo off
%
pt=1

%
% 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=0;    %  initial value for z
             clc
             home,disp([blanks(30)]),
             disp(['Initial tableau ' blanks(10)])

           if minmax==1
	      tc=-c;
              [tc z
               a b]
           else
              [c z
                a b]
           end

            if pt > 20,
                 disp(['Hit return to continue'])
                 pause
            else
                 pause(pt)
            end
%  make cannonical form

if minmax==1
  z=-c(bas)*b;
else
  z=c(bas)*b;
end

clc
for i=1:m,
        c=c-c(bas(i))*a(i,:);                  
end
            clc
                  home,disp([blanks(30)]),
             %disp(['Price out the cost vector ' blanks(10) ])
             disp(['standard or canonical form of this problem' ])
 
          if minmax==1
	      tc=-c;
              tz=-z;
              [tc tz
               a b]
           else
              [c z
                a b]
           end

    
            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:'])
	
		if minmax==1
		      c0(bas)*b
	        else
 	            -c0(bas)*b
 	        end
           

            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)])
		 if minmax==1
		    tc=-c;
                    tz=-z;
        	     [tc tz
        	      a b]
       		 else
       		     [c z
        	      a b]
        	 end

                           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;
       
            if pt > 20,
                 disp(['Hit return to continue'])
                 pause
            else
                 pause(pt)
            end
%
                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;
               
                end
%  Now for the objective row update
                ratio=c(t)/ap(i,t);
                c=c-ap(i,:)*ratio;
                z=z-bp(i)*ratio;
            
%  Now for the pivot row update
                a(i,:)=ap(i,:)/ap(i,t);
                b(i)=bp(i)/ap(i,t);
                  clc
                  home,disp([blanks(30)]),
               disp(['pivot= a(' int2str(i)+1 ',' int2str(t) ')' blanks(10)])
              disp(['after pivoting                 '])
          
		if minmax==1
		    tc=-c;
 	            tz=-z;
		    [tc tz
 	              a b]
  	      else
 	            [c z
   	             a b]
   	      end

            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
