<<>>
<<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.