PROGRAM Hilbert;

{Uses the decompose and solve procedures from the linear equations
 solver (lineq.pas) to invert matrices.  Uses a Hilbert matrix of
 rank 1..5 as an example.  Hilbert matrices are ill-conditioned and
 difficult to invert accurately.  This program inverts a Hilbert
 matrix, and then it inverts the inverse.}

CONST
    SIZE = 5;

TYPE
    vector = ARRAY [1..SIZE]         OF real;
    matrix = ARRAY [1..SIZE, 1..SIZE] OF real;
    
    why = (SINGULARMATRIX, ZEROROW, NOCONVERGENCE);

VAR
    ps   : ARRAY [1..SIZE] OF integer;  {global pivot index array}
    H    : matrix;                     {Hilbert matrix}
    Hinv : matrix;                     {Hilbert matrix inverse}
    Z    : matrix;                     {Identity matrix}
    i, j : integer;

PROCEDURE singular(w : why);

    BEGIN
        CASE w OF
            SINGULARMATRIX : writeln('Singular matrix in decompose.');
            ZEROROW        : writeln('Matrix with zero row in decompose.');
            NOCONVERGENCE  : writeln('No convergence in improve.');
        END;
    END;
    
FUNCTION absvalue(x : real) : real;

    BEGIN
        IF x < 0 THEN x := -x;
        absvalue := x;
    END;

PROCEDURE decompose(n : integer; VAR A, LU : matrix);

    {Computes triangular matrices L and U and permutation matrix P
     so that LU = PA.  Stores L-I and U in LU.  Vector ps contains
     permuted row indices.  Note that A and LU are often passed the
     same matrix.}

    VAR
        scales : vector;
        i, j, k, pivotindex : integer;
        normrow, pivot, size, biggest, mult : real;

    BEGIN

        {Initialize ps, LU, and scales.}
        FOR i := 1 TO n DO BEGIN  {rows}
            ps[i] := i;
            normrow := 0;

            FOR j := 1 TO n DO BEGIN  {columns}
                LU[i,j] := A[i,j];

                {Find the largest row element.}
                IF normrow < absvalue(LU[i,j]) THEN normrow := absvalue(LU[i,j]);
            END;

            {Set the scaling factor for row equilibration.}
            IF normrow <> 0 THEN scales[i] := 1/normrow
            ELSE BEGIN
                scales[i] := 0;
                singular(ZEROROW);
            END;
        END;

        {Gaussian elimination with partial pivoting.}
        FOR k := 1 TO n-1 DO BEGIN  {pivot row k}
            pivotindex := 0;
            biggest := 0;

            {Go down rows from row k.}
            FOR i := k TO n DO BEGIN

                {Divide by the largest row element.}
                size := absvalue(LU[ps[i], k])*scales[ps[i]];

                IF biggest < size THEN BEGIN
                    biggest    := size;  {biggest scales column element}
                    pivotindex := i;     {row index of this element}
                END;
            END;

            IF biggest = 0 THEN singular(SINGULARMATRIX)
            ELSE BEGIN
                IF pivotindex <> k THEN BEGIN

                    {Exchange rows.}
                    j := ps[k];
                    ps[k] := ps[pivotindex];
                    ps[pivotindex] := j;
                END;

                pivot := LU[ps[k], k];  {pivot element}

                {Go down rest of rows.}
                FOR i := k+1 TO n DO BEGIN
                    mult := LU[ps[i], k]/pivot;
                    LU[ps[i], k] := mult;

                    IF mult <> 0 THEN BEGIN

                        {Inner loop.  Only column subscript varies.}
                        FOR j := k+1 TO n DO BEGIN
                            LU[ps[i], j] := LU[ps[i], j]
                                                - mult*LU[ps[k], j];
                        END;
                    END;
                END;
            END;
        END;

        {Check bottom right element of permuted matrix.}
        IF LU[ps[n], n] = 0 THEN singular(SINGULARMATRIX);
    END;

PROCEDURE solve(n : integer; VAR LU : matrix; VAR b, x : vector);

    {Solves Ax = b using LU from decompose.}

    VAR
        i, j : integer;
        dot  : real;

    BEGIN

        {Ly = b : solve for y.}
        FOR i := 1 TO n DO BEGIN
            dot := 0;
            FOR j := 1 TO i-1 DO BEGIN
                dot := dot + LU[ps[i], j]*x[j];
            END;
            x[i] := b[ps[i]] - dot;
        END;

        {Ux = y : solve for x.}
        FOR i := n DOWNTO 1 DO BEGIN
            dot := 0;
            FOR j := i+1 TO n DO BEGIN
                dot := dot + LU[ps[i], j]*x[j];
            END;
            x[i] := (x[i] - dot)/LU[ps[i], i];
        END;
    END;

PROCEDURE invert(n : integer; VAR A, Ainv : matrix);

    {Compute Ainv = inverse(A).  Note that A and Ainv are often
     passed the same matrix.}

    VAR
        LU :   matrix;
        b, x : vector;
        i, j : integer;

    BEGIN
        decompose(n, A, LU);

        FOR j := 1 TO n DO BEGIN
            FOR i := 1 TO n DO BEGIN
                IF i = j THEN b[i] := 1
                         ELSE b[i] := 0;
            END;

            solve(n, LU, b, x);

            FOR i := 1 TO n DO Ainv[i,j] := x[i];
        END;
    END;
    
PROCEDURE multiply(n : integer; VAR A, B, P : matrix);

    {Compute P = A x B.}
    
    VAR
        i, j, k : integer;
        sum : real;
        
    BEGIN
        sum := 0.0;
        
        FOR i := 1 TO n DO BEGIN
            FOR j := 1 TO n DO BEGIN
                sum := 0.0;
                
                FOR k := 1 TO n DO BEGIN
                    sum := sum + A[i,k]*B[k,j];
                END;
                
                P[i,j] := sum;
            END;
        END;
    END;

PROCEDURE printmatrix(VAR M : matrix);

    VAR
        i, j : integer;

    BEGIN
        writeln;
        FOR i := 1 TO SIZE DO BEGIN
            FOR j := 1 TO SIZE DO write(M[i,j]:15:6);
            writeln;
        END;
        writeln;
    END;

BEGIN
    {Compute the Hilbert matrix.}
    FOR i := 1 TO SIZE DO BEGIN
        FOR j := 1 TO SIZE DO BEGIN
            H[i,j] := 1.0/(i + j - 1);
        END;
    END;

    writeln;
    writeln('Hilbert matrix:');
    printmatrix(H);

    {Invert the Hilbert matrix.}
    invert(SIZE, H, Hinv);

    writeln('Hilbert matrix inverted:');
    printmatrix(Hinv);
    
    {Multiply the Hilbert matrix by its inverse.}
    multiply(SIZE, H, Hinv, Z);
    
    writeln('Hilbert matrix multiplied by its inverse:');
    printmatrix(Z);

    {Invert the inverse.}
    invert(SIZE, Hinv, H);

    writeln('Inverse matrix inverted:');
    printmatrix(H);
END.