Programmer's ROM - The Computer Language Library/programmersrom.iso/ada/math/matrix.ada

Text File | 1988-05-03 | 17.2 KB | 532 lines

------------ SIMTEL20 Ada Software Repository Proloque --------------
--                                   -*
-- Unit name    : MATRIX_PACKAGE
-- Version       : 1.0
-- Author    : Dr. Roger Lee
--        : Naval Air Developement Center
--        : Advanced Software Technology Division
--        : Warminster Pa. 18974
-- DDN Address    : LEE at NADC
-- Copyright    : (c) 1984 Roger Lee
-- Date created    : Feb 84
-- Release date    : 7 Jun 85
-- Last update     : Feb 84
-- Machine/System Compiled/Run on : VMS VAX11/780 Telesoft Ada compiler
--                                     -*
-----------------------------------------------------------------------
--                                     -*
-- Keywords    :
----------------: MATRIX,VECTOR
--
-- ABSTRACT    :
----------------: MATRIX_PACKAGE is a general purpose matrix package.
-- It defines data types VECTOR and MATRIX, and contains functions
-- to perform general matrix algebra operations.  It provides for addition,
-- subtraction, and multiplication of VECTORS, MATRICES and SCALARS.
-- It also provides for matrix inversion and vector dot product.
--
--                                      -*
---------------- Revision history --------------------------------------
--                                      -*
-- DATE        VERSION AUTHOR            HISTORY
-- 2/84          1.0    Roger Lee               Initial Release
--                                      -*
---------------- Distribution and Copyright ----------------------------
--                                      -*
-- This prologue must be included in all copies of this software.
--
-- This software is copyright by the author.
--
-- This software is released to the Ada community.
-- This software is released to the Public Domain (note:
--   software released to the Public Domain is not subject
--   to copyright protection).
-- Restrictions on use or distribution:  NONE
--                                      -*
--------------------- Disclaimer ---------------------------------------
--                                      -*
-- This software and its documentation are provided "AS IS" and
-- without any expressed or implied warranties whatsoever.
-- No warranties as to performance, merchantability, or fitness
-- for a particular purpose exist.
--
-- Because of the diverity of conditions and hardware under
-- which this software may be used, no warranty of fitness for
-- a particular purpose is offered.  The user is advised to
-- test the software thoroughly before relying on it.  The user
-- must assume the entire risk and liability of using this
-- software.
--
-- In no event shall any person or organization of people  be
-- held responsible for any direct, indirect, consequential
-- or inconsequential damages or lost profits.
--                                       -*
-------------------- END-PROLOGUE ---------------------------------------
--
--
--
package MATRIX_PACKAGE is
--> PROLOGUE
-->
-->**************************************************************************
--> This package is a general purpose matrix package. It defines data types
--> VECTOR and MATRIX, and contains functions to perform general matrix
--> algebra operations.
-->**************************************************************************
-->

    type VECTOR is array(INTEGER range<>) of FLOAT ;
    type MATRIX is array(INTEGER range<>,INTEGER range <>) of FLOAT;
    INCOMPARABLE_DIMENSION :exception; -- the dimension of matrices
                                           -- or vectors to be operated are
                       -- incomparable
    SINGULAR : exception;    -- matrix to be inverted is singular
    function TRANSPOSE(A : MATRIX) return MATRIX ; -- transpose of matrix
        function TRANSPOSE(A : VECTOR) return VECTOR ; -- transpose of vector
    function "+" (A : VECTOR; B : VECTOR) return VECTOR ; -- sum of vector
    function "+" (A : MATRIX; B : MATRIX) return MATRIX ; -- sum of matrix
    function "-" (A : VECTOR; B : VECTOR) return VECTOR ;
                           -- difference of vector
    function "-" (A : MATRIX; B : MATRIX) return MATRIX ;
                           -- difference of matrix
    function "*" (A : FLOAT;  B : VECTOR) return VECTOR ;
                           -- scalar, vector multiplication
    function "*" (A : VECTOR; B : FLOAT)  return VECTOR ;
                           -- vector, scalar multiplication
    function "*" (A : VECTOR; B : VECTOR) return FLOAT ;
                           -- inner(dot) product of two vectors
    function "*" (A : MATRIX; B : VECTOR) return VECTOR ;
                           -- matrix,column vector multiplication
    function "*" (A : VECTOR; B : MATRIX) return VECTOR ;
                           -- row vector,matrix multiplication
    function "*" (A : FLOAT;  B : MATRIX) return MATRIX ;
                           -- scalar, matrix multiplication
    function "*" (A : MATRIX; B : FLOAT)  return MATRIX ;
                           -- matrix, scalar multiplication
    function "*" (A : MATRIX; B : MATRIX) return MATRIX ;
                           -- matrix, matrix multiplication
    function "**"(A : MATRIX; P : INTEGER) return MATRIX;
                           -- square matrix raised to integer
                           -- power
