diff options
Diffstat (limited to 'gcc-4.7/gcc/ada/s-gearop.adb')
| -rw-r--r-- | gcc-4.7/gcc/ada/s-gearop.adb | 948 |
1 files changed, 0 insertions, 948 deletions
diff --git a/gcc-4.7/gcc/ada/s-gearop.adb b/gcc-4.7/gcc/ada/s-gearop.adb deleted file mode 100644 index db18a7ebe..000000000 --- a/gcc-4.7/gcc/ada/s-gearop.adb +++ /dev/null @@ -1,948 +0,0 @@ ------------------------------------------------------------------------------- --- -- --- 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-2012, 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. -- --- -- ------------------------------------------------------------------------------- - -with Ada.Numerics; use Ada.Numerics; - -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); - - -------------- - -- Diagonal -- - -------------- - - function Diagonal (A : Matrix) return Vector is - N : constant Natural := Natural'Min (A'Length (1), A'Length (2)); - R : Vector (A'First (1) .. A'First (1) + N - 1); - - begin - for J in 0 .. N - 1 loop - R (R'First + J) := A (A'First (1) + J, A'First (2) + J); - end loop; - - return R; - end Diagonal; - - -------------------------- - -- Square_Matrix_Length -- - -------------------------- - - 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; - - --------------------- - -- Back_Substitute -- - --------------------- - - procedure Back_Substitute (M, N : in out Matrix) is - pragma Assert (M'First (1) = N'First (1) - and then - M'Last (1) = N'Last (1)); - - procedure Sub_Row - (M : in out Matrix; - Target : Integer; - Source : Integer; - Factor : Scalar); - -- Elementary row operation that subtracts Factor * M (Source, <>) from - -- M (Target, <>) - - procedure Sub_Row - (M : in out Matrix; - Target : Integer; - Source : Integer; - Factor : Scalar) - is - begin - for J in M'Range (2) loop - M (Target, J) := M (Target, J) - Factor * M (Source, J); - end loop; - end Sub_Row; - - -- Local declarations - - Max_Col : Integer := M'Last (2); - - -- Start of processing for Back_Substitute - - begin - Do_Rows : for Row in reverse M'Range (1) loop - Find_Non_Zero : for Col in reverse M'First (2) .. Max_Col loop - if Is_Non_Zero (M (Row, Col)) then - - -- Found first non-zero element, so subtract a multiple of this - -- element from all higher rows, to reduce all other elements - -- in this column to zero. - - declare - -- We can't use a for loop, as we'd need to iterate to - -- Row - 1, but that expression will overflow if M'First - -- equals Integer'First, which is true for aggregates - -- without explicit bounds.. - - J : Integer := M'First (1); - - begin - while J < Row loop - Sub_Row (N, J, Row, (M (J, Col) / M (Row, Col))); - Sub_Row (M, J, Row, (M (J, Col) / M (Row, Col))); - J := J + 1; - end loop; - end; - - -- Avoid potential overflow in the subtraction below - - exit Do_Rows when Col = M'First (2); - - Max_Col := Col - 1; - - exit Find_Non_Zero; - end if; - end loop Find_Non_Zero; - end loop Do_Rows; - end Back_Substitute; - - ----------------------- - -- Forward_Eliminate -- - ----------------------- - - procedure Forward_Eliminate - (M : in out Matrix; - N : in out Matrix; - Det : out Scalar) - is - pragma Assert (M'First (1) = N'First (1) - and then - M'Last (1) = N'Last (1)); - - -- The following are variations of the elementary matrix row operations: - -- row switching, row multiplication and row addition. Because in this - -- algorithm the addition factor is always a negated value, we chose to - -- use row subtraction instead. Similarly, instead of multiplying by - -- a reciprocal, we divide. - - procedure Sub_Row - (M : in out Matrix; - Target : Integer; - Source : Integer; - Factor : Scalar); - -- Subtrace Factor * M (Source, <>) from M (Target, <>) - - procedure Divide_Row - (M, N : in out Matrix; - Row : Integer; - Scale : Scalar); - -- Divide M (Row) and N (Row) by Scale, and update Det - - procedure Switch_Row - (M, N : in out Matrix; - Row_1 : Integer; - Row_2 : Integer); - -- Exchange M (Row_1) and N (Row_1) with M (Row_2) and N (Row_2), - -- negating Det in the process. - - ------------- - -- Sub_Row -- - ------------- - - procedure Sub_Row - (M : in out Matrix; - Target : Integer; - Source : Integer; - Factor : Scalar) - is - begin - for J in M'Range (2) loop - M (Target, J) := M (Target, J) - Factor * M (Source, J); - end loop; - end Sub_Row; - - ---------------- - -- Divide_Row -- - ---------------- - - procedure Divide_Row - (M, N : in out Matrix; - Row : Integer; - Scale : Scalar) - is - begin - Det := Det * Scale; - - for J in M'Range (2) loop - M (Row, J) := M (Row, J) / Scale; - end loop; - - for J in N'Range (2) loop - N (Row - M'First (1) + N'First (1), J) := - N (Row - M'First (1) + N'First (1), J) / Scale; - end loop; - end Divide_Row; - - ---------------- - -- Switch_Row -- - ---------------- - - procedure Switch_Row - (M, N : in out Matrix; - Row_1 : Integer; - Row_2 : Integer) - is - procedure Swap (X, Y : in out Scalar); - -- Exchange the values of X and Y - - procedure Swap (X, Y : in out Scalar) is - T : constant Scalar := X; - begin - X := Y; - Y := T; - end Swap; - - -- Start of processing for Switch_Row - - begin - if Row_1 /= Row_2 then - Det := Zero - Det; - - for J in M'Range (2) loop - Swap (M (Row_1, J), M (Row_2, J)); - end loop; - - for J in N'Range (2) loop - Swap (N (Row_1 - M'First (1) + N'First (1), J), - N (Row_2 - M'First (1) + N'First (1), J)); - end loop; - end if; - end Switch_Row; - - -- Local declarations - - Row : Integer := M'First (1); - - -- Start of processing for Forward_Eliminate - - begin - Det := One; - - for J in M'Range (2) loop - declare - Max_Row : Integer := Row; - Max_Abs : Real'Base := 0.0; - - begin - -- Find best pivot in column J, starting in row Row - - for K in Row .. M'Last (1) loop - declare - New_Abs : constant Real'Base := abs M (K, J); - begin - if Max_Abs < New_Abs then - Max_Abs := New_Abs; - Max_Row := K; - end if; - end; - end loop; - - if Max_Abs > 0.0 then - Switch_Row (M, N, Row, Max_Row); - - -- The temporaries below are necessary to force a copy of the - -- value and avoid improper aliasing. - - declare - Scale : constant Scalar := M (Row, J); - begin - Divide_Row (M, N, Row, Scale); - end; - - for U in Row + 1 .. M'Last (1) loop - declare - Factor : constant Scalar := M (U, J); - begin - Sub_Row (N, U, Row, Factor); - Sub_Row (M, U, Row, Factor); - end; - end loop; - - exit when Row >= M'Last (1); - - Row := Row + 1; - - else - -- Set zero (note that we do not have literals) - - Det := Zero; - end if; - end; - end loop; - end Forward_Eliminate; - - ------------------- - -- 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; - - ------------- - -- L2_Norm -- - ------------- - - function L2_Norm (X : X_Vector) return Result_Real'Base is - Sum : Result_Real'Base := 0.0; - - begin - for J in X'Range loop - Sum := Sum + Result_Real'Base (abs X (J))**2; - end loop; - - return Sqrt (Sum); - end L2_Norm; - - ---------------------------------- - -- 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; - - ---------- - -- Sqrt -- - ---------- - - function Sqrt (X : Real'Base) return Real'Base is - Root, Next : Real'Base; - - begin - -- Be defensive: any comparisons with NaN values will yield False. - - if not (X > 0.0) then - if X = 0.0 then - return X; - else - raise Argument_Error; - end if; - - elsif X > Real'Base'Last then - - -- X is infinity, which is its own square root - - return X; - end if; - - -- Compute an initial estimate based on: - - -- X = M * R**E and Sqrt (X) = Sqrt (M) * R**(E / 2.0), - - -- where M is the mantissa, R is the radix and E the exponent. - - -- By ignoring the mantissa and ignoring the case of an odd - -- exponent, we get a final error that is at most R. In other words, - -- the result has about a single bit precision. - - Root := Real'Base (Real'Machine_Radix) ** (Real'Exponent (X) / 2); - - -- Because of the poor initial estimate, use the Babylonian method of - -- computing the square root, as it is stable for all inputs. Every step - -- will roughly double the precision of the result. Just a few steps - -- suffice in most cases. Eight iterations should give about 2**8 bits - -- of precision. - - for J in 1 .. 8 loop - Next := (Root + X / Root) / 2.0; - exit when Root = Next; - Root := Next; - end loop; - - return Root; - end Sqrt; - - --------------------------- - -- 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_Solution -- - ---------------------------- - - function Matrix_Vector_Solution (A : Matrix; X : Vector) return Vector is - N : constant Natural := A'Length (1); - MA : Matrix := A; - MX : Matrix (A'Range (1), 1 .. 1); - R : Vector (A'Range (2)); - Det : Scalar; - - begin - if A'Length (2) /= N then - raise Constraint_Error with "matrix is not square"; - end if; - - if X'Length /= N then - raise Constraint_Error with "incompatible vector length"; - end if; - - for J in 0 .. MX'Length (1) - 1 loop - MX (MX'First (1) + J, 1) := X (X'First + J); - end loop; - - Forward_Eliminate (MA, MX, Det); - Back_Substitute (MA, MX); - - for J in 0 .. R'Length - 1 loop - R (R'First + J) := MX (MX'First (1) + J, 1); - end loop; - - return R; - end Matrix_Vector_Solution; - - ---------------------------- - -- Matrix_Matrix_Solution -- - ---------------------------- - - function Matrix_Matrix_Solution (A, X : Matrix) return Matrix is - N : constant Natural := A'Length (1); - MA : Matrix (A'Range (2), A'Range (2)); - MB : Matrix (A'Range (2), X'Range (2)); - Det : Scalar; - - begin - if A'Length (2) /= N then - raise Constraint_Error with "matrix is not square"; - end if; - - if X'Length (1) /= N then - raise Constraint_Error with "matrices have unequal number of rows"; - end if; - - for J in 0 .. A'Length (1) - 1 loop - for K in MA'Range (2) loop - MA (MA'First (1) + J, K) := A (A'First (1) + J, K); - end loop; - - for K in MB'Range (2) loop - MB (MB'First (1) + J, K) := X (X'First (1) + J, K); - end loop; - end loop; - - Forward_Eliminate (MA, MB, Det); - Back_Substitute (MA, MB); - - return MB; - end Matrix_Matrix_Solution; - - --------------------------- - -- 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; - - ----------------- - -- Swap_Column -- - ----------------- - - procedure Swap_Column (A : in out Matrix; Left, Right : Integer) is - Temp : Scalar; - begin - for J in A'Range (1) loop - Temp := A (J, Left); - A (J, Left) := A (J, Right); - A (J, Right) := Temp; - end loop; - end Swap_Column; - - --------------- - -- 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; |