DataMuseum.dk

top - metrics - download

    Length: 4261 (0x10a5)
    Types: TextFile
    Names: »B«

└─⟦180fe333a⟧ Bits:30000405 8mm tape, Rational 1000, SW CATALOG, 10_20_0
└─⟦180fe333a⟧ Bits:30000537 8mm tape, Rational 1000, SW Catalog 10_20_0
    └─⟦5cb1d1d7f⟧ »DATA« 
        └─⟦3b1ee7bd8⟧ 
            └─⟦this⟧ 

TextFile

--
-- Version: @(#)cauchfl.ada 1.1  Date: 6/3/84
--

with Text_Io;
use Text_Io;
procedure Cauchy is
    --
    --  This test of floating point accuracy based on computing the inverses
    --  of Cauchy matricies.  These are N x N matricies for which the i, jth
    --  entry is 1 / (i + j - 1).  The inverse is computed using determinants.
    --  As N increases, the determinant rapidly approaches zero.  The inverse
    --  is computed exactly and then checked by multiplying it by the original
    --  matrix.
    --
    --     Gerry Fisher
    --     Computer Sciences Corporation
    --     May 27, 1984

    type Real is digits 6;

    type Matrix is array (Positive range <>, Positive range <>) of Real;

    Trials : constant := 5;
    Failed : Boolean := False;

    function Cofactor (A : Matrix; I, J : Positive) return Matrix is
        B : Matrix (A'First (1) .. A'Last (1) - 1,
                    A'First (2) .. A'Last (2) - 1);
        X : Real;
    begin
        for P in A'Range (1) loop
            for Q in A'Range (2) loop
                X := A (P, Q);

                if P < I and then Q < J then
                    B (P, Q) := X;
                elsif P < I and then Q > J then
                    B (P, Q - 1) := X;
                elsif P > I and then Q < J then
                    B (P - 1, Q) := X;
                elsif P > I and then Q > J then
                    B (P - 1, Q - 1) := X;
                end if;
            end loop;
        end loop;

        return B;
    end Cofactor;

    function Det (A : Matrix) return Real is
        D : Real;
        K : Integer;
    begin
        if A'Length = 1 then
            D := A (A'First (1), A'First (2));
        else
            D := 0.0;
            K := 1;

            for J in A'Range (2) loop
                D := D + Real (K) * A (A'First (1), J) *
                            Det (Cofactor (A, A'First (1), J));
                K := -K;
            end loop;
        end if;

        return D;
    end Det;

    function Init (N : Positive) return Matrix is
        B : Matrix (1 .. N, 1 .. N);
    begin
        for I in B'Range (1) loop
            for J in B'Range (2) loop
                B (I, J) := 1.0 / Real (I + J - 1);
            end loop;
        end loop;

        return B;
    end Init;

    function Inverse (A : Matrix) return Matrix is
        B : Matrix (A'Range (1), A'Range (2));
        D : Real := Det (A);
        E : Real;
    begin
        if A'Length = 1 then
            return (1 .. 1 => (1 .. 1 => 1.0 / D));
        end if;

        for I in B'Range (1) loop
            for J in B'Range (2) loop
                B (I, J) := Real ((-1) ** (I + J)) *
                               (Det (Cofactor (A, I, J)) / D);
            end loop;
        end loop;

        -- Now check the inverse

        for I in A'Range loop
            for J in A'Range loop
                E := 0.0;

                for K in A'Range loop
                    E := E + A (I, K) * B (K, J);
                end loop;

                if (I = J and then E /= 1.0) or else
                   (I /= J and then E /= 0.0) then
                    raise Program_Error;
                end if;
            end loop;
        end loop;

        return B;
    end Inverse;


begin
    Put_Line ("*** TEST Inversion of Cauchy Matricies.");

    for N in 1 .. Trials loop
        begin
            declare
                A : constant Matrix := Init (N);
                B : constant Matrix := Inverse (A);
            begin
                Put_Line ("*** REMARK: The Cauchy Matrix of size" &
                          Integer'Image (N) & " successfully inverted.");
            end;
        exception
            when Program_Error =>
                Put_Line ("*** REMARK: The Cauchy Matrix of size" &
                          Integer'Image (N) & " not successfully inverted.");
            when Numeric_Error =>
                Put_Line ("*** REMARK: The Cauchy Matrix of size" &
                          Integer'Image (N) & " appears singular.");
            when others =>
                Put_Line ("*** REMARK: Unexpected exception raised.");
                raise;
        end;
    end loop;

    Put_Line ("*** FINISHED Matrix Inversion Test.");

end Cauchy;