--
--> CALLS
-->       NONE
-->
end MATRIX_PACKAGE;


---
package body MATRIX_PACKAGE is

---
    function TRANSPOSE(A : MATRIX) return MATRIX is
        B : MATRIX(A'FIRST(2)..A'LAST(2),A'FIRST(1)..A'LAST(1)) ;
--        ******************************************************************
--      This function performs the tranpose of input matrix A
--    ******************************************************************
    begin
        for I in A'RANGE(2) loop
            for J in A'RANGE(1) loop
                B(I,J) := A(J,I) ;
            end loop ;
        end loop ;
        return B ;
    end TRANSPOSE;
---
    function TRANSPOSE(A : VECTOR) return VECTOR is
--    *****************************************************************
--      This function returns the transpose of a vector. In programming
--      a vector is always stored as one-dimensional array. Therefore,
--        there is no difference between row vector and column vector.
--        Thus, this function just returns the input vector (do nothing).
--    *****************************************************************
        begin
            return A;
        end TRANSPOSE;
---
    function "+" (A : VECTOR; B : VECTOR) return VECTOR is
        C : VECTOR(A'FIRST..A'LAST) ;
--    **************************************************************
--      This function performs the addition of vector A and vector B
--       resulting in a vector. Comparability of dimensions is checked.
--    **************************************************************
    begin
        if A'FIRST /= B'FIRST or A'LAST /= B'LAST then
            raise INCOMPARABLE_DIMENSION;
        end if ;
        for I in A'RANGE loop
            C(I) := A(I)+B(I) ;
        end loop ;
        return C ;
    end "+";
---
    function "+" (A : MATRIX; B : MATRIX) return MATRIX is
        C : MATRIX(A'FIRST(1)..A'LAST(1),A'FIRST(2)..A'LAST(2)) ;
--    *******************************************************************
--      This function performs the addition of matrix A and matrix B
--      resulting in a matrix. Comparability of dimensions is checked.
--    *******************************************************************
    begin
        if (A'FIRST(1) /= B'FIRST(1) or A'LAST(1) /= B'LAST(1))
          or (A'FIRST(2) /= B'FIRST(2) or A'LAST(2) /= B'LAST(2)) then
            raise INCOMPARABLE_DIMENSION;
        end if ;
        for I in A'RANGE(1) loop
            for J in A'RANGE(2) loop
                C(I,J) := A(I,J)+B(I,J) ;
            end loop ;
        end loop ;
        return C ;
    end "+";
---
    function "-" (A : VECTOR; B : VECTOR) return VECTOR is
        C : VECTOR(A'FIRST..A'LAST) ;
--    ******************************************************************
--      This function performs the subtraction of vector B from vector A
--      resulting in a vector. Comparability of dimensions is checked.
--    ******************************************************************
    begin
        if A'FIRST /= B'FIRST or A'LAST /= B'LAST then
            raise INCOMPARABLE_DIMENSION;
        end if ;
        for I in A'RANGE loop
            C(I) := A(I)-B(I) ;
        end loop ;
        return C ;
    end "-";
---
    function "-" (A : MATRIX; B : MATRIX) return MATRIX is
        C : MATRIX(A'FIRST(1)..A'LAST(1),A'FIRST(2)..A'LAST(2)) ;
--    ******************************************************************
--      This function performs the subtraction of matrix B from matrix A
--      resulting in a matrix. Comparability of dimensions is checked.
--    ******************************************************************
    begin
        if (A'FIRST(1) /= B'FIRST(1) or A'LAST(1) /= B'LAST(1))
          or (A'FIRST(2) /= B'FIRST(2) or A'LAST(2) /= B'LAST(2)) then
            raise INCOMPARABLE_DIMENSION;
        end if ;
        for I in A'RANGE(1) loop
            for J in A'RANGE(2) loop
                C(I,J) := A(I,J)-B(I,J) ;
            end loop ;
        end loop ;
        return C ;
    end "-";
---
    function "*" (A:FLOAT; B:VECTOR) return VECTOR is
        C: VECTOR(B'FIRST..B'LAST);
--    ******************************************************************
--      This function performs the scalar multiplication of a floating
--      number A and a vector B resulting in a vector.
--    ******************************************************************
    begin
        for I in B'RANGE loop
            C(I):=A*B(I);
        end loop;
        return C ;
    end "*";

---
    function "*" (A:VECTOR; B:FLOAT) return VECTOR is
    begin
--    ********************************************************************
--      This function performs the scalar multiplication of a vector A and
--      a floating number B resulting in a vector.
--    ********************************************************************
        return B*A;
    end "*";

---
    function "*" (A : VECTOR; B : VECTOR) return FLOAT is
        S :FLOAT:=0.0;
