LinearSystemImpl.mesa
Copyright Ó 1985, 1986, 1987, 1992 by Xerox Corporation. All rights reserved.
Last edited by Maureen Stone 27-Oct-81 11:16:26
Written by Michael Plass, 8-Oct-81
Tim Diebert May 21, 1985 5:51:30 pm PDT
Stone, October 16, 1985 11:22:42 am PDT
Christian LeCocq March 13, 1987 5:56:05 pm PST
Maureen Stone, March 24, 1987 11:23:59 am PST
DIRECTORY LinearSystem;
LinearSystemImpl: CEDAR 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
ai: RowN;
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};
ai ¬ A[bestk];
A[bestk] ¬ A[i];
A[i] ¬ ai;
{t:REAL ¬ b[i]; b[i] ¬ b[bestk]; b[bestk] ¬ t};
FOR k:
INTEGER
IN (i..n)
DO
ak: RowN ¬ A[k];
IF ak[i]#0
THEN
BEGIN
FOR j:INTEGER IN [i..n) DO
aij: REAL ← ai[j];
IF aij#0.0 THEN ak[j] ← ak[j] - aij*r
ENDLOOP;
r:REAL;
TRUSTED {
akjp: LONG POINTER TO REAL ¬ @ak[i];
limitKJ: LONG POINTER TO REAL ¬ @ak[n-1];
aijp: LONG POINTER TO REAL ¬ @ai[i];
r ¬ akjp/aijp; -- Singular A causes Real.RealException = divide by zero
DO
IF aijp#0.0 THEN akjp ¬ akjp - aijp*r;
IF akjp = limitKJ THEN EXIT;
aijp ¬ aijp + SIZE[REAL];
akjp ¬ akjp + SIZE[REAL];
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: [0..n)
IN [0..n)
DO
bestk: [0..n) ¬ i;
FOR k: [0..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:(0..n)
IN (i..n)
DO
r:REAL = A[k][i]/A[i][i]; -- Singular A causes divide by zero
FOR j: [0..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: [0..n)
DECREASING
IN [0..n)
DO
xi:REAL ¬ b[i];
FOR j: [0..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: [0..n)
IN [0..n)
DO
bestk: [0..n) ¬ i;
FOR k: [0..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:(0..n)
IN (i..n)
DO
r:REAL = A[k][i]/A[i][i]; -- Singular A causes divide by zero
FOR j: [0..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: [0..n)
DECREASING
IN [0..n)
DO
xi:REAL ¬ b[i];
FOR j: [0..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: [0..n)
IN [0..n)
DO
bestk: [0..n) ¬ i;
FOR k: [0..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:(0..n)
IN (i..n)
DO
r:REAL = A[k][i]/A[i][i]; -- Singular A causes divide by zero
FOR j: [0..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: [0..n)
DECREASING
IN [0..n)
DO
xi:REAL ¬ b[i];
FOR j: [0..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: [0..n)
IN [0..n)
DO
bestk: [0..n) ¬ i;
FOR k: [0..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:(0..n)
IN (i..n)
DO
r:REAL = A[k][i]/A[i][i]; -- Singular A causes divide by zero
FOR j: [0..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: [0..n)
DECREASING
IN [0..n)
DO
xi:REAL ¬ b[i];
FOR j: [0..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: [0..n)
IN [0..n)
DO
bestk: [0..n) ¬ i;
FOR k: [0..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:(0..n)
IN (i..n)
DO
r:REAL = A[k][i]/A[i][i]; -- Singular A causes divide by zero
FOR j: [0..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: [0..n)
DECREASING
IN [0..n)
DO
xi:REAL ¬ b[i];
FOR j: [0..n)
IN (i..n)
DO
xi ¬ xi - A[i][j]*x[j];
ENDLOOP;
x[i] ¬ xi / A[i][i]
ENDLOOP
END
END;
InvalidMatrix: PUBLIC SIGNAL = CODE;
InvalidOperation: PUBLIC SIGNAL = CODE;
ValidMatrix:
PROCEDURE [a: MatrixN]
RETURNS [nrows,ncols:
INTEGER] = {
IF a=NIL OR a.nrows=0 THEN SIGNAL InvalidMatrix;
nrows ¬ a.nrows;
ncols ¬ a[0].ncols;
FOR i:
INTEGER
IN [0..nrows)
DO
IF a[i].ncols#ncols THEN SIGNAL InvalidMatrix;
ENDLOOP;
RETURN[nrows,ncols];
};
Sign: PROC[i,j: INTEGER] RETURNS[REAL] = {RETURN[IF (i+j) MOD 2 = 0 THEN 1 ELSE -1]};
MakeAij:
PROC[a, Aij: MatrixN, i,j:
INTEGER] = {
--assume Aij is well formed
n,m: INTEGER ¬ 0; --row index, column index for new matrix
FOR row:
INTEGER
IN [0..a.nrows)
DO
IF row=i THEN LOOP; --row index for original matrix
m ¬ 0;
FOR col:
INTEGER
IN [0..a[0].ncols)
DO
--column index for new matrix
IF col=j THEN LOOP;
Aij[n][m] ¬ a[row][col]; --column 0 of old matrix supplies the aij values
m ¬ m+1;
ENDLOOP;
n ¬ n+1;
ENDLOOP;
};
Invert:
PUBLIC PROCEDURE [a: MatrixN]
RETURNS [ai: MatrixN] = {
nrows,ncols: INTEGER;
det: REAL;
Aij: MatrixN;
[nrows,ncols] ¬ ValidMatrix[a];
IF nrows#ncols THEN SIGNAL InvalidOperation;
det ¬ Determinant[a];
ai ¬ Create[nrows,ncols];
Aij ¬ Create[nrows-1,ncols-1];
FOR i:
INTEGER
IN [0..nrows)
DO
FOR j:
INTEGER
IN [0..ncols)
DO
MakeAij[a,Aij,i,j];
ai[j][i] ¬ Sign[i,j]*Determinant[Aij]/det;
ENDLOOP;
ENDLOOP;
RETURN[ai];
};
Determinant:
PUBLIC PROC[a: MatrixN]
RETURNS [det:
REAL] = {
nrows,ncols: INTEGER;
[nrows,ncols] ¬ ValidMatrix[a];
IF nrows#ncols THEN SIGNAL InvalidOperation;
IF nrows=1 THEN det ¬ a[0][0]
ELSE IF nrows=2 THEN det ¬ a[0][0]*a[1][1]-a[0][1]*a[1][0]
ELSE {
i,j: INTEGER;
Aij: MatrixN ¬ Create[nrows-1,ncols-1];
det ¬ 0;
j ¬ 0; --always use column 0 for now
FOR i
IN [0..nrows)
DO
MakeAij[a,Aij,i,j];
det ¬ det + a[i][j]*Sign[i,j]*Determinant[Aij];
ENDLOOP;
};
RETURN[det];
};
Create:
PUBLIC
PROC [nrows, ncols:
INTEGER]
RETURNS [a: MatrixN] = {
a ¬ NEW[MatrixSeq[nrows]];
FOR i: INTEGER IN [0..nrows) DO a[i] ¬ NEW[VecSeq[ncols]]; ENDLOOP;
};
Copy:
PUBLIC
PROC [a: MatrixN]
RETURNS[MatrixN] = {
nrows,ncols: INTEGER;
new: MatrixN;
[nrows,ncols] ¬ ValidMatrix[a];
new ¬ Create[nrows, ncols];
FOR i:
INTEGER
IN [0..nrows)
DO
FOR j: INTEGER IN [0..ncols) DO new[i][j] ¬ a[i][j]; ENDLOOP;
ENDLOOP;
RETURN[new];
};
Transpose:
PUBLIC
PROCEDURE [a: MatrixN]
RETURNS [transpose: MatrixN] = {
nrows,ncols: INTEGER;
[nrows,ncols] ¬ ValidMatrix[a];
transpose ¬ Create[ncols,nrows];
FOR i:
INTEGER
IN [0..nrows)
DO
FOR j:
INTEGER
IN [0..ncols)
DO
transpose[j][i] ¬ a[i][j];
ENDLOOP;
ENDLOOP;
RETURN[transpose];
};
Multiply:
PUBLIC PROCEDURE [a: MatrixN, b: MatrixN]
RETURNS [c: MatrixN] = {
nra,nca: INTEGER;
nrb,ncb: INTEGER;
[nra,nca] ¬ ValidMatrix[a];
[nrb,ncb] ¬ ValidMatrix[b];
IF nca#nrb THEN SIGNAL InvalidOperation;
c ¬ Create[nra,ncb];
FOR j:
INTEGER
IN [0..nra)
DO
FOR i:
INTEGER
IN [0..ncb)
DO
c[j][i] ¬ 0;
FOR k:
INTEGER
IN [0..nca)
DO
c[j][i] ¬ c[j][i] + a[j][k]*b[k][i];
ENDLOOP;
ENDLOOP;
ENDLOOP;
};
MultiplyVec:
PUBLIC
PROC[a: MatrixN, v: ColumnN]
RETURNS [c: RowN] = {
nrows,ncols,nels: INTEGER;
[nrows,ncols] ¬ ValidMatrix[a];
IF v=NIL THEN SIGNAL InvalidMatrix;
nels ¬ v.ncols; --actually nrows for this case. historical
IF nels#ncols THEN SIGNAL InvalidOperation;
c ¬ NEW[VecSeq[nrows]];
FOR i:
INTEGER
IN [0..nrows)
DO
c[i] ¬ 0;
FOR j:
INTEGER
IN [0..nels)
DO
c[i] ¬ c[i]+a[i][j]*v[j];
ENDLOOP;
ENDLOOP;
};
END.