[Cyan]<Tioga>TableTool>TableauImpl.Mesa!1

-- TableauImpl.mesa

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

-- Last Edited by: Beach, January 12, 1984 1:35 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, champion: 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, champion, pec] ← MaximalPivot[v];

SELECT TRUE FROM

pec = ToPivotAtUWouldMakeVPositive OR

pec = VIsManifestlyUnbounded =>

{Pivot[champion.i, u.j]; KillColVar[champion]; 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

MaximalPivot: PROC[v: Variable] RETURNS [u: Variable, champion: Variable ← NIL, pec: PivotEffectCode] =

{[u, champion, pec] ← FindPivot[v];

WHILE pec = ToPivotAtUWouldIncreaseVButNotMakeItPositive

DO Pivot[champion.i, u.j]; [u, champion, pec] ← FindPivot[v] ENDLOOP};

FindPivot: PROC[v: Variable] RETURNS [u: Variable, champion: Variable ← NIL, 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 ← - 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;

champion ← y[champ];

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, champion, pec] ← MaximalPivot[v];

SELECT TRUE FROM

pec = ToPivotAtUWouldMakeVPositive => {v.restricted ← TRUE; Pivot[champion.i, u.j]};

pec = VIsMaximizedAtZero => {KillAllCols[champion, PropEQ]};

pec = VIsMaximizedAndNegative => {satisfiable ← FALSE};

pec = VIsManifestlyUnbounded => {v.restricted ← TRUE; Pivot[v.i, u.j]};

ENDCASE => ERROR}};

Maximize: PROC[v:Variable] RETURNS [unbounded: BOOLEAN ← FALSE] = {

u, champion: Variable;

pec: PivotEffectCode;

IF v.i # 0 THEN MakeVOwnARow[v]; -- try to make v a row

IF v.i = 0 THEN RETURN[NOT v.restricted];

[u, champion, pec] ← FindPivot[v];

DO SELECT TRUE FROM

pec = ToPivotAtUWouldMakeVPositive OR

pec = ToPivotAtUWouldIncreaseVButNotMakeItPositive => {

Pivot[champion.i, u.j]; [u, champion, pec] ← FindPivot[v]};

ENDCASE => EXIT; ENDLOOP;

unbounded ← pec = VIsManifestlyUnbounded;

};

MakeVOwnARow: PROC [v:Variable] = {

champ: INT;

score: REAL;

sgn: REAL;

i: INT ← 1;

WHILE i < n AND matrix[i][v.j] = 0 DO i ← i + 1 ENDLOOP;

SELECT TRUE FROM

i = n => RETURN;

matrix[i][v.j] > 0 => sgn ← +1;

matrix[i][v.j] < 0 => sgn ← -1

ENDCASE => ERROR;

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;

champion ← y[champ];

Pivot[champ.i, v.j];

};