
(*   INPUT	: The associated datafile for dual simplex method is 
		  called "DualplexDatafile".

		  FIRST NUMBER in "DualplexDatafile" represents # of
		  VARIABLES in Linear  Program problem(LP).

		  SECOND NUMBER represents # of CONSTRINTS in LP.
		
	          First number in each of rest of rows in
		  "DualplexDatafile" represents the Right Hand Side(RHS)
		  vector b.  Initially objective value is 0.
  
		  THIRD set of numbers represents the cost matrix C of
		  then object function. Cost matrix strats on second #.

		  FOURTH set of numbers represents the constraint MATRIX A
		  Matrix A also starts on the second number.

    Algorithm	: The dual simplex method solves the LP problem in the
		  following form:

                     minimize(or maximize)   z=cTx+c^Tx^
		     s.t.
					     Ax + Ix^ = b
					     x, x^ >= 0

		  Maximum # of variables, and maximum # of constraints are
		  set to maxvar =100 and maxconstraint=50 respectively.
		  One can modify them according to their needs. But
		  maxconstraint is <= maxvar.

		  EPSis small ral number such that if for ant real
		  number a, |a|< EPS, then a =0.0
		  Initially, EPS = 0.0001

		  INF is maximal real number available in the floating
		  point number system that is used.
		  Initially, INF = 999.00

    OUTPUT	: Outputs are
		  1.Check if a LP is feasible and optimal solution exists.
		  2.Determine the optimal basic variables and their
		    corresponding values.
		  3.Determine the optimal value of the objective function.
								     *)



	        
program Dualplex (input, output, DualplexDatafile,DualplexOutfile);

const   INF 	= 999.00;
	maxvar	= 100;
	maxconstraint	= 50;
	EPS	= 0.0001;

type	CHARFILE	= file of char;
	ARRMN		= array[0..maxconstraint,0..maxvar] of real;
	ARRN		= array[0..maxvar] of integer;
	ARRMIN		= array[0..maxvar-maxconstraint] of integer;

var	DualplexDatafile : CHARFILE;
	DualplexOutfile  : CHARFILE;
	N,M, NMINM	 : integer;
	FOPT		 : integer;
        A		 : ARRMN;
	U		 : ARRN;
	NOFEAS, NOSOL	 : boolean;
	Nextint		 : integer;
	Nextreal	 : real;

procedure Infile(var	N,M,NMINM : integer;
		 var	Nextint	  : integer;
		 var	Nextreal  : real;
		 var	A	  : ARRMN);

var	row,column : integer;

begin
  reset(DualplexDatafile);
  readln(DualplexDatafile,Nextint);
  N := Nextint;
  readln(DualplexDatafile,Nextint);
  M := Nextint;
  NMINM := N-M;
  for row := 0 to M do
  begin
    if row = 0 then
    begin
      for column := 0 to N do
      begin
        read(DualplexDatafile,Nextreal);
	A[row,column] := Nextreal;
      end;
      readln(DualplexDatafile);
    end
    else
    begin
      for column := 0 to NMINM do
      begin
        read(DualplexDatafile, Nextreal);
	A[row,column] := Nextreal;
      end;
      readln(DualplexDatafile);
    end;
  end;
end;

procedure DSIMPLEX(
       M,N,FOPT    :integer;
   var A           :ARRMN;
   var U           :ARRN;
       EPS,INF     :real;
   var NOFEAS,NOSOL:boolean);

   var I,J,K,K1,K2,K3,K4,L,W:integer;
       MIN,XM,XS            :real;
       B,StoP               :boolean;
       Z,Z1                 :ARRMIN;
