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; 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; END. „LinearSystemImpl.mesa Last edited by Maureen Stone 27-Oct-81 11:16:26 Written by Michael Plass, 8-Oct-81 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 Ê ‡˜Jšœ™Jšœ/™/Jšœ#™#J˜JšÏk œ˜Jšœœœ˜0J˜Jšœœ˜J˜š Ïnœœ œœœœ˜QšœÏc%˜+šœœœ˜Jšœœ˜šœœœ˜Jš œœœ œœ œ ˜2Jšœ˜—Jš œ œœœ œ ˜/Jšœœ(˜/šœœœ˜šœœ œ˜JšœœœœŸ8˜Ršœœœ˜Jšœ œ œ˜Jšœ˜—J˜Jšœ˜Jšœ˜——Jšœ˜—Jšœ-™-Jšœœ ˜Jš œœœœ œ˜-š œœ œœ˜%Jšœœ˜šœœœ˜Jšœ œ ˜Jšœ˜—Jšœ œ˜Jšœ˜ ——Jšœ˜J˜—š žœœ œœœ˜EJš˜šœœ˜ JšœŸ%˜+šœ œ˜J˜šœ œ˜Jš œœœ œœ œ ˜2Jšœ˜—Jš œ œœœ œŸ"˜RJšœœ(˜/šœ œ˜JšœœœœŸ#˜=šœ œ˜Jšœ œ œ˜Jšœ˜—J˜Jš˜—Jšœ˜—Jšœ-™-šœ œœ˜$Jšœœ˜šœ œ˜Jšœ œ ˜Jšœ˜—Jšœ œ˜Jšœ˜ —Jš˜—Jšœ˜J˜—š žœœ œœœ˜EJš˜šœœ˜ JšœŸ%˜+šœ œ˜J˜šœ œ˜Jš œœœ œœ œ ˜2Jšœ˜—Jš œ œœœ œŸ"˜RJšœœ(˜/šœ œ˜JšœœœœŸ#˜=šœ œ˜Jšœ œ œ˜Jšœ˜—J˜Jš˜—Jšœ˜—Jšœ-™-šœ œœ˜$Jšœœ˜šœ œ˜Jšœ œ ˜Jšœ˜—Jšœ œ˜Jšœ˜ —Jš˜—Jšœ˜J˜—š žœœ œœœ˜EJš˜šœœ˜ JšœŸ%˜+šœ œ˜J˜šœ œ˜Jš œœœ œœ œ ˜2Jšœ˜—Jš œ œœœ œŸ"˜RJšœœ(˜/šœ œ˜JšœœœœŸ#˜=šœ œ˜Jšœ œ œ˜Jšœ˜—J˜Jš˜—Jšœ˜—Jšœ-™-šœ œœ˜$Jšœœ˜šœ œ˜Jšœ œ ˜Jšœ˜—Jšœ œ˜Jšœ˜ —Jš˜—Jšœ˜J˜—š žœœ œœœ˜EJš˜šœœ˜ JšœŸ%˜+šœ œ˜J˜šœ œ˜Jš œœœ œœ œ ˜2Jšœ˜—Jš œ œœœ œŸ"˜RJšœœ(˜/šœ œ˜JšœœœœŸ#˜=šœ œ˜Jšœ œ œ˜Jšœ˜—J˜Jš˜—Jšœ˜—Jšœ-™-šœ œœ˜$Jšœœ˜šœ œ˜Jšœ œ ˜Jšœ˜—Jšœ œ˜Jšœ˜ —Jš˜—Jšœ˜J˜—š žœœ œœœ˜EJš˜šœœ˜ JšœŸ%˜+šœ œ˜J˜šœ œ˜Jš œœœ œœ œ ˜2Jšœ˜—Jš œ œœœ œŸ"˜RJšœœ(˜/šœ œ˜JšœœœœŸ#˜=šœ œ˜Jšœ œ œ˜Jšœ˜—J˜Jš˜—Jšœ˜—Jšœ-™-šœ œœ˜$Jšœœ˜šœ œ˜Jšœ œ ˜Jšœ˜—Jšœ œ˜Jšœ˜ —Jš˜—Jšœ˜J˜—Jšœ˜J˜—…—þ