-- LinearSolverImpl.mesa
-- Last Edited by: Gnelson, December 11, 1983 11:48 pm

DIRECTORY LinearSolver, Real, Rope, IO;

LinearSolverImpl: PROGRAM IMPORTS IO EXPORTS LinearSolver = BEGIN OPEN LinearSolver;

RowMax: INT = 100;
ColMax: INT = 100;
Row: TYPE = [0 .. RowMax);
Col: TYPE = [0 .. ColMax);
n: [0 .. RowMax];
m: [0 .. ColMax];
a: REF ARRAY Row OF ARRAY Col OF REAL
← NEW[ARRAY Row OF ARRAY Col OF REAL];

-- m[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 100 by 100.

x: REF ARRAY Col OF Unknown ← NEW[ARRAY Col OF Unknown];
y: REF ARRAY Row OF Unknown ← NEW[ARRAY Row OF Unknown];
unity: PUBLIC Unknown;

-- 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.

satisfiable: BOOL;

Init: PUBLIC PROC =
BEGIN n ← 0; m ← 0; unity ← NewUnknown["unity"]; satisfiable ← TRUE END;

NewUnknown: PUBLIC PROC [name: REF] RETURNS [Unknown] =
BEGIN
v: Unknown ← NEW[UnknownRec];
Debug["NewUnknown"];
v.name ← name;
v.i ← -1;
v.j ← m;
FOR i: INT IN [0 .. n) DO a[i][m] ← 0 ENDLOOP;
x[m] ← v;
m ← m + 1;
RETURN [v]
END;

-- 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.

Assert: PUBLIC PROC [f: Constraint] =
-- Asserts f = 0, in two steps. Step 1: create a new unknown c and assert c = f.
-- This is done by creating a new row owned by c that represents f.
-- Step 2: Assert c = 0. This is done by pivoting c into a column and deleting the column
BEGIN
c: Unknown ← NEW[UnknownRec];
Debug["Assert"];
c.i ← n;
c.j ← -1;
y[n] ← c;
FOR j: INT IN [0 .. m) DO a[n][j] ← 0 ENDLOOP;
WHILE f # NIL DO
{g: REF LinearMonomial ← NARROW[f.first];
AddToRow[n, g.coefficient, g.unknown]};
f ← f.rest
ENDLOOP;
n ← n + 1;
Kill[c]
END;

AddToRow: PROC[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)
BEGIN
SELECT TRUE FROM
v.j # -1 => a[n][v.j] ← a[n][v.j] + c;
v.i # -1 => FOR j: INT IN [0 .. m) DO a[n][j] ← a[n][j] + a[v.i][j] * c ENDLOOP;
ENDCASE => ERROR -- Client called Assert[f] where f contains an inactive variable
END;

Kill: PROC[c: Unknown] =
-- Add the constraint c = 0, by pivoting c into a column and deleting the column
-- If there is no place to pivot, then the constraints are unsatisfiable.
BEGIN
j: INT ← 0;
epsilon: REAL = .00001;
Debug["Kill"];
WHILE ABS[a[c.i][j]] < epsilon OR x[j] = unity DO j ← j + 1 ENDLOOP;
-- Pivot can't be too tiny, and can't be in unit column
IF ABS[a[c.i][j]] < epsilon THEN
BEGIN
IF ABS[a[c.i][unity.j]] > epsilon THEN satisfiable ← FALSE;
RETURN
END;
Pivot[c.i, j];
-- Now c owns a column instead of a row. Delete the column.
FOR i: INT IN [0 .. n) DO a[i][c.j] ← a[i][m - 1] ENDLOOP;
x[c.j] ← x[m-1];
x[c.j].j ← c.j;
m ← m - 1;
c.j ← -1
END;

Satisfiable: PUBLIC PROC[] RETURNS [BOOL] = {RETURN [satisfiable]};

Solution: PUBLIC PROC[v: Unknown] RETURNS [REAL] =
BEGIN
SELECT TRUE FROM
v = unity => RETURN[1];
v.j # -1 AND v # unity => RETURN[0];
v.i # -1 => RETURN[a[v.i][unity.j]]
ENDCASE => ERROR
END;

Pivot: PROC[k: Row, l: Col] =
BEGIN
pp: REAL ← 1.0 / a[k][l]; -- compute reciprocal of pivot
Debug["Pivot"];
a[k][l] ← 1.0;
{t: Unknown ← y[k]; y[k] ← x[l]; x[l] ← t; y[k].i ← k; y[k].j ← -1; x[l].i ← -1; x[l].j ← l};
-- swap the row and column headers
FOR j: INT IN [0 .. m) DO a[k][j] ← - a[k][j] * pp ENDLOOP;
-- divide everything in the pivot row by minus the pivot
FOR i: INT IN [0 .. n) DO
IF i # k THEN
  FOR j: INT IN [0 .. m) DO
  IF j # l THEN
  -- i.e., for each a[i,j] outside the pivot row and column:
  a[i][j] ← a[i][j] + a[i][l] * a[k][j];
  -- note that a[k, j] was already divided by minus the pivot
  ENDLOOP
  ENDLOOP;
  -- Finally process pivot column:
  FOR i: INT IN [0 .. n) DO IF i # k THEN a[i][l] ← -a[i][l] * pp ENDLOOP
END;

-- The rest of the code in this module is only for debugging.

in, out: IO.STREAM; -- for debugging
debug: BOOL ← FALSE;

StartDebug: PROC =
{[in, out] ← IO.CreateViewerStreams["Debug LinearSolver"]; debug ← TRUE};

Debug: PROC[r: Rope.ROPE] = {
IF debug THEN {
IO.Put[out, IO.rope[r]];
IO.Put[out, IO.rope["\n\n"]];
Print[];
IO.Put[out, IO.rope["\n"]]}};

Print: PROC = {
IO.PutF[out, "%12g", IO.rope[""]];
FOR j : INT IN [0 .. m) DO Out[x[j]] ENDLOOP;
Newline[];
FOR i : INT IN [0 .. n) DO
Out[y[i]]; FOR j: INT IN [0 .. m) DO OutR[a[i][j]] ENDLOOP; Newline[]
ENDLOOP };

Out: PROC[r: Unknown] = {
SELECT TRUE FROM
r.name = NIL => IO.PutF[out, "%12g", IO.rope["**"]];
TRUE => IO.PutF[out, "%12g", IO.refAny[r.name]];
ENDCASE => ERROR};

OutR: PROC[r: REAL] = {IO.PutF[out, "%12g", IO.real[r]]};

Newline: PROC = {IO.Put[out, IO.rope["\n"]]};


END.