LinearSystemImpl.mesa
Last edited by Maureen Stone 27-Oct-81 11:16:26
Written by Michael Plass, 8-Oct-81
DIRECTORY LinearSystem;
LinearSystemImpl: PROGRAM EXPORTS LinearSystem =
BEGIN OPEN LinearSystem;
SolveN: PUBLIC PROCEDURE [A:MatrixN, b:ColumnN, n: INTEGER] RETURNS [x:ColumnN] =
BEGIN -- solve Ax=b by Gaussian Elimination
FOR i:INTEGER IN [0..n) DO
bestk:INTEGER ← i;
FOR k:INTEGER IN [i..n) DO
IF ABS[A[k][i]] > ABS[A[bestk][i]] THEN bestk ← k;
ENDLOOP;
{t:RowN ← A[i]; A[i] ← A[bestk]; A[bestk] ← t};
{t:REAL ← b[i]; b[i] ← b[bestk]; b[bestk] ← t};
FOR k:INTEGER IN (i..n) DO
IF A[k][i]#0 THEN BEGIN
r:REAL = A[k][i]/A[i][i]; -- Singular A causes Real.RealException = divide by zero
FOR j:INTEGER IN [i..n) DO
A[k][j] ← A[k][j] - A[i][j]*r
ENDLOOP;
b[k] ← b[k] - b[i]*r
END;
ENDLOOP;
ENDLOOP;
Now A is upper-triangular, so back substitute
x ← NEW[VecSeq[n]];
FOR i:INTEGER IN [0..n) DO x[i] ← 0; ENDLOOP;
FOR i:INTEGER DECREASING IN [0..n) DO
xi:REAL ← b[i];
FOR j:INTEGER IN (i..n) DO
xi ← xi - A[i][j]*x[j];
ENDLOOP;
x[i] ← xi / A[i][i]
ENDLOOP
END;
Solve2: PUBLIC PROCEDURE [A:Matrix2, b:Column2] RETURNS [x:Column2] =
BEGIN
n:NAT = 2;
BEGIN -- solve Ax=b by Gaussian Elimination
FOR i:[1..n] IN [1..n] DO
bestk:[1..n] ← i;
FOR k:[1..n] IN [i..n] DO
IF ABS[A[k][i]] > ABS[A[bestk][i]] THEN bestk ← k;
ENDLOOP;
{t:Row2 ← A[i]; A[i] ← A[bestk]; A[bestk] ← t}; -- sorry about the dependence on n
{t:REAL ← b[i]; b[i] ← b[bestk]; b[bestk] ← t};
FOR k:(1..n] IN (i..n] DO
r:REAL = A[k][i]/A[i][i]; -- Singular A causes divide by zero
FOR j:[1..n] IN [i..n] DO
A[k][j] ← A[k][j] - A[i][j]*r
ENDLOOP;
b[k] ← b[k] - b[i]*r
ENDLOOP
ENDLOOP;
Now A is upper-triangular, so back substitute
FOR i:[1..n] DECREASING IN [1..n] DO
xi:REAL ← b[i];
FOR j:[1..n] IN (i..n] DO
xi ← xi - A[i][j]*x[j];
ENDLOOP;
x[i] ← xi / A[i][i]
ENDLOOP
END
END;
Solve3: PUBLIC PROCEDURE [A:Matrix3, b:Column3] RETURNS [x:Column3] =
BEGIN
n:NAT = 3;
BEGIN -- solve Ax=b by Gaussian Elimination
FOR i:[1..n] IN [1..n] DO
bestk:[1..n] ← i;
FOR k:[1..n] IN [i..n] DO
IF ABS[A[k][i]] > ABS[A[bestk][i]] THEN bestk ← k;
ENDLOOP;
{t:Row3 ← A[i]; A[i] ← A[bestk]; A[bestk] ← t}; -- sorry about the dependence on n
{t:REAL ← b[i]; b[i] ← b[bestk]; b[bestk] ← t};
FOR k:(1..n] IN (i..n] DO
r:REAL = A[k][i]/A[i][i]; -- Singular A causes divide by zero
FOR j:[1..n] IN [i..n] DO
A[k][j] ← A[k][j] - A[i][j]*r
ENDLOOP;
b[k] ← b[k] - b[i]*r
ENDLOOP
ENDLOOP;
Now A is upper-triangular, so back substitute
FOR i:[1..n] DECREASING IN [1..n] DO
xi:REAL ← b[i];
FOR j:[1..n] IN (i..n] DO
xi ← xi - A[i][j]*x[j];
ENDLOOP;
x[i] ← xi / A[i][i]
ENDLOOP
END
END;
Solve4: PUBLIC PROCEDURE [A:Matrix4, b:Column4] RETURNS [x:Column4] =
BEGIN
n:NAT = 4;
BEGIN -- solve Ax=b by Gaussian Elimination
FOR i:[1..n] IN [1..n] DO
bestk:[1..n] ← i;
FOR k:[1..n] IN [i..n] DO
IF ABS[A[k][i]] > ABS[A[bestk][i]] THEN bestk ← k;
ENDLOOP;
{t:Row4 ← A[i]; A[i] ← A[bestk]; A[bestk] ← t}; -- sorry about the dependence on n
{t:REAL ← b[i]; b[i] ← b[bestk]; b[bestk] ← t};
FOR k:(1..n] IN (i..n] DO
r:REAL = A[k][i]/A[i][i]; -- Singular A causes divide by zero
FOR j:[1..n] IN [i..n] DO
A[k][j] ← A[k][j] - A[i][j]*r
ENDLOOP;
b[k] ← b[k] - b[i]*r
ENDLOOP
ENDLOOP;
Now A is upper-triangular, so back substitute
FOR i:[1..n] DECREASING IN [1..n] DO
xi:REAL ← b[i];
FOR j:[1..n] IN (i..n] DO
xi ← xi - A[i][j]*x[j];
ENDLOOP;
x[i] ← xi / A[i][i]
ENDLOOP
END
END;
Solve5: PUBLIC PROCEDURE [A:Matrix5, b:Column5] RETURNS [x:Column5] =
BEGIN
n:NAT = 5;
BEGIN -- solve Ax=b by Gaussian Elimination
FOR i:[1..n] IN [1..n] DO
bestk:[1..n] ← i;
FOR k:[1..n] IN [i..n] DO
IF ABS[A[k][i]] > ABS[A[bestk][i]] THEN bestk ← k;
ENDLOOP;
{t:Row5 ← A[i]; A[i] ← A[bestk]; A[bestk] ← t}; -- sorry about the dependence on n
{t:REAL ← b[i]; b[i] ← b[bestk]; b[bestk] ← t};
FOR k:(1..n] IN (i..n] DO
r:REAL = A[k][i]/A[i][i]; -- Singular A causes divide by zero
FOR j:[1..n] IN [i..n] DO
A[k][j] ← A[k][j] - A[i][j]*r
ENDLOOP;
b[k] ← b[k] - b[i]*r
ENDLOOP
ENDLOOP;
Now A is upper-triangular, so back substitute
FOR i:[1..n] DECREASING IN [1..n] DO
xi:REAL ← b[i];
FOR j:[1..n] IN (i..n] DO
xi ← xi - A[i][j]*x[j];
ENDLOOP;
x[i] ← xi / A[i][i]
ENDLOOP
END
END;
Solve6: PUBLIC PROCEDURE [A:Matrix6, b:Column6] RETURNS [x:Column6] =
BEGIN
n:NAT = 6;
BEGIN -- solve Ax=b by Gaussian Elimination
FOR i:[1..n] IN [1..n] DO
bestk:[1..n] ← i;
FOR k:[1..n] IN [i..n] DO
IF ABS[A[k][i]] > ABS[A[bestk][i]] THEN bestk ← k;
ENDLOOP;
{t:Row6 ← A[i]; A[i] ← A[bestk]; A[bestk] ← t}; -- sorry about the dependence on n
{t:REAL ← b[i]; b[i] ← b[bestk]; b[bestk] ← t};
FOR k:(1..n] IN (i..n] DO
r:REAL = A[k][i]/A[i][i]; -- Singular A causes divide by zero
FOR j:[1..n] IN [i..n] DO
A[k][j] ← A[k][j] - A[i][j]*r
ENDLOOP;
b[k] ← b[k] - b[i]*r
ENDLOOP
ENDLOOP;
Now A is upper-triangular, so back substitute
FOR i:[1..n] DECREASING IN [1..n] DO
xi:REAL ← b[i];
FOR j:[1..n] IN (i..n] DO
xi ← xi - A[i][j]*x[j];
ENDLOOP;
x[i] ← xi / A[i][i]
ENDLOOP
END
END;
END.