527 lines
15 KiB
Ada
527 lines
15 KiB
Ada
------------------------------------------------------------------------------
|
|
-- --
|
|
-- GNAT RUN-TIME COMPONENTS --
|
|
-- --
|
|
-- S Y S T E M . G E N E R I C _ A R R A Y _ O P E R A T I O N S --
|
|
-- --
|
|
-- B o d y --
|
|
-- --
|
|
-- Copyright (C) 2006-2009, Free Software Foundation, Inc. --
|
|
-- --
|
|
-- GNAT is free software; you can redistribute it and/or modify it under --
|
|
-- terms of the GNU General Public License as published by the Free Soft- --
|
|
-- ware Foundation; either version 3, or (at your option) any later ver- --
|
|
-- sion. GNAT is distributed in the hope that it will be useful, but WITH- --
|
|
-- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY --
|
|
-- or FITNESS FOR A PARTICULAR PURPOSE. --
|
|
-- --
|
|
-- As a special exception under Section 7 of GPL version 3, you are granted --
|
|
-- additional permissions described in the GCC Runtime Library Exception, --
|
|
-- version 3.1, as published by the Free Software Foundation. --
|
|
-- --
|
|
-- You should have received a copy of the GNU General Public License and --
|
|
-- a copy of the GCC Runtime Library Exception along with this program; --
|
|
-- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see --
|
|
-- <http://www.gnu.org/licenses/>. --
|
|
-- --
|
|
-- GNAT was originally developed by the GNAT team at New York University. --
|
|
-- Extensive contributions were provided by Ada Core Technologies Inc. --
|
|
-- --
|
|
------------------------------------------------------------------------------
|
|
|
|
package body System.Generic_Array_Operations is
|
|
|
|
-- The local function Check_Unit_Last computes the index
|
|
-- of the last element returned by Unit_Vector or Unit_Matrix.
|
|
-- A separate function is needed to allow raising Constraint_Error
|
|
-- before declaring the function result variable. The result variable
|
|
-- needs to be declared first, to allow front-end inlining.
|
|
|
|
function Check_Unit_Last
|
|
(Index : Integer;
|
|
Order : Positive;
|
|
First : Integer) return Integer;
|
|
pragma Inline_Always (Check_Unit_Last);
|
|
|
|
function Square_Matrix_Length (A : Matrix) return Natural is
|
|
begin
|
|
if A'Length (1) /= A'Length (2) then
|
|
raise Constraint_Error with "matrix is not square";
|
|
end if;
|
|
|
|
return A'Length (1);
|
|
end Square_Matrix_Length;
|
|
|
|
---------------------
|
|
-- Check_Unit_Last --
|
|
---------------------
|
|
|
|
function Check_Unit_Last
|
|
(Index : Integer;
|
|
Order : Positive;
|
|
First : Integer) return Integer is
|
|
begin
|
|
-- Order the tests carefully to avoid overflow
|
|
|
|
if Index < First
|
|
or else First > Integer'Last - Order + 1
|
|
or else Index > First + (Order - 1)
|
|
then
|
|
raise Constraint_Error;
|
|
end if;
|
|
|
|
return First + (Order - 1);
|
|
end Check_Unit_Last;
|
|
|
|
-------------------
|
|
-- Inner_Product --
|
|
-------------------
|
|
|
|
function Inner_Product
|
|
(Left : Left_Vector;
|
|
Right : Right_Vector)
|
|
return Result_Scalar
|
|
is
|
|
R : Result_Scalar := Zero;
|
|
|
|
begin
|
|
if Left'Length /= Right'Length then
|
|
raise Constraint_Error with
|
|
"vectors are of different length in inner product";
|
|
end if;
|
|
|
|
for J in Left'Range loop
|
|
R := R + Left (J) * Right (J - Left'First + Right'First);
|
|
end loop;
|
|
|
|
return R;
|
|
end Inner_Product;
|
|
|
|
----------------------------------
|
|
-- Matrix_Elementwise_Operation --
|
|
----------------------------------
|
|
|
|
function Matrix_Elementwise_Operation (X : X_Matrix) return Result_Matrix is
|
|
R : Result_Matrix (X'Range (1), X'Range (2));
|
|
|
|
begin
|
|
for J in R'Range (1) loop
|
|
for K in R'Range (2) loop
|
|
R (J, K) := Operation (X (J, K));
|
|
end loop;
|
|
end loop;
|
|
|
|
return R;
|
|
end Matrix_Elementwise_Operation;
|
|
|
|
----------------------------------
|
|
-- Vector_Elementwise_Operation --
|
|
----------------------------------
|
|
|
|
function Vector_Elementwise_Operation (X : X_Vector) return Result_Vector is
|
|
R : Result_Vector (X'Range);
|
|
|
|
begin
|
|
for J in R'Range loop
|
|
R (J) := Operation (X (J));
|
|
end loop;
|
|
|
|
return R;
|
|
end Vector_Elementwise_Operation;
|
|
|
|
-----------------------------------------
|
|
-- Matrix_Matrix_Elementwise_Operation --
|
|
-----------------------------------------
|
|
|
|
function Matrix_Matrix_Elementwise_Operation
|
|
(Left : Left_Matrix;
|
|
Right : Right_Matrix)
|
|
return Result_Matrix
|
|
is
|
|
R : Result_Matrix (Left'Range (1), Left'Range (2));
|
|
begin
|
|
if Left'Length (1) /= Right'Length (1)
|
|
or else Left'Length (2) /= Right'Length (2)
|
|
then
|
|
raise Constraint_Error with
|
|
"matrices are of different dimension in elementwise operation";
|
|
end if;
|
|
|
|
for J in R'Range (1) loop
|
|
for K in R'Range (2) loop
|
|
R (J, K) :=
|
|
Operation
|
|
(Left (J, K),
|
|
Right
|
|
(J - R'First (1) + Right'First (1),
|
|
K - R'First (2) + Right'First (2)));
|
|
end loop;
|
|
end loop;
|
|
|
|
return R;
|
|
end Matrix_Matrix_Elementwise_Operation;
|
|
|
|
------------------------------------------------
|
|
-- Matrix_Matrix_Scalar_Elementwise_Operation --
|
|
------------------------------------------------
|
|
|
|
function Matrix_Matrix_Scalar_Elementwise_Operation
|
|
(X : X_Matrix;
|
|
Y : Y_Matrix;
|
|
Z : Z_Scalar) return Result_Matrix
|
|
is
|
|
R : Result_Matrix (X'Range (1), X'Range (2));
|
|
|
|
begin
|
|
if X'Length (1) /= Y'Length (1)
|
|
or else X'Length (2) /= Y'Length (2)
|
|
then
|
|
raise Constraint_Error with
|
|
"matrices are of different dimension in elementwise operation";
|
|
end if;
|
|
|
|
for J in R'Range (1) loop
|
|
for K in R'Range (2) loop
|
|
R (J, K) :=
|
|
Operation
|
|
(X (J, K),
|
|
Y (J - R'First (1) + Y'First (1),
|
|
K - R'First (2) + Y'First (2)),
|
|
Z);
|
|
end loop;
|
|
end loop;
|
|
|
|
return R;
|
|
end Matrix_Matrix_Scalar_Elementwise_Operation;
|
|
|
|
-----------------------------------------
|
|
-- Vector_Vector_Elementwise_Operation --
|
|
-----------------------------------------
|
|
|
|
function Vector_Vector_Elementwise_Operation
|
|
(Left : Left_Vector;
|
|
Right : Right_Vector) return Result_Vector
|
|
is
|
|
R : Result_Vector (Left'Range);
|
|
|
|
begin
|
|
if Left'Length /= Right'Length then
|
|
raise Constraint_Error with
|
|
"vectors are of different length in elementwise operation";
|
|
end if;
|
|
|
|
for J in R'Range loop
|
|
R (J) := Operation (Left (J), Right (J - R'First + Right'First));
|
|
end loop;
|
|
|
|
return R;
|
|
end Vector_Vector_Elementwise_Operation;
|
|
|
|
------------------------------------------------
|
|
-- Vector_Vector_Scalar_Elementwise_Operation --
|
|
------------------------------------------------
|
|
|
|
function Vector_Vector_Scalar_Elementwise_Operation
|
|
(X : X_Vector;
|
|
Y : Y_Vector;
|
|
Z : Z_Scalar) return Result_Vector
|
|
is
|
|
R : Result_Vector (X'Range);
|
|
|
|
begin
|
|
if X'Length /= Y'Length then
|
|
raise Constraint_Error with
|
|
"vectors are of different length in elementwise operation";
|
|
end if;
|
|
|
|
for J in R'Range loop
|
|
R (J) := Operation (X (J), Y (J - X'First + Y'First), Z);
|
|
end loop;
|
|
|
|
return R;
|
|
end Vector_Vector_Scalar_Elementwise_Operation;
|
|
|
|
-----------------------------------------
|
|
-- Matrix_Scalar_Elementwise_Operation --
|
|
-----------------------------------------
|
|
|
|
function Matrix_Scalar_Elementwise_Operation
|
|
(Left : Left_Matrix;
|
|
Right : Right_Scalar) return Result_Matrix
|
|
is
|
|
R : Result_Matrix (Left'Range (1), Left'Range (2));
|
|
|
|
begin
|
|
for J in R'Range (1) loop
|
|
for K in R'Range (2) loop
|
|
R (J, K) := Operation (Left (J, K), Right);
|
|
end loop;
|
|
end loop;
|
|
|
|
return R;
|
|
end Matrix_Scalar_Elementwise_Operation;
|
|
|
|
-----------------------------------------
|
|
-- Vector_Scalar_Elementwise_Operation --
|
|
-----------------------------------------
|
|
|
|
function Vector_Scalar_Elementwise_Operation
|
|
(Left : Left_Vector;
|
|
Right : Right_Scalar) return Result_Vector
|
|
is
|
|
R : Result_Vector (Left'Range);
|
|
|
|
begin
|
|
for J in R'Range loop
|
|
R (J) := Operation (Left (J), Right);
|
|
end loop;
|
|
|
|
return R;
|
|
end Vector_Scalar_Elementwise_Operation;
|
|
|
|
-----------------------------------------
|
|
-- Scalar_Matrix_Elementwise_Operation --
|
|
-----------------------------------------
|
|
|
|
function Scalar_Matrix_Elementwise_Operation
|
|
(Left : Left_Scalar;
|
|
Right : Right_Matrix) return Result_Matrix
|
|
is
|
|
R : Result_Matrix (Right'Range (1), Right'Range (2));
|
|
|
|
begin
|
|
for J in R'Range (1) loop
|
|
for K in R'Range (2) loop
|
|
R (J, K) := Operation (Left, Right (J, K));
|
|
end loop;
|
|
end loop;
|
|
|
|
return R;
|
|
end Scalar_Matrix_Elementwise_Operation;
|
|
|
|
-----------------------------------------
|
|
-- Scalar_Vector_Elementwise_Operation --
|
|
-----------------------------------------
|
|
|
|
function Scalar_Vector_Elementwise_Operation
|
|
(Left : Left_Scalar;
|
|
Right : Right_Vector) return Result_Vector
|
|
is
|
|
R : Result_Vector (Right'Range);
|
|
|
|
begin
|
|
for J in R'Range loop
|
|
R (J) := Operation (Left, Right (J));
|
|
end loop;
|
|
|
|
return R;
|
|
end Scalar_Vector_Elementwise_Operation;
|
|
|
|
---------------------------
|
|
-- Matrix_Matrix_Product --
|
|
---------------------------
|
|
|
|
function Matrix_Matrix_Product
|
|
(Left : Left_Matrix;
|
|
Right : Right_Matrix) return Result_Matrix
|
|
is
|
|
R : Result_Matrix (Left'Range (1), Right'Range (2));
|
|
|
|
begin
|
|
if Left'Length (2) /= Right'Length (1) then
|
|
raise Constraint_Error with
|
|
"incompatible dimensions in matrix multiplication";
|
|
end if;
|
|
|
|
for J in R'Range (1) loop
|
|
for K in R'Range (2) loop
|
|
declare
|
|
S : Result_Scalar := Zero;
|
|
begin
|
|
for M in Left'Range (2) loop
|
|
S := S + Left (J, M)
|
|
* Right (M - Left'First (2) + Right'First (1), K);
|
|
end loop;
|
|
|
|
R (J, K) := S;
|
|
end;
|
|
end loop;
|
|
end loop;
|
|
|
|
return R;
|
|
end Matrix_Matrix_Product;
|
|
|
|
---------------------------
|
|
-- Matrix_Vector_Product --
|
|
---------------------------
|
|
|
|
function Matrix_Vector_Product
|
|
(Left : Matrix;
|
|
Right : Right_Vector) return Result_Vector
|
|
is
|
|
R : Result_Vector (Left'Range (1));
|
|
|
|
begin
|
|
if Left'Length (2) /= Right'Length then
|
|
raise Constraint_Error with
|
|
"incompatible dimensions in matrix-vector multiplication";
|
|
end if;
|
|
|
|
for J in Left'Range (1) loop
|
|
declare
|
|
S : Result_Scalar := Zero;
|
|
begin
|
|
for K in Left'Range (2) loop
|
|
S := S + Left (J, K) * Right (K - Left'First (2) + Right'First);
|
|
end loop;
|
|
|
|
R (J) := S;
|
|
end;
|
|
end loop;
|
|
|
|
return R;
|
|
end Matrix_Vector_Product;
|
|
|
|
-------------------
|
|
-- Outer_Product --
|
|
-------------------
|
|
|
|
function Outer_Product
|
|
(Left : Left_Vector;
|
|
Right : Right_Vector) return Matrix
|
|
is
|
|
R : Matrix (Left'Range, Right'Range);
|
|
|
|
begin
|
|
for J in R'Range (1) loop
|
|
for K in R'Range (2) loop
|
|
R (J, K) := Left (J) * Right (K);
|
|
end loop;
|
|
end loop;
|
|
|
|
return R;
|
|
end Outer_Product;
|
|
|
|
---------------
|
|
-- Transpose --
|
|
---------------
|
|
|
|
procedure Transpose (A : Matrix; R : out Matrix) is
|
|
begin
|
|
for J in R'Range (1) loop
|
|
for K in R'Range (2) loop
|
|
R (J, K) := A (K - R'First (2) + A'First (1),
|
|
J - R'First (1) + A'First (2));
|
|
end loop;
|
|
end loop;
|
|
end Transpose;
|
|
|
|
-------------------------------
|
|
-- Update_Matrix_With_Matrix --
|
|
-------------------------------
|
|
|
|
procedure Update_Matrix_With_Matrix (X : in out X_Matrix; Y : Y_Matrix) is
|
|
begin
|
|
if X'Length (1) /= Y'Length (1)
|
|
or else X'Length (2) /= Y'Length (2)
|
|
then
|
|
raise Constraint_Error with
|
|
"matrices are of different dimension in update operation";
|
|
end if;
|
|
|
|
for J in X'Range (1) loop
|
|
for K in X'Range (2) loop
|
|
Update (X (J, K), Y (J - X'First (1) + Y'First (1),
|
|
K - X'First (2) + Y'First (2)));
|
|
end loop;
|
|
end loop;
|
|
end Update_Matrix_With_Matrix;
|
|
|
|
-------------------------------
|
|
-- Update_Vector_With_Vector --
|
|
-------------------------------
|
|
|
|
procedure Update_Vector_With_Vector (X : in out X_Vector; Y : Y_Vector) is
|
|
begin
|
|
if X'Length /= Y'Length then
|
|
raise Constraint_Error with
|
|
"vectors are of different length in update operation";
|
|
end if;
|
|
|
|
for J in X'Range loop
|
|
Update (X (J), Y (J - X'First + Y'First));
|
|
end loop;
|
|
end Update_Vector_With_Vector;
|
|
|
|
-----------------
|
|
-- Unit_Matrix --
|
|
-----------------
|
|
|
|
function Unit_Matrix
|
|
(Order : Positive;
|
|
First_1 : Integer := 1;
|
|
First_2 : Integer := 1) return Matrix
|
|
is
|
|
R : Matrix (First_1 .. Check_Unit_Last (First_1, Order, First_1),
|
|
First_2 .. Check_Unit_Last (First_2, Order, First_2));
|
|
|
|
begin
|
|
R := (others => (others => Zero));
|
|
|
|
for J in 0 .. Order - 1 loop
|
|
R (First_1 + J, First_2 + J) := One;
|
|
end loop;
|
|
|
|
return R;
|
|
end Unit_Matrix;
|
|
|
|
-----------------
|
|
-- Unit_Vector --
|
|
-----------------
|
|
|
|
function Unit_Vector
|
|
(Index : Integer;
|
|
Order : Positive;
|
|
First : Integer := 1) return Vector
|
|
is
|
|
R : Vector (First .. Check_Unit_Last (Index, Order, First));
|
|
begin
|
|
R := (others => Zero);
|
|
R (Index) := One;
|
|
return R;
|
|
end Unit_Vector;
|
|
|
|
---------------------------
|
|
-- Vector_Matrix_Product --
|
|
---------------------------
|
|
|
|
function Vector_Matrix_Product
|
|
(Left : Left_Vector;
|
|
Right : Matrix) return Result_Vector
|
|
is
|
|
R : Result_Vector (Right'Range (2));
|
|
|
|
begin
|
|
if Left'Length /= Right'Length (2) then
|
|
raise Constraint_Error with
|
|
"incompatible dimensions in vector-matrix multiplication";
|
|
end if;
|
|
|
|
for J in Right'Range (2) loop
|
|
declare
|
|
S : Result_Scalar := Zero;
|
|
|
|
begin
|
|
for K in Right'Range (1) loop
|
|
S := S + Left (J - Right'First (1) + Left'First) * Right (K, J);
|
|
end loop;
|
|
|
|
R (J) := S;
|
|
end;
|
|
end loop;
|
|
|
|
return R;
|
|
end Vector_Matrix_Product;
|
|
|
|
end System.Generic_Array_Operations;
|