--    *******************************************************************
--      This function performs the inner (dot) product of two vectors A
--      and B resulting in a floating number.
--      Comparability of dimensions is checked.
--    *******************************************************************
    begin
        if A'FIRST /= B'FIRST or A'LAST /= B'LAST then
            raise INCOMPARABLE_DIMENSION;
        end if ;
        for I in A'RANGE loop
            S := S+A(I)*B(I) ;
        end loop ;
        return S ;
    end "*";
---
    function "*" (A:MATRIX; B:VECTOR) return VECTOR is
        C:VECTOR(A'FIRST(1)..A'LAST(1));
        SUM:FLOAT;
--    **********************************************************************
--      This function performs the multiplication of a matrix A and a column
--      vector B resulting in a column vector.
--      Comparability of dimensions is checked.
--    **********************************************************************
    begin
        if A'FIRST(2)/=B'FIRST or A'LAST(2) /= B'LAST then
            raise INCOMPARABLE_DIMENSION;
        end if ;
        for I in A'RANGE(1) loop
            SUM := 0.0 ;
            for K in A'RANGE(2) loop
                SUM := SUM+A(I,K)*B(K);
            end loop;
            C(I):=SUM;
        end loop ;
        return C ;
    end "*";
---
    function "*" (A:VECTOR; B:MATRIX) return VECTOR is
        C:VECTOR(B'FIRST(2)..B'LAST(2));
        SUM:FLOAT;
--    ********************************************************************
--      This function performs the multiplication of a row vector A and a
--      matrix B resulting in a row vector.
--      Comparability of dimensions is checked.
--    ********************************************************************
    begin
        if A'FIRST/=B'FIRST(1) or A'LAST/=B'LAST(1) then
            raise INCOMPARABLE_DIMENSION;
        end if ;
        for J in B'RANGE(2) loop
            SUM := 0.0 ;
            for K in A'RANGE loop
                SUM := SUM+A(K)*B(K,J);
            end loop;
            C(J):=SUM;
        end loop ;
        return C ;
    end "*";
---
    function "*" (A:FLOAT; B:MATRIX) return MATRIX is
        C:MATRIX(B'FIRST(1)..B'LAST(1),B'FIRST(2)..B'LAST(2));
--    ********************************************************************
--      This function performs the scalar multipliction of a matrix B by
--      a floating number A resulting in a matrix.
--    ********************************************************************
    begin
        for I in B'RANGE(1) loop
            for J in B'RANGE(2) loop
                C(I,J) := A*B(I,J);
            end loop ;
        end loop ;
        return C ;
    end "*";
---
    function "*" (A:MATRIX; B:FLOAT) return MATRIX is
        C:MATRIX(A'FIRST(1)..A'LAST(1),A'FIRST(2)..A'LAST(2));
--    *****************************************************************
--      This function performs the scalar multipliction of a matrix A
--      by a floating number B resulting in a matrix.
--    *****************************************************************
    begin
        return B*A ;
    end "*";
---
    function "*" (A:MATRIX; B:MATRIX) return MATRIX is
        C:MATRIX(A'FIRST(1)..A'LAST(1),B'FIRST(2)..B'LAST(2));
        SUM: FLOAT;
--    ********************************************************************
--      This function performs the multiplication of matrix A and matrix B
--      resulting in a matrix. Comparability of dimensions is checked.
--    ********************************************************************
    begin
        if A'FIRST(2)/=B'FIRST(1) or A'LAST(2)/=B'LAST(1) then
            raise INCOMPARABLE_DIMENSION;
        end if ;
        for I in A'RANGE(1) loop
            for J in B'RANGE(2) loop
                SUM := 0.0 ;
                for K in A'RANGE(2) loop
                    SUM := SUM+A(I,K)*B(K,J);
                end loop;
                C(I,J) := SUM;
            end loop ;
        end loop ;
        return C ;
    end "*";
