From 554fd8c5195424bdbcabf5de30fdc183aba391bd Mon Sep 17 00:00:00 2001 From: upstream source tree Date: Sun, 15 Mar 2015 20:14:05 -0400 Subject: obtained gcc-4.6.4.tar.bz2 from upstream website; verified gcc-4.6.4.tar.bz2.sig; imported gcc-4.6.4 source tree from verified upstream tarball. downloading a git-generated archive based on the 'upstream' tag should provide you with a source tree that is binary identical to the one extracted from the above tarball. if you have obtained the source via the command 'git clone', however, do note that line-endings of files in your working directory might differ from line-endings of the respective files in the upstream repository. --- gcc/ada/a-ngrear.adb | 784 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 784 insertions(+) create mode 100644 gcc/ada/a-ngrear.adb (limited to 'gcc/ada/a-ngrear.adb') diff --git a/gcc/ada/a-ngrear.adb b/gcc/ada/a-ngrear.adb new file mode 100644 index 000000000..5c8a00924 --- /dev/null +++ b/gcc/ada/a-ngrear.adb @@ -0,0 +1,784 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- ADA.NUMERICS.GENERIC_REAL_ARRAYS -- +-- -- +-- 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 -- +-- . -- +-- -- +-- GNAT was originally developed by the GNAT team at New York University. -- +-- Extensive contributions were provided by Ada Core Technologies Inc. -- +-- -- +------------------------------------------------------------------------------ + +with System; use System; +with System.Generic_Real_BLAS; +with System.Generic_Real_LAPACK; +with System.Generic_Array_Operations; use System.Generic_Array_Operations; + +package body Ada.Numerics.Generic_Real_Arrays is + + -- Operations involving inner products use BLAS library implementations. + -- This allows larger matrices and vectors to be computed efficiently, + -- taking into account memory hierarchy issues and vector instructions + -- that vary widely between machines. + + -- Operations that are defined in terms of operations on the type Real, + -- such as addition, subtraction and scaling, are computed in the canonical + -- way looping over all elements. + + -- Operations for solving linear systems and computing determinant, + -- eigenvalues, eigensystem and inverse, are implemented using the + -- LAPACK library. + + package BLAS is + new Generic_Real_BLAS (Real'Base, Real_Vector, Real_Matrix); + + package LAPACK is + new Generic_Real_LAPACK (Real'Base, Real_Vector, Real_Matrix); + + use BLAS, LAPACK; + + -- Procedure versions of functions returning unconstrained values. + -- This allows for inlining the function wrapper. + + procedure Eigenvalues (A : Real_Matrix; Values : out Real_Vector); + procedure Inverse (A : Real_Matrix; R : out Real_Matrix); + procedure Solve (A : Real_Matrix; X : Real_Vector; B : out Real_Vector); + procedure Solve (A : Real_Matrix; X : Real_Matrix; B : out Real_Matrix); + + procedure Transpose is new + Generic_Array_Operations.Transpose + (Scalar => Real'Base, + Matrix => Real_Matrix); + + -- Helper function that raises a Constraint_Error is the argument is + -- not a square matrix, and otherwise returns its length. + + function Length is new Square_Matrix_Length (Real'Base, Real_Matrix); + + -- Instantiating the following subprograms directly would lead to + -- name clashes, so use a local package. + + package Instantiations is + + function "+" is new + Vector_Elementwise_Operation + (X_Scalar => Real'Base, + Result_Scalar => Real'Base, + X_Vector => Real_Vector, + Result_Vector => Real_Vector, + Operation => "+"); + + function "+" is new + Matrix_Elementwise_Operation + (X_Scalar => Real'Base, + Result_Scalar => Real'Base, + X_Matrix => Real_Matrix, + Result_Matrix => Real_Matrix, + Operation => "+"); + + function "+" is new + Vector_Vector_Elementwise_Operation + (Left_Scalar => Real'Base, + Right_Scalar => Real'Base, + Result_Scalar => Real'Base, + Left_Vector => Real_Vector, + Right_Vector => Real_Vector, + Result_Vector => Real_Vector, + Operation => "+"); + + function "+" is new + Matrix_Matrix_Elementwise_Operation + (Left_Scalar => Real'Base, + Right_Scalar => Real'Base, + Result_Scalar => Real'Base, + Left_Matrix => Real_Matrix, + Right_Matrix => Real_Matrix, + Result_Matrix => Real_Matrix, + Operation => "+"); + + function "-" is new + Vector_Elementwise_Operation + (X_Scalar => Real'Base, + Result_Scalar => Real'Base, + X_Vector => Real_Vector, + Result_Vector => Real_Vector, + Operation => "-"); + + function "-" is new + Matrix_Elementwise_Operation + (X_Scalar => Real'Base, + Result_Scalar => Real'Base, + X_Matrix => Real_Matrix, + Result_Matrix => Real_Matrix, + Operation => "-"); + + function "-" is new + Vector_Vector_Elementwise_Operation + (Left_Scalar => Real'Base, + Right_Scalar => Real'Base, + Result_Scalar => Real'Base, + Left_Vector => Real_Vector, + Right_Vector => Real_Vector, + Result_Vector => Real_Vector, + Operation => "-"); + + function "-" is new + Matrix_Matrix_Elementwise_Operation + (Left_Scalar => Real'Base, + Right_Scalar => Real'Base, + Result_Scalar => Real'Base, + Left_Matrix => Real_Matrix, + Right_Matrix => Real_Matrix, + Result_Matrix => Real_Matrix, + Operation => "-"); + + function "*" is new + Scalar_Vector_Elementwise_Operation + (Left_Scalar => Real'Base, + Right_Scalar => Real'Base, + Result_Scalar => Real'Base, + Right_Vector => Real_Vector, + Result_Vector => Real_Vector, + Operation => "*"); + + function "*" is new + Scalar_Matrix_Elementwise_Operation + (Left_Scalar => Real'Base, + Right_Scalar => Real'Base, + Result_Scalar => Real'Base, + Right_Matrix => Real_Matrix, + Result_Matrix => Real_Matrix, + Operation => "*"); + + function "*" is new + Vector_Scalar_Elementwise_Operation + (Left_Scalar => Real'Base, + Right_Scalar => Real'Base, + Result_Scalar => Real'Base, + Left_Vector => Real_Vector, + Result_Vector => Real_Vector, + Operation => "*"); + + function "*" is new + Matrix_Scalar_Elementwise_Operation + (Left_Scalar => Real'Base, + Right_Scalar => Real'Base, + Result_Scalar => Real'Base, + Left_Matrix => Real_Matrix, + Result_Matrix => Real_Matrix, + Operation => "*"); + + function "*" is new + Outer_Product + (Left_Scalar => Real'Base, + Right_Scalar => Real'Base, + Result_Scalar => Real'Base, + Left_Vector => Real_Vector, + Right_Vector => Real_Vector, + Matrix => Real_Matrix); + + function "/" is new + Vector_Scalar_Elementwise_Operation + (Left_Scalar => Real'Base, + Right_Scalar => Real'Base, + Result_Scalar => Real'Base, + Left_Vector => Real_Vector, + Result_Vector => Real_Vector, + Operation => "/"); + + function "/" is new + Matrix_Scalar_Elementwise_Operation + (Left_Scalar => Real'Base, + Right_Scalar => Real'Base, + Result_Scalar => Real'Base, + Left_Matrix => Real_Matrix, + Result_Matrix => Real_Matrix, + Operation => "/"); + + function "abs" is new + Vector_Elementwise_Operation + (X_Scalar => Real'Base, + Result_Scalar => Real'Base, + X_Vector => Real_Vector, + Result_Vector => Real_Vector, + Operation => "abs"); + + function "abs" is new + Matrix_Elementwise_Operation + (X_Scalar => Real'Base, + Result_Scalar => Real'Base, + X_Matrix => Real_Matrix, + Result_Matrix => Real_Matrix, + Operation => "abs"); + + function Unit_Matrix is new + Generic_Array_Operations.Unit_Matrix + (Scalar => Real'Base, + Matrix => Real_Matrix, + Zero => 0.0, + One => 1.0); + + function Unit_Vector is new + Generic_Array_Operations.Unit_Vector + (Scalar => Real'Base, + Vector => Real_Vector, + Zero => 0.0, + One => 1.0); + + end Instantiations; + + --------- + -- "+" -- + --------- + + function "+" (Right : Real_Vector) return Real_Vector + renames Instantiations."+"; + + function "+" (Right : Real_Matrix) return Real_Matrix + renames Instantiations."+"; + + function "+" (Left, Right : Real_Vector) return Real_Vector + renames Instantiations."+"; + + function "+" (Left, Right : Real_Matrix) return Real_Matrix + renames Instantiations."+"; + + --------- + -- "-" -- + --------- + + function "-" (Right : Real_Vector) return Real_Vector + renames Instantiations."-"; + + function "-" (Right : Real_Matrix) return Real_Matrix + renames Instantiations."-"; + + function "-" (Left, Right : Real_Vector) return Real_Vector + renames Instantiations."-"; + + function "-" (Left, Right : Real_Matrix) return Real_Matrix + renames Instantiations."-"; + + --------- + -- "*" -- + --------- + + -- Scalar multiplication + + function "*" (Left : Real'Base; Right : Real_Vector) return Real_Vector + renames Instantiations."*"; + + function "*" (Left : Real_Vector; Right : Real'Base) return Real_Vector + renames Instantiations."*"; + + function "*" (Left : Real'Base; Right : Real_Matrix) return Real_Matrix + renames Instantiations."*"; + + function "*" (Left : Real_Matrix; Right : Real'Base) return Real_Matrix + renames Instantiations."*"; + + -- Vector multiplication + + function "*" (Left, Right : Real_Vector) return Real'Base is + begin + if Left'Length /= Right'Length then + raise Constraint_Error with + "vectors are of different length in inner product"; + end if; + + return dot (Left'Length, X => Left, Y => Right); + end "*"; + + function "*" (Left, Right : Real_Vector) return Real_Matrix + renames Instantiations."*"; + + function "*" + (Left : Real_Vector; + Right : Real_Matrix) return Real_Vector + is + R : Real_Vector (Right'Range (2)); + + begin + if Left'Length /= Right'Length (1) then + raise Constraint_Error with + "incompatible dimensions in vector-matrix multiplication"; + end if; + + gemv (Trans => No_Trans'Access, + M => Right'Length (2), + N => Right'Length (1), + A => Right, + Ld_A => Right'Length (2), + X => Left, + Y => R); + + return R; + end "*"; + + function "*" + (Left : Real_Matrix; + Right : Real_Vector) return Real_Vector + is + R : Real_Vector (Left'Range (1)); + + begin + if Left'Length (2) /= Right'Length then + raise Constraint_Error with + "incompatible dimensions in matrix-vector multiplication"; + end if; + + gemv (Trans => Trans'Access, + M => Left'Length (2), + N => Left'Length (1), + A => Left, + Ld_A => Left'Length (2), + X => Right, + Y => R); + + return R; + end "*"; + + -- Matrix Multiplication + + function "*" (Left, Right : Real_Matrix) return Real_Matrix is + R : Real_Matrix (Left'Range (1), Right'Range (2)); + + begin + if Left'Length (2) /= Right'Length (1) then + raise Constraint_Error with + "incompatible dimensions in matrix-matrix multiplication"; + end if; + + gemm (Trans_A => No_Trans'Access, + Trans_B => No_Trans'Access, + M => Right'Length (2), + N => Left'Length (1), + K => Right'Length (1), + A => Right, + Ld_A => Right'Length (2), + B => Left, + Ld_B => Left'Length (2), + C => R, + Ld_C => R'Length (2)); + + return R; + end "*"; + + --------- + -- "/" -- + --------- + + function "/" (Left : Real_Vector; Right : Real'Base) return Real_Vector + renames Instantiations."/"; + + function "/" (Left : Real_Matrix; Right : Real'Base) return Real_Matrix + renames Instantiations."/"; + + ----------- + -- "abs" -- + ----------- + + function "abs" (Right : Real_Vector) return Real'Base is + begin + return nrm2 (Right'Length, Right); + end "abs"; + + function "abs" (Right : Real_Vector) return Real_Vector + renames Instantiations."abs"; + + function "abs" (Right : Real_Matrix) return Real_Matrix + renames Instantiations."abs"; + + ----------------- + -- Determinant -- + ----------------- + + function Determinant (A : Real_Matrix) return Real'Base is + N : constant Integer := Length (A); + LU : Real_Matrix (1 .. N, 1 .. N) := A; + Piv : Integer_Vector (1 .. N); + Info : aliased Integer := -1; + Det : Real := 1.0; + + begin + getrf (M => N, + N => N, + A => LU, + Ld_A => N, + I_Piv => Piv, + Info => Info'Access); + + if Info /= 0 then + raise Constraint_Error with "ill-conditioned matrix"; + end if; + + for J in 1 .. N loop + Det := (if Piv (J) /= J then -Det * LU (J, J) else Det * LU (J, J)); + end loop; + + return Det; + end Determinant; + + ----------------- + -- Eigensystem -- + ----------------- + + procedure Eigensystem + (A : Real_Matrix; + Values : out Real_Vector; + Vectors : out Real_Matrix) + is + N : constant Natural := Length (A); + Tau : Real_Vector (1 .. N); + L_Work : Real_Vector (1 .. 1); + Info : aliased Integer; + + E : Real_Vector (1 .. N); + pragma Warnings (Off, E); + + begin + if Values'Length /= N then + raise Constraint_Error with "wrong length for output vector"; + end if; + + if N = 0 then + return; + end if; + + -- Initialize working matrix and check for symmetric input matrix + + Transpose (A, Vectors); + + if A /= Vectors then + raise Argument_Error with "matrix not symmetric"; + end if; + + -- Compute size of additional working space + + sytrd (Uplo => Lower'Access, + N => N, + A => Vectors, + Ld_A => N, + D => Values, + E => E, + Tau => Tau, + Work => L_Work, + L_Work => -1, + Info => Info'Access); + + declare + Work : Real_Vector (1 .. Integer'Max (Integer (L_Work (1)), 2 * N)); + pragma Warnings (Off, Work); + + Comp_Z : aliased constant Character := 'V'; + + begin + -- Reduce matrix to tridiagonal form + + sytrd (Uplo => Lower'Access, + N => N, + A => Vectors, + Ld_A => A'Length (1), + D => Values, + E => E, + Tau => Tau, + Work => Work, + L_Work => Work'Length, + Info => Info'Access); + + if Info /= 0 then + raise Program_Error; + end if; + + -- Generate the real orthogonal matrix determined by sytrd + + orgtr (Uplo => Lower'Access, + N => N, + A => Vectors, + Ld_A => N, + Tau => Tau, + Work => Work, + L_Work => Work'Length, + Info => Info'Access); + + if Info /= 0 then + raise Program_Error; + end if; + + -- Compute all eigenvalues and eigenvectors using QR algorithm + + steqr (Comp_Z => Comp_Z'Access, + N => N, + D => Values, + E => E, + Z => Vectors, + Ld_Z => N, + Work => Work, + Info => Info'Access); + + if Info /= 0 then + raise Constraint_Error with + "eigensystem computation failed to converge"; + end if; + end; + end Eigensystem; + + ----------------- + -- Eigenvalues -- + ----------------- + + procedure Eigenvalues + (A : Real_Matrix; + Values : out Real_Vector) + is + N : constant Natural := Length (A); + L_Work : Real_Vector (1 .. 1); + Info : aliased Integer; + + B : Real_Matrix (1 .. N, 1 .. N); + Tau : Real_Vector (1 .. N); + E : Real_Vector (1 .. N); + pragma Warnings (Off, B); + pragma Warnings (Off, Tau); + pragma Warnings (Off, E); + + begin + if Values'Length /= N then + raise Constraint_Error with "wrong length for output vector"; + end if; + + if N = 0 then + return; + end if; + + -- Initialize working matrix and check for symmetric input matrix + + Transpose (A, B); + + if A /= B then + raise Argument_Error with "matrix not symmetric"; + end if; + + -- Find size of work area + + sytrd (Uplo => Lower'Access, + N => N, + A => B, + Ld_A => N, + D => Values, + E => E, + Tau => Tau, + Work => L_Work, + L_Work => -1, + Info => Info'Access); + + declare + Work : Real_Vector (1 .. Integer'Min (Integer (L_Work (1)), 4 * N)); + pragma Warnings (Off, Work); + + begin + -- Reduce matrix to tridiagonal form + + sytrd (Uplo => Lower'Access, + N => N, + A => B, + Ld_A => A'Length (1), + D => Values, + E => E, + Tau => Tau, + Work => Work, + L_Work => Work'Length, + Info => Info'Access); + + if Info /= 0 then + raise Constraint_Error; + end if; + + -- Compute all eigenvalues using QR algorithm + + sterf (N => N, + D => Values, + E => E, + Info => Info'Access); + + if Info /= 0 then + raise Constraint_Error with + "eigenvalues computation failed to converge"; + end if; + end; + end Eigenvalues; + + function Eigenvalues (A : Real_Matrix) return Real_Vector is + R : Real_Vector (A'Range (1)); + begin + Eigenvalues (A, R); + return R; + end Eigenvalues; + + ------------- + -- Inverse -- + ------------- + + procedure Inverse (A : Real_Matrix; R : out Real_Matrix) is + N : constant Integer := Length (A); + Piv : Integer_Vector (1 .. N); + L_Work : Real_Vector (1 .. 1); + Info : aliased Integer := -1; + + begin + -- All computations are done using column-major order, but this works + -- out fine, because Transpose (Inverse (Transpose (A))) = Inverse (A). + + R := A; + + -- Compute LU decomposition + + getrf (M => N, + N => N, + A => R, + Ld_A => N, + I_Piv => Piv, + Info => Info'Access); + + if Info /= 0 then + raise Constraint_Error with "inverting singular matrix"; + end if; + + -- Determine size of work area + + getri (N => N, + A => R, + Ld_A => N, + I_Piv => Piv, + Work => L_Work, + L_Work => -1, + Info => Info'Access); + + if Info /= 0 then + raise Constraint_Error; + end if; + + declare + Work : Real_Vector (1 .. Integer (L_Work (1))); + pragma Warnings (Off, Work); + + begin + -- Compute inverse from LU decomposition + + getri (N => N, + A => R, + Ld_A => N, + I_Piv => Piv, + Work => Work, + L_Work => Work'Length, + Info => Info'Access); + + if Info /= 0 then + raise Constraint_Error with "inverting singular matrix"; + end if; + + -- ??? Should iterate with gerfs, based on implementation advice + end; + end Inverse; + + function Inverse (A : Real_Matrix) return Real_Matrix is + R : Real_Matrix (A'Range (2), A'Range (1)); + begin + Inverse (A, R); + return R; + end Inverse; + + ----------- + -- Solve -- + ----------- + + procedure Solve (A : Real_Matrix; X : Real_Vector; B : out Real_Vector) is + begin + if Length (A) /= X'Length then + raise Constraint_Error with + "incompatible matrix and vector dimensions"; + end if; + + -- ??? Should solve directly, is faster and more accurate + + B := Inverse (A) * X; + end Solve; + + procedure Solve (A : Real_Matrix; X : Real_Matrix; B : out Real_Matrix) is + begin + if Length (A) /= X'Length (1) then + raise Constraint_Error with "incompatible matrix dimensions"; + end if; + + -- ??? Should solve directly, is faster and more accurate + + B := Inverse (A) * X; + end Solve; + + function Solve (A : Real_Matrix; X : Real_Vector) return Real_Vector is + B : Real_Vector (A'Range (2)); + begin + Solve (A, X, B); + return B; + end Solve; + + function Solve (A, X : Real_Matrix) return Real_Matrix is + B : Real_Matrix (A'Range (2), X'Range (2)); + begin + Solve (A, X, B); + return B; + end Solve; + + --------------- + -- Transpose -- + --------------- + + function Transpose (X : Real_Matrix) return Real_Matrix is + R : Real_Matrix (X'Range (2), X'Range (1)); + begin + Transpose (X, R); + + return R; + end Transpose; + + ----------------- + -- Unit_Matrix -- + ----------------- + + function Unit_Matrix + (Order : Positive; + First_1 : Integer := 1; + First_2 : Integer := 1) return Real_Matrix + renames Instantiations.Unit_Matrix; + + ----------------- + -- Unit_Vector -- + ----------------- + + function Unit_Vector + (Index : Integer; + Order : Positive; + First : Integer := 1) return Real_Vector + renames Instantiations.Unit_Vector; + +end Ada.Numerics.Generic_Real_Arrays; -- cgit v1.2.3