----------------------------------------------------------------------
     -- Title       :  LineAK
     -- Version     :  1.0
     -- Author      :  David Wallace Croft, CompuServe [76600,102]
     -- Compiler    :  Ada
     -- Unit Type   :  package body
     -- Copyright   :  1994 David Wallace Croft.  All rights reserved.
     -- Description :  Linear Programming & Game Theory
     ----------------------------------------------------------------------

     with ConsAK ; use ConsAK ;
     with FloaAK ; use FloaAK ;
     with InteAK ; use InteAK ;
     with Text_IO; use Text_IO;

     package body LineAK is
     ----------------------------------------------------------------------
     ----------------------------------------------------------------------

     procedure Demo is
     ----------------------------------------------------------------------
     begin
       Simplex_Revised_Demo;
       Pause;
     end Demo;

     procedure Simplex_Revised (
       Degenerate :    out boolean;
       X          :    out Vector_Type;
       Y          :    out Vector_Type;
       Z          :    out float;
       A          : in     Matrix_Type;
       B          : in     Vector_Type;
       C          : in     Vector_Type;
       Show       : in     boolean := false;
       Near_0     : in     float   := 1.0e-10 ) is
     ----------------------------------------------------------------------
     -- Assumes A non-degenerate.
     --
     -- Franklin, Joel.  "Methods of Mathematical Economics:  Linear and
     --   Non-linear Programming, Fixed-Point Theorem".  Springer-Verlag:
     --   NY, 1980, pp93-.
     ----------------------------------------------------------------------
       M         : constant positive := A'last ( 1 );
       N         : constant positive := A'last ( 2 );
       I         :          Matrix_Type ( 1..M, 1..M );
       A1        :          Matrix_Type ( 1..M, 1..( N + M ) )
	 := ( others => ( others => 0.0 ) );
       C1        :          Vector_Type ( 1..( N + M ) );
       Tableau   :          Matrix_Type ( 0..M, 0..M )
	 := ( others => ( others => 0.0 ) );
       U         :          Matrix_Type ( 1..M, 1..M );
       Basis     :          Vector_Integer_Type ( 1..M );
       Old_Basis :          Vector_Integer_Type ( 1..M );
       X_Temp    :          Vector_Type ( X'range );
       Degenerate_Temp : boolean;
     begin
       Degenerate := false;
       if Show then Put_Line ( "Phase 1" ); end if;
     -- A1 := [ A, I ]
       A1 := Concatenate ( A, Identity ( I ) );
       C1 := Concatenate ( Fill ( C ), Fill ( B, 1.0 ) );
     -- T_i0 = B_i
       for I in 1..M loop
	 Tableau ( I, 0 ) := B ( I );
       end loop;
     -- U_ij = Delta_ij, initially the identity matrix.
       U := Identity ( U );
       for I in 1..M loop
	 for J in 1..M loop
	   Tableau ( I, J ) := U ( I, J );
	 end loop;
       end loop;
     -- Z0 = Sum ( B_I )
       for I in B'range loop
	 Tableau ( 0, 0 ) := Tableau ( 0, 0 ) + B ( I );
       end loop;
     -- Y_i = 1
       for I in 1..M loop
	 Tableau ( 0, I ) := 1.0;
       end loop;
       for index in Basis'range loop
	 Basis ( index ) := N + index;
       end loop;
       loop
	 if Show then
	   Put_Line ( "Extended Tableau" );
	   Put_Line ( Tableau );
	   Put_Line ( "Basis" );
	   Put_Line ( Basis );
	   Pause;
	 end if;
	 Old_Basis := Basis;
	 Simplex_Revised_2 (
	   Degenerate_Temp, Tableau, Basis, A1, C1, Show, Near_0 );
	 if Degenerate_Temp then
	   Degenerate := true;
	   return;
	 end if;
	 exit when Basis = Old_Basis;
       end loop;
     if Show then Put_Line ( "Phase 2" ); end if;
     -- recalculate Z(0)
       Tableau ( 0, 0 ) := 0.0;
       for index in Basis'range loop
	 begin
	   Tableau ( 0, 0 ) := Tableau ( 0, 0 )
	     + Tableau ( 0, index ) * C ( Basis ( index ) );
	 exception
	   when constraint_error =>
	     if Show then
	       Put_Line ( "Degenerate:  B is a linear combination of fewer"
		 & " than m columns of A" );
	       Put_Line ( "where m is the number of rows in A." );
	       Pause;
	     end if;
	     Degenerate := true;
	     return;
	 end;
       end loop;
     -- recalculate Y vector
       for J in 1..M loop
	 Tableau ( 0, J ) := 0.0;
	 for I in Basis'range loop
	   Tableau ( 0, J ) := Tableau ( 0, J )
	     + Tableau ( I, J ) * C ( Basis ( I ) );
	   if abs ( Tableau ( 0, J ) ) < Near_0 then
	     Tableau ( 0, J ) := 0.0;
	   end if;
	 end loop;
       end loop;
       loop
	 Old_Basis := Basis;
	 if Show then
	   Put_Line ( "Extended Tableau" );
	   Put_Line ( Tableau );
	   Put_Line ( "Basis" );
	   Put_Line ( Basis );
	   Pause;
	 end if;
	 Simplex_Revised_2 (
	   Degenerate_Temp, Tableau, Basis, A, C, Show, Near_0 );
	 if Degenerate_Temp then
	   Degenerate := true;
	   return;
	 end if;
	 exit when Basis = Old_Basis;
       end loop;
       X_Temp := ( others => 0.0 );
       for index in Basis'range loop
	 begin
	   X_Temp ( Basis ( index ) ) := Tableau ( index, 0 );
	 exception
	   when constraint_error =>
	     if Show then
	       Put_Line ( "Degenerate:  no feasible solution exists." );
	       Pause;
	     end if;
	     Degenerate := true;
	     return;
	 end;
       end loop;
       X := X_Temp;
       for index in 1..M loop
	 Y ( index ) := Tableau ( 0, index );
       end loop;
       Z := Dot ( C, X_Temp );
     end Simplex_Revised;

     procedure Simplex_Revised_2 (
       Degenerate :    out boolean;
       Tableau    : in out Matrix_Type;
       Basis      : in out Vector_Integer_Type;
       A          : in     Matrix_Type;
       C          : in     Vector_Type;
       Show       : in     boolean := false;
       Near_0     : in     float   := 1.0e-10 ) is
     ----------------------------------------------------------------------
     -- p94
     --       0    1     2 ...   M
     --    -------------------------
     --  0 |  Z0 |  Y1  Y2 ...  YM |   Extended Tableau Matrix format
     --    -------------------------
     --  1 | T10 | U11 U12 ... U1M |
     --    |     |                 |
     --  2 | T20 | U21 U22 ... U2M |
     --    |     |                 |
     --  : |  :  |  :   :  ...  :  |
     --    |     |                 |
     --  M | TM0 | UM1 UM2 ... UMM |
     --    -------------------------
     ----------------------------------------------------------------------
       Case_2a   :          boolean;
       M         : constant positive                := Basis'last;
       R         :          natural                 := 0;
       S         :          natural                 := 0;
       T_S       :          Vector_Type ( 1..M );
       Theta     :          Vector_Type ( 0..M );
       Z         :          Vector_Type ( A'range ( 2 ) );
       Non_Basic :          boolean;
       Old_Row_R :          Vector_Type ( Tableau'range ( 2 ) );
     begin
       Degenerate := false;
       for J in A'range ( 2 ) loop
     -- Find a non-basic index j.
	 Non_Basic := true;
	 for index in Basis'range loop
	   Non_Basic := Non_Basic and ( J /= Basis ( index ) );
	   exit when not Non_Basic;
	 end loop;
	 if Non_Basic then
     -- Check if the cost can be reduced.
     -- Z ( J ) := Y_t * A_J;
	   Z ( J ) := 0.0;
	   for I in 1..M loop
	     Z ( J ) := Z ( J ) + Tableau ( 0, I ) * A ( I, J );
	   end loop;
	   if Z ( J ) > C ( J ) then
     -- If so, move column J into the basis.
	     S := J;
	     exit;
	   end if;
	 end if;
       end loop;
     -- Otherwise, X is already optimal.
       if S = 0 then
	 return;
       end if;
       if Show then
	 Put_Line ( "Pivot column =" & positive'image ( S ) );
       end if;
     -- Compute the pivot column of the simplex tableau.
       Case_2a := true;
       for I in 1..M loop
	 T_S ( I ) := 0.0;
	 for K in 1..M loop
	   T_S ( I ) := T_S ( I ) + Tableau ( I, K ) * A ( K, S );
	 end loop;
     -- Check to see if costs can be driven to minus infinity.
	 Case_2a := Case_2a and ( T_S ( I ) < 0.0 );
       end loop;
       if Case_2a then
	 if Show then
	   Put_Line ( "Degenerate:  cost can go to negative infinity." );
	   Pause;
	 end if;
	 Degenerate := true;
	 return;
       end if;
     -- Choose the pivot row R.
       R := 0;
       for I in T_S'range loop
	 if T_S ( I ) > 0.0 then
	   if R = 0 then
	     R := I;
	   else
	     if Tableau ( I, 0 ) / T_S ( I )
	       < Tableau ( R, 0 ) / T_S ( R ) then
		 R := I;
	     end if;
	   end if;
	 end if;
       end loop;

     -- Replace the Basis with S.
       if Show then
	 Put_Line ( "Pivot row =" & positive'image ( R ) );
       end if;
       Basis ( R ) := S;
     -- Setup the theta vector
       Theta ( 0 ) := ( Z ( S ) - C ( S ) ) / T_S ( R );
       for I in 1..Tableau'last loop
	 Theta ( I ) := T_S ( I ) / T_S ( R );
       end loop;
     -- Compute new rows
       for Col in Tableau'range ( 2 ) loop
	 Old_Row_R ( Col ) := Tableau ( R, Col );
       end loop;
       for Row in Tableau'range ( 1 ) loop
	 for Col in Tableau'range ( 2 ) loop
	   if Row = R then
	     Tableau ( R, Col ) := Old_Row_R ( Col ) / T_S ( R );
	   else
	     Tableau ( Row, Col ) := Tableau ( Row, Col )
	       - Theta ( Row ) * Old_Row_R ( Col );
	   end if;
	   if abs ( Tableau ( Row, Col ) ) < Near_0 then
	     Tableau ( Row, Col ) := 0.0;
	   end if;
	 end loop;
       end loop;
     end Simplex_Revised_2;

     procedure Simplex_Revised_Demo is
     ----------------------------------------------------------------------
       B_Last, C_Last : positive;
       Degenerate : boolean;
     begin
       Put_Line ( "Simplex_Revised_Demo (C) 1994 David Wallace Croft" );
       New_Line;
       Put_Line (
	 "Given matrix A, requirement vector B, and cost vector C," );
       Put_Line ( "find the solution vector X and shadow vector Y" );
       Put_Line (
	 "for the primal A*X >= B, X >= 0, Transpose (C) * X = minimum" );
       Put_Line ( "and the dual Transpose (Y) * A <= Transpose (C)," );
       Put_Line ( "Transpose (Y) * B = maximum; also show the" );
       Put_Line ( "cost Z = Transpose (C) * X = Transpose (Y) * B." );
       Put_Line ( "This assumes the problem is non-degenerate." );
       New_Line;
       B_Last := Ask_Pos ( "How many rows are in your vector B? ", 2, 1 );
       C_Last := Ask_Pos ( "How many rows are in your vector C? ", 3, 1 );
       declare
	 A : Matrix_Type ( 1..B_Last, 1..C_Last );
	 B : Vector_Type ( 1..B_Last );
	 C : Vector_Type ( 1..C_Last );
	 X : Vector_Type ( 1..C_Last );
	 Y : Vector_Type ( 1..B_Last );
	 Z : float;
       begin
	 Ask ( A, "Please enter the matrix A." );
	 Ask ( B, "Please enter the requirement vector B." );
	 Ask ( C, "Please enter the cost vector C." );
	 Simplex_Revised ( Degenerate, X, Y, Z, A, B, C, true );
	 if not Degenerate then
	   Put_Line ( "Your solution vector X is" );
	   Put_Line ( X, 2, 4 );
	   Put_Line ( "Your shadow vector Y is" );
	   Put_Line ( Y, 2, 4 );
	   Put_Line ( "Your cost Z is" );
	   Put_Line ( Z, 2, 4 );
	 else
	   Put_Line (
	     "The problem is degenerate; no solution is given using the" );
	   Put_Line ( "Revised Simplex Algorithm." );
	 end if;
       end;
     end Simplex_Revised_Demo;

     ----------------------------------------------------------------------
     ----------------------------------------------------------------------
     end LineAK;