-- TableauImpl.mesa
-- Last Edited by: Gnelson, November 21, 1983 4:59 pm

DIRECTORY Tableau, Real;

TableauImpl: PROGRAM IMPORTS Tableau EXPORTS Tableau = {
OPEN Tableau;

matrix: PUBLIC REF ARRAY[1 .. nrows] OF ARRAY[1 .. ncols] OF REAL ←
  NEW[ARRAY[1 .. nrows] OF ARRAY[1 .. ncols] OF REAL];

n, m: PUBLIC INT;

dcol: PUBLIC INT;

unit: PUBLIC Variable ← NewVariable[];

Init: PUBLIC PROC = {n ← 0; m ← 1; unit.i ← 0; unit.j ← 1; satisfiable ← TRUE};

satisfiable: PUBLIC BOOL;

NewVariable: PUBLIC PROC RETURNS [Variable] = 
  {v: Variable ← NEW[VariableRec];
   m ← m + 1;
   v.i ← 0;
   v.j ← m;
   x[m] ← v;
   FOR i: INT IN [1 .. n] DO matrix[i][m] ← 0 ENDLOOP;
   RETURN [v]};

Pivot: PUBLIC PROC[k, l: Variable] = 
  {k: INT ← rowvar.i;
    l: INT ← rowvar.j;
    -- We will pivot at the matrix element matrix[k, l].
    i, j: INT; -- These are temporaries in the loops below.
   {temp:Variable ← y[k]; y[k] ← x[l]; x[l] ← temp};  -- swap the row and column headers
   {pp: REAL ← 1.0 / matrix[k][l]; matrix[k][l] ← 1.0;
     FOR j IN [1..m] DO 
       matrix[k][j] ← matrix[k][j] * pp -- divide everything in pivot row by the pivot element:
     ENDLOOP;
     FOR i IN [1..n] DO 
       IF i # k THEN 
	    FOR j IN [1..m] DO 
	      IF j # l THEN 
	      -- for each matrix[i,j] outside the pivot row and column:
	      matrix[i][j] ← matrix[i][j] - matrix[i][l] * matrix[k][j]; 
	        -- note that matrix[k,j] was already * pp.
       ENDLOOP 
     ENDLOOP;
     -- Finally process pivot column:
     FOR i IN [1..n] 
       DO IF i # k THEN matrix[i][l] ← -matrix[i][l] * pp 
     ENDLOOP}};

AssertZ: PROC[l: LinearForm, PropEQ: PROC[u, v: Variable]] =
  {IF NOT satisfiable THEN ERROR;  
   {temp: Variable ← NewRowVar[];
    WHILE l # NIL DO Add[to: temp, what: l.first]; l ← l.rest ENDLOOP;
    -- now temp is a row variable that represents the linear form l
    Kill[temp, PropEQ]}};

NewRowVar: PROC[] RETURNS [v: Variable] =
  {v ← NEW[VariableRec];
   n ← n + 1;
   v.i ← n;
   v.j ← 0;
   y[n] ← v;
   FOR j: INT IN [1 .. m] DO matrix[n][j] ← 0 ENDLOOP;
   RETURN[v]};

Add: PROC[to: Variable, what: LinearMonomial] =
 -- to is assumed to be a row variable.
 {IF what.var.i = 0 
  THEN matrix[to.i][what.var.j] ← matrix[to.i][what.var.j] + what.coef
  ELSE FOR j: INT IN [1 .. m] 
           DO matrix[to.i][j] ← matrix[to.i][j] + what.coef * matrix[what.var.i][j]
           ENDLOOP};

Kill: PROC[v: Variable, PropEQ: PROC[u, v: Variable]] =
  {u: Variable;
   pec: PivotEffectCode;
   IF matrix[v.i][unit.j] < 0 
   THEN FOR j: INT IN [1 .. m] DO matrix[v.i][j] ← - matrix[v.i][j] ENDLOOP;
     -- Asserting -v = 0 is the same as asserting v = 0.
     -- Now the sample value of v is non-positive.  Thus maximizing v will
     -- move it towards zero.
   [u, pec] ← Maximize[v];
   SELECT TRUE FROM
      pec = ToPivotAtUWouldMakeVPositive OR pec = VIsManifestlyUnbounded =>  
         {Pivot[v.i, u.j]; KillColVar[v]; CheckRowEqualities[PropEQ]};
      pec = VIsMaximizedAtZero => {KillAllCols[v, PropEQ]; CheckRowEqualities[PropEQ]};
      pec = VIsMaximizedAndNegative => satisfiable ← FALSE;
   ENDCASE => ERROR};

KillColVar: PROC[v: Variable, PropEQ: PROC[u, v: Variable]] =
  {IF v.j # m 
   THEN {FOR i: INT IN [1 .. n] DO matrix[i][v.j] ← matrix[i][m] ENDLOOP;
            x[v.j] ← x[m]};
   m ← m - 1};

KillAllCols: PROC[v:Variable, PropEQ: PROC[u, v: Variable]] =
 -- kill all col vars that have a non-zero entry in v's row.
 {j: INT ← 1;
  WHILE j <= m DO
    IF j # unit.j THEN KillColVar[x[j], PropEQ];
    j ← j + 1
  ENDLOOP};

PivotEffectCode: TYPE = 
  {ToPivotAtUWouldMakeVPositive, VIsMaximizedAtZero, VIsMaximizedAndNegative,
   ToPivotAtUWouldIncreaseVButNotMakeItPositive, VIsManifestlyUnbounded};
   -- The last code is misleadingly named: it may "increase" v by zero.
   -- To guarantee termination, some cycle-avoding 
 
Maximize: PROC[v: Variable] RETURNS [u: Variable, pec: PivotEffectCode] =
  {[u, pec] ← FindPivot[v];
   WHILE pec = ToPivotAtUWouldIncreaseVButNotMakeItPositive
   DO Pivot[v.i, u.j]; [u, pec] ← FindPivot[v] ENDLOOP};

FindPivot: PROC[v: Variable] RETURNS [u: Variable, pec: PivotEffectCode] =
  {j: INT ← 0;
   sgn: REAL;
   FOR jj: INT IN [1 .. m] 
   DO IF (matrix[v.i][jj] # 0 AND NOT y[i].restricted) OR matrix[v.i][jj] > 0 
       THEN {j ← jj; sgn ← - matrix[v.i][jj]}
   ENDLOOP;   
   IF j = 0 THEN 
     SELECT TRUE FROM
       matrix[v.i][unit.j] < 0 => {pec ← VIsMaximizedAndNegative; RETURN};
       matrix[v.i][unit.j] = 0 => {pec ← VIsMaximizedAtZero; RETURN}
     ENDCASE => ERROR;
   u ← x[j]; -- we will pivot in column j
  {champ: INT;
   score: REAL;
   champ ← v.i;
   score ← Real.PlusInfinity;
   FOR ii: INT IN [1..n] 
   DO IF ii # v.i AND y[ii].restricted AND sgn * matrix[ii][j] < 0
       THEN IF ABS[matrix[ii][unit.j] / matrix[ii][j]] < score
               THEN {champ ← ii; score ← ABS[matrix[ii][unit.j] / matrix[ii][j]]}
   ENDLOOP;
   IF champ = i THEN {pec ← VIsManifestlyUnbounded; RETURN};
   {newSampleValueOfV: REAL = 
     matrix[v.i][unit.j] - matrix[v.i][j] * matrix[i][unit.j] / matrix[i][j];
    SELECT TRUE FROM
      newSampleValueOfV <= 0 => pec ← ToPivotAtUWouldIncreaseVButNotMakeItPositive;
      newSampleValueOfV > 0  => pec ← ToPivotAtUWouldMakeVPositive
    ENDCASE => ERROR}}};

Restrict: PROC[v:Variable, PropEQ: PROC[u, v: Variable]] =
  {IF v.restricted THEN RETURN;
   MakeVOwnARow[v];
   -- Now v owns a row, or else it owns an empty column.
   IF v.i = 0 THEN {v.restricted ← TRUE; RETURN};
   IF matrix[v.i][unit.j] < 0 THEN {v.restricted ← TRUE; RETURN};
   {u: Variable; pec: PivotEffectCode;
    [u, pec] ← Maximize[v];
    SELECT TRUE FROM
      pec = ToPivotAtUWouldMakeVPositive => {v.restricted ← TRUE; Pivot[v.i, u.j]};
      pec = VIsMaximizedAtZero => {KillAllCols[v, PropEQ]};
      pec = VIsMaximizedAndNegative => {satisfiable ← FALSE};
      pec = VIsManifestlyUnbounded => {v.restricted ← TRUE; Pivot[v.i, u.j]};
    ENDCASE => ERROR}};
     
 }.