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