<> <> <> <> DIRECTORY Convert, IO, LinearSolver, Real, Rope, ViewerIO; LinearSolverImpl: CEDAR PROGRAM IMPORTS Convert, IO, Rope, ViewerIO EXPORTS LinearSolver = BEGIN OPEN LinearSolver; <> <> <<>> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <<];>> <<>> doesntOwnARow: INT ~ -1; doesntOwnACol: INT ~ -1; <> <<>> <<(and i: i in [0, n): y(i) = (sum j: j in [0, m): a(i, j) * x(j))).>> <<>> <> <> <> <> <<>> <<(E j: j in [0, m): x(j) = unity).>> <<>> <> <> <> <<>> <<(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)>> <<>> <> <<>> <> <> <> <<>> <> 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] = { <> <> <> 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] = { <= 0, in two steps. Step 1: create a new unknown v, a new slack variable s (s >= 0), and assert v = c + s. >> <> <> 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] }; <> 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] = { <> <> <> 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 { <> 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 { <> 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] = { <> 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 }; <> <> 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]; <> 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]]]; }; <> 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"]]}; END. <> < test) , MakeVOwnARow could be called by Restrict when v already owns a row (inserted a test to return if v.j = doesntOwnACol), Maximize>> <<>> <> <> <<>> <<>> <<>> <> <> <<>> <> <> <<>> <> <> <<>> <> <> <> <> <> <> <<>> <<>> <<>> <<>> <> <> <> <> <<>> <> <> <<>>