=
BEGIN
OPEN LinearSolver;
RowMax: INT = 100;
Row: TYPE = [0 .. RowMax);
Col: TYPE = [0 .. ColMax);
Tableau: TYPE = REF TableauRec;
TableauRec: TYPE = RECORD [
n: [0 .. RowMax],
m: [0 .. ColMax],
maxN: [0 .. RowMax],
maxM: [0 .. ColMax],
a: REF ARRAY Row OF ARRAY Col OF REAL,
a[i][j] is only relevant if i is in [0, n) and j is in [0, m). Thus we declare a matrix of n rows by m columns, which can grow as large as RowMax by ColMax.
x: REF ARRAY Col OF Unknown,
y: REF ARRAY Row OF Unknown,
unity: Unknown,
satisfiable: BOOL,
nameCount: NAT
];
doesntOwnARow: INT ~ -1;
x, y, a, n, and m are together called the *Tableau*. They represent the n constraints
(and i: i in [0, n): y(i) = (sum j: j in [0, m): a(i, j) * x(j))).
Every active unknown appears once in either the x or y array, and does not appear in both.
The unknowns in x are called column variables, the unknowns in y are called row variables.
The module maintains the invariant that the variable unity is one of the column variables.
That is,
(E j: j in [0, m): x(j) = unity).
The module also maintains the inverse of the arrays x and y, using the space in the records
that represent the unknowns. The i field is the inverse of the y array, and the j field is
the inverse of the x array. The unused inverse field is set to -1. That is:
(A i: i in [0, n): x(i).i = i and x(i).j = -1)
(A j: j in [0, m): x(j).j = j and x(j).i = -1)
Notice that one solution to the constraints is
unity = 1,
x(j) = 0 for j in [0, n) and x(j) # unity
y(i) = a[i, unity.j] for i in [0, m)
This is the solution returned by Solution.
signalUnsatisfiable: BOOL ← TRUE;
Unsatisfiable: ERROR = CODE;
C:
PUBLIC
PROCEDURE [x: LinearMonomial]
RETURNS [r:
REF] = {
r ← NEW[LinearMonomial ← x];
};
AssertZero:
PUBLIC
PROCEDURE [t: Tableau, c: Constraint, propEQ: PropEQProc ←
NIL] = {
Asserts c = 0, in two steps. Step 1: create a new unknown v and assert v = c.
This is done by creating a new row owned by v that represents c.
Step 2: Assert v = 0. This is done by pivoting v into a column and deleting the column
IF NOT t.satisfiable THEN ERROR;
{v: Unknown ← NewRowUnknown[t];
IF debug THEN Debug[t, IO.PutFR["AssertZero c = %g", IO.refAny[c]]];
WHILE c #
NIL
DO
g: REF LinearMonomial ← NARROW[c.first];
AddToRow[t, t.n-1, g.coefficient, g.unknown];
c ← c.rest
ENDLOOP;
KillRowVar[t, v, propEQ] }
};
AssertZero2:
PUBLIC
PROCEDURE [t: Tableau, c1:
REAL, v1: Unknown, c2:
REAL, v2: Unknown, propEQ: PropEQProc ←
NIL] = {
v: Unknown ← NewRowUnknown[t];
IF debug THEN Debug[t, IO.PutFR["AssertZero2 c = %g", IO.refAny[ LIST[C[[c1, v1]], C[[c2, v2]]]]]];
AddToRow[t, t.n-1, c1, v1];
AddToRow[t, t.n-1, c2, v2];
KillRowVar[t, v, propEQ]
};
AssertZero3:
PUBLIC
PROCEDURE [t: Tableau, c1:
REAL, v1: Unknown, c2:
REAL, v2: Unknown, c3:
REAL, v3: Unknown, propEQ: PropEQProc ←
NIL] = {
v: Unknown ← NewRowUnknown[t];
IF debug THEN Debug[t, IO.PutFR["AssertZero3 c = %g", IO.refAny[ LIST[C[[c1, v1]], C[[c2, v2]], C[[c3, v3]]]]]];
AddToRow[t, t.n-1, c1, v1];
AddToRow[t, t.n-1, c2, v2];
AddToRow[t, t.n-1, c3, v3];
KillRowVar[t, v, propEQ]
};
AssertGEZero:
PUBLIC
PROCEDURE [t: Tableau, c: Constraint, propEQ: PropEQProc ←
NIL] = {
Asserts c >= 0, in two steps. Step 1: create a new unknown v, a new slack variable s (s >= 0), and assert v = c + s.
This is done by creating a new row owned by v that represents c.
Step 2: Assert v = 0. This is done by pivoting v into a column and deleting the column
IF NOT t.satisfiable THEN ERROR;
{v: Unknown ← NewRowUnknown[t];
slack: Unknown ← NewUnknown[t, c];
IF debug THEN Debug[t, IO.PutFR["AssertGEZero c = %g", IO.refAny[c]]];
Restrict[t, slack];
WHILE c #
NIL
DO
g: REF LinearMonomial ← NARROW[c.first];
AddToRow[t, t.n-1, g.coefficient, g.unknown];
c ← c.rest
ENDLOOP;
AddToRow[t, t.n-1, -1.0, slack];
KillRowVar[t, v, propEQ] }
};
AssertGEZero2:
PUBLIC
PROCEDURE [t: Tableau, c1:
REAL, v1: Unknown, c2:
REAL, v2: Unknown, propEQ: PropEQProc ←
NIL] = {
v: Unknown ← NewRowUnknown[t];
slack: Unknown ← NewUnknown[t, InventName[t, "slack"]];
IF debug THEN Debug[t, IO.PutFR["AssertGEZero2 c = %g", IO.refAny[ LIST[C[[c1, v1]], C[[c2, v2]]]]]];
Restrict[t, slack];
AddToRow[t, t.n-1, c1, v1];
AddToRow[t, t.n-1, c2, v2];
AddToRow[t, t.n-1, -1.0, slack];
KillRowVar[t, v, propEQ] ;
};
AssertGEZero3:
PUBLIC
PROCEDURE [t: Tableau, c1:
REAL, v1: Unknown, c2:
REAL, v2: Unknown, c3:
REAL, v3: Unknown, propEQ: PropEQProc ←
NIL] = {
v: Unknown ← NewRowUnknown[t];
slack: Unknown ← NewUnknown[t, InventName[t, "slack"]];
IF debug THEN Debug[t, IO.PutFR["AssertGEZero3 c = %g", IO.refAny[ LIST[C[[c1, v1]], C[[c2, v2]], C[[c3, v3]]]]]];
Restrict[t, slack];
AddToRow[t, t.n-1, c1, v1];
AddToRow[t, t.n-1, c2, v2];
AddToRow[t, t.n-1, c3, v3];
AddToRow[t, t.n-1, -1.0, slack];
KillRowVar[t, v, propEQ]
};
NewUnknown:
PUBLIC
PROCEDURE [t: Tableau, name:
REF]
RETURNS [Unknown] = {
v: Unknown ← NEW[UnknownRec];
IF debug THEN Debug[t, IO.PutFR["NewUnknown %g", IO.refAny[name]]];
v.name ← name;
v.i ← doesntOwnARow;
v.j ← t.m;
FOR i: INT IN [0 .. t.n) DO t.a[i][t.m] ← 0 ENDLOOP;
t.x[t.m] ← v;
t.m ← t.m + 1;
t.maxM ← MAX[t.m, t.maxM];
RETURN [v]
};
NewUnknown is easily seen to maintain the module invariant: since the new column is zeroed out, none of the linear combinations represented by the rows are changed.
NewRowUnknown:
PROCEDURE [t: Tableau]
RETURNS [v: Unknown] = {
v ← NEW[UnknownRec];
v.name ← InventName[t, "rVar"];
v.i ← t.n;
v.j ← doesntOwnACol;
FOR j: INT IN [0 .. t.m) DO t.a[t.n][j] ← 0 ENDLOOP;
t.y[t.n] ← v;
t.n ← t.n + 1;
t.maxN ← MAX[t.n, t.maxN];
RETURN[v];
};
AddToRow:
PROCEDURE [t: Tableau, n: Row, c:
REAL, v: Unknown] = {
adds c*v to the linear combination represented by the nth row,
either by incrementing one element of the row (if v owns a column)
or by doing a row operation (if v owns a row)
SELECT
TRUE
FROM
v.j # doesntOwnACol => t.a[n][v.j] ← t.a[n][v.j] + c;
v.i # doesntOwnARow =>
FOR j: INT IN [0 .. t.m) DO t.a[n][j] ← t.a[n][j] + t.a[v.i][j] * c ENDLOOP;
ENDCASE => ERROR -- Client called AssertZero[f] where f contains an inactive variable
};
Satisfiable: PUBLIC PROCEDURE [t: Tableau] RETURNS [BOOL] = {RETURN [t.satisfiable]};
Solution:
PUBLIC
PROCEDURE [t: Tableau, v: Unknown]
RETURNS [
REAL] =
BEGIN
SELECT TRUE FROM
v = t.unity => RETURN[1];
v.j # doesntOwnACol AND v # t.unity => RETURN[0];
v.i # doesntOwnARow => RETURN[t.a[v.i][t.unity.j]]
ENDCASE => ERROR
END;
KillRowVar:
PROCEDURE [t: Tableau, v: Unknown, propEQ: PropEQProc ←
NIL] =
{u, champion: Unknown;
pec: PivotEffectCode;
IF debug THEN Debug[t, IO.PutFR["Kill row variable %g", IO.refAny[v.name]]];
IF
ABS[t.a[v.i][t.unity.j]] < 0.0001
THEN {
We may have a redundant constraint, silly boy!
allZeros: BOOLEAN ← TRUE;
FOR j:
INT
IN [0 .. t.m)
DO
allZeros ← allZeros AND ABS[t.a[v.i][j]] < 0.0001;
ENDLOOP;
IF allZeros
THEN {
Back up tableau one row
t.n ← t.n-1; v← NIL;
RETURN;
};
};
IF t.a[v.i][t.unity.j] > 0
THEN
FOR j: INT IN [0 .. t.m) DO t.a[v.i][j] ← - t.a[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[t, v];
SELECT
TRUE
FROM
pec = ToPivotAtUWouldMakeVPositive OR
pec = VIsManifestlyUnbounded =>
{Pivot[t, champion.i, u.j]; KillColVar[t, champion, propEQ]; CheckRowEqualities[t, propEQ]};
pec = VIsMaximizedAtZero OR
pec = VIsManifestlyMaximized =>
{KillAllCols[t, v, propEQ]; CheckRowEqualities[t, propEQ]};
pec = VIsMaximizedAndNegative =>
IF signalUnsatisfiable THEN ERROR Unsatisfiable ELSE t.satisfiable ← FALSE;
ENDCASE => ERROR;
};
CheckRowEqualities:
PROCEDURE [t: Tableau, propEQ: PropEQProc] = {
};
KillColVar:
PROCEDURE [t: Tableau, v: Unknown, propEQ: PropEQProc ←
NIL] = {
IF debug THEN Debug[t, IO.PutFR["Kill column %g", IO.refAny[v.name]]];
t.m ← t.m - 1;
IF v.j # t.m
THEN {FOR i: INT IN [0 .. t.n) DO t.a[i][v.j] ← t.a[i][t.m] ENDLOOP; t.x[v.j] ← t.x[t.m]; t.x[t.m].j ← v.j}};
KillAllCols:
PROCEDURE [t: Tableau, v:Unknown, propEQ: PropEQProc ←
NIL] = {
kill all col vars that have a non-zero entry in v's row.
IF debug THEN Debug[t, IO.PutFR["Kill all columns for %g", IO.refAny[v.name]]];
FOR j:
INT
IN [0..t.m)
DO
IF j # t.unity.j AND t.a[v.i][j] # 0 THEN KillColVar[t, t.x[j], propEQ] ENDLOOP;
};
PivotEffectCode:
TYPE = {
ToPivotAtUWouldMakeVPositive,
VIsMaximizedAtZero,
VIsMaximizedAndNegative,
ToPivotAtUWouldIncreaseVButNotMakeItPositive,
VIsManifestlyUnbounded,
VIsManifestlyMaximized };
The last code is misleadingly named: it may "increase" v by zero.
To guarantee termination, some cycle-avoding
MaximalPivot:
PROCEDURE [t: Tableau, v: Unknown]
RETURNS [u, champion: Unknown ←
NIL, pec: PivotEffectCode] = {
[u, champion, pec] ← FindPivot[t, v];
WHILE pec = ToPivotAtUWouldIncreaseVButNotMakeItPositive
DO
Pivot[t, champion.i, u.j];
[u, champion, pec] ← FindPivot[t, v];
ENDLOOP;
};
FindPivot:
PROCEDURE [t: Tableau, v: Unknown]
RETURNS [u, champion: Unknown ←
NIL, pec: PivotEffectCode] = {
SGN:
PROCEDURE [x:
REAL]
RETURNS [
REAL] =
INLINE {
RETURN [IF x < 0 THEN -1.0 ELSE +1.0]
};
j: INT ← doesntOwnACol;
sgn: REAL;
FOR jj:
INT
IN [0 .. t.m)
DO
IF jj # t.unity.j
THEN
IF (t.a[v.i][jj] # 0
AND
NOT t.x[jj].restricted)
OR t.a[v.i][jj] > 0
THEN {j ← jj; sgn ← SGN[t.a[v.i][jj]]}
ENDLOOP;
IF j = doesntOwnACol
THEN
SELECT
TRUE
FROM
t.a[v.i][t.unity.j] < 0 => {pec ← VIsMaximizedAndNegative; RETURN};
t.a[v.i][t.unity.j] = 0 => {pec ← VIsMaximizedAtZero; RETURN};
t.a[v.i][t.unity.j] > 0 => {pec ← VIsManifestlyMaximized; RETURN};
ENDCASE => ERROR;
u ← t.x[j]; -- we will pivot in column j
{champ: INT;
score: REAL;
champ ← v.i;
score ← Real.LargestNumber;
FOR ii:
INT
IN [0..t.n)
DO
IF ii # v.i
AND t.y[ii].restricted
AND sgn * t.a[ii][j] < 0
THEN
IF
ABS[t.a[ii][t.unity.j] / t.a[ii][j]] < score
THEN {champ ← ii; score ← ABS[t.a[ii][t.unity.j] / t.a[ii][j]]}
ENDLOOP;
champion ← t.y[champ];
IF champ = v.i THEN {pec ← VIsManifestlyUnbounded; RETURN};
{newSampleValueOfV: REAL = t.a[champ][t.unity.j] - t.a[v.i][t.unity.j] * t.a[champ][u.j] / t.a[v.i][u.j];
SELECT
TRUE
FROM
newSampleValueOfV <= 0 => pec ← ToPivotAtUWouldIncreaseVButNotMakeItPositive;
newSampleValueOfV > 0 => pec ← ToPivotAtUWouldMakeVPositive
ENDCASE => ERROR}}};
Pivot:
PROCEDURE [t: Tableau, k: Row, l: Col] = {
pp: REAL ← 1.0 / t.a[k][l]; -- compute reciprocal of pivot
IF debug THEN Debug[t, IO.PutFR["Pivot row %g and column %g", IO.int[k], IO.int[l]]];
t.a[k][l] ← -1.0; -- to compensate for dividing by minus the pivot
-- swap the row and column headers
{z: Unknown ← t.y[k];t.y[k] ← t.x[l]; t.x[l] ← z;
t.y[k].i ← k; t.y[k].j ← doesntOwnACol; t.x[l].i ← doesntOwnARow; t.x[l].j ← l};
FOR j: INT IN [0 .. t.m) DO t.a[k][j] ← - t.a[k][j] * pp ENDLOOP;
-- divide everything in the pivot row by minus the pivot
FOR i:
INT
IN [0 .. t.n)
DO
IF i # k
THEN
FOR j:
INT
IN [0 .. t.m)
DO
IF j # l
THEN
-- i.e., for each t.a[i,j] outside the pivot row and column:
t.a[i][j] ← t.a[i][j] + t.a[i][l] * t.a[k][j];
-- note that t.a[k, j] was already divided by minus the pivot
ENDLOOP
ENDLOOP;
-- Finally process pivot column:
FOR i: INT IN [0 .. t.n) DO IF i # k THEN t.a[i][l] ← t.a[i][l] * pp ENDLOOP
};
Restrict:
PUBLIC PROCEDURE [t: Tableau, v: Unknown, propEQ: PropEQProc ←
NIL] = {
IF debug THEN Debug[t, IO.PutFR["Restrict variable %g", IO.refAny[v.name]]];
IF v.restricted THEN RETURN;
MakeVOwnARow[t, v];
Now v owns a row, or else it owns an empty column.
IF v.i = doesntOwnARow THEN {v.restricted ← TRUE; RETURN};
IF t.a[v.i][t.unity.j] < 0 THEN {v.restricted ← TRUE; RETURN};
{u, champion: Unknown;
pec: PivotEffectCode;
[u, champion, pec] ← MaximalPivot[t, v];
SELECT
TRUE
FROM
pec = ToPivotAtUWouldMakeVPositive => {v.restricted ← TRUE; Pivot[t, champion.i, u.j]};
pec = VIsMaximizedAtZero OR
pec = VIsManifestlyMaximized => {KillAllCols[t, champion, propEQ]};
pec = VIsMaximizedAndNegative => {IF signalUnsatisfiable THEN ERROR Unsatisfiable ELSE t.satisfiable ← FALSE};
pec = VIsManifestlyUnbounded => {v.restricted ← TRUE; Pivot[t, v.i, u.j]};
ENDCASE => ERROR}};
Maximize:
PUBLIC PROCEDURE [t: Tableau, v: Unknown]
RETURNS [unbounded:
BOOLEAN ←
FALSE] = {
u, champion: Unknown;
pec: PivotEffectCode;
IF debug THEN Debug[t, IO.PutFR["Maximize variable %g", IO.refAny[v.name]]];
MakeVOwnARow[t, v];
IF v.i = doesntOwnARow THEN RETURN[NOT v.restricted];
[u, champion, pec] ← FindPivot[t, v];
DO
SELECT
TRUE
FROM
pec = ToPivotAtUWouldMakeVPositive OR
pec = ToPivotAtUWouldIncreaseVButNotMakeItPositive =>
{Pivot[t, champion.i, u.j]; [u, champion, pec] ← FindPivot[t, v]};
ENDCASE => EXIT; ENDLOOP;
unbounded ← pec = VIsManifestlyUnbounded;
};
MakeVOwnARow:
PROCEDURE [t: Tableau, v: Unknown] = {
champ: INT;
score: REAL;
sgn: REAL;
i: INT ← 0;
IF v.j = doesntOwnACol THEN RETURN;
WHILE i < t.n AND t.a[i][v.j] = 0.0 DO i ← i + 1 ENDLOOP;
SELECT
TRUE
FROM
i = t.n => RETURN;
t.a[i][v.j] > 0 => sgn ← +1;
t.a[i][v.j] < 0 => sgn ← -1
ENDCASE => ERROR;
champ ← i;
score ← Real.LargestNumber;
FOR ii:
INT
IN [0..t.n)
DO
IF ii # i
AND t.y[ii].restricted
AND sgn * t.a[ii][v.j] < 0
THEN
IF
ABS[t.a[ii][t.unity.j] / t.a[ii][v.j]] < score
THEN {champ ← ii; score ← ABS[t.a[ii][t.unity.j] / t.a[ii][v.j]]}
ENDLOOP;
Pivot[t, champ, v.j];
};
Init:
PUBLIC
PROCEDURE
RETURNS [t: Tableau] = {
t ← NEW[TableauRec];
t.a ← NEW[ARRAY Row OF ARRAY Col OF REAL];
t.x ← NEW[ARRAY Col OF Unknown];
t.y ← NEW[ARRAY Row OF Unknown];
t.n ← t.m ← 0;
t.maxN ← t.maxM ← 0;
t.unity ← NewUnknown[t, "unity"];
t.satisfiable ← TRUE;
t.nameCount ← 0;
};
InventName:
PROC [t: Tableau, r: Rope.
ROPE]
RETURNS [Rope.
ROPE] ~ {
t.nameCount ← t.nameCount + 1;
RETURN [r.Concat[Convert.RopeFromInt[t.nameCount]]];
};
The rest of the code in this module is only for debugging.
in, out: IO.STREAM; -- for debugging
debug: BOOL ← FALSE;
StartDebug:
PROCEDURE =
{[in, out] ← ViewerIO.CreateViewerStreams["Debug LinearSolver"]; debug ← TRUE};
Debug:
PROCEDURE
[t: Tableau, r: Rope.
ROPE] = {
IF debug
THEN {
IO.Put[out, IO.rope[r]];
Newline[];
Print[t];
IO.PutF[out, "%12gSatisfiable = %g\n", IO.rope[""], IO.bool[t.satisfiable]];
Newline[]; Newline[]}};
Print:
PROCEDURE [t: Tableau] = {
IO.PutF[out, "%12g", IO.rope[""]];
FOR j : INT IN [0 .. t.m) DO Out[t.x[j]] ENDLOOP;
Newline[];
FOR i :
INT
IN [0 .. t.n)
DO Out[t.y[i]]; FOR j: INT IN [0 .. t.m) DO OutR[t.a[i][j]] ENDLOOP; Newline[]
ENDLOOP };
Out:
PROCEDURE [r: Unknown] = {
SELECT
TRUE
FROM
r.name = NIL => IO.PutF[out, "%12g", IO.rope["**"]];
r.restricted => IO.PutF[out, "%11g>", IO.refAny[r.name]];
NOT r.restricted => IO.PutF[out, "%12g", IO.refAny[r.name]];
ENDCASE => ERROR};
OutR: PROCEDURE [r: REAL] = {IO.PutF[out, "%12g", IO.real[r]]};
Newline: PROCEDURE = {IO.Put[out, IO.rope["\n"]]};