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; 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 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; 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; 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; 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; 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; 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; 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.  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 {t:RowN _ A[i]; A[i] _ A[bestk]; A[bestk] _ t}; FOR j:INTEGER IN [i..n) DO aij: REAL _ ai[j]; IF aij#0.0 THEN ak[j] _ ak[j] - aij*r ENDLOOP; 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 Κ—•NewlineDelimiter –(cedarcode) style™šœ™Icodešœ ΟeœC™NJšœ/™/Jšœ#™#K™'K™'K™.K™-—K˜KšΟk œ˜K˜Kš Οnœž œžœžœžœ˜OK˜š Ÿœžœž œžœžœžœ˜QšžœΟc%˜+Kšœ ˜ šžœžœžœž˜Kšœžœ˜šžœžœžœž˜Kš žœžœžœ žœžœ žœ ˜2Kšžœ˜—Jš œ žœžœžœ žœ ™/Kšœžœ˜Kšžœ žœ˜Kšžœ ˜ Kšœžœ(˜/šžœžœžœž˜Kšœ žœ˜šžœ žœž˜šžœžœžœž™Jšœžœ ™Jšžœ žœ™%Jšžœ™—Kšœžœ˜šžœ˜ Kš œžœžœžœžœ ˜$Kš œ žœžœžœžœ ˜)Kš œžœžœžœžœ ˜$Kšœ 8˜Išž˜Kšžœ žœ˜*Kšžœžœžœ˜Kšœžœžœ˜Kšœžœžœ˜Kšžœ˜—Kšœ˜—K˜Kšžœ˜Kšžœ˜——Kšžœ˜—Jšœ-™-Kšœžœ ˜Kš žœžœžœžœ žœ˜-š žœžœž œžœž˜%Kšœžœ˜šžœžœžœž˜Kšœ žœ ˜Kšžœ˜—Kšœ žœ˜Kšžœ˜ ——Kšžœ˜K˜—š Ÿœžœž œžœžœ˜EKšž˜šœžœ˜ Kšžœ %˜+šžœ žœž˜K˜šžœ žœž˜Kš žœžœžœ žœžœ žœ ˜2Kšžœ˜—Kš œ žœžœžœ žœ "˜RKšœžœ(˜/šžœ žœž˜Kšœžœžœžœ #˜=šžœ žœž˜Kšžœ žœ žœ˜Kšžœ˜—K˜Kšž˜—Kšžœ˜—Jšœ-™-šžœ ž œžœž˜%Kšœžœ˜šžœ žœž˜Kšœ žœ ˜Kšžœ˜—Kšœ žœ˜Kšžœ˜ —Kšž˜—Kšžœ˜K˜—š Ÿœžœž œžœžœ˜EKšž˜šœžœ˜ Kšžœ %˜+šžœ žœž˜K˜šžœ žœž˜Kš žœžœžœ žœžœ žœ ˜2Kšžœ˜—Kš œ žœžœžœ žœ "˜RKšœžœ(˜/šžœ žœž˜Kšœžœžœžœ #˜=šžœ žœž˜Kšžœ žœ žœ˜Kšžœ˜—K˜Kšž˜—Kšžœ˜—Jšœ-™-šžœ ž œžœž˜%Kšœžœ˜šžœ žœž˜Kšœ žœ ˜Kšžœ˜—Kšœ žœ˜Kšžœ˜ —Kšž˜—Kšžœ˜K˜—š Ÿœžœž œžœžœ˜EKšž˜šœžœ˜ Kšžœ %˜+šžœ žœž˜K˜šžœ žœž˜Kš žœžœžœ žœžœ žœ ˜2Kšžœ˜—Kš œ žœžœžœ žœ "˜RKšœžœ(˜/šžœ žœž˜Kšœžœžœžœ #˜=šžœ žœž˜Kšžœ žœ žœ˜Kšžœ˜—K˜Kšž˜—Kšžœ˜—Jšœ-™-šžœ ž œžœž˜%Kšœžœ˜šžœ žœž˜Kšœ žœ ˜Kšžœ˜—Kšœ žœ˜Kšžœ˜ —Kšž˜—Kšžœ˜K˜—š Ÿœžœž œžœžœ˜EKšž˜šœžœ˜ Kšžœ %˜+šžœ žœž˜K˜šžœ žœž˜Kš žœžœžœ žœžœ žœ ˜2Kšžœ˜—Kš œ žœžœžœ žœ "˜RKšœžœ(˜/šžœ žœž˜Kšœžœžœžœ #˜=šžœ žœž˜Kšžœ žœ žœ˜Kšžœ˜—K˜Kšž˜—Kšžœ˜—Jšœ-™-šžœ ž œžœž˜%Kšœžœ˜šžœ žœž˜Kšœ žœ ˜Kšžœ˜—Kšœ žœ˜Kšžœ˜ —Kšž˜—Kšžœ˜K˜—š Ÿœžœž œžœžœ˜EKšž˜šœžœ˜ Kšžœ %˜+šžœ žœž˜K˜šžœ žœž˜Kš žœžœžœ žœžœ žœ ˜2Kšžœ˜—Kš œ žœžœžœ žœ "˜RKšœžœ(˜/šžœ žœž˜Kšœžœžœžœ #˜=šžœ žœž˜Kšžœ žœ žœ˜Kšžœ˜—K˜Kšž˜—Kšžœ˜—Jšœ-™-šžœ ž œžœž˜%Kšœžœ˜šžœ žœž˜Kšœ žœ ˜Kšžœ˜—Kšœ žœ˜Kšžœ˜ —Kšž˜—Kšžœ˜K˜—Kšœž œžœ˜$Kšœž œžœ˜'K˜šŸ œž œžœžœ˜FKš žœžœžœ žœžœ˜0K˜K˜šžœžœžœ ž˜Kšžœžœžœ˜.Kšžœ˜—Kšžœ˜K˜—KšŸœžœžœžœžœžœžœžœžœžœ˜UšŸœžœžœ ˜LKšœžœ (˜:šžœžœžœž˜#Kšžœžœžœ ˜3K˜š žœžœžœžœ ˜DKšžœžœžœ˜Kšœ 0˜IK˜Kšžœ˜—K˜Kšžœ˜K˜——šŸœžœžœ˜?Kšœ žœ˜Kšœžœ˜ Kšœ ˜ K˜Kšžœ žœžœ˜,K˜K˜K˜šžœžœžœ ž˜šžœžœžœ ž˜K˜K˜*Kšžœ˜—Kšžœ˜—Kšžœ˜ K˜—šŸ œž œ žœžœ˜