---
    function "**" (A : MATRIX; P : INTEGER)    return MATRIX is
        B,C : MATRIX(A'FIRST(1)..A'LAST(1), A'FIRST(1)..A'LAST(1));
        I_PIVOT,J_PIVOT : INTEGER range A'FIRST(1)..A'LAST(1);
        BIG_ENTRY, TEMP, EPSILON : FLOAT ;
        L, M : array(A'RANGE(1)) of INTEGER ;
--    *******************************************************************
--      This function performs the square matrix operation of " matrix A
--      raise to integer power P ".  When  P is negative , say P = -N ,
--      A**(-N) = (inverse(A))**N , that is, the inverse of A raise to
--      power N .  In this case, matrix A must be non-singular.
--      Exceptions will be raised if the matrix A is not a square matrix,
--      or if matrix A is singular.
--    *******************************************************************
    begin
        if A'FIRST(1)/=A'FIRST(2) or A'LAST(1)/=A'LAST(2) then
                    -- if not a square matrix
            raise INCOMPARABLE_DIMENSION ;
        end if;

        if P=0 then
            --& B = identity matrix

                  for I in A'RANGE(1) loop
            for J in A'RANGE(1) loop
                if I /= J then
                    B(I,J) := 0.0;
                else
                B(I,J) := 1.0;
                end if;
            end loop;
            end loop;
            return B;
        end if ;

        B := A ;

        if P>0 then
                --& B = A multiplied itself for P times

            for I in 1..P-1 loop
                B := B*A ;
            end loop ;
            return B ;
        end if;

    -- P is negative, find inverse first

        -- initiate the row and column interchange information

        for K in B'RANGE(1) loop
            L(K) := K ;    -- row interchage information
            M(K) := K ;    -- column interchange information
        end loop;

            -- major loop for inverse

        for K in B'RANGE(1) loop

            -- & search for row and column index I_PIVOT, J_PIVOT
            -- & both in (K .. B'LAST(1) ) for maximum B(I,J)
            -- & in absolute value :BIG_ENTRY

            BIG_ENTRY := 0.0 ;
                   --
             -- check matrix singularity
                 --
            for I in K..B'LAST(1) loop
                for J in K..B'LAST(1) loop
                    if abs(B(I,J)) > abs(BIG_ENTRY) then
                        BIG_ENTRY := B(I,J) ;
                        I_PIVOT := I ;
                        J_PIVOT := J ;
                    end if;
                end loop;
            end loop;
            if K = B'FIRST(1) then
                if BIG_ENTRY = 0.0 then
                raise SINGULAR;
                else
                EPSILON := FLOAT(A'LENGTH(1))*abs(BIG_ENTRY)
                                 *0.000001;
                end if;
            else
                if abs(BIG_ENTRY) < EPSILON then
                raise SINGULAR ;
                end if;
            end if;

            -- interchange row and column

                --& interchange K-th and I_PIVOT-th rows
            if I_PIVOT/=K then
                for J in B'RANGE(1) loop
                    TEMP := B(I_PIVOT,J);
                    B(I_PIVOT,J) := B(K,J) ;
                    B(K,J) := TEMP ;
                end loop;
                L(K) := I_PIVOT ;
            end if;
                --& interchange K-th and J_PIVOT-th columns
            if J_PIVOT/=K then
                for I in B'RANGE(1) loop
                    TEMP := B(I,J_PIVOT) ;
                    B(I,J_PIVOT) := B(I,K) ;
                    B(I,K) := TEMP ;
                end loop ;
                M(K) := J_PIVOT ;
            end if ;

            --& divide K-th column by minus pivot (-BIG_ENTRY)

            for I in B'RANGE(1) loop
                if I/=K then
                    B(I,K) := B(I,K)/(-BIG_ENTRY) ;
                end if;
            end loop ;

            -- reduce matrix row by row

            for I in B'RANGE(1) loop
                if I/=K then
                    for J in B'RANGE(1) loop
                         if J/=K then
                          B(I,J):=B(I,J)+B(I,K)*B(K,J);
                         end if ;
                    end loop ;
                end if ;
            end loop ;

            --& divide K-th row by pivot

            for J in B'RANGE(1) loop
                if J/=K then
                    B(K,J) := B(K,J)/BIG_ENTRY ;
                end if ;
            end loop ;
            B(K,K) := 1.0/BIG_ENTRY ;

        end loop ;    -- end of major inverse loop

        -- final column and row interchange to obtain
        -- inverse of A, i.e. A**(-1)

        for K in reverse B'RANGE(1) loop
                -- column interchage
            J_PIVOT := L(K) ;
            if J_PIVOT/=K then
                    --& intechange B(I,J_PIVOT) and B(I,K) for each row I
                for I in B'RANGE(1) loop
                    TEMP := B(I,J_PIVOT) ;
                    B(I,J_PIVOT) := B(I,K) ;
                    B(I,K) := TEMP ;
                end loop ;
            end if ;
                -- row interchage
            I_PIVOT := M(K) ;
            if I_PIVOT/=K then
                --& INTECHANGE B(I_PIVOT,J) and B(K,J) for each column J
                for J in B'RANGE(1) loop
                    TEMP := B(I_PIVOT,J) ;
                    B(I_PIVOT,J) := B(K,J) ;
                    B(K,J) := TEMP ;
                end loop ;
            end if ;
        end loop ;

        -- inverse of A is obtained and stored in B
        -- now ready to handle the negative power

                -- & C = B**(-P)
        if P=-1 then
            return B ;
        end if ;

        C := B ;
        for I in P+1..-1 loop
            C:= C*B ;
        end loop ;
        return C;
    end "**" ;
---
end MATRIX_PACKAGE;