begin
   NOFEAS:=false;  NOSOL:=false;
   K4:=N-M;
   A[0,0]:=0.0;
   for I:=0 to K4 do begin
      XS:=0.0;
      for J:=1 to M do XS:=XS+A[J,I]*A[0,K4+J];
      A[0,I]:=XS-A[0,I]
   end;
   for I:=K4+1 to N do
      for J:=1 to M do
         if I = K4+J then A[J,I]:=1.0
         else A[J,I]:=0.0;
   I:=0;
   while (not NOFEAS) and (I < K4) do begin
      I:=I+1;  XS:=A[0,I];
      NOFEAS:=(abs(XS) > EPS) and (XS*FOPT < 0);
      if not NOFEAS then U[M+I]:=I
   end;
   if not NOFEAS then begin
      for I:=1 to M do U[I]:=K4+I;
      StoP:=false;
      repeat  (* until StoP *)
         MIN:=0.0;  B:=true;  I:=0;
         repeat  (* until StoP or (I >= M) *)
            I:=I+1;  J:=M;  XS:=A[I,0];
            if XS < -EPS then begin
               StoP:=true;
               while (J < N) and StoP do begin
                  J:=J+1;  W:=U[J];
                  StoP:=A[I,W] >= -EPS
               end;
               if StoP then NOSOL:=true
               else begin
                  B:=false;
                  if XS-MIN < -EPS then begin
                     MIN:=XS;  L:=I
                  end
               end  (* else: not StoP *)
            end  (* if XS < -EPS *)
         until StoP or (I >= M);
         if not StoP then begin
            if B then begin NOSOL:=false;  StoP:=true end
            else begin
               MIN:=INF;
               for J:=1 to K4 do Z1[J]:=M+J;
               for I:=0 to M do
                  if (I <> 1) and (not B) then begin
                     K:=0;
                     for J:=1 to K4 do Z[J]:=Z1[J];
                     K3:=1;
                     for J:=M+1 to N do
                        if J = Z[K3] then begin
                           K3:=K3+1;  W:=U[J];  XS:=A[L,W];
                           if XS < -EPS then begin
                              XS:=abs(A[I,W]/XS);  XM:=XS-MIN;
                              if abs(XM) < EPS then begin
                                 K:=K+1;  Z1[K]:=J;
                                 B:=false
                              end
                              else
                                 if XM < 0.0 then begin
                                    MIN:=XS;  K1:=J;  K2:=W;
                                    Z1[1]:=1;  K:=1;
                                    for W:=2 to K4 do Z1[W]:=0;
                                    B:=true
                                 end
                           end  (* if XS < -EPS *)
                        end  (* if J = Z[K3], for J *)
                  end;  (* if I <> 1 and (not B), for I *)
               MIN:=1.0/A[L,K2];
               U[K1]:=U[L];
               if L = 0 then I:=1 else I:=0;
               repeat
                  XS:=A[I,K2]*MIN;
                  A[I,0]:=A[I,0]-A[L,0]*XS;
                  for J:=M+1 to N do begin
                     W:=U[J];
                     A[I,W]:=A[I,W]-A[L,W]*XS
                  end;
                  if I = L-1 then I:=I+2 else I:=I+1
               until I > M;
               for J:=M+1 to N do begin
                  W:=U[J];
                  A[L,W]:=A[L,W]*MIN
               end;
               A[L,0]:=A[L,0]*MIN;
               for I:=0 to M do
                  if I = 1 then A[I,K2]:=1.0
                  else A[I,K2]:=0.0;
               U[L]:=K2
            end  (* else: not B *)
         end  (* if not StoP *)
      until StoP
   end  (* if not NOFEAS *)
end;  (* DSIMPLEX *)


procedure Outfile(NOFEAS,NOSOL : boolean;
		  M	       : integer;
		  A	       : ARRMN;
		  U	       : ARRN);

var counter : integer;

begin
  rewrite(DualplexOutfile);
  writeln (DualplexOutfile,'   NOFEAS = ',NOFEAS,'     NOSOL  =  ',NOSOL);
  if (NOFEAS = false) and (NOSOL = false) then
  begin
    writeln (DualplexOutfile,'  Optimal Value is  ',A[0,0]:9:2);
    writeln (DualplexOutfile,'       Basic-Var   Value  ');
    for counter := 1 to M do
    begin
      write(DualplexOutfile,U[counter]);
      writeln(DualplexOutfile,'     ',A[counter,0]:9:2);
    end;
  end;
end;
   
begin (* main *)
  Infile(N,M,NMINM,Nextint,Nextreal,A);
  DSIMPLEX(M,N,FOPT,A,U,EPS,INF,NOFEAS,NOSOL);
  Outfile(NOFEAS,NOSOL,M,A,U);
end.
