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 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; 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; 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; 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; 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; 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; 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; 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. LinearSystemImpl.mesa Copyright c 1985 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 Now A is upper-triangular, so back substitute Now A is upper-triangular, so back substitute Now A is upper-triangular, so back substitute Now A is upper-triangular, so back substitute Now A is upper-triangular, so back substitute Now A is upper-triangular, so back substitute Κ˜šœ™Icodešœ Οmœ1™