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