From 815f44d0cd91db97ee2793bcdf007f498f78f8aa Mon Sep 17 00:00:00 2001 From: Geert Bosch Date: Fri, 6 Apr 2007 11:23:23 +0200 Subject: i-fortra.ads: Add Double_Complex type. 2007-04-06 Geert Bosch Robert Dewar * i-fortra.ads: Add Double_Complex type. * impunit.adb: (Is_Known_Unit): New function Add Gnat.Byte_Swapping Add GNAT.SHA1 Add new Ada 2005 units Ada.Numerics.Generic_Complex_Arrays, Ada.Numerics.Generic_Real_Arrays, Ada.Numerics.Complex_Arrays, Ada.Numerics.Real_Arrays, Ada.Numerics.Long_Complex_Arrays, Ada.Numerics.Long_Long_Complex_Arrays, Ada.Numerics.Long_Long_Real_Arrays and Ada.Numerics.Long_Real_Arrays * impunit.ads (Is_Known_Unit): New function * a-ngcoar.adb, a-ngcoar.ads, a-ngrear.adb, a-ngrear.ads, a-nlcoar.ads, a-nllcar.ads, a-nllrar.ads, a-nlrear.ads, a-nucoar.ads, a-nurear.ads, g-bytswa.adb, g-bytswa-x86.adb, g-bytswa.ads, g-sha1.adb, g-sha1.ads, i-forbla.ads, i-forlap.ads, s-gearop.adb, s-gearop.ads, s-gecobl.adb, s-gecobl.ads, s-gecola.adb, s-gecola.ads, s-gerebl.adb, s-gerebl.ads, s-gerela.adb, s-gerela.ads: New files. * Makefile.rtl: Add g-bytswa, g-sha1, a-fzteio and a-izteio * a-fzteio.ads, a-izteio.ads: New Ada 2005 run-time units. From-SVN: r123579 --- gcc/ada/Makefile.rtl | 5 +- gcc/ada/a-fzteio.ads | 19 + gcc/ada/a-izteio.ads | 19 + gcc/ada/a-ngcoar.adb | 1501 ++++++++++++++++++++++++++++++++++++++++++++++++++ gcc/ada/a-ngcoar.ads | 281 ++++++++++ gcc/ada/a-ngrear.adb | 779 ++++++++++++++++++++++++++ gcc/ada/a-ngrear.ads | 141 +++++ gcc/ada/a-nlcoar.ads | 23 + gcc/ada/a-nllcar.ads | 24 + gcc/ada/a-nllrar.ads | 21 + gcc/ada/a-nlrear.ads | 21 + gcc/ada/a-nucoar.ads | 23 + gcc/ada/a-nurear.ads | 21 + gcc/ada/g-sha1.adb | 379 +++++++++++++ gcc/ada/g-sha1.ads | 116 ++++ gcc/ada/i-forbla.ads | 251 +++++++++ gcc/ada/i-forlap.ads | 416 ++++++++++++++ gcc/ada/i-fortra.ads | 7 +- gcc/ada/impunit.adb | 89 ++- gcc/ada/impunit.ads | 8 +- gcc/ada/s-gearop.adb | 518 +++++++++++++++++ gcc/ada/s-gearop.ads | 398 +++++++++++++ gcc/ada/s-gecobl.adb | 352 ++++++++++++ gcc/ada/s-gecobl.ads | 104 ++++ gcc/ada/s-gecola.adb | 495 +++++++++++++++++ gcc/ada/s-gecola.ads | 133 +++++ gcc/ada/s-gerebl.adb | 313 +++++++++++ gcc/ada/s-gerebl.ads | 98 ++++ gcc/ada/s-gerela.adb | 566 +++++++++++++++++++ gcc/ada/s-gerela.ads | 130 +++++ 30 files changed, 7246 insertions(+), 5 deletions(-) create mode 100755 gcc/ada/a-fzteio.ads create mode 100755 gcc/ada/a-izteio.ads create mode 100644 gcc/ada/a-ngcoar.adb create mode 100644 gcc/ada/a-ngcoar.ads create mode 100644 gcc/ada/a-ngrear.adb create mode 100644 gcc/ada/a-ngrear.ads create mode 100644 gcc/ada/a-nlcoar.ads create mode 100644 gcc/ada/a-nllcar.ads create mode 100644 gcc/ada/a-nllrar.ads create mode 100644 gcc/ada/a-nlrear.ads create mode 100644 gcc/ada/a-nucoar.ads create mode 100644 gcc/ada/a-nurear.ads create mode 100644 gcc/ada/g-sha1.adb create mode 100644 gcc/ada/g-sha1.ads create mode 100644 gcc/ada/i-forbla.ads create mode 100644 gcc/ada/i-forlap.ads create mode 100644 gcc/ada/s-gearop.adb create mode 100644 gcc/ada/s-gearop.ads create mode 100644 gcc/ada/s-gecobl.adb create mode 100644 gcc/ada/s-gecobl.ads create mode 100644 gcc/ada/s-gecola.adb create mode 100644 gcc/ada/s-gecola.ads create mode 100644 gcc/ada/s-gerebl.adb create mode 100644 gcc/ada/s-gerebl.ads create mode 100644 gcc/ada/s-gerela.adb create mode 100644 gcc/ada/s-gerela.ads diff --git a/gcc/ada/Makefile.rtl b/gcc/ada/Makefile.rtl index 64c3ef4..aed7b61 100644 --- a/gcc/ada/Makefile.rtl +++ b/gcc/ada/Makefile.rtl @@ -26,7 +26,6 @@ # Objects needed only for tasking GNATRTL_TASKING_OBJS= \ - a-diroro$(objext) \ a-dispat$(objext) \ a-dynpri$(objext) \ a-interr$(objext) \ @@ -135,9 +134,11 @@ GNATRTL_NONTASKING_OBJS= \ a-finali$(objext) \ a-flteio$(objext) \ a-fwteio$(objext) \ + a-fzteio$(objext) \ a-inteio$(objext) \ a-ioexce$(objext) \ a-iwteio$(objext) \ + a-izteio$(objext) \ a-lcteio$(objext) \ a-lfteio$(objext) \ a-llctio$(objext) \ @@ -313,6 +314,7 @@ GNATRTL_NONTASKING_OBJS= \ g-bubsor$(objext) \ g-busora$(objext) \ g-busorg$(objext) \ + g-bytswa$(objext) \ g-calend$(objext) \ g-casuti$(objext) \ g-catiio$(objext) \ @@ -350,6 +352,7 @@ GNATRTL_NONTASKING_OBJS= \ g-regexp$(objext) \ g-regpat$(objext) \ g-sestin$(objext) \ + g-sha1$(objext) \ g-soccon$(objext) \ g-socket$(objext) \ g-socthi$(objext) \ diff --git a/gcc/ada/a-fzteio.ads b/gcc/ada/a-fzteio.ads new file mode 100755 index 0000000..81bf7b2 --- /dev/null +++ b/gcc/ada/a-fzteio.ads @@ -0,0 +1,19 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- A D A . F L O A T _ W I D E _ W I D E _ T E X T _ I O -- +-- -- +-- S p e c -- +-- -- +-- This specification is derived from the Ada Reference Manual for use with -- +-- GNAT. In accordance with the copyright of that document, you can freely -- +-- copy and modify this specification, provided that if you redistribute a -- +-- modified version, any changes that you have made are clearly indicated. -- +-- -- +------------------------------------------------------------------------------ + +with Ada.Wide_Wide_Text_IO; + +package Ada.Float_Wide_Wide_Text_IO is + new Ada.Wide_Wide_Text_IO.Float_IO (Float); diff --git a/gcc/ada/a-izteio.ads b/gcc/ada/a-izteio.ads new file mode 100755 index 0000000..8eb5466 --- /dev/null +++ b/gcc/ada/a-izteio.ads @@ -0,0 +1,19 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- A D A . I N T E G E R _ W I D E _ W I D E _ T E X T _ I O -- +-- -- +-- S p e c -- +-- -- +-- This specification is derived from the Ada Reference Manual for use with -- +-- GNAT. In accordance with the copyright of that document, you can freely -- +-- copy and modify this specification, provided that if you redistribute a -- +-- modified version, any changes that you have made are clearly indicated. -- +-- -- +------------------------------------------------------------------------------ + +with Ada.Wide_Wide_Text_IO; + +package Ada.Integer_Wide_Wide_Text_IO is + new Ada.Wide_Wide_Text_IO.Integer_IO (Integer); diff --git a/gcc/ada/a-ngcoar.adb b/gcc/ada/a-ngcoar.adb new file mode 100644 index 0000000..7d5a51e --- /dev/null +++ b/gcc/ada/a-ngcoar.adb @@ -0,0 +1,1501 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT COMPILER COMPONENTS -- +-- -- +-- ADA.NUMERICS.GENERIC_COMPLEX_ARRAYS -- +-- -- +-- B o d y -- +-- -- +-- Copyright (C) 2006, 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 2, 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. See the GNU General Public License -- +-- for more details. You should have received a copy of the GNU General -- +-- Public License distributed with GNAT; see file COPYING. If not, write -- +-- to the Free Software Foundation, 51 Franklin Street, Fifth Floor, -- +-- Boston, MA 02110-1301, USA. -- +-- -- +-- As a special exception, if other files instantiate generics from this -- +-- unit, or you link this unit with other files to produce an executable, -- +-- this unit does not by itself cause the resulting executable to be -- +-- covered by the GNU General Public License. This exception does not -- +-- however invalidate any other reasons why the executable file might be -- +-- covered by the GNU Public License. -- +-- -- +-- GNAT was originally developed by the GNAT team at New York University. -- +-- Extensive contributions were provided by Ada Core Technologies Inc. -- +-- -- +------------------------------------------------------------------------------ + +with System.Generic_Array_Operations; use System.Generic_Array_Operations; +with System.Generic_Complex_BLAS; +with System.Generic_Complex_LAPACK; + +package body Ada.Numerics.Generic_Complex_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. + + type BLAS_Real_Vector is array (Integer range <>) of Real; + + package BLAS is new System.Generic_Complex_BLAS + (Real => Real, + Complex_Types => Complex_Types, + Complex_Vector => Complex_Vector, + Complex_Matrix => Complex_Matrix); + + package LAPACK is new System.Generic_Complex_LAPACK + (Real => Real, + Real_Vector => BLAS_Real_Vector, + Complex_Types => Complex_Types, + Complex_Vector => Complex_Vector, + Complex_Matrix => Complex_Matrix); + + use BLAS, LAPACK; + + -- Procedure versions of functions returning unconstrained values. + -- This allows for inlining the function wrapper. + + procedure Eigenvalues + (A : Complex_Matrix; + Values : out Real_Vector); + + procedure Inverse + (A : Complex_Matrix; + R : out Complex_Matrix); + + procedure Solve + (A : Complex_Matrix; + X : Complex_Vector; + B : out Complex_Vector); + + procedure Solve + (A : Complex_Matrix; + X : Complex_Matrix; + B : out Complex_Matrix); + + procedure Transpose is new System.Generic_Array_Operations.Transpose + (Scalar => Complex, + Matrix => Complex_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 (Complex, Complex_Matrix); + + -- Instantiating the following subprograms directly would lead to + -- name clashes, so use a local package. + + package Instantiations is + + --------- + -- "*" -- + --------- + + function "*" is new Vector_Scalar_Elementwise_Operation + (Left_Scalar => Complex, + Right_Scalar => Complex, + Result_Scalar => Complex, + Left_Vector => Complex_Vector, + Result_Vector => Complex_Vector, + Operation => "*"); + + function "*" is new Vector_Scalar_Elementwise_Operation + (Left_Scalar => Complex, + Right_Scalar => Real'Base, + Result_Scalar => Complex, + Left_Vector => Complex_Vector, + Result_Vector => Complex_Vector, + Operation => "*"); + + function "*" is new Scalar_Vector_Elementwise_Operation + (Left_Scalar => Complex, + Right_Scalar => Complex, + Result_Scalar => Complex, + Right_Vector => Complex_Vector, + Result_Vector => Complex_Vector, + Operation => "*"); + + function "*" is new Scalar_Vector_Elementwise_Operation + (Left_Scalar => Real'Base, + Right_Scalar => Complex, + Result_Scalar => Complex, + Right_Vector => Complex_Vector, + Result_Vector => Complex_Vector, + Operation => "*"); + + function "*" is new Inner_Product + (Left_Scalar => Complex, + Right_Scalar => Real'Base, + Result_Scalar => Complex, + Left_Vector => Complex_Vector, + Right_Vector => Real_Vector, + Zero => (0.0, 0.0)); + + function "*" is new Inner_Product + (Left_Scalar => Real'Base, + Right_Scalar => Complex, + Result_Scalar => Complex, + Left_Vector => Real_Vector, + Right_Vector => Complex_Vector, + Zero => (0.0, 0.0)); + + function "*" is new Outer_Product + (Left_Scalar => Complex, + Right_Scalar => Complex, + Result_Scalar => Complex, + Left_Vector => Complex_Vector, + Right_Vector => Complex_Vector, + Matrix => Complex_Matrix); + + function "*" is new Outer_Product + (Left_Scalar => Real'Base, + Right_Scalar => Complex, + Result_Scalar => Complex, + Left_Vector => Real_Vector, + Right_Vector => Complex_Vector, + Matrix => Complex_Matrix); + + function "*" is new Outer_Product + (Left_Scalar => Complex, + Right_Scalar => Real'Base, + Result_Scalar => Complex, + Left_Vector => Complex_Vector, + Right_Vector => Real_Vector, + Matrix => Complex_Matrix); + + function "*" is new Matrix_Scalar_Elementwise_Operation + (Left_Scalar => Complex, + Right_Scalar => Complex, + Result_Scalar => Complex, + Left_Matrix => Complex_Matrix, + Result_Matrix => Complex_Matrix, + Operation => "*"); + + function "*" is new Matrix_Scalar_Elementwise_Operation + (Left_Scalar => Complex, + Right_Scalar => Real'Base, + Result_Scalar => Complex, + Left_Matrix => Complex_Matrix, + Result_Matrix => Complex_Matrix, + Operation => "*"); + + function "*" is new Scalar_Matrix_Elementwise_Operation + (Left_Scalar => Complex, + Right_Scalar => Complex, + Result_Scalar => Complex, + Right_Matrix => Complex_Matrix, + Result_Matrix => Complex_Matrix, + Operation => "*"); + + function "*" is new Scalar_Matrix_Elementwise_Operation + (Left_Scalar => Real'Base, + Right_Scalar => Complex, + Result_Scalar => Complex, + Right_Matrix => Complex_Matrix, + Result_Matrix => Complex_Matrix, + Operation => "*"); + + function "*" is new Matrix_Vector_Product + (Left_Scalar => Real'Base, + Right_Scalar => Complex, + Result_Scalar => Complex, + Matrix => Real_Matrix, + Right_Vector => Complex_Vector, + Result_Vector => Complex_Vector, + Zero => (0.0, 0.0)); + + function "*" is new Matrix_Vector_Product + (Left_Scalar => Complex, + Right_Scalar => Real'Base, + Result_Scalar => Complex, + Matrix => Complex_Matrix, + Right_Vector => Real_Vector, + Result_Vector => Complex_Vector, + Zero => (0.0, 0.0)); + + function "*" is new Vector_Matrix_Product + (Left_Scalar => Real'Base, + Right_Scalar => Complex, + Result_Scalar => Complex, + Left_Vector => Real_Vector, + Matrix => Complex_Matrix, + Result_Vector => Complex_Vector, + Zero => (0.0, 0.0)); + + function "*" is new Vector_Matrix_Product + (Left_Scalar => Complex, + Right_Scalar => Real'Base, + Result_Scalar => Complex, + Left_Vector => Complex_Vector, + Matrix => Real_Matrix, + Result_Vector => Complex_Vector, + Zero => (0.0, 0.0)); + + function "*" is new Matrix_Matrix_Product + (Left_Scalar => Real'Base, + Right_Scalar => Complex, + Result_Scalar => Complex, + Left_Matrix => Real_Matrix, + Right_Matrix => Complex_Matrix, + Result_Matrix => Complex_Matrix, + Zero => (0.0, 0.0)); + + function "*" is new Matrix_Matrix_Product + (Left_Scalar => Complex, + Right_Scalar => Real'Base, + Result_Scalar => Complex, + Left_Matrix => Complex_Matrix, + Right_Matrix => Real_Matrix, + Result_Matrix => Complex_Matrix, + Zero => (0.0, 0.0)); + + --------- + -- "+" -- + --------- + + function "+" is new Vector_Elementwise_Operation + (X_Scalar => Complex, + Result_Scalar => Complex, + X_Vector => Complex_Vector, + Result_Vector => Complex_Vector, + Operation => "+"); + + function "+" is new Vector_Vector_Elementwise_Operation + (Left_Scalar => Complex, + Right_Scalar => Complex, + Result_Scalar => Complex, + Left_Vector => Complex_Vector, + Right_Vector => Complex_Vector, + Result_Vector => Complex_Vector, + Operation => "+"); + + function "+" is new Vector_Vector_Elementwise_Operation + (Left_Scalar => Real'Base, + Right_Scalar => Complex, + Result_Scalar => Complex, + Left_Vector => Real_Vector, + Right_Vector => Complex_Vector, + Result_Vector => Complex_Vector, + Operation => "+"); + + function "+" is new Vector_Vector_Elementwise_Operation + (Left_Scalar => Complex, + Right_Scalar => Real'Base, + Result_Scalar => Complex, + Left_Vector => Complex_Vector, + Right_Vector => Real_Vector, + Result_Vector => Complex_Vector, + Operation => "+"); + + function "+" is new Matrix_Elementwise_Operation + (X_Scalar => Complex, + Result_Scalar => Complex, + X_Matrix => Complex_Matrix, + Result_Matrix => Complex_Matrix, + Operation => "+"); + + function "+" is new Matrix_Matrix_Elementwise_Operation + (Left_Scalar => Complex, + Right_Scalar => Complex, + Result_Scalar => Complex, + Left_Matrix => Complex_Matrix, + Right_Matrix => Complex_Matrix, + Result_Matrix => Complex_Matrix, + Operation => "+"); + + function "+" is new Matrix_Matrix_Elementwise_Operation + (Left_Scalar => Real'Base, + Right_Scalar => Complex, + Result_Scalar => Complex, + Left_Matrix => Real_Matrix, + Right_Matrix => Complex_Matrix, + Result_Matrix => Complex_Matrix, + Operation => "+"); + + function "+" is new Matrix_Matrix_Elementwise_Operation + (Left_Scalar => Complex, + Right_Scalar => Real'Base, + Result_Scalar => Complex, + Left_Matrix => Complex_Matrix, + Right_Matrix => Real_Matrix, + Result_Matrix => Complex_Matrix, + Operation => "+"); + + --------- + -- "-" -- + --------- + + function "-" is new Vector_Elementwise_Operation + (X_Scalar => Complex, + Result_Scalar => Complex, + X_Vector => Complex_Vector, + Result_Vector => Complex_Vector, + Operation => "-"); + + function "-" is new Vector_Vector_Elementwise_Operation + (Left_Scalar => Complex, + Right_Scalar => Complex, + Result_Scalar => Complex, + Left_Vector => Complex_Vector, + Right_Vector => Complex_Vector, + Result_Vector => Complex_Vector, + Operation => "-"); + + function "-" is new Vector_Vector_Elementwise_Operation + (Left_Scalar => Real'Base, + Right_Scalar => Complex, + Result_Scalar => Complex, + Left_Vector => Real_Vector, + Right_Vector => Complex_Vector, + Result_Vector => Complex_Vector, + Operation => "-"); + + function "-" is new Vector_Vector_Elementwise_Operation + (Left_Scalar => Complex, + Right_Scalar => Real'Base, + Result_Scalar => Complex, + Left_Vector => Complex_Vector, + Right_Vector => Real_Vector, + Result_Vector => Complex_Vector, + Operation => "-"); + + function "-" is new Matrix_Elementwise_Operation + (X_Scalar => Complex, + Result_Scalar => Complex, + X_Matrix => Complex_Matrix, + Result_Matrix => Complex_Matrix, + Operation => "-"); + + function "-" is new Matrix_Matrix_Elementwise_Operation + (Left_Scalar => Complex, + Right_Scalar => Complex, + Result_Scalar => Complex, + Left_Matrix => Complex_Matrix, + Right_Matrix => Complex_Matrix, + Result_Matrix => Complex_Matrix, + Operation => "-"); + + function "-" is new Matrix_Matrix_Elementwise_Operation + (Left_Scalar => Real'Base, + Right_Scalar => Complex, + Result_Scalar => Complex, + Left_Matrix => Real_Matrix, + Right_Matrix => Complex_Matrix, + Result_Matrix => Complex_Matrix, + Operation => "-"); + + function "-" is new Matrix_Matrix_Elementwise_Operation + (Left_Scalar => Complex, + Right_Scalar => Real'Base, + Result_Scalar => Complex, + Left_Matrix => Complex_Matrix, + Right_Matrix => Real_Matrix, + Result_Matrix => Complex_Matrix, + Operation => "-"); + + --------- + -- "/" -- + --------- + + function "/" is new Vector_Scalar_Elementwise_Operation + (Left_Scalar => Complex, + Right_Scalar => Complex, + Result_Scalar => Complex, + Left_Vector => Complex_Vector, + Result_Vector => Complex_Vector, + Operation => "/"); + + function "/" is new Vector_Scalar_Elementwise_Operation + (Left_Scalar => Complex, + Right_Scalar => Real'Base, + Result_Scalar => Complex, + Left_Vector => Complex_Vector, + Result_Vector => Complex_Vector, + Operation => "/"); + + function "/" is new Matrix_Scalar_Elementwise_Operation + (Left_Scalar => Complex, + Right_Scalar => Complex, + Result_Scalar => Complex, + Left_Matrix => Complex_Matrix, + Result_Matrix => Complex_Matrix, + Operation => "/"); + + function "/" is new Matrix_Scalar_Elementwise_Operation + (Left_Scalar => Complex, + Right_Scalar => Real'Base, + Result_Scalar => Complex, + Left_Matrix => Complex_Matrix, + Result_Matrix => Complex_Matrix, + Operation => "/"); + + -------------- + -- Argument -- + -------------- + + function Argument is new Vector_Elementwise_Operation + (X_Scalar => Complex, + Result_Scalar => Real'Base, + X_Vector => Complex_Vector, + Result_Vector => Real_Vector, + Operation => Argument); + + function Argument is new Vector_Scalar_Elementwise_Operation + (Left_Scalar => Complex, + Right_Scalar => Real'Base, + Result_Scalar => Real'Base, + Left_Vector => Complex_Vector, + Result_Vector => Real_Vector, + Operation => Argument); + + function Argument is new Matrix_Elementwise_Operation + (X_Scalar => Complex, + Result_Scalar => Real'Base, + X_Matrix => Complex_Matrix, + Result_Matrix => Real_Matrix, + Operation => Argument); + + function Argument is new Matrix_Scalar_Elementwise_Operation + (Left_Scalar => Complex, + Right_Scalar => Real'Base, + Result_Scalar => Real'Base, + Left_Matrix => Complex_Matrix, + Result_Matrix => Real_Matrix, + Operation => Argument); + + ---------------------------- + -- Compose_From_Cartesian -- + ---------------------------- + + function Compose_From_Cartesian is new Vector_Elementwise_Operation + (X_Scalar => Real'Base, + Result_Scalar => Complex, + X_Vector => Real_Vector, + Result_Vector => Complex_Vector, + Operation => Compose_From_Cartesian); + + function Compose_From_Cartesian is + new Vector_Vector_Elementwise_Operation + (Left_Scalar => Real'Base, + Right_Scalar => Real'Base, + Result_Scalar => Complex, + Left_Vector => Real_Vector, + Right_Vector => Real_Vector, + Result_Vector => Complex_Vector, + Operation => Compose_From_Cartesian); + + function Compose_From_Cartesian is new Matrix_Elementwise_Operation + (X_Scalar => Real'Base, + Result_Scalar => Complex, + X_Matrix => Real_Matrix, + Result_Matrix => Complex_Matrix, + Operation => Compose_From_Cartesian); + + function Compose_From_Cartesian is + new Matrix_Matrix_Elementwise_Operation + (Left_Scalar => Real'Base, + Right_Scalar => Real'Base, + Result_Scalar => Complex, + Left_Matrix => Real_Matrix, + Right_Matrix => Real_Matrix, + Result_Matrix => Complex_Matrix, + Operation => Compose_From_Cartesian); + + ------------------------ + -- Compose_From_Polar -- + ------------------------ + + function Compose_From_Polar is + new Vector_Vector_Elementwise_Operation + (Left_Scalar => Real'Base, + Right_Scalar => Real'Base, + Result_Scalar => Complex, + Left_Vector => Real_Vector, + Right_Vector => Real_Vector, + Result_Vector => Complex_Vector, + Operation => Compose_From_Polar); + + function Compose_From_Polar is + new Vector_Vector_Scalar_Elementwise_Operation + (X_Scalar => Real'Base, + Y_Scalar => Real'Base, + Z_Scalar => Real'Base, + Result_Scalar => Complex, + X_Vector => Real_Vector, + Y_Vector => Real_Vector, + Result_Vector => Complex_Vector, + Operation => Compose_From_Polar); + + function Compose_From_Polar is + new Matrix_Matrix_Elementwise_Operation + (Left_Scalar => Real'Base, + Right_Scalar => Real'Base, + Result_Scalar => Complex, + Left_Matrix => Real_Matrix, + Right_Matrix => Real_Matrix, + Result_Matrix => Complex_Matrix, + Operation => Compose_From_Polar); + + function Compose_From_Polar is + new Matrix_Matrix_Scalar_Elementwise_Operation + (X_Scalar => Real'Base, + Y_Scalar => Real'Base, + Z_Scalar => Real'Base, + Result_Scalar => Complex, + X_Matrix => Real_Matrix, + Y_Matrix => Real_Matrix, + Result_Matrix => Complex_Matrix, + Operation => Compose_From_Polar); + + --------------- + -- Conjugate -- + --------------- + + function Conjugate is new Vector_Elementwise_Operation + (X_Scalar => Complex, + Result_Scalar => Complex, + X_Vector => Complex_Vector, + Result_Vector => Complex_Vector, + Operation => Conjugate); + + function Conjugate is new Matrix_Elementwise_Operation + (X_Scalar => Complex, + Result_Scalar => Complex, + X_Matrix => Complex_Matrix, + Result_Matrix => Complex_Matrix, + Operation => Conjugate); + + -------- + -- Im -- + -------- + + function Im is new Vector_Elementwise_Operation + (X_Scalar => Complex, + Result_Scalar => Real'Base, + X_Vector => Complex_Vector, + Result_Vector => Real_Vector, + Operation => Im); + + function Im is new Matrix_Elementwise_Operation + (X_Scalar => Complex, + Result_Scalar => Real'Base, + X_Matrix => Complex_Matrix, + Result_Matrix => Real_Matrix, + Operation => Im); + + ------------- + -- Modulus -- + ------------- + + function Modulus is new Vector_Elementwise_Operation + (X_Scalar => Complex, + Result_Scalar => Real'Base, + X_Vector => Complex_Vector, + Result_Vector => Real_Vector, + Operation => Modulus); + + function Modulus is new Matrix_Elementwise_Operation + (X_Scalar => Complex, + Result_Scalar => Real'Base, + X_Matrix => Complex_Matrix, + Result_Matrix => Real_Matrix, + Operation => Modulus); + + -------- + -- Re -- + -------- + + function Re is new Vector_Elementwise_Operation + (X_Scalar => Complex, + Result_Scalar => Real'Base, + X_Vector => Complex_Vector, + Result_Vector => Real_Vector, + Operation => Re); + + function Re is new Matrix_Elementwise_Operation + (X_Scalar => Complex, + Result_Scalar => Real'Base, + X_Matrix => Complex_Matrix, + Result_Matrix => Real_Matrix, + Operation => Re); + + ------------ + -- Set_Im -- + ------------ + + procedure Set_Im is new Update_Vector_With_Vector + (X_Scalar => Complex, + Y_Scalar => Real'Base, + X_Vector => Complex_Vector, + Y_Vector => Real_Vector, + Update => Set_Im); + + procedure Set_Im is new Update_Matrix_With_Matrix + (X_Scalar => Complex, + Y_Scalar => Real'Base, + X_Matrix => Complex_Matrix, + Y_Matrix => Real_Matrix, + Update => Set_Im); + + ------------ + -- Set_Re -- + ------------ + + procedure Set_Re is new Update_Vector_With_Vector + (X_Scalar => Complex, + Y_Scalar => Real'Base, + X_Vector => Complex_Vector, + Y_Vector => Real_Vector, + Update => Set_Re); + + procedure Set_Re is new Update_Matrix_With_Matrix + (X_Scalar => Complex, + Y_Scalar => Real'Base, + X_Matrix => Complex_Matrix, + Y_Matrix => Real_Matrix, + Update => Set_Re); + + ----------------- + -- Unit_Matrix -- + ----------------- + + function Unit_Matrix is new System.Generic_Array_Operations.Unit_Matrix + (Scalar => Complex, + Matrix => Complex_Matrix, + Zero => (0.0, 0.0), + One => (1.0, 0.0)); + + function Unit_Vector is new System.Generic_Array_Operations.Unit_Vector + (Scalar => Complex, + Vector => Complex_Vector, + Zero => (0.0, 0.0), + One => (1.0, 0.0)); + + end Instantiations; + + --------- + -- "*" -- + --------- + + function "*" + (Left : Complex_Vector; + Right : Complex_Vector) return Complex + 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 : Real_Vector; + Right : Complex_Vector) return Complex + renames Instantiations."*"; + + function "*" + (Left : Complex_Vector; + Right : Real_Vector) return Complex + renames Instantiations."*"; + + function "*" + (Left : Complex; + Right : Complex_Vector) return Complex_Vector + renames Instantiations."*"; + + function "*" + (Left : Complex_Vector; + Right : Complex) return Complex_Vector + renames Instantiations."*"; + + function "*" + (Left : Real'Base; + Right : Complex_Vector) return Complex_Vector + renames Instantiations."*"; + + function "*" + (Left : Complex_Vector; + Right : Real'Base) return Complex_Vector + renames Instantiations."*"; + + function "*" + (Left : Complex_Matrix; + Right : Complex_Matrix) + return Complex_Matrix + is + R : Complex_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 multipication"; + 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 : Complex_Vector; + Right : Complex_Vector) return Complex_Matrix + renames Instantiations."*"; + + function "*" + (Left : Complex_Vector; + Right : Complex_Matrix) return Complex_Vector + is + R : Complex_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 : Complex_Matrix; + Right : Complex_Vector) return Complex_Vector + is + R : Complex_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 "*"; + + function "*" + (Left : Real_Matrix; + Right : Complex_Matrix) return Complex_Matrix + renames Instantiations."*"; + + function "*" + (Left : Complex_Matrix; + Right : Real_Matrix) return Complex_Matrix + renames Instantiations."*"; + + function "*" + (Left : Real_Vector; + Right : Complex_Vector) return Complex_Matrix + renames Instantiations."*"; + + function "*" + (Left : Complex_Vector; + Right : Real_Vector) return Complex_Matrix + renames Instantiations."*"; + + function "*" + (Left : Real_Vector; + Right : Complex_Matrix) return Complex_Vector + renames Instantiations."*"; + + function "*" + (Left : Complex_Vector; + Right : Real_Matrix) return Complex_Vector + renames Instantiations."*"; + + function "*" + (Left : Real_Matrix; + Right : Complex_Vector) return Complex_Vector + renames Instantiations."*"; + + function "*" + (Left : Complex_Matrix; + Right : Real_Vector) return Complex_Vector + renames Instantiations."*"; + + function "*" + (Left : Complex; + Right : Complex_Matrix) return Complex_Matrix + renames Instantiations."*"; + + function "*" + (Left : Complex_Matrix; + Right : Complex) return Complex_Matrix + renames Instantiations."*"; + + function "*" + (Left : Real'Base; + Right : Complex_Matrix) return Complex_Matrix + renames Instantiations."*"; + + function "*" + (Left : Complex_Matrix; + Right : Real'Base) return Complex_Matrix + renames Instantiations."*"; + + --------- + -- "+" -- + --------- + + function "+" (Right : Complex_Vector) return Complex_Vector + renames Instantiations."+"; + + function "+" + (Left : Complex_Vector; + Right : Complex_Vector) return Complex_Vector + renames Instantiations."+"; + + function "+" + (Left : Real_Vector; + Right : Complex_Vector) return Complex_Vector + renames Instantiations."+"; + + function "+" + (Left : Complex_Vector; + Right : Real_Vector) return Complex_Vector + renames Instantiations."+"; + + function "+" (Right : Complex_Matrix) return Complex_Matrix + renames Instantiations."+"; + + function "+" + (Left : Complex_Matrix; + Right : Complex_Matrix) return Complex_Matrix + renames Instantiations."+"; + + function "+" + (Left : Real_Matrix; + Right : Complex_Matrix) return Complex_Matrix + renames Instantiations."+"; + + function "+" + (Left : Complex_Matrix; + Right : Real_Matrix) return Complex_Matrix + renames Instantiations."+"; + + --------- + -- "-" -- + --------- + + function "-" + (Right : Complex_Vector) return Complex_Vector + renames Instantiations."-"; + + function "-" + (Left : Complex_Vector; + Right : Complex_Vector) return Complex_Vector + renames Instantiations."-"; + + function "-" + (Left : Real_Vector; + Right : Complex_Vector) return Complex_Vector + renames Instantiations."-"; + + function "-" + (Left : Complex_Vector; + Right : Real_Vector) return Complex_Vector + renames Instantiations."-"; + + function "-" (Right : Complex_Matrix) return Complex_Matrix + renames Instantiations."-"; + + function "-" + (Left : Complex_Matrix; + Right : Complex_Matrix) return Complex_Matrix + renames Instantiations."-"; + + function "-" + (Left : Real_Matrix; + Right : Complex_Matrix) return Complex_Matrix + renames Instantiations."-"; + + function "-" + (Left : Complex_Matrix; + Right : Real_Matrix) return Complex_Matrix + renames Instantiations."-"; + + --------- + -- "/" -- + --------- + + function "/" + (Left : Complex_Vector; + Right : Complex) return Complex_Vector + renames Instantiations."/"; + + function "/" + (Left : Complex_Vector; + Right : Real'Base) return Complex_Vector + renames Instantiations."/"; + + function "/" + (Left : Complex_Matrix; + Right : Complex) return Complex_Matrix + renames Instantiations."/"; + + function "/" + (Left : Complex_Matrix; + Right : Real'Base) return Complex_Matrix + renames Instantiations."/"; + + ----------- + -- "abs" -- + ----------- + + function "abs" (Right : Complex_Vector) return Complex is + begin + return (nrm2 (Right'Length, Right), 0.0); + end "abs"; + + -------------- + -- Argument -- + -------------- + + function Argument (X : Complex_Vector) return Real_Vector + renames Instantiations.Argument; + + function Argument + (X : Complex_Vector; + Cycle : Real'Base) return Real_Vector + renames Instantiations.Argument; + + function Argument (X : Complex_Matrix) return Real_Matrix + renames Instantiations.Argument; + + function Argument + (X : Complex_Matrix; + Cycle : Real'Base) return Real_Matrix + renames Instantiations.Argument; + + ---------------------------- + -- Compose_From_Cartesian -- + ---------------------------- + + function Compose_From_Cartesian (Re : Real_Vector) return Complex_Vector + renames Instantiations.Compose_From_Cartesian; + + function Compose_From_Cartesian + (Re : Real_Vector; + Im : Real_Vector) return Complex_Vector + renames Instantiations.Compose_From_Cartesian; + + function Compose_From_Cartesian (Re : Real_Matrix) return Complex_Matrix + renames Instantiations.Compose_From_Cartesian; + + function Compose_From_Cartesian + (Re : Real_Matrix; + Im : Real_Matrix) return Complex_Matrix + renames Instantiations.Compose_From_Cartesian; + + ------------------------ + -- Compose_From_Polar -- + ------------------------ + + function Compose_From_Polar + (Modulus : Real_Vector; + Argument : Real_Vector) return Complex_Vector + renames Instantiations.Compose_From_Polar; + + function Compose_From_Polar + (Modulus : Real_Vector; + Argument : Real_Vector; + Cycle : Real'Base) return Complex_Vector + renames Instantiations.Compose_From_Polar; + + function Compose_From_Polar + (Modulus : Real_Matrix; + Argument : Real_Matrix) return Complex_Matrix + renames Instantiations.Compose_From_Polar; + + function Compose_From_Polar + (Modulus : Real_Matrix; + Argument : Real_Matrix; + Cycle : Real'Base) return Complex_Matrix + renames Instantiations.Compose_From_Polar; + + --------------- + -- Conjugate -- + --------------- + + function Conjugate (X : Complex_Vector) return Complex_Vector + renames Instantiations.Conjugate; + + function Conjugate (X : Complex_Matrix) return Complex_Matrix + renames Instantiations.Conjugate; + + ----------------- + -- Determinant -- + ----------------- + + function Determinant (A : Complex_Matrix) return Complex is + N : constant Integer := Length (A); + LU : Complex_Matrix (1 .. N, 1 .. N) := A; + Piv : Integer_Vector (1 .. N); + Info : aliased Integer := -1; + Neg : Boolean; + Det : Complex; + + begin + if N = 0 then + return (0.0, 0.0); + end if; + + getrf (N, N, LU, N, Piv, Info'Access); + + if Info /= 0 then + raise Constraint_Error with "ill-conditioned matrix"; + end if; + + Det := LU (1, 1); + Neg := Piv (1) /= 1; + + for J in 2 .. N loop + Det := Det * LU (J, J); + Neg := Neg xor (Piv (J) /= J); + end loop; + + if Neg then + return -Det; + + else + return Det; + end if; + end Determinant; + + ----------------- + -- Eigensystem -- + ----------------- + + procedure Eigensystem + (A : in Complex_Matrix; + Values : out Real_Vector; + Vectors : out Complex_Matrix) + is + Job_Z : aliased Character := 'V'; + Rng : aliased Character := 'A'; + Uplo : aliased Character := 'U'; + + N : constant Natural := Length (A); + W : BLAS_Real_Vector (Values'Range); + M : Integer; + B : Complex_Matrix (1 .. N, 1 .. N); + L_Work : Complex_Vector (1 .. 1); + LR_Work : BLAS_Real_Vector (1 .. 1); + LI_Work : Integer_Vector (1 .. 1); + I_Supp_Z : Integer_Vector (1 .. 2 * N); + Info : aliased Integer; + + begin + if Values'Length /= N then + raise Constraint_Error with "wrong length for output vector"; + end if; + + if Vectors'First (1) /= A'First (1) + or else Vectors'Last (1) /= A'Last (1) + or else Vectors'First (2) /= A'First (2) + or else Vectors'Last (2) /= A'Last (2) + then + raise Constraint_Error with "wrong dimensions for output matrix"; + end if; + + if N = 0 then + return; + end if; + + -- Check for hermitian matrix ??? + -- Copy only required triangle ??? + + B := A; + + -- Find size of work area + + heevr + (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N, + M => M, + W => W, + Z => Vectors, + Ld_Z => N, + I_Supp_Z => I_Supp_Z, + Work => L_Work, + L_Work => -1, + R_Work => LR_Work, + LR_Work => -1, + I_Work => LI_Work, + LI_Work => -1, + Info => Info'Access); + + if Info /= 0 then + raise Constraint_Error; + end if; + + declare + Work : Complex_Vector (1 .. Integer (L_Work (1).Re)); + R_Work : BLAS_Real_Vector (1 .. Integer (LR_Work (1))); + I_Work : Integer_Vector (1 .. LI_Work (1)); + + begin + heevr + (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N, + M => M, + W => W, + Z => Vectors, + Ld_Z => N, + I_Supp_Z => I_Supp_Z, + Work => Work, + L_Work => Work'Length, + R_Work => R_Work, + LR_Work => LR_Work'Length, + I_Work => I_Work, + LI_Work => LI_Work'Length, + Info => Info'Access); + + if Info /= 0 then + raise Constraint_Error with "inverting non-Hermetian matrix"; + end if; + + for J in Values'Range loop + Values (J) := W (J); + end loop; + end; + end Eigensystem; + + ----------------- + -- Eigenvalues -- + ----------------- + + procedure Eigenvalues + (A : Complex_Matrix; + Values : out Real_Vector) + is + Job_Z : aliased Character := 'N'; + Rng : aliased Character := 'A'; + Uplo : aliased Character := 'U'; + N : constant Natural := Length (A); + B : Complex_Matrix (1 .. N, 1 .. N) := A; + Z : Complex_Matrix (1 .. 1, 1 .. 1); + W : BLAS_Real_Vector (Values'Range); + L_Work : Complex_Vector (1 .. 1); + LR_Work : BLAS_Real_Vector (1 .. 1); + LI_Work : Integer_Vector (1 .. 1); + I_Supp_Z : Integer_Vector (1 .. 2 * N); + M : Integer; + Info : aliased Integer; + + begin + if Values'Length /= N then + raise Constraint_Error with "wrong length for output vector"; + end if; + + if N = 0 then + return; + end if; + + -- Check for hermitian matrix ??? + + -- Find size of work area + + heevr (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N, + M => M, + W => W, + Z => Z, + Ld_Z => 1, + I_Supp_Z => I_Supp_Z, + Work => L_Work, L_Work => -1, + R_Work => LR_Work, LR_Work => -1, + I_Work => LI_Work, LI_Work => -1, + Info => Info'Access); + + if Info /= 0 then + raise Constraint_Error; + end if; + + declare + Work : Complex_Vector (1 .. Integer (L_Work (1).Re)); + R_Work : BLAS_Real_Vector (1 .. Integer (LR_Work (1))); + I_Work : Integer_Vector (1 .. LI_Work (1)); + begin + heevr (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N, + M => M, + W => W, + Z => Z, + Ld_Z => 1, + I_Supp_Z => I_Supp_Z, + Work => Work, L_Work => Work'Length, + R_Work => R_Work, LR_Work => R_Work'Length, + I_Work => I_Work, LI_Work => I_Work'Length, + Info => Info'Access); + + if Info /= 0 then + raise Constraint_Error with "inverting singular matrix"; + end if; + + for J in Values'Range loop + Values (J) := W (J); + end loop; + end; + end Eigenvalues; + + function Eigenvalues (A : Complex_Matrix) return Real_Vector is + R : Real_Vector (A'Range (1)); + begin + Eigenvalues (A, R); + return R; + end Eigenvalues; + + -------- + -- Im -- + -------- + + function Im (X : Complex_Vector) return Real_Vector + renames Instantiations.Im; + + function Im (X : Complex_Matrix) return Real_Matrix + renames Instantiations.Im; + + ------------- + -- Inverse -- + ------------- + + procedure Inverse (A : Complex_Matrix; R : out Complex_Matrix) is + N : constant Integer := Length (A); + Piv : Integer_Vector (1 .. N); + L_Work : Complex_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 : Complex_Vector (1 .. Integer (L_Work (1).Re)); + + 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 : Complex_Matrix) return Complex_Matrix is + R : Complex_Matrix (A'Range (2), A'Range (1)); + begin + Inverse (A, R); + return R; + end Inverse; + + ------------- + -- Modulus -- + ------------- + + function Modulus (X : Complex_Vector) return Real_Vector + renames Instantiations.Modulus; + + function Modulus (X : Complex_Matrix) return Real_Matrix + renames Instantiations.Modulus; + + -------- + -- Re -- + -------- + + function Re (X : Complex_Vector) return Real_Vector + renames Instantiations.Re; + + function Re (X : Complex_Matrix) return Real_Matrix + renames Instantiations.Re; + + ------------ + -- Set_Im -- + ------------ + + procedure Set_Im + (X : in out Complex_Matrix; + Im : Real_Matrix) + renames Instantiations.Set_Im; + + procedure Set_Im + (X : in out Complex_Vector; + Im : Real_Vector) + renames Instantiations.Set_Im; + + ------------ + -- Set_Re -- + ------------ + + procedure Set_Re + (X : in out Complex_Matrix; + Re : Real_Matrix) + renames Instantiations.Set_Re; + + procedure Set_Re + (X : in out Complex_Vector; + Re : Real_Vector) + renames Instantiations.Set_Re; + + ----------- + -- Solve -- + ----------- + + procedure Solve + (A : Complex_Matrix; + X : Complex_Vector; + B : out Complex_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 : Complex_Matrix; + X : Complex_Matrix; + B : out Complex_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 : Complex_Matrix; + X : Complex_Vector) return Complex_Vector + is + B : Complex_Vector (A'Range (2)); + begin + Solve (A, X, B); + return B; + end Solve; + + function Solve (A, X : Complex_Matrix) return Complex_Matrix is + B : Complex_Matrix (A'Range (2), X'Range (2)); + begin + Solve (A, X, B); + return B; + end Solve; + + --------------- + -- Transpose -- + --------------- + + function Transpose + (X : Complex_Matrix) return Complex_Matrix + is + R : Complex_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 Complex_Matrix + renames Instantiations.Unit_Matrix; + + ----------------- + -- Unit_Vector -- + ----------------- + + function Unit_Vector + (Index : Integer; + Order : Positive; + First : Integer := 1) return Complex_Vector + renames Instantiations.Unit_Vector; + +end Ada.Numerics.Generic_Complex_Arrays; diff --git a/gcc/ada/a-ngcoar.ads b/gcc/ada/a-ngcoar.ads new file mode 100644 index 0000000..f94721c --- /dev/null +++ b/gcc/ada/a-ngcoar.ads @@ -0,0 +1,281 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- ADA.NUMERICS.GENERIC_COMPLEX_ARRAYS -- +-- -- +-- S p e c -- +-- -- +-- This specification is derived from the Ada Reference Manual for use with -- +-- GNAT. In accordance with the copyright of that document, you can freely -- +-- copy and modify this specification, provided that if you redistribute a -- +-- modified version, any changes that you have made are clearly indicated. -- +-- -- +------------------------------------------------------------------------------ + +with Ada.Numerics.Generic_Real_Arrays, Ada.Numerics.Generic_Complex_Types; + +generic + with package Real_Arrays is new Ada.Numerics.Generic_Real_Arrays (<>); + use Real_Arrays; + with package Complex_Types is new Ada.Numerics.Generic_Complex_Types (Real); + use Complex_Types; +package Ada.Numerics.Generic_Complex_Arrays is + pragma Pure (Generic_Complex_Arrays); + + -- Types + + type Complex_Vector is array (Integer range <>) of Complex; + type Complex_Matrix is + array (Integer range <>, Integer range <>) of Complex; + + -- Subprograms for Complex_Vector types + -- Complex_Vector selection, conversion and composition operations + + function Re (X : Complex_Vector) return Real_Vector; + function Im (X : Complex_Vector) return Real_Vector; + + procedure Set_Re (X : in out Complex_Vector; Re : in Real_Vector); + procedure Set_Im (X : in out Complex_Vector; Im : in Real_Vector); + + function Compose_From_Cartesian + (Re : Real_Vector) return Complex_Vector; + function Compose_From_Cartesian + (Re, Im : Real_Vector) return Complex_Vector; + + function Modulus (X : Complex_Vector) return Real_Vector; + function "abs" (Right : Complex_Vector) return Real_Vector renames Modulus; + function Argument (X : Complex_Vector) return Real_Vector; + + function Argument + (X : Complex_Vector; + Cycle : Real'Base) return Real_Vector; + + function Compose_From_Polar + (Modulus, Argument : Real_Vector) return Complex_Vector; + + function Compose_From_Polar + (Modulus, Argument : Real_Vector; + Cycle : Real'Base) return Complex_Vector; + + -- Complex_Vector arithmetic operations + + function "+" (Right : Complex_Vector) return Complex_Vector; + function "-" (Right : Complex_Vector) return Complex_Vector; + function Conjugate (X : Complex_Vector) return Complex_Vector; + function "+" (Left, Right : Complex_Vector) return Complex_Vector; + function "-" (Left, Right : Complex_Vector) return Complex_Vector; + function "*" (Left, Right : Complex_Vector) return Complex; + function "abs" (Right : Complex_Vector) return Complex; + + -- Mixed Real_Vector and Complex_Vector arithmetic operations + + function "+" + (Left : Real_Vector; + Right : Complex_Vector) return Complex_Vector; + + function "+" + (Left : Complex_Vector; + Right : Real_Vector) return Complex_Vector; + + function "-" + (Left : Real_Vector; + Right : Complex_Vector) return Complex_Vector; + + function "-" + (Left : Complex_Vector; + Right : Real_Vector) return Complex_Vector; + + function "*" (Left : Real_Vector; Right : Complex_Vector) return Complex; + function "*" (Left : Complex_Vector; Right : Real_Vector) return Complex; + + -- Complex_Vector scaling operations + + function "*" + (Left : Complex; + Right : Complex_Vector) return Complex_Vector; + + function "*" + (Left : Complex_Vector; + Right : Complex) return Complex_Vector; + + function "/" + (Left : Complex_Vector; + Right : Complex) return Complex_Vector; + + function "*" + (Left : Real'Base; + Right : Complex_Vector) return Complex_Vector; + + function "*" + (Left : Complex_Vector; + Right : Real'Base) return Complex_Vector; + + function "/" + (Left : Complex_Vector; + Right : Real'Base) return Complex_Vector; + + -- Other Complex_Vector operations + + function Unit_Vector + (Index : Integer; + Order : Positive; + First : Integer := 1) return Complex_Vector; + + -- Subprograms for Complex_Matrix types + + -- Complex_Matrix selection, conversion and composition operations + + function Re (X : Complex_Matrix) return Real_Matrix; + function Im (X : Complex_Matrix) return Real_Matrix; + + procedure Set_Re (X : in out Complex_Matrix; Re : in Real_Matrix); + procedure Set_Im (X : in out Complex_Matrix; Im : in Real_Matrix); + + function Compose_From_Cartesian (Re : Real_Matrix) return Complex_Matrix; + + function Compose_From_Cartesian + (Re, Im : Real_Matrix) return Complex_Matrix; + + function Modulus (X : Complex_Matrix) return Real_Matrix; + function "abs" (Right : Complex_Matrix) return Real_Matrix renames Modulus; + + function Argument (X : Complex_Matrix) return Real_Matrix; + + function Argument + (X : Complex_Matrix; + Cycle : Real'Base) return Real_Matrix; + + function Compose_From_Polar + (Modulus, Argument : Real_Matrix) return Complex_Matrix; + + function Compose_From_Polar + (Modulus : Real_Matrix; + Argument : Real_Matrix; + Cycle : Real'Base) return Complex_Matrix; + + -- Complex_Matrix arithmetic operations + + function "+" (Right : Complex_Matrix) return Complex_Matrix; + function "-" (Right : Complex_Matrix) return Complex_Matrix; + + function Conjugate (X : Complex_Matrix) return Complex_Matrix; + function Transpose (X : Complex_Matrix) return Complex_Matrix; + + function "+" (Left, Right : Complex_Matrix) return Complex_Matrix; + function "-" (Left, Right : Complex_Matrix) return Complex_Matrix; + function "*" (Left, Right : Complex_Matrix) return Complex_Matrix; + function "*" (Left, Right : Complex_Vector) return Complex_Matrix; + + function "*" + (Left : Complex_Vector; + Right : Complex_Matrix) return Complex_Vector; + + function "*" + (Left : Complex_Matrix; + Right : Complex_Vector) return Complex_Vector; + + -- Mixed Real_Matrix and Complex_Matrix arithmetic operations + + function "+" + (Left : Real_Matrix; + Right : Complex_Matrix) return Complex_Matrix; + + function "+" + (Left : Complex_Matrix; + Right : Real_Matrix) return Complex_Matrix; + + function "-" + (Left : Real_Matrix; + Right : Complex_Matrix) return Complex_Matrix; + + function "-" + (Left : Complex_Matrix; + Right : Real_Matrix) return Complex_Matrix; + + function "*" + (Left : Real_Matrix; + Right : Complex_Matrix) return Complex_Matrix; + + function "*" + (Left : Complex_Matrix; + Right : Real_Matrix) return Complex_Matrix; + + function "*" + (Left : Real_Vector; + Right : Complex_Vector) return Complex_Matrix; + + function "*" + (Left : Complex_Vector; + Right : Real_Vector) return Complex_Matrix; + + function "*" + (Left : Real_Vector; + Right : Complex_Matrix) return Complex_Vector; + + function "*" + (Left : Complex_Vector; + Right : Real_Matrix) return Complex_Vector; + + function "*" + (Left : Real_Matrix; + Right : Complex_Vector) return Complex_Vector; + + function "*" + (Left : Complex_Matrix; + Right : Real_Vector) return Complex_Vector; + + -- Complex_Matrix scaling operations + + function "*" + (Left : Complex; + Right : Complex_Matrix) return Complex_Matrix; + + function "*" + (Left : Complex_Matrix; + Right : Complex) return Complex_Matrix; + + function "/" + (Left : Complex_Matrix; + Right : Complex) return Complex_Matrix; + + function "*" + (Left : Real'Base; + Right : Complex_Matrix) return Complex_Matrix; + + function "*" + (Left : Complex_Matrix; + Right : Real'Base) return Complex_Matrix; + + function "/" + (Left : Complex_Matrix; + Right : Real'Base) return Complex_Matrix; + + -- Complex_Matrix inversion and related operations + + function Solve + (A : Complex_Matrix; + X : Complex_Vector) return Complex_Vector; + + function Solve (A, X : Complex_Matrix) return Complex_Matrix; + + function Inverse (A : Complex_Matrix) return Complex_Matrix; + + function Determinant (A : Complex_Matrix) return Complex; + + -- Eigenvalues and vectors of a Hermitian matrix + + function Eigenvalues (A : Complex_Matrix) return Real_Vector; + + procedure Eigensystem + (A : in Complex_Matrix; + Values : out Real_Vector; + Vectors : out Complex_Matrix); + + -- Other Complex_Matrix operations + + function Unit_Matrix + (Order : Positive; + First_1, First_2 : Integer := 1) return Complex_Matrix; + +end Ada.Numerics.Generic_Complex_Arrays; diff --git a/gcc/ada/a-ngrear.adb b/gcc/ada/a-ngrear.adb new file mode 100644 index 0000000..2ff5d01 --- /dev/null +++ b/gcc/ada/a-ngrear.adb @@ -0,0 +1,779 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- ADA.NUMERICS.GENERIC_REAL_ARRAYS -- +-- -- +-- B o d y -- +-- -- +-- Copyright (C) 2006, 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 2, 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. See the GNU General Public License -- +-- for more details. You should have received a copy of the GNU General -- +-- Public License distributed with GNAT; see file COPYING. If not, write -- +-- to the Free Software Foundation, 51 Franklin Street, Fifth Floor, -- +-- Boston, MA 02110-1301, USA. -- +-- -- +-- As a special exception, if other files instantiate generics from this -- +-- unit, or you link this unit with other files to produce an executable, -- +-- this unit does not by itself cause the resulting executable to be -- +-- covered by the GNU General Public License. This exception does not -- +-- however invalidate any other reasons why the executable file might be -- +-- covered by the GNU Public License. -- +-- -- +-- 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 multipication"; + 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 + if Piv (J) /= J then + Det := -Det * LU (J, J); + else + Det := Det * LU (J, J); + end if; + 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); + E : Real_Vector (1 .. N); + Tau : Real_Vector (1 .. N); + L_Work : Real_Vector (1 .. 1); + Info : aliased Integer; + + 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)); + 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); + B : Real_Matrix (1 .. N, 1 .. N); + E : Real_Vector (1 .. N); + Tau : Real_Vector (1 .. N); + L_Work : Real_Vector (1 .. 1); + Info : aliased Integer; + + 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)); + + 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))); + 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; diff --git a/gcc/ada/a-ngrear.ads b/gcc/ada/a-ngrear.ads new file mode 100644 index 0000000..61e1f34 --- /dev/null +++ b/gcc/ada/a-ngrear.ads @@ -0,0 +1,141 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- ADA.NUMERICS.GENERIC_REAL_ARRAYS -- +-- -- +-- S p e c -- +-- -- +-- Copyright (C) 2006, Free Software Foundation, Inc. -- +-- -- +-- This specification is derived from the Ada Reference Manual for use with -- +-- GNAT. The copyright notice above, and the license provisions that follow -- +-- apply solely to the contents of the part following the private keyword. -- +-- -- +-- 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 2, 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. See the GNU General Public License -- +-- for more details. You should have received a copy of the GNU General -- +-- Public License distributed with GNAT; see file COPYING. If not, write -- +-- to the Free Software Foundation, 51 Franklin Street, Fifth Floor, -- +-- Boston, MA 02110-1301, USA. -- +-- -- +-- As a special exception, if other files instantiate generics from this -- +-- unit, or you link this unit with other files to produce an executable, -- +-- this unit does not by itself cause the resulting executable to be -- +-- covered by the GNU General Public License. This exception does not -- +-- however invalidate any other reasons why the executable file might be -- +-- covered by the GNU Public License. -- +-- -- +-- GNAT was originally developed by the GNAT team at New York University. -- +-- Extensive contributions were provided by Ada Core Technologies Inc. -- +-- -- +------------------------------------------------------------------------------ + +generic + type Real is digits <>; +package Ada.Numerics.Generic_Real_Arrays is + pragma Pure (Generic_Real_Arrays); + + -- Types + + type Real_Vector is array (Integer range <>) of Real'Base; + type Real_Matrix is array (Integer range <>, Integer range <>) of Real'Base; + + -- Subprograms for Real_Vector types + + -- Real_Vector arithmetic operations + + function "+" (Right : Real_Vector) return Real_Vector; + function "-" (Right : Real_Vector) return Real_Vector; + function "abs" (Right : Real_Vector) return Real_Vector; + + function "+" (Left, Right : Real_Vector) return Real_Vector; + function "-" (Left, Right : Real_Vector) return Real_Vector; + + function "*" (Left, Right : Real_Vector) return Real'Base; + + function "abs" (Right : Real_Vector) return Real'Base; + + -- Real_Vector scaling operations + + function "*" (Left : Real'Base; Right : Real_Vector) return Real_Vector; + function "*" (Left : Real_Vector; Right : Real'Base) return Real_Vector; + function "/" (Left : Real_Vector; Right : Real'Base) return Real_Vector; + + -- Other Real_Vector operations + + function Unit_Vector + (Index : Integer; + Order : Positive; + First : Integer := 1) return Real_Vector; + + -- Subprograms for Real_Matrix types + + -- Real_Matrix arithmetic operations + + function "+" (Right : Real_Matrix) return Real_Matrix; + function "-" (Right : Real_Matrix) return Real_Matrix; + function "abs" (Right : Real_Matrix) return Real_Matrix; + function Transpose (X : Real_Matrix) return Real_Matrix; + + function "+" (Left, Right : Real_Matrix) return Real_Matrix; + function "-" (Left, Right : Real_Matrix) return Real_Matrix; + function "*" (Left, Right : Real_Matrix) return Real_Matrix; + + function "*" (Left, Right : Real_Vector) return Real_Matrix; + + function "*" (Left : Real_Vector; Right : Real_Matrix) return Real_Vector; + function "*" (Left : Real_Matrix; Right : Real_Vector) return Real_Vector; + + -- Real_Matrix scaling operations + + function "*" (Left : Real'Base; Right : Real_Matrix) return Real_Matrix; + function "*" (Left : Real_Matrix; Right : Real'Base) return Real_Matrix; + function "/" (Left : Real_Matrix; Right : Real'Base) return Real_Matrix; + + -- Real_Matrix inversion and related operations + + function Solve (A : Real_Matrix; X : Real_Vector) return Real_Vector; + function Solve (A, X : Real_Matrix) return Real_Matrix; + function Inverse (A : Real_Matrix) return Real_Matrix; + function Determinant (A : Real_Matrix) return Real'Base; + + -- Eigenvalues and vectors of a real symmetric matrix + + function Eigenvalues (A : Real_Matrix) return Real_Vector; + + procedure Eigensystem + (A : Real_Matrix; + Values : out Real_Vector; + Vectors : out Real_Matrix); + + -- Other Real_Matrix operations + + function Unit_Matrix + (Order : Positive; + First_1 : Integer := 1; + First_2 : Integer := 1) return Real_Matrix; + +private + -- The following operations are either relatively simple compared to the + -- expense of returning unconstrained arrays, or are just function wrappers + -- calling procedures implementing the actual operation. By having the + -- front end always inline these, the expense of the unconstrained returns + -- can be avoided. + + pragma Inline_Always ("+"); + pragma Inline_Always ("-"); + pragma Inline_Always ("*"); + pragma Inline_Always ("/"); + pragma Inline_Always ("abs"); + pragma Inline_Always (Eigenvalues); + pragma Inline_Always (Inverse); + pragma Inline_Always (Solve); + pragma Inline_Always (Transpose); + pragma Inline_Always (Unit_Matrix); + pragma Inline_Always (Unit_Vector); +end Ada.Numerics.Generic_Real_Arrays; diff --git a/gcc/ada/a-nlcoar.ads b/gcc/ada/a-nlcoar.ads new file mode 100644 index 0000000..35e97a5 --- /dev/null +++ b/gcc/ada/a-nlcoar.ads @@ -0,0 +1,23 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- ADA.NUMERICS.LONG_COMPLEX_ARRAYS -- +-- -- +-- S p e c -- +-- -- +-- This specification is derived from the Ada Reference Manual for use with -- +-- GNAT. In accordance with the copyright of that document, you can freely -- +-- copy and modify this specification, provided that if you redistribute a -- +-- modified version, any changes that you have made are clearly indicated. -- +-- -- +------------------------------------------------------------------------------ + +with Ada.Numerics.Generic_Complex_Arrays; +with Ada.Numerics.Long_Real_Arrays; +with Ada.Numerics.Long_Complex_Types; + +package Ada.Numerics.Long_Complex_Arrays is new + Ada.Numerics.Generic_Complex_Arrays (Long_Real_Arrays, Long_Complex_Types); + +pragma Pure (Long_Complex_Arrays); diff --git a/gcc/ada/a-nllcar.ads b/gcc/ada/a-nllcar.ads new file mode 100644 index 0000000..48fd91a --- /dev/null +++ b/gcc/ada/a-nllcar.ads @@ -0,0 +1,24 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- ADA.NUMERICS.LONG_LONG_COMPLEX_ARRAYS -- +-- -- +-- S p e c -- +-- -- +-- This specification is derived from the Ada Reference Manual for use with -- +-- GNAT. In accordance with the copyright of that document, you can freely -- +-- copy and modify this specification, provided that if you redistribute a -- +-- modified version, any changes that you have made are clearly indicated. -- +-- -- +------------------------------------------------------------------------------ + +with Ada.Numerics.Generic_Complex_Arrays; +with Ada.Numerics.Long_Long_Real_Arrays; +with Ada.Numerics.Long_Long_Complex_Types; + +package Ada.Numerics.Long_Long_Complex_Arrays is + new Ada.Numerics.Generic_Complex_Arrays (Long_Long_Real_Arrays, + Long_Long_Complex_Types); + +pragma Pure (Long_Long_Complex_Arrays); diff --git a/gcc/ada/a-nllrar.ads b/gcc/ada/a-nllrar.ads new file mode 100644 index 0000000..8a1713b --- /dev/null +++ b/gcc/ada/a-nllrar.ads @@ -0,0 +1,21 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- ADA.NUMERICS.LONG_LONG_REAL_ARRAYS -- +-- -- +-- S p e c -- +-- -- +-- This specification is derived from the Ada Reference Manual for use with -- +-- GNAT. In accordance with the copyright of that document, you can freely -- +-- copy and modify this specification, provided that if you redistribute a -- +-- modified version, any changes that you have made are clearly indicated. -- +-- -- +------------------------------------------------------------------------------ + +with Ada.Numerics.Generic_Real_Arrays; + +package Ada.Numerics.Long_Long_Real_Arrays is + new Ada.Numerics.Generic_Real_Arrays (Long_Long_Float); + +pragma Pure (Long_Long_Real_Arrays); diff --git a/gcc/ada/a-nlrear.ads b/gcc/ada/a-nlrear.ads new file mode 100644 index 0000000..252c5d2 --- /dev/null +++ b/gcc/ada/a-nlrear.ads @@ -0,0 +1,21 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- ADA.NUMERICS.LONG_REAL_ARRAYS -- +-- -- +-- S p e c -- +-- -- +-- This specification is derived from the Ada Reference Manual for use with -- +-- GNAT. In accordance with the copyright of that document, you can freely -- +-- copy and modify this specification, provided that if you redistribute a -- +-- modified version, any changes that you have made are clearly indicated. -- +-- -- +------------------------------------------------------------------------------ + +with Ada.Numerics.Generic_Real_Arrays; + +package Ada.Numerics.Long_Real_Arrays is + new Ada.Numerics.Generic_Real_Arrays (Long_Float); + +pragma Pure (Long_Real_Arrays); diff --git a/gcc/ada/a-nucoar.ads b/gcc/ada/a-nucoar.ads new file mode 100644 index 0000000..938e70a --- /dev/null +++ b/gcc/ada/a-nucoar.ads @@ -0,0 +1,23 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- ADA.NUMERICS.COMPLEX_ARRAYS -- +-- -- +-- S p e c -- +-- -- +-- This specification is derived from the Ada Reference Manual for use with -- +-- GNAT. In accordance with the copyright of that document, you can freely -- +-- copy and modify this specification, provided that if you redistribute a -- +-- modified version, any changes that you have made are clearly indicated. -- +-- -- +------------------------------------------------------------------------------ + +with Ada.Numerics.Generic_Complex_Arrays; +with Ada.Numerics.Real_Arrays; +with Ada.Numerics.Complex_Types; + +package Ada.Numerics.Complex_Arrays is + new Ada.Numerics.Generic_Complex_Arrays (Real_Arrays, Complex_Types); + +pragma Pure (Complex_Arrays); diff --git a/gcc/ada/a-nurear.ads b/gcc/ada/a-nurear.ads new file mode 100644 index 0000000..5141a16 --- /dev/null +++ b/gcc/ada/a-nurear.ads @@ -0,0 +1,21 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- ADA.NUMERICS.REAL_ARRAYS -- +-- -- +-- S p e c -- +-- -- +-- This specification is derived from the Ada Reference Manual for use with -- +-- GNAT. In accordance with the copyright of that document, you can freely -- +-- copy and modify this specification, provided that if you redistribute a -- +-- modified version, any changes that you have made are clearly indicated. -- +-- -- +------------------------------------------------------------------------------ + +with Ada.Numerics.Generic_Real_Arrays; + +package Ada.Numerics.Real_Arrays is + new Ada.Numerics.Generic_Real_Arrays (Float); + +pragma Pure (Real_Arrays); diff --git a/gcc/ada/g-sha1.adb b/gcc/ada/g-sha1.adb new file mode 100644 index 0000000..72b1924 --- /dev/null +++ b/gcc/ada/g-sha1.adb @@ -0,0 +1,379 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT LIBRARY COMPONENTS -- +-- -- +-- G N A T . S H A 1 -- +-- -- +-- B o d y -- +-- -- +-- Copyright (C) 2002-2006, AdaCore -- +-- -- +-- 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 2, 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. See the GNU General Public License -- +-- for more details. You should have received a copy of the GNU General -- +-- Public License distributed with GNAT; see file COPYING. If not, write -- +-- to the Free Software Foundation, 51 Franklin Street, Fifth Floor, -- +-- Boston, MA 02110-1301, USA. -- +-- -- +-- As a special exception, if other files instantiate generics from this -- +-- unit, or you link this unit with other files to produce an executable, -- +-- this unit does not by itself cause the resulting executable to be -- +-- covered by the GNU General Public License. This exception does not -- +-- however invalidate any other reasons why the executable file might be -- +-- covered by the GNU Public License. -- +-- -- +-- GNAT was originally developed by the GNAT team at New York University. -- +-- Extensive contributions were provided by Ada Core Technologies Inc. -- +-- -- +------------------------------------------------------------------------------ + +-- Note: the code for this unit is derived from GNAT.MD5 + +with Ada.Unchecked_Conversion; + +package body GNAT.SHA1 is + + use Interfaces; + + Padding : constant String := + (1 => Character'Val (16#80#), 2 .. 64 => ASCII.NUL); + + Hex_Digit : constant array (Unsigned_32 range 0 .. 15) of Character := + ('0', '1', '2', '3', '4', '5', '6', '7', + '8', '9', 'a', 'b', 'c', 'd', 'e', 'f'); + -- Look-up table for each hex digit of the Message-Digest. + -- Used by function Digest (Context). + + type Sixteen_Words is array (Natural range 0 .. 15) + of Interfaces.Unsigned_32; + -- Sixteen 32-bit words, converted from block of 64 characters. + -- Used in procedure Decode and Transform. + + procedure Decode (Block : String; X : out Sixteen_Words); + -- Convert a String of 64 characters into 16 32-bit numbers + + -- The following functions are the four elementary components of each + -- of the four round groups (0 .. 19, 20 .. 39, 40 .. 59, and 60 .. 79) + -- defined in RFC 3174. + + function F0 (B, C, D : Unsigned_32) return Unsigned_32; + pragma Inline (F0); + + function F1 (B, C, D : Unsigned_32) return Unsigned_32; + pragma Inline (F1); + + function F2 (B, C, D : Unsigned_32) return Unsigned_32; + pragma Inline (F2); + + function F3 (B, C, D : Unsigned_32) return Unsigned_32; + pragma Inline (F3); + + procedure Transform (Ctx : in out Context; Block : String); + -- Process one block of 64 characters + + ------------ + -- Decode -- + ------------ + + procedure Decode (Block : String; X : out Sixteen_Words) is + Cur : Positive := Block'First; + + begin + pragma Assert (Block'Length = 64); + + for Index in X'Range loop + X (Index) := + Unsigned_32 (Character'Pos (Block (Cur + 3))) + + Shift_Left (Unsigned_32 (Character'Pos (Block (Cur + 2))), 8) + + Shift_Left (Unsigned_32 (Character'Pos (Block (Cur + 1))), 16) + + Shift_Left (Unsigned_32 (Character'Pos (Block (Cur))), 24); + Cur := Cur + 4; + end loop; + end Decode; + + ------------ + -- Digest -- + ------------ + + function Digest (C : Context) return Message_Digest is + Result : Message_Digest; + + Cur : Natural := 1; + -- Index in Result where the next character will be placed + + Last_Block : String (1 .. 64); + + C1 : Context := C; + + procedure Convert (X : Unsigned_32); + -- Put the contribution of one of the five H words of the Context in + -- Result. Increments Cur. + + ------------- + -- Convert -- + ------------- + + procedure Convert (X : Unsigned_32) is + Y : Unsigned_32 := X; + begin + for J in 1 .. 8 loop + Y := Rotate_Left (Y, 4); + Result (Cur) := Hex_Digit (Y and Unsigned_32'(16#0F#)); + Cur := Cur + 1; + end loop; + end Convert; + + -- Start of processing for Digest + + begin + -- Process characters in the context buffer, if any + + pragma Assert (C.Last /= C.Buffer'Last); + Last_Block (1 .. C.Last) := C.Buffer (1 .. C.Last); + + if C.Last > 55 then + Last_Block (C.Last + 1 .. 64) := Padding (1 .. 64 - C.Last); + Transform (C1, Last_Block); + Last_Block := (others => ASCII.NUL); + + else + Last_Block (C.Last + 1 .. 56) := Padding (1 .. 56 - C.Last); + end if; + + -- Add the input length (as stored in the context) as 8 characters + + Last_Block (57 .. 64) := (others => ASCII.NUL); + + declare + L : Unsigned_64 := Unsigned_64 (C.Length) * 8; + Idx : Positive := 64; + begin + while L > 0 loop + Last_Block (Idx) := Character'Val (L and 16#Ff#); + L := Shift_Right (L, 8); + Idx := Idx - 1; + end loop; + end; + + Transform (C1, Last_Block); + + Convert (C1.H (0)); + Convert (C1.H (1)); + Convert (C1.H (2)); + Convert (C1.H (3)); + Convert (C1.H (4)); + return Result; + end Digest; + + function Digest (S : String) return Message_Digest is + C : Context; + begin + Update (C, S); + return Digest (C); + end Digest; + + function Digest + (A : Ada.Streams.Stream_Element_Array) return Message_Digest + is + C : Context; + begin + Update (C, A); + return Digest (C); + end Digest; + + -------- + -- F0 -- + -------- + + function F0 + (B, C, D : Interfaces.Unsigned_32) return Interfaces.Unsigned_32 + is + begin + return (B and C) or ((not B) and D); + end F0; + + -------- + -- F1 -- + -------- + + function F1 + (B, C, D : Interfaces.Unsigned_32) return Interfaces.Unsigned_32 + is + begin + return B xor C xor D; + end F1; + + -------- + -- F2 -- + -------- + + function F2 + (B, C, D : Interfaces.Unsigned_32) return Interfaces.Unsigned_32 + is + begin + return (B and C) or (B and D) or (C and D); + end F2; + + -------- + -- F3 -- + -------- + + function F3 + (B, C, D : Interfaces.Unsigned_32) return Interfaces.Unsigned_32 + renames F1; + + --------------- + -- Transform -- + --------------- + + procedure Transform + (Ctx : in out Context; + Block : String) + is + W : array (0 .. 79) of Interfaces.Unsigned_32; + + A, B, C, D, E, Temp : Interfaces.Unsigned_32; + + begin + pragma Assert (Block'Length = 64); + + -- a. Divide data block into sixteen words + + Decode (Block, Sixteen_Words (W (0 .. 15))); + + -- b. Prepare working block of 80 words + + for T in 16 .. 79 loop + + -- W(t) = S^1(W(t-3) XOR W(t-8) XOR W(t-14) XOR W(t-16)) + + W (T) := Rotate_Left + (W (T - 3) xor W (T - 8) xor W (T - 14) xor W (T - 16), 1); + + end loop; + + -- c. Set up transformation variables + + A := Ctx.H (0); + B := Ctx.H (1); + C := Ctx.H (2); + D := Ctx.H (3); + E := Ctx.H (4); + + -- d. For each of the 80 rounds, compute: + + -- TEMP = S^5(A) + f(t;B,C,D) + E + W(t) + K(t); + -- E = D; D = C; C = S^30(B); B = A; A = TEMP; + + for T in 0 .. 19 loop + Temp := Rotate_Left (A, 5) + F0 (B, C, D) + E + W (T) + 16#5A827999#; + E := D; D := C; C := Rotate_Left (B, 30); B := A; A := Temp; + end loop; + + for T in 20 .. 39 loop + Temp := Rotate_Left (A, 5) + F1 (B, C, D) + E + W (T) + 16#6ED9EBA1#; + E := D; D := C; C := Rotate_Left (B, 30); B := A; A := Temp; + end loop; + + for T in 40 .. 59 loop + Temp := Rotate_Left (A, 5) + F2 (B, C, D) + E + W (T) + 16#8F1BBCDC#; + E := D; D := C; C := Rotate_Left (B, 30); B := A; A := Temp; + end loop; + + for T in 60 .. 79 loop + Temp := Rotate_Left (A, 5) + F3 (B, C, D) + E + W (T) + 16#CA62C1D6#; + E := D; D := C; C := Rotate_Left (B, 30); B := A; A := Temp; + end loop; + + -- e. Update context: + -- H0 = H0 + A, H1 = H1 + B, H2 = H2 + C, H3 = H3 + D, H4 = H4 + E + + Ctx.H (0) := Ctx.H (0) + A; + Ctx.H (1) := Ctx.H (1) + B; + Ctx.H (2) := Ctx.H (2) + C; + Ctx.H (3) := Ctx.H (3) + D; + Ctx.H (4) := Ctx.H (4) + E; + end Transform; + + ------------ + -- Update -- + ------------ + + procedure Update + (C : in out Context; + Input : String) + is + Inp : constant String := C.Buffer (1 .. C.Last) & Input; + Cur : Positive := Inp'First; + + begin + C.Length := C.Length + Input'Length; + + while Cur + 63 <= Inp'Last loop + Transform (C, Inp (Cur .. Cur + 63)); + Cur := Cur + 64; + end loop; + + C.Last := Inp'Last - Cur + 1; + C.Buffer (1 .. C.Last) := Inp (Cur .. Inp'Last); + end Update; + + procedure Update + (C : in out Context; + Input : Ada.Streams.Stream_Element_Array) + is + subtype Stream_Array is Ada.Streams.Stream_Element_Array (Input'Range); + subtype Stream_String is + String (1 + Integer (Input'First) .. 1 + Integer (Input'Last)); + + function To_String is new Ada.Unchecked_Conversion + (Stream_Array, Stream_String); + + String_Input : constant String := To_String (Input); + begin + Update (C, String_Input); + end Update; + + ----------------- + -- Wide_Digest -- + ----------------- + + function Wide_Digest (W : Wide_String) return Message_Digest is + C : Context; + begin + Wide_Update (C, W); + return Digest (C); + end Wide_Digest; + + ----------------- + -- Wide_Update -- + ----------------- + + procedure Wide_Update + (C : in out Context; + Input : Wide_String) + is + String_Input : String (1 .. 2 * Input'Length); + Cur : Positive := 1; + + begin + for Index in Input'Range loop + String_Input (Cur) := + Character'Val + (Unsigned_32 (Wide_Character'Pos (Input (Index))) and 16#FF#); + Cur := Cur + 1; + String_Input (Cur) := + Character'Val + (Shift_Right (Unsigned_32 (Wide_Character'Pos (Input (Index))), 8) + and 16#FF#); + Cur := Cur + 1; + end loop; + + Update (C, String_Input); + end Wide_Update; + +end GNAT.SHA1; diff --git a/gcc/ada/g-sha1.ads b/gcc/ada/g-sha1.ads new file mode 100644 index 0000000..36e2e25 --- /dev/null +++ b/gcc/ada/g-sha1.ads @@ -0,0 +1,116 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT LIBRARY COMPONENTS -- +-- -- +-- G N A T . S H A 1 -- +-- -- +-- S p e c -- +-- -- +-- Copyright (C) 2002-2006, AdaCore -- +-- -- +-- 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 2, 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. See the GNU General Public License -- +-- for more details. You should have received a copy of the GNU General -- +-- Public License distributed with GNAT; see file COPYING. If not, write -- +-- to the Free Software Foundation, 51 Franklin Street, Fifth Floor, -- +-- Boston, MA 02110-1301, USA. -- +-- -- +-- As a special exception, if other files instantiate generics from this -- +-- unit, or you link this unit with other files to produce an executable, -- +-- this unit does not by itself cause the resulting executable to be -- +-- covered by the GNU General Public License. This exception does not -- +-- however invalidate any other reasons why the executable file might be -- +-- covered by the GNU Public License. -- +-- -- +-- GNAT was originally developed by the GNAT team at New York University. -- +-- Extensive contributions were provided by Ada Core Technologies Inc. -- +-- -- +------------------------------------------------------------------------------ + +-- This package implements the US Secure Hash Algorithm 1 (SHA1) as described +-- in RFC 3174. The complete text of RFC 3174 can be found at: + +-- http://www.ietf.org/rfc/rfc3174.txt + +-- Note: the code for this unit is derived from GNAT.MD5 + +with Ada.Streams; +with Interfaces; + +package GNAT.SHA1 is + + type Context is private; + -- This type holds the five-word (20 byte) buffer H, as described in + -- RFC 3174 (6.1). Its initial value is Initial_Context below. + + Initial_Context : constant Context; + -- Initial value of a Context object. May be used to reinitialize + -- a Context value by simple assignment of this value to the object. + + procedure Update + (C : in out Context; + Input : String); + procedure Wide_Update + (C : in out Context; + Input : Wide_String); + procedure Update + (C : in out Context; + Input : Ada.Streams.Stream_Element_Array); + -- Modify the Context C. If C has the initial value Initial_Context, + -- then, after a call to one of these procedures, Digest (C) will return + -- the Message-Digest of Input. + -- + -- These procedures may be called successively with the same context and + -- different inputs, and these several successive calls will produce + -- the same final context as a call with the concatenation of the inputs. + + subtype Message_Digest is String (1 .. 40); + -- The string type returned by function Digest + + function Digest (C : Context) return Message_Digest; + -- Extracts the Message-Digest from a context. This function should be + -- used after one or several calls to Update. + + function Digest (S : String) return Message_Digest; + function Wide_Digest (W : Wide_String) return Message_Digest; + function Digest + (A : Ada.Streams.Stream_Element_Array) return Message_Digest; + -- These functions are equivalent to the corresponding Update (or + -- Wide_Update) on a default initialized Context, followed by Digest + -- on the resulting Context. + +private + + -- Magic numbers + + Initial_H0 : constant := 16#67452301#; + Initial_H1 : constant := 16#EFCDAB89#; + Initial_H2 : constant := 16#98BADCFE#; + Initial_H3 : constant := 16#10325476#; + Initial_H4 : constant := 16#C3D2E1F0#; + + type H_Type is array (0 .. 4) of Interfaces.Unsigned_32; + + Initial_H : constant H_Type := + (0 => Initial_H0, + 1 => Initial_H1, + 2 => Initial_H2, + 3 => Initial_H3, + 4 => Initial_H4); + + type Context is record + H : H_Type := Initial_H; + Buffer : String (1 .. 64) := (others => ASCII.NUL); + Last : Natural := 0; + Length : Natural := 0; + end record; + + Initial_Context : constant Context := + (H => Initial_H, + Buffer => (others => ASCII.NUL), Last => 0, Length => 0); + +end GNAT.SHA1; diff --git a/gcc/ada/i-forbla.ads b/gcc/ada/i-forbla.ads new file mode 100644 index 0000000..db093dd --- /dev/null +++ b/gcc/ada/i-forbla.ads @@ -0,0 +1,251 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- INTERFACES.FORTRAN.BLAS -- +-- -- +-- S p e c -- +-- -- +-- Copyright (C) 2006, 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 2, 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. See the GNU General Public License -- +-- for more details. You should have received a copy of the GNU General -- +-- Public License distributed with GNAT; see file COPYING. If not, write -- +-- to the Free Software Foundation, 51 Franklin Street, Fifth Floor, -- +-- Boston, MA 02110-1301, USA. -- +-- -- +-- As a special exception, if other files instantiate generics from this -- +-- unit, or you link this unit with other files to produce an executable, -- +-- this unit does not by itself cause the resulting executable to be -- +-- covered by the GNU General Public License. This exception does not -- +-- however invalidate any other reasons why the executable file might be -- +-- covered by the GNU Public License. -- +-- -- +-- GNAT was originally developed by the GNAT team at New York University. -- +-- Extensive contributions were provided by Ada Core Technologies Inc. -- +-- -- +------------------------------------------------------------------------------ + +-- Comment required if non-RM package ??? + +package Interfaces.Fortran.BLAS is + pragma Pure; + + No_Trans : aliased constant Character := 'N'; + Trans : aliased constant Character := 'T'; + Conj_Trans : aliased constant Character := 'C'; + + -- Vector types + + type Real_Vector is array (Integer range <>) of Real; + + type Complex_Vector is array (Integer range <>) of Complex; + + type Double_Precision_Vector is array (Integer range <>) + of Double_Precision; + + type Double_Complex_Vector is array (Integer range <>) of Double_Complex; + + -- Matrix types + + type Real_Matrix is array (Integer range <>, Integer range <>) + of Real; + + type Double_Precision_Matrix is array (Integer range <>, Integer range <>) + of Double_Precision; + + type Complex_Matrix is array (Integer range <>, Integer range <>) + of Complex; + + type Double_Complex_Matrix is array (Integer range <>, Integer range <>) + of Double_Complex; + + -- BLAS Level 1 + + function sdot + (N : Positive; + X : Real_Vector; + Inc_X : Integer := 1; + Y : Real_Vector; + Inc_Y : Integer := 1) return Real; + + function ddot + (N : Positive; + X : Double_Precision_Vector; + Inc_X : Integer := 1; + Y : Double_Precision_Vector; + Inc_Y : Integer := 1) return Double_Precision; + + function cdot + (N : Positive; + X : Complex_Vector; + Inc_X : Integer := 1; + Y : Complex_Vector; + Inc_Y : Integer := 1) return Complex; + + function zdot + (N : Positive; + X : Double_Complex_Vector; + Inc_X : Integer := 1; + Y : Double_Complex_Vector; + Inc_Y : Integer := 1) return Double_Complex; + + function snrm2 + (N : Natural; + X : Real_Vector; + Inc_X : Integer := 1) return Real; + + function dnrm2 + (N : Natural; + X : Double_Precision_Vector; + Inc_X : Integer := 1) return Double_Precision; + + function scnrm2 + (N : Natural; + X : Complex_Vector; + Inc_X : Integer := 1) return Real; + + function dznrm2 + (N : Natural; + X : Double_Complex_Vector; + Inc_X : Integer := 1) return Double_Precision; + + -- BLAS Level 2 + + procedure sgemv + (Trans : access constant Character; + M : Natural := 0; + N : Natural := 0; + Alpha : Real := 1.0; + A : Real_Matrix; + Ld_A : Positive; + X : Real_Vector; + Inc_X : Integer := 1; -- must be non-zero + Beta : Real := 0.0; + Y : in out Real_Vector; + Inc_Y : Integer := 1); -- must be non-zero + + procedure dgemv + (Trans : access constant Character; + M : Natural := 0; + N : Natural := 0; + Alpha : Double_Precision := 1.0; + A : Double_Precision_Matrix; + Ld_A : Positive; + X : Double_Precision_Vector; + Inc_X : Integer := 1; -- must be non-zero + Beta : Double_Precision := 0.0; + Y : in out Double_Precision_Vector; + Inc_Y : Integer := 1); -- must be non-zero + + procedure cgemv + (Trans : access constant Character; + M : Natural := 0; + N : Natural := 0; + Alpha : Complex := (1.0, 1.0); + A : Complex_Matrix; + Ld_A : Positive; + X : Complex_Vector; + Inc_X : Integer := 1; -- must be non-zero + Beta : Complex := (0.0, 0.0); + Y : in out Complex_Vector; + Inc_Y : Integer := 1); -- must be non-zero + + procedure zgemv + (Trans : access constant Character; + M : Natural := 0; + N : Natural := 0; + Alpha : Double_Complex := (1.0, 1.0); + A : Double_Complex_Matrix; + Ld_A : Positive; + X : Double_Complex_Vector; + Inc_X : Integer := 1; -- must be non-zero + Beta : Double_Complex := (0.0, 0.0); + Y : in out Double_Complex_Vector; + Inc_Y : Integer := 1); -- must be non-zero + + -- BLAS Level 3 + + procedure sgemm + (Trans_A : access constant Character; + Trans_B : access constant Character; + M : Positive; + N : Positive; + K : Positive; + Alpha : Real := 1.0; + A : Real_Matrix; + Ld_A : Integer; + B : Real_Matrix; + Ld_B : Integer; + Beta : Real := 0.0; + C : in out Real_Matrix; + Ld_C : Integer); + + procedure dgemm + (Trans_A : access constant Character; + Trans_B : access constant Character; + M : Positive; + N : Positive; + K : Positive; + Alpha : Double_Precision := 1.0; + A : Double_Precision_Matrix; + Ld_A : Integer; + B : Double_Precision_Matrix; + Ld_B : Integer; + Beta : Double_Precision := 0.0; + C : in out Double_Precision_Matrix; + Ld_C : Integer); + + procedure cgemm + (Trans_A : access constant Character; + Trans_B : access constant Character; + M : Positive; + N : Positive; + K : Positive; + Alpha : Complex := (1.0, 1.0); + A : Complex_Matrix; + Ld_A : Integer; + B : Complex_Matrix; + Ld_B : Integer; + Beta : Complex := (0.0, 0.0); + C : in out Complex_Matrix; + Ld_C : Integer); + + procedure zgemm + (Trans_A : access constant Character; + Trans_B : access constant Character; + M : Positive; + N : Positive; + K : Positive; + Alpha : Double_Complex := (1.0, 1.0); + A : Double_Complex_Matrix; + Ld_A : Integer; + B : Double_Complex_Matrix; + Ld_B : Integer; + Beta : Double_Complex := (0.0, 0.0); + C : in out Double_Complex_Matrix; + Ld_C : Integer); + +private + pragma Import (Fortran, cdot, "cdot_"); + pragma Import (Fortran, cgemm, "cgemm_"); + pragma Import (Fortran, cgemv, "cgemv_"); + pragma Import (Fortran, ddot, "ddot_"); + pragma Import (Fortran, dgemm, "dgemm_"); + pragma Import (Fortran, dgemv, "dgemv_"); + pragma Import (Fortran, dnrm2, "dnrm2_"); + pragma Import (Fortran, dznrm2, "dznrm2_"); + pragma Import (Fortran, scnrm2, "scnrm2_"); + pragma Import (Fortran, sdot, "sdot_"); + pragma Import (Fortran, sgemm, "sgemm_"); + pragma Import (Fortran, sgemv, "sgemv_"); + pragma Import (Fortran, snrm2, "snrm2_"); + pragma Import (Fortran, zdot, "zdot_"); + pragma Import (Fortran, zgemm, "zgemm_"); + pragma Import (Fortran, zgemv, "zgemv_"); +end Interfaces.Fortran.BLAS; diff --git a/gcc/ada/i-forlap.ads b/gcc/ada/i-forlap.ads new file mode 100644 index 0000000..9faacce --- /dev/null +++ b/gcc/ada/i-forlap.ads @@ -0,0 +1,416 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- INTERFACES.FORTRAN.LAPACK -- +-- -- +-- S p e c -- +-- -- +-- Copyright (C) 2006, 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 2, 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. See the GNU General Public License -- +-- for more details. You should have received a copy of the GNU General -- +-- Public License distributed with GNAT; see file COPYING. If not, write -- +-- to the Free Software Foundation, 51 Franklin Street, Fifth Floor, -- +-- Boston, MA 02110-1301, USA. -- +-- -- +-- As a special exception, if other files instantiate generics from this -- +-- unit, or you link this unit with other files to produce an executable, -- +-- this unit does not by itself cause the resulting executable to be -- +-- covered by the GNU General Public License. This exception does not -- +-- however invalidate any other reasons why the executable file might be -- +-- covered by the GNU Public License. -- +-- -- +-- GNAT was originally developed by the GNAT team at New York University. -- +-- Extensive contributions were provided by Ada Core Technologies Inc. -- +-- -- +------------------------------------------------------------------------------ + +-- Package comment required if non-RM package ??? + +with Interfaces.Fortran.BLAS; +package Interfaces.Fortran.LAPACK is + pragma Pure; + + type Integer_Vector is array (Integer range <>) of Integer; + + Upper : aliased constant Character := 'U'; + Lower : aliased constant Character := 'L'; + + subtype Real_Vector is BLAS.Real_Vector; + subtype Real_Matrix is BLAS.Real_Matrix; + subtype Double_Precision_Vector is BLAS.Double_Precision_Vector; + subtype Double_Precision_Matrix is BLAS.Double_Precision_Matrix; + subtype Complex_Vector is BLAS.Complex_Vector; + subtype Complex_Matrix is BLAS.Complex_Matrix; + subtype Double_Complex_Vector is BLAS.Double_Complex_Vector; + subtype Double_Complex_Matrix is BLAS.Double_Complex_Matrix; + + -- LAPACK Computational Routines + + -- gerfs Refines the solution of a system of linear equations with + -- a general matrix and estimates its error + -- getrf Computes LU factorization of a general m-by-n matrix + -- getri Computes inverse of an LU-factored general matrix + -- square matrix, with multiple right-hand sides + -- getrs Solves a system of linear equations with an LU-factored + -- square matrix, with multiple right-hand sides + -- hetrd Reduces a complex Hermitian matrix to tridiagonal form + -- heevr Computes selected eigenvalues and, optionally, eigenvectors of + -- a Hermitian matrix using the Relatively Robust Representations + -- orgtr Generates the real orthogonal matrix Q determined by sytrd + -- steqr Computes all eigenvalues and eigenvectors of a symmetric or + -- Hermitian matrix reduced to tridiagonal form (QR algorithm) + -- sterf Computes all eigenvalues of a real symmetric + -- tridiagonal matrix using QR algorithm + -- sytrd Reduces a real symmetric matrix to tridiagonal form + + procedure sgetrf + (M : Natural; + N : Natural; + A : in out Real_Matrix; + Ld_A : Positive; + I_Piv : out Integer_Vector; + Info : access Integer); + + procedure dgetrf + (M : Natural; + N : Natural; + A : in out Double_Precision_Matrix; + Ld_A : Positive; + I_Piv : out Integer_Vector; + Info : access Integer); + + procedure cgetrf + (M : Natural; + N : Natural; + A : in out Complex_Matrix; + Ld_A : Positive; + I_Piv : out Integer_Vector; + Info : access Integer); + + procedure zgetrf + (M : Natural; + N : Natural; + A : in out Double_Complex_Matrix; + Ld_A : Positive; + I_Piv : out Integer_Vector; + Info : access Integer); + + procedure sgetri + (N : Natural; + A : in out Real_Matrix; + Ld_A : Positive; + I_Piv : Integer_Vector; + Work : in out Real_Vector; + L_Work : Integer; + Info : access Integer); + + procedure dgetri + (N : Natural; + A : in out Double_Precision_Matrix; + Ld_A : Positive; + I_Piv : Integer_Vector; + Work : in out Double_Precision_Vector; + L_Work : Integer; + Info : access Integer); + + procedure cgetri + (N : Natural; + A : in out Complex_Matrix; + Ld_A : Positive; + I_Piv : Integer_Vector; + Work : in out Complex_Vector; + L_Work : Integer; + Info : access Integer); + + procedure zgetri + (N : Natural; + A : in out Double_Complex_Matrix; + Ld_A : Positive; + I_Piv : Integer_Vector; + Work : in out Double_Complex_Vector; + L_Work : Integer; + Info : access Integer); + + procedure sgetrs + (Trans : access constant Character; + N : Natural; + N_Rhs : Natural; + A : Real_Matrix; + Ld_A : Positive; + I_Piv : Integer_Vector; + B : in out Real_Matrix; + Ld_B : Positive; + Info : access Integer); + + procedure dgetrs + (Trans : access constant Character; + N : Natural; + N_Rhs : Natural; + A : Double_Precision_Matrix; + Ld_A : Positive; + I_Piv : Integer_Vector; + B : in out Double_Precision_Matrix; + Ld_B : Positive; + Info : access Integer); + + procedure cgetrs + (Trans : access constant Character; + N : Natural; + N_Rhs : Natural; + A : Complex_Matrix; + Ld_A : Positive; + I_Piv : Integer_Vector; + B : in out Complex_Matrix; + Ld_B : Positive; + Info : access Integer); + + procedure zgetrs + (Trans : access constant Character; + N : Natural; + N_Rhs : Natural; + A : Double_Complex_Matrix; + Ld_A : Positive; + I_Piv : Integer_Vector; + B : in out Double_Complex_Matrix; + Ld_B : Positive; + Info : access Integer); + + procedure cheevr + (Job_Z : access constant Character; + Rng : access constant Character; + Uplo : access constant Character; + N : Natural; + A : in out Complex_Matrix; + Ld_A : Positive; + Vl, Vu : Real := 0.0; + Il, Iu : Integer := 1; + Abs_Tol : Real := 0.0; + M : out Integer; + W : out Real_Vector; + Z : out Complex_Matrix; + Ld_Z : Positive; + I_Supp_Z : out Integer_Vector; + Work : out Complex_Vector; + L_Work : Integer; + R_Work : out Real_Vector; + LR_Work : Integer; + I_Work : out Integer_Vector; + LI_Work : Integer; + Info : access Integer); + + procedure zheevr + (Job_Z : access constant Character; + Rng : access constant Character; + Uplo : access constant Character; + N : Natural; + A : in out Double_Complex_Matrix; + Ld_A : Positive; + Vl, Vu : Double_Precision := 0.0; + Il, Iu : Integer := 1; + Abs_Tol : Double_Precision := 0.0; + M : out Integer; + W : out Double_Precision_Vector; + Z : out Double_Complex_Matrix; + Ld_Z : Positive; + I_Supp_Z : out Integer_Vector; + Work : out Double_Complex_Vector; + L_Work : Integer; + R_Work : out Double_Precision_Vector; + LR_Work : Integer; + I_Work : out Integer_Vector; + LI_Work : Integer; + Info : access Integer); + + procedure chetrd + (Uplo : access constant Character; + N : Natural; + A : in out Complex_Matrix; + Ld_A : Positive; + D : out Real_Vector; + E : out Real_Vector; + Tau : out Complex_Vector; + Work : out Complex_Vector; + L_Work : Integer; + Info : access Integer); + + procedure zhetrd + (Uplo : access constant Character; + N : Natural; + A : in out Double_Complex_Matrix; + Ld_A : Positive; + D : out Double_Precision_Vector; + E : out Double_Precision_Vector; + Tau : out Double_Complex_Vector; + Work : out Double_Complex_Vector; + L_Work : Integer; + Info : access Integer); + + procedure ssytrd + (Uplo : access constant Character; + N : Natural; + A : in out Real_Matrix; + Ld_A : Positive; + D : out Real_Vector; + E : out Real_Vector; + Tau : out Real_Vector; + Work : out Real_Vector; + L_Work : Integer; + Info : access Integer); + + procedure dsytrd + (Uplo : access constant Character; + N : Natural; + A : in out Double_Precision_Matrix; + Ld_A : Positive; + D : out Double_Precision_Vector; + E : out Double_Precision_Vector; + Tau : out Double_Precision_Vector; + Work : out Double_Precision_Vector; + L_Work : Integer; + Info : access Integer); + + procedure ssterf + (N : Natural; + D : in out Real_Vector; + E : in out Real_Vector; + Info : access Integer); + + procedure dsterf + (N : Natural; + D : in out Double_Precision_Vector; + E : in out Double_Precision_Vector; + Info : access Integer); + + procedure sorgtr + (Uplo : access constant Character; + N : Natural; + A : in out Real_Matrix; + Ld_A : Positive; + Tau : in Real_Vector; + Work : out Real_Vector; + L_Work : Integer; + Info : access Integer); + + procedure dorgtr + (Uplo : access constant Character; + N : Natural; + A : in out Double_Precision_Matrix; + Ld_A : Positive; + Tau : in Double_Precision_Vector; + Work : out Double_Precision_Vector; + L_Work : Integer; + Info : access Integer); + + procedure sstebz + (Rng : access constant Character; + Order : access constant Character; + N : in Natural; + Vl, Vu : in Real := 0.0; + Il, Iu : in Integer := 1; + Abs_Tol : in Real := 0.0; + D : in Real_Vector; + E : in Real_Vector; + M : out Natural; + N_Split : out Natural; + W : out Real_Vector; + I_Block : out Integer_Vector; + I_Split : out Integer_Vector; + Work : out Real_Vector; + I_Work : out Integer_Vector; + Info : access Integer); + + procedure dstebz + (Rng : access constant Character; + Order : access constant Character; + N : in Natural; + Vl, Vu : in Double_Precision := 0.0; + Il, Iu : in Integer := 1; + Abs_Tol : in Double_Precision := 0.0; + D : in Double_Precision_Vector; + E : in Double_Precision_Vector; + M : out Natural; + N_Split : out Natural; + W : out Double_Precision_Vector; + I_Block : out Integer_Vector; + I_Split : out Integer_Vector; + Work : out Double_Precision_Vector; + I_Work : out Integer_Vector; + Info : access Integer); + + procedure ssteqr + (Comp_Z : access constant Character; + N : Natural; + D : in out Real_Vector; + E : in out Real_Vector; + Z : in out Real_Matrix; + Ld_Z : Positive; + Work : out Real_Vector; + Info : access Integer); + + procedure dsteqr + (Comp_Z : access constant Character; + N : Natural; + D : in out Double_Precision_Vector; + E : in out Double_Precision_Vector; + Z : in out Double_Precision_Matrix; + Ld_Z : Positive; + Work : out Double_Precision_Vector; + Info : access Integer); + + procedure csteqr + (Comp_Z : access constant Character; + N : Natural; + D : in out Real_Vector; + E : in out Real_Vector; + Z : in out Complex_Matrix; + Ld_Z : Positive; + Work : out Real_Vector; + Info : access Integer); + + procedure zsteqr + (Comp_Z : access constant Character; + N : Natural; + D : in out Double_Precision_Vector; + E : in out Double_Precision_Vector; + Z : in out Double_Complex_Matrix; + Ld_Z : Positive; + Work : out Double_Precision_Vector; + Info : access Integer); + +private + pragma Import (Fortran, csteqr, "csteqr_"); + pragma Import (Fortran, cgetrf, "cgetrf_"); + pragma Import (Fortran, cgetri, "cgetri_"); + pragma Import (Fortran, cgetrs, "cgetrs_"); + pragma Import (Fortran, cheevr, "cheevr_"); + pragma Import (Fortran, chetrd, "chetrd_"); + pragma Import (Fortran, dgetrf, "dgetrf_"); + pragma Import (Fortran, dgetri, "dgetri_"); + pragma Import (Fortran, dgetrs, "dgetrs_"); + pragma Import (Fortran, dsytrd, "dsytrd_"); + pragma Import (Fortran, dstebz, "dstebz_"); + pragma Import (Fortran, dsterf, "dsterf_"); + pragma Import (Fortran, dorgtr, "dorgtr_"); + pragma Import (Fortran, dsteqr, "dsteqr_"); + pragma Import (Fortran, sgetrf, "sgetrf_"); + pragma Import (Fortran, sgetri, "sgetri_"); + pragma Import (Fortran, sgetrs, "sgetrs_"); + pragma Import (Fortran, sorgtr, "sorgtr_"); + pragma Import (Fortran, sstebz, "sstebz_"); + pragma Import (Fortran, ssterf, "ssterf_"); + pragma Import (Fortran, ssteqr, "ssteqr_"); + pragma Import (Fortran, ssytrd, "ssytrd_"); + pragma Import (Fortran, zgetrf, "zgetrf_"); + pragma Import (Fortran, zgetri, "zgetri_"); + pragma Import (Fortran, zgetrs, "zgetrs_"); + pragma Import (Fortran, zheevr, "zheevr_"); + pragma Import (Fortran, zhetrd, "zhetrd_"); + pragma Import (Fortran, zsteqr, "zsteqr_"); +end Interfaces.Fortran.LAPACK; diff --git a/gcc/ada/i-fortra.ads b/gcc/ada/i-fortra.ads index da330d2..992eb28 100644 --- a/gcc/ada/i-fortra.ads +++ b/gcc/ada/i-fortra.ads @@ -6,7 +6,7 @@ -- -- -- S p e c -- -- -- --- This specification is adapted from the Ada Reference Manual for use with -- +-- This specification is derived from the Ada Reference Manual for use with -- -- GNAT. In accordance with the copyright of that document, you can freely -- -- copy and modify this specification, provided that if you redistribute a -- -- modified version, any changes that you have made are clearly indicated. -- @@ -35,8 +35,13 @@ package Interfaces.Fortran is package Single_Precision_Complex_Types is new Ada.Numerics.Generic_Complex_Types (Real); + package Double_Precision_Complex_Types is + new Ada.Numerics.Generic_Complex_Types (Double_Precision); + type Complex is new Single_Precision_Complex_Types.Complex; + type Double_Complex is new Double_Precision_Complex_Types.Complex; + subtype Imaginary is Single_Precision_Complex_Types.Imaginary; i : Imaginary renames Single_Precision_Complex_Types.i; j : Imaginary renames Single_Precision_Complex_Types.j; diff --git a/gcc/ada/impunit.adb b/gcc/ada/impunit.adb index e849797..ff5e88b 100644 --- a/gcc/ada/impunit.adb +++ b/gcc/ada/impunit.adb @@ -24,8 +24,12 @@ -- -- ------------------------------------------------------------------------------ -with Lib; use Lib; -with Namet; use Namet; +with Atree; use Atree; +with Sinfo; use Sinfo; +with Fname.UF; use Fname.UF; +with Lib; use Lib; +with Namet; use Namet; +with Uname; use Uname; package body Impunit is @@ -207,6 +211,7 @@ package body Impunit is "g-bubsor", -- GNAT.Bubble_Sort "g-busora", -- GNAT.Bubble_Sort_A "g-busorg", -- GNAT.Bubble_Sort_G + "g-bytswa", -- Gnat.Byte_Swapping "g-calend", -- GNAT.Calendar "g-casuti", -- GNAT.Case_Util "g-catiio", -- GNAT.Calendar.Time_IO @@ -246,6 +251,7 @@ package body Impunit is "g-regpat", -- GNAT.Regpat "g-semaph", -- GNAT.Semaphores "g-sestin", -- GNAT.Secondary_Stack_Info + "g-sha1 ", -- GNAT.SHA1 "g-signal", -- GNAT.Signals "g-socket", -- GNAT.Sockets "g-souinf", -- GNAT.Source_Info @@ -359,6 +365,10 @@ package body Impunit is "a-dispat", -- Ada.Dispatching "a-envvar", -- Ada.Environment_Variables "a-rttiev", -- Ada.Real_Time.Timing_Events + "a-ngcoar", -- Ada.Numerics.Generic_Complex_Arrays + "a-ngrear", -- Ada.Numerics.Generic_Real_Arrays + "a-nucoar", -- Ada.Numerics.Complex_Arrays + "a-nurear", -- Ada.Numerics.Real_Arrays "a-stboha", -- Ada.Strings.Bounded.Hash "a-stfiha", -- Ada.Strings.Fixed.Hash "a-strhas", -- Ada.Strings.Hash @@ -401,6 +411,10 @@ package body Impunit is "a-llctio", -- Ada.Long_Long_Complex_Text_IO "a-llfzti", -- Ada.Long_Long_Float_Wide_Wide_Text_IO "a-llizti", -- Ada.Long_Long_Integer_Wide_Wide_Text_IO + "a-nlcoar", -- Ada.Numerics.Long_Complex_Arrays + "a-nllcar", -- Ada.Numerics.Long_Long_Complex_Arrays + "a-nllrar", -- Ada.Numerics.Long_Long_Real_Arrays + "a-nlrear", -- Ada.Numerics.Long_Real_Arrays "a-scteio", -- Ada.Short_Complex_Text_IO "a-sfztio", -- Ada.Short_Float_Wide_Wide_Text_IO "a-siztio", -- Ada.Short_Integer_Wide_Wide_Text_IO @@ -536,4 +550,75 @@ package body Impunit is return Implementation_Unit; end Get_Kind_Of_Unit; + ------------------- + -- Is_Known_Unit -- + ------------------- + + function Is_Known_Unit (Nam : Node_Id) return Boolean is + Unam : Unit_Name_Type; + Fnam : File_Name_Type; + + begin + -- If selector is not an identifier (e.g. it is a character literal or + -- some junk from a previous error), then definitely not a known unit. + + if Nkind (Selector_Name (Nam)) /= N_Identifier then + return False; + end if; + + -- Otherwise get corresponding file name + + Unam := Get_Unit_Name (Nam); + Fnam := Get_File_Name (Unam, Subunit => False); + Get_Name_String (Fnam); + + -- Remove extension from file name + + if Name_Buffer (Name_Len - 3 .. Name_Len) = ".adb" then + Name_Len := Name_Len - 4; + else + return False; + end if; + + -- Pad name to 8 characters + + while Name_Len < 8 loop + Name_Len := Name_Len + 1; + Name_Buffer (Name_Len) := ' '; + end loop; + + -- If length more than 8, definitely not a match + + if Name_Len /= 8 then + return False; + end if; + + -- If length is 8, search our tables + + for J in Non_Imp_File_Names_95'Range loop + if Name_Buffer (1 .. 8) = Non_Imp_File_Names_95 (J) then + return True; + end if; + end loop; + + for J in Non_Imp_File_Names_05'Range loop + if Name_Buffer (1 .. 8) = Non_Imp_File_Names_05 (J) then + return True; + end if; + end loop; + + -- If not found, not known + + return False; + + -- A safety guard, if we get an exception during this processing then it + -- is most likely the result of a previous error, or a peculiar case we + -- have not thought of. Since this routine is only used for error message + -- refinement, we will just return False. + + exception + when others => + return False; + end Is_Known_Unit; + end Impunit; diff --git a/gcc/ada/impunit.ads b/gcc/ada/impunit.ads index f14ac7c..f2ffc4e 100644 --- a/gcc/ada/impunit.ads +++ b/gcc/ada/impunit.ads @@ -6,7 +6,7 @@ -- -- -- S p e c -- -- -- --- Copyright (C) 2000-2005, Free Software Foundation, Inc. -- +-- Copyright (C) 2000-2006, 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- -- @@ -58,4 +58,10 @@ package Impunit is -- Given the unit number of a unit, this function determines the type -- of the unit, as defined above. + function Is_Known_Unit (Nam : Node_Id) return Boolean; + -- Nam is the possible name of a child unit, represented as a selected + -- component node. This function determines whether the name matches + -- one of the known library units, and if so, returns True. If the name + -- does not match any known library unit, False is returned. + end Impunit; diff --git a/gcc/ada/s-gearop.adb b/gcc/ada/s-gearop.adb new file mode 100644 index 0000000..b6a3c79 --- /dev/null +++ b/gcc/ada/s-gearop.adb @@ -0,0 +1,518 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- SYSTEM.GENERIC_ARRAY_OPERATIONS -- +-- -- +-- B o d y -- +-- -- +-- Copyright (C) 2006, 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 2, 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. See the GNU General Public License -- +-- for more details. You should have received a copy of the GNU General -- +-- Public License distributed with GNAT; see file COPYING. If not, write -- +-- to the Free Software Foundation, 51 Franklin Street, Fifth Floor, -- +-- Boston, MA 02110-1301, USA. -- +-- -- +-- As a special exception, if other files instantiate generics from this -- +-- unit, or you link this unit with other files to produce an executable, -- +-- this unit does not by itself cause the resulting executable to be -- +-- covered by the GNU General Public License. This exception does not -- +-- however invalidate any other reasons why the executable file might be -- +-- covered by the GNU Public License. -- +-- -- +-- 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, K)); + 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, K), 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)); + 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), 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 (J - R'First (1) + A'First (1), + K - R'First (2) + 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; diff --git a/gcc/ada/s-gearop.ads b/gcc/ada/s-gearop.ads new file mode 100644 index 0000000..b922871 --- /dev/null +++ b/gcc/ada/s-gearop.ads @@ -0,0 +1,398 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- SYSTEM.GENERIC_ARRAY_OPERATIONS -- +-- -- +-- S p e c -- +-- -- +-- Copyright (C) 2006, 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 2, 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. See the GNU General Public License -- +-- for more details. You should have received a copy of the GNU General -- +-- Public License distributed with GNAT; see file COPYING. If not, write -- +-- to the Free Software Foundation, 51 Franklin Street, Fifth Floor, -- +-- Boston, MA 02110-1301, USA. -- +-- -- +-- As a special exception, if other files instantiate generics from this -- +-- unit, or you link this unit with other files to produce an executable, -- +-- this unit does not by itself cause the resulting executable to be -- +-- covered by the GNU General Public License. This exception does not -- +-- however invalidate any other reasons why the executable file might be -- +-- covered by the GNU Public License. -- +-- -- +-- GNAT was originally developed by the GNAT team at New York University. -- +-- Extensive contributions were provided by Ada Core Technologies Inc. -- +-- -- +------------------------------------------------------------------------------ + +package System.Generic_Array_Operations is +pragma Pure (Generic_Array_Operations); + + -------------------------- + -- Square_Matrix_Length -- + -------------------------- + + generic + type Scalar is private; + type Matrix is array (Integer range <>, Integer range <>) of Scalar; + function Square_Matrix_Length (A : Matrix) return Natural; + -- If A is non-square, raise Constraint_Error, else return its dimension + + ---------------------------------- + -- Vector_Elementwise_Operation -- + ---------------------------------- + + generic + type X_Scalar is private; + type Result_Scalar is private; + type X_Vector is array (Integer range <>) of X_Scalar; + type Result_Vector is array (Integer range <>) of Result_Scalar; + with function Operation (X : X_Scalar) return Result_Scalar; + function Vector_Elementwise_Operation (X : X_Vector) return Result_Vector; + + ---------------------------------- + -- Matrix_Elementwise_Operation -- + ---------------------------------- + + generic + type X_Scalar is private; + type Result_Scalar is private; + type X_Matrix is array (Integer range <>, Integer range <>) of X_Scalar; + type Result_Matrix is array (Integer range <>, Integer range <>) + of Result_Scalar; + with function Operation (X : X_Scalar) return Result_Scalar; + function Matrix_Elementwise_Operation (X : X_Matrix) return Result_Matrix; + + ----------------------------------------- + -- Vector_Vector_Elementwise_Operation -- + ----------------------------------------- + + generic + type Left_Scalar is private; + type Right_Scalar is private; + type Result_Scalar is private; + type Left_Vector is array (Integer range <>) of Left_Scalar; + type Right_Vector is array (Integer range <>) of Right_Scalar; + type Result_Vector is array (Integer range <>) of Result_Scalar; + with function Operation + (Left : Left_Scalar; + Right : Right_Scalar) return Result_Scalar; + function Vector_Vector_Elementwise_Operation + (Left : Left_Vector; + Right : Right_Vector) return Result_Vector; + + ------------------------------------------------ + -- Vector_Vector_Scalar_Elementwise_Operation -- + ------------------------------------------------ + + generic + type X_Scalar is private; + type Y_Scalar is private; + type Z_Scalar is private; + type Result_Scalar is private; + type X_Vector is array (Integer range <>) of X_Scalar; + type Y_Vector is array (Integer range <>) of Y_Scalar; + type Result_Vector is array (Integer range <>) of Result_Scalar; + with function Operation + (X : X_Scalar; + Y : Y_Scalar; + Z : Z_Scalar) return Result_Scalar; + function Vector_Vector_Scalar_Elementwise_Operation + (X : X_Vector; + Y : Y_Vector; + Z : Z_Scalar) return Result_Vector; + + ----------------------------------------- + -- Matrix_Matrix_Elementwise_Operation -- + ----------------------------------------- + + generic + type Left_Scalar is private; + type Right_Scalar is private; + type Result_Scalar is private; + type Left_Matrix is array (Integer range <>, Integer range <>) + of Left_Scalar; + type Right_Matrix is array (Integer range <>, Integer range <>) + of Right_Scalar; + type Result_Matrix is array (Integer range <>, Integer range <>) + of Result_Scalar; + with function Operation + (Left : Left_Scalar; + Right : Right_Scalar) return Result_Scalar; + function Matrix_Matrix_Elementwise_Operation + (Left : Left_Matrix; + Right : Right_Matrix) return Result_Matrix; + + ------------------------------------------------ + -- Matrix_Matrix_Scalar_Elementwise_Operation -- + ------------------------------------------------ + + generic + type X_Scalar is private; + type Y_Scalar is private; + type Z_Scalar is private; + type Result_Scalar is private; + type X_Matrix is array (Integer range <>, Integer range <>) of X_Scalar; + type Y_Matrix is array (Integer range <>, Integer range <>) of Y_Scalar; + type Result_Matrix is array (Integer range <>, Integer range <>) + of Result_Scalar; + with function Operation + (X : X_Scalar; + Y : Y_Scalar; + Z : Z_Scalar) return Result_Scalar; + function Matrix_Matrix_Scalar_Elementwise_Operation + (X : X_Matrix; + Y : Y_Matrix; + Z : Z_Scalar) return Result_Matrix; + + ----------------------------------------- + -- Vector_Scalar_Elementwise_Operation -- + ----------------------------------------- + + generic + type Left_Scalar is private; + type Right_Scalar is private; + type Result_Scalar is private; + type Left_Vector is array (Integer range <>) of Left_Scalar; + type Result_Vector is array (Integer range <>) of Result_Scalar; + with function Operation + (Left : Left_Scalar; + Right : Right_Scalar) return Result_Scalar; + function Vector_Scalar_Elementwise_Operation + (Left : Left_Vector; + Right : Right_Scalar) return Result_Vector; + + ----------------------------------------- + -- Matrix_Scalar_Elementwise_Operation -- + ----------------------------------------- + + generic + type Left_Scalar is private; + type Right_Scalar is private; + type Result_Scalar is private; + type Left_Matrix is array (Integer range <>, Integer range <>) + of Left_Scalar; + type Result_Matrix is array (Integer range <>, Integer range <>) + of Result_Scalar; + with function Operation + (Left : Left_Scalar; + Right : Right_Scalar) return Result_Scalar; + function Matrix_Scalar_Elementwise_Operation + (Left : Left_Matrix; + Right : Right_Scalar) return Result_Matrix; + + ----------------------------------------- + -- Scalar_Vector_Elementwise_Operation -- + ----------------------------------------- + + generic + type Left_Scalar is private; + type Right_Scalar is private; + type Result_Scalar is private; + type Right_Vector is array (Integer range <>) of Right_Scalar; + type Result_Vector is array (Integer range <>) of Result_Scalar; + with function Operation + (Left : Left_Scalar; + Right : Right_Scalar) return Result_Scalar; + function Scalar_Vector_Elementwise_Operation + (Left : Left_Scalar; + Right : Right_Vector) return Result_Vector; + + ----------------------------------------- + -- Scalar_Matrix_Elementwise_Operation -- + ----------------------------------------- + + generic + type Left_Scalar is private; + type Right_Scalar is private; + type Result_Scalar is private; + type Right_Matrix is array (Integer range <>, Integer range <>) + of Right_Scalar; + type Result_Matrix is array (Integer range <>, Integer range <>) + of Result_Scalar; + with function Operation + (Left : Left_Scalar; + Right : Right_Scalar) return Result_Scalar; + function Scalar_Matrix_Elementwise_Operation + (Left : Left_Scalar; + Right : Right_Matrix) return Result_Matrix; + + ------------------- + -- Inner_Product -- + ------------------- + + generic + type Left_Scalar is private; + type Right_Scalar is private; + type Result_Scalar is private; + type Left_Vector is array (Integer range <>) of Left_Scalar; + type Right_Vector is array (Integer range <>) of Right_Scalar; + Zero : Result_Scalar; + with function "*" + (Left : Left_Scalar; + Right : Right_Scalar) return Result_Scalar is <>; + with function "+" + (Left : Result_Scalar; + Right : Result_Scalar) return Result_Scalar is <>; + function Inner_Product + (Left : Left_Vector; + Right : Right_Vector) return Result_Scalar; + + ------------------- + -- Outer_Product -- + ------------------- + + generic + type Left_Scalar is private; + type Right_Scalar is private; + type Result_Scalar is private; + type Left_Vector is array (Integer range <>) of Left_Scalar; + type Right_Vector is array (Integer range <>) of Right_Scalar; + type Matrix is array (Integer range <>, Integer range <>) + of Result_Scalar; + with function "*" + (Left : Left_Scalar; + Right : Right_Scalar) return Result_Scalar is <>; + function Outer_Product + (Left : Left_Vector; + Right : Right_Vector) return Matrix; + + --------------------------- + -- Matrix_Vector_Product -- + --------------------------- + + generic + type Left_Scalar is private; + type Right_Scalar is private; + type Result_Scalar is private; + type Matrix is array (Integer range <>, Integer range <>) + of Left_Scalar; + type Right_Vector is array (Integer range <>) of Right_Scalar; + type Result_Vector is array (Integer range <>) of Result_Scalar; + Zero : Result_Scalar; + with function "*" + (Left : Left_Scalar; + Right : Right_Scalar) return Result_Scalar is <>; + with function "+" + (Left : Result_Scalar; + Right : Result_Scalar) return Result_Scalar is <>; + function Matrix_Vector_Product + (Left : Matrix; + Right : Right_Vector) return Result_Vector; + + --------------------------- + -- Vector_Matrix_Product -- + --------------------------- + + generic + type Left_Scalar is private; + type Right_Scalar is private; + type Result_Scalar is private; + type Left_Vector is array (Integer range <>) of Left_Scalar; + type Matrix is array (Integer range <>, Integer range <>) + of Right_Scalar; + type Result_Vector is array (Integer range <>) of Result_Scalar; + Zero : Result_Scalar; + with function "*" + (Left : Left_Scalar; + Right : Right_Scalar) return Result_Scalar is <>; + with function "+" + (Left : Result_Scalar; + Right : Result_Scalar) return Result_Scalar is <>; + function Vector_Matrix_Product + (Left : Left_Vector; + Right : Matrix) return Result_Vector; + + --------------------------- + -- Matrix_Matrix_Product -- + --------------------------- + + generic + type Left_Scalar is private; + type Right_Scalar is private; + type Result_Scalar is private; + type Left_Matrix is array (Integer range <>, Integer range <>) + of Left_Scalar; + type Right_Matrix is array (Integer range <>, Integer range <>) + of Right_Scalar; + type Result_Matrix is array (Integer range <>, Integer range <>) + of Result_Scalar; + Zero : Result_Scalar; + with function "*" + (Left : Left_Scalar; + Right : Right_Scalar) return Result_Scalar is <>; + with function "+" + (Left : Result_Scalar; + Right : Result_Scalar) return Result_Scalar is <>; + function Matrix_Matrix_Product + (Left : Left_Matrix; + Right : Right_Matrix) return Result_Matrix; + + --------------- + -- Transpose -- + --------------- + + generic + type Scalar is private; + type Matrix is array (Integer range <>, Integer range <>) of Scalar; + procedure Transpose (A : Matrix; R : out Matrix); + + ------------------------------- + -- Update_Vector_With_Vector -- + ------------------------------- + + generic + type X_Scalar is private; + type Y_Scalar is private; + type X_Vector is array (Integer range <>) of X_Scalar; + type Y_Vector is array (Integer range <>) of Y_Scalar; + with procedure Update (X : in out X_Scalar; Y : Y_Scalar); + procedure Update_Vector_With_Vector (X : in out X_Vector; Y : Y_Vector); + + ------------------------------- + -- Update_Matrix_With_Matrix -- + ------------------------------- + + generic + type X_Scalar is private; + type Y_Scalar is private; + type X_Matrix is array (Integer range <>, Integer range <>) of X_Scalar; + type Y_Matrix is array (Integer range <>, Integer range <>) of Y_Scalar; + with procedure Update (X : in out X_Scalar; Y : Y_Scalar); + procedure Update_Matrix_With_Matrix (X : in out X_Matrix; Y : Y_Matrix); + + ----------------- + -- Unit_Matrix -- + ----------------- + + generic + type Scalar is private; + type Matrix is array (Integer range <>, Integer range <>) of Scalar; + Zero : Scalar; + One : Scalar; + function Unit_Matrix + (Order : Positive; + First_1 : Integer := 1; + First_2 : Integer := 1) return Matrix; + + ----------------- + -- Unit_Vector -- + ----------------- + + generic + type Scalar is private; + type Vector is array (Integer range <>) of Scalar; + Zero : Scalar; + One : Scalar; + function Unit_Vector + (Index : Integer; + Order : Positive; + First : Integer := 1) return Vector; + +end System.Generic_Array_Operations; diff --git a/gcc/ada/s-gecobl.adb b/gcc/ada/s-gecobl.adb new file mode 100644 index 0000000..db657c3 --- /dev/null +++ b/gcc/ada/s-gecobl.adb @@ -0,0 +1,352 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- SYSTEM.GENERIC_COMPLEX_BLAS -- +-- -- +-- B o d y -- +-- -- +-- Copyright (C) 2006, 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 2, 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. See the GNU General Public License -- +-- for more details. You should have received a copy of the GNU General -- +-- Public License distributed with GNAT; see file COPYING. If not, write -- +-- to the Free Software Foundation, 51 Franklin Street, Fifth Floor, -- +-- Boston, MA 02110-1301, USA. -- +-- -- +-- As a special exception, if other files instantiate generics from this -- +-- unit, or you link this unit with other files to produce an executable, -- +-- this unit does not by itself cause the resulting executable to be -- +-- covered by the GNU General Public License. This exception does not -- +-- however invalidate any other reasons why the executable file might be -- +-- covered by the GNU Public License. -- +-- -- +-- GNAT was originally developed by the GNAT team at New York University. -- +-- Extensive contributions were provided by Ada Core Technologies Inc. -- +-- -- +------------------------------------------------------------------------------ + +with Ada.Unchecked_Conversion; use Ada; +with Interfaces; use Interfaces; +with Interfaces.Fortran; use Interfaces.Fortran; +with Interfaces.Fortran.BLAS; use Interfaces.Fortran.BLAS; +with System.Generic_Array_Operations; use System.Generic_Array_Operations; + +package body System.Generic_Complex_BLAS is + + Is_Single : constant Boolean := + Real'Machine_Mantissa = Fortran.Real'Machine_Mantissa + and then Fortran.Real (Real'First) = Fortran.Real'First + and then Fortran.Real (Real'Last) = Fortran.Real'Last; + + Is_Double : constant Boolean := + Real'Machine_Mantissa = Double_Precision'Machine_Mantissa + and then + Double_Precision (Real'First) = Double_Precision'First + and then + Double_Precision (Real'Last) = Double_Precision'Last; + + subtype Complex is Complex_Types.Complex; + + -- Local subprograms + + function To_Double_Precision (X : Real) return Double_Precision; + pragma Inline (To_Double_Precision); + + function To_Double_Complex (X : Complex) return Double_Complex; + pragma Inline (To_Double_Complex); + + function To_Complex (X : Double_Complex) return Complex; + function To_Complex (X : Fortran.Complex) return Complex; + pragma Inline (To_Complex); + + function To_Fortran (X : Complex) return Fortran.Complex; + pragma Inline (To_Fortran); + + -- Instantiations + + function To_Double_Complex is new + Vector_Elementwise_Operation + (X_Scalar => Complex_Types.Complex, + Result_Scalar => Fortran.Double_Complex, + X_Vector => Complex_Vector, + Result_Vector => BLAS.Double_Complex_Vector, + Operation => To_Double_Complex); + + function To_Complex is new + Vector_Elementwise_Operation + (X_Scalar => Fortran.Double_Complex, + Result_Scalar => Complex, + X_Vector => BLAS.Double_Complex_Vector, + Result_Vector => Complex_Vector, + Operation => To_Complex); + + function To_Double_Complex is new + Matrix_Elementwise_Operation + (X_Scalar => Complex, + Result_Scalar => Double_Complex, + X_Matrix => Complex_Matrix, + Result_Matrix => BLAS.Double_Complex_Matrix, + Operation => To_Double_Complex); + + function To_Complex is new + Matrix_Elementwise_Operation + (X_Scalar => Double_Complex, + Result_Scalar => Complex, + X_Matrix => BLAS.Double_Complex_Matrix, + Result_Matrix => Complex_Matrix, + Operation => To_Complex); + + function To_Double_Precision (X : Real) return Double_Precision is + begin + return Double_Precision (X); + end To_Double_Precision; + + function To_Double_Complex (X : Complex) return Double_Complex is + begin + return (To_Double_Precision (X.Re), To_Double_Precision (X.Im)); + end To_Double_Complex; + + function To_Complex (X : Double_Complex) return Complex is + begin + return (Real (X.Re), Real (X.Im)); + end To_Complex; + + function To_Complex (X : Fortran.Complex) return Complex is + begin + return (Real (X.Re), Real (X.Im)); + end To_Complex; + + function To_Fortran (X : Complex) return Fortran.Complex is + begin + return (Fortran.Real (X.Re), Fortran.Real (X.Im)); + end To_Fortran; + + --------- + -- dot -- + --------- + + function dot + (N : Positive; + X : Complex_Vector; + Inc_X : Integer := 1; + Y : Complex_Vector; + Inc_Y : Integer := 1) return Complex + is + begin + if Is_Single then + declare + type X_Ptr is access all BLAS.Complex_Vector (X'Range); + type Y_Ptr is access all BLAS.Complex_Vector (Y'Range); + function Conv_X is new Unchecked_Conversion (Address, X_Ptr); + function Conv_Y is new Unchecked_Conversion (Address, Y_Ptr); + begin + return To_Complex (BLAS.cdot (N, Conv_X (X'Address).all, Inc_X, + Conv_Y (Y'Address).all, Inc_Y)); + end; + + elsif Is_Double then + declare + type X_Ptr is access all BLAS.Double_Complex_Vector (X'Range); + type Y_Ptr is access all BLAS.Double_Complex_Vector (Y'Range); + function Conv_X is new Unchecked_Conversion (Address, X_Ptr); + function Conv_Y is new Unchecked_Conversion (Address, Y_Ptr); + begin + return To_Complex (BLAS.zdot (N, Conv_X (X'Address).all, Inc_X, + Conv_Y (Y'Address).all, Inc_Y)); + end; + + else + return To_Complex (BLAS.zdot (N, To_Double_Complex (X), Inc_X, + To_Double_Complex (Y), Inc_Y)); + end if; + end dot; + + ---------- + -- gemm -- + ---------- + + procedure gemm + (Trans_A : access constant Character; + Trans_B : access constant Character; + M : Positive; + N : Positive; + K : Positive; + Alpha : Complex := (1.0, 1.0); + A : Complex_Matrix; + Ld_A : Integer; + B : Complex_Matrix; + Ld_B : Integer; + Beta : Complex := (0.0, 0.0); + C : in out Complex_Matrix; + Ld_C : Integer) + is + begin + if Is_Single then + declare + subtype A_Type is BLAS.Complex_Matrix (A'Range (1), A'Range (2)); + subtype B_Type is BLAS.Complex_Matrix (B'Range (1), B'Range (2)); + type C_Ptr is + access all BLAS.Complex_Matrix (C'Range (1), C'Range (2)); + function Conv_A is + new Unchecked_Conversion (Complex_Matrix, A_Type); + function Conv_B is + new Unchecked_Conversion (Complex_Matrix, B_Type); + function Conv_C is + new Unchecked_Conversion (Address, C_Ptr); + begin + BLAS.cgemm (Trans_A, Trans_B, M, N, K, To_Fortran (Alpha), + Conv_A (A), Ld_A, Conv_B (B), Ld_B, To_Fortran (Beta), + Conv_C (C'Address).all, Ld_C); + end; + + elsif Is_Double then + declare + subtype A_Type is + BLAS.Double_Complex_Matrix (A'Range (1), A'Range (2)); + subtype B_Type is + BLAS.Double_Complex_Matrix (B'Range (1), B'Range (2)); + type C_Ptr is access all + BLAS.Double_Complex_Matrix (C'Range (1), C'Range (2)); + function Conv_A is + new Unchecked_Conversion (Complex_Matrix, A_Type); + function Conv_B is + new Unchecked_Conversion (Complex_Matrix, B_Type); + function Conv_C is new Unchecked_Conversion (Address, C_Ptr); + begin + BLAS.zgemm (Trans_A, Trans_B, M, N, K, To_Double_Complex (Alpha), + Conv_A (A), Ld_A, Conv_B (B), Ld_B, + To_Double_Complex (Beta), + Conv_C (C'Address).all, Ld_C); + end; + + else + declare + DP_C : BLAS.Double_Complex_Matrix (C'Range (1), C'Range (2)); + begin + if Beta.Re /= 0.0 or else Beta.Im /= 0.0 then + DP_C := To_Double_Complex (C); + end if; + + BLAS.zgemm (Trans_A, Trans_B, M, N, K, To_Double_Complex (Alpha), + To_Double_Complex (A), Ld_A, + To_Double_Complex (B), Ld_B, To_Double_Complex (Beta), + DP_C, Ld_C); + + C := To_Complex (DP_C); + end; + end if; + end gemm; + + ---------- + -- gemv -- + ---------- + + procedure gemv + (Trans : access constant Character; + M : Natural := 0; + N : Natural := 0; + Alpha : Complex := (1.0, 1.0); + A : Complex_Matrix; + Ld_A : Positive; + X : Complex_Vector; + Inc_X : Integer := 1; + Beta : Complex := (0.0, 0.0); + Y : in out Complex_Vector; + Inc_Y : Integer := 1) + is + begin + if Is_Single then + declare + subtype A_Type is BLAS.Complex_Matrix (A'Range (1), A'Range (2)); + subtype X_Type is BLAS.Complex_Vector (X'Range); + type Y_Ptr is access all BLAS.Complex_Vector (Y'Range); + function Conv_A is + new Unchecked_Conversion (Complex_Matrix, A_Type); + function Conv_X is + new Unchecked_Conversion (Complex_Vector, X_Type); + function Conv_Y is + new Unchecked_Conversion (Address, Y_Ptr); + begin + BLAS.cgemv (Trans, M, N, To_Fortran (Alpha), + Conv_A (A), Ld_A, Conv_X (X), Inc_X, To_Fortran (Beta), + Conv_Y (Y'Address).all, Inc_Y); + end; + + elsif Is_Double then + declare + subtype A_Type is + BLAS.Double_Complex_Matrix (A'Range (1), A'Range (2)); + subtype X_Type is + BLAS.Double_Complex_Vector (X'Range); + type Y_Ptr is access all BLAS.Double_Complex_Vector (Y'Range); + function Conv_A is + new Unchecked_Conversion (Complex_Matrix, A_Type); + function Conv_X is + new Unchecked_Conversion (Complex_Vector, X_Type); + function Conv_Y is + new Unchecked_Conversion (Address, Y_Ptr); + begin + BLAS.zgemv (Trans, M, N, To_Double_Complex (Alpha), + Conv_A (A), Ld_A, Conv_X (X), Inc_X, + To_Double_Complex (Beta), + Conv_Y (Y'Address).all, Inc_Y); + end; + + else + declare + DP_Y : BLAS.Double_Complex_Vector (Y'Range); + begin + if Beta.Re /= 0.0 or else Beta.Im /= 0.0 then + DP_Y := To_Double_Complex (Y); + end if; + + BLAS.zgemv (Trans, M, N, To_Double_Complex (Alpha), + To_Double_Complex (A), Ld_A, + To_Double_Complex (X), Inc_X, To_Double_Complex (Beta), + DP_Y, Inc_Y); + + Y := To_Complex (DP_Y); + end; + end if; + end gemv; + + ---------- + -- nrm2 -- + ---------- + + function nrm2 + (N : Natural; + X : Complex_Vector; + Inc_X : Integer := 1) return Real + is + begin + if Is_Single then + declare + subtype X_Type is BLAS.Complex_Vector (X'Range); + function Conv_X is + new Unchecked_Conversion (Complex_Vector, X_Type); + begin + return Real (BLAS.scnrm2 (N, Conv_X (X), Inc_X)); + end; + + elsif Is_Double then + declare + subtype X_Type is BLAS.Double_Complex_Vector (X'Range); + function Conv_X is + new Unchecked_Conversion (Complex_Vector, X_Type); + begin + return Real (BLAS.dznrm2 (N, Conv_X (X), Inc_X)); + end; + + else + return Real (BLAS.dznrm2 (N, To_Double_Complex (X), Inc_X)); + end if; + end nrm2; + +end System.Generic_Complex_BLAS; diff --git a/gcc/ada/s-gecobl.ads b/gcc/ada/s-gecobl.ads new file mode 100644 index 0000000..9d9d216 --- /dev/null +++ b/gcc/ada/s-gecobl.ads @@ -0,0 +1,104 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- SYSTEM.GENERIC_COMPLEX_BLAS -- +-- -- +-- S p e c -- +-- -- +-- Copyright (C) 2006, 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 2, 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. See the GNU General Public License -- +-- for more details. You should have received a copy of the GNU General -- +-- Public License distributed with GNAT; see file COPYING. If not, write -- +-- to the Free Software Foundation, 51 Franklin Street, Fifth Floor, -- +-- Boston, MA 02110-1301, USA. -- +-- -- +-- As a special exception, if other files instantiate generics from this -- +-- unit, or you link this unit with other files to produce an executable, -- +-- this unit does not by itself cause the resulting executable to be -- +-- covered by the GNU General Public License. This exception does not -- +-- however invalidate any other reasons why the executable file might be -- +-- covered by the GNU Public License. -- +-- -- +-- GNAT was originally developed by the GNAT team at New York University. -- +-- Extensive contributions were provided by Ada Core Technologies Inc. -- +-- -- +------------------------------------------------------------------------------ + +-- Package comment required ??? + +with Ada.Numerics.Generic_Complex_Types; + +generic + type Real is digits <>; + with package Complex_Types is new Ada.Numerics.Generic_Complex_Types (Real); + use Complex_Types; + + type Complex_Vector is array (Integer range <>) of Complex; + type Complex_Matrix is array (Integer range <>, Integer range <>) + of Complex; +package System.Generic_Complex_BLAS is + pragma Pure; + + -- Although BLAS support is only available for IEEE single and double + -- compatible floating-point types, this unit will accept any type + -- and apply conversions as necessary, with possible loss of + -- precision and range. + + No_Trans : aliased constant Character := 'N'; + Trans : aliased constant Character := 'T'; + Conj_Trans : aliased constant Character := 'C'; + + -- BLAS Level 1 Subprograms and Types + + function dot + (N : Positive; + X : Complex_Vector; + Inc_X : Integer := 1; + Y : Complex_Vector; + Inc_Y : Integer := 1) return Complex; + + function nrm2 + (N : Natural; + X : Complex_Vector; + Inc_X : Integer := 1) return Real; + + procedure gemv + (Trans : access constant Character; + M : Natural := 0; + N : Natural := 0; + Alpha : Complex := (1.0, 1.0); + A : Complex_Matrix; + Ld_A : Positive; + X : Complex_Vector; + Inc_X : Integer := 1; -- must be non-zero + Beta : Complex := (0.0, 0.0); + Y : in out Complex_Vector; + Inc_Y : Integer := 1); -- must be non-zero + + -- BLAS Level 3 + + -- gemm s, d, c, z Matrix-matrix product of general matrices + + procedure gemm + (Trans_A : access constant Character; + Trans_B : access constant Character; + M : Positive; + N : Positive; + K : Positive; + Alpha : Complex := (1.0, 1.0); + A : Complex_Matrix; + Ld_A : Integer; + B : Complex_Matrix; + Ld_B : Integer; + Beta : Complex := (0.0, 0.0); + C : in out Complex_Matrix; + Ld_C : Integer); + +end System.Generic_Complex_BLAS; diff --git a/gcc/ada/s-gecola.adb b/gcc/ada/s-gecola.adb new file mode 100644 index 0000000..ef3ea1c --- /dev/null +++ b/gcc/ada/s-gecola.adb @@ -0,0 +1,495 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- SYSTEM.GENERIC_COMPLEX_LAPACK -- +-- -- +-- B o d y -- +-- -- +-- Copyright (C) 2006, 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 2, 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. See the GNU General Public License -- +-- for more details. You should have received a copy of the GNU General -- +-- Public License distributed with GNAT; see file COPYING. If not, write -- +-- to the Free Software Foundation, 51 Franklin Street, Fifth Floor, -- +-- Boston, MA 02110-1301, USA. -- +-- -- +-- As a special exception, if other files instantiate generics from this -- +-- unit, or you link this unit with other files to produce an executable, -- +-- this unit does not by itself cause the resulting executable to be -- +-- covered by the GNU General Public License. This exception does not -- +-- however invalidate any other reasons why the executable file might be -- +-- covered by the GNU Public License. -- +-- -- +-- GNAT was originally developed by the GNAT team at New York University. -- +-- Extensive contributions were provided by Ada Core Technologies Inc. -- +-- -- +------------------------------------------------------------------------------ + +with Ada.Unchecked_Conversion; use Ada; +with Interfaces; use Interfaces; +with Interfaces.Fortran; use Interfaces.Fortran; +with Interfaces.Fortran.BLAS; use Interfaces.Fortran.BLAS; +with Interfaces.Fortran.LAPACK; use Interfaces.Fortran.LAPACK; +with System.Generic_Array_Operations; use System.Generic_Array_Operations; + +package body System.Generic_Complex_LAPACK is + + Is_Single : constant Boolean := + Real'Machine_Mantissa = Fortran.Real'Machine_Mantissa + and then Fortran.Real (Real'First) = Fortran.Real'First + and then Fortran.Real (Real'Last) = Fortran.Real'Last; + + Is_Double : constant Boolean := + Real'Machine_Mantissa = Double_Precision'Machine_Mantissa + and then + Double_Precision (Real'First) = Double_Precision'First + and then + Double_Precision (Real'Last) = Double_Precision'Last; + + subtype Complex is Complex_Types.Complex; + + -- Local subprograms + + function To_Double_Precision (X : Real) return Double_Precision; + pragma Inline (To_Double_Precision); + + function To_Real (X : Double_Precision) return Real; + pragma Inline (To_Real); + + function To_Double_Complex (X : Complex) return Double_Complex; + pragma Inline (To_Double_Complex); + + function To_Complex (X : Double_Complex) return Complex; + pragma Inline (To_Complex); + + -- Instantiations + + function To_Double_Precision is new + Vector_Elementwise_Operation + (X_Scalar => Real, + Result_Scalar => Double_Precision, + X_Vector => Real_Vector, + Result_Vector => Double_Precision_Vector, + Operation => To_Double_Precision); + + function To_Real is new + Vector_Elementwise_Operation + (X_Scalar => Double_Precision, + Result_Scalar => Real, + X_Vector => Double_Precision_Vector, + Result_Vector => Real_Vector, + Operation => To_Real); + + function To_Double_Complex is new + Matrix_Elementwise_Operation + (X_Scalar => Complex, + Result_Scalar => Double_Complex, + X_Matrix => Complex_Matrix, + Result_Matrix => Double_Complex_Matrix, + Operation => To_Double_Complex); + + function To_Complex is new + Matrix_Elementwise_Operation + (X_Scalar => Double_Complex, + Result_Scalar => Complex, + X_Matrix => Double_Complex_Matrix, + Result_Matrix => Complex_Matrix, + Operation => To_Complex); + + function To_Double_Precision (X : Real) return Double_Precision is + begin + return Double_Precision (X); + end To_Double_Precision; + + function To_Real (X : Double_Precision) return Real is + begin + return Real (X); + end To_Real; + + function To_Double_Complex (X : Complex) return Double_Complex is + begin + return (To_Double_Precision (X.Re), To_Double_Precision (X.Im)); + end To_Double_Complex; + + function To_Complex (X : Double_Complex) return Complex is + begin + return (Real (X.Re), Real (X.Im)); + end To_Complex; + + ----------- + -- getrf -- + ----------- + + procedure getrf + (M : Natural; + N : Natural; + A : in out Complex_Matrix; + Ld_A : Positive; + I_Piv : out Integer_Vector; + Info : access Integer) + is + begin + if Is_Single then + declare + type A_Ptr is + access all BLAS.Complex_Matrix (A'Range (1), A'Range (2)); + function Conv_A is new Unchecked_Conversion (Address, A_Ptr); + begin + cgetrf (M, N, Conv_A (A'Address).all, Ld_A, + LAPACK.Integer_Vector (I_Piv), Info); + end; + + elsif Is_Double then + declare + type A_Ptr is + access all Double_Complex_Matrix (A'Range (1), A'Range (2)); + function Conv_A is new Unchecked_Conversion (Address, A_Ptr); + begin + zgetrf (M, N, Conv_A (A'Address).all, Ld_A, + LAPACK.Integer_Vector (I_Piv), Info); + end; + + else + declare + DP_A : Double_Complex_Matrix (A'Range (1), A'Range (2)); + begin + DP_A := To_Double_Complex (A); + zgetrf (M, N, DP_A, Ld_A, LAPACK.Integer_Vector (I_Piv), Info); + A := To_Complex (DP_A); + end; + end if; + end getrf; + + ----------- + -- getri -- + ----------- + + procedure getri + (N : Natural; + A : in out Complex_Matrix; + Ld_A : Positive; + I_Piv : Integer_Vector; + Work : in out Complex_Vector; + L_Work : Integer; + Info : access Integer) + is + begin + if Is_Single then + declare + type A_Ptr is + access all BLAS.Complex_Matrix (A'Range (1), A'Range (2)); + type Work_Ptr is + access all BLAS.Complex_Vector (Work'Range); + function Conv_A is new Unchecked_Conversion (Address, A_Ptr); + function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr); + begin + cgetri (N, Conv_A (A'Address).all, Ld_A, + LAPACK.Integer_Vector (I_Piv), + Conv_Work (Work'Address).all, L_Work, + Info); + end; + + elsif Is_Double then + declare + type A_Ptr is + access all Double_Complex_Matrix (A'Range (1), A'Range (2)); + type Work_Ptr is + access all Double_Complex_Vector (Work'Range); + function Conv_A is new Unchecked_Conversion (Address, A_Ptr); + function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr); + begin + zgetri (N, Conv_A (A'Address).all, Ld_A, + LAPACK.Integer_Vector (I_Piv), + Conv_Work (Work'Address).all, L_Work, + Info); + end; + + else + declare + DP_A : Double_Complex_Matrix (A'Range (1), A'Range (2)); + DP_Work : Double_Complex_Vector (Work'Range); + begin + DP_A := To_Double_Complex (A); + zgetri (N, DP_A, Ld_A, LAPACK.Integer_Vector (I_Piv), + DP_Work, L_Work, Info); + A := To_Complex (DP_A); + Work (1) := To_Complex (DP_Work (1)); + end; + end if; + end getri; + + ----------- + -- getrs -- + ----------- + + procedure getrs + (Trans : access constant Character; + N : Natural; + N_Rhs : Natural; + A : Complex_Matrix; + Ld_A : Positive; + I_Piv : Integer_Vector; + B : in out Complex_Matrix; + Ld_B : Positive; + Info : access Integer) + is + begin + if Is_Single then + declare + subtype A_Type is BLAS.Complex_Matrix (A'Range (1), A'Range (2)); + type B_Ptr is + access all BLAS.Complex_Matrix (B'Range (1), B'Range (2)); + function Conv_A is + new Unchecked_Conversion (Complex_Matrix, A_Type); + function Conv_B is new Unchecked_Conversion (Address, B_Ptr); + begin + cgetrs (Trans, N, N_Rhs, + Conv_A (A), Ld_A, + LAPACK.Integer_Vector (I_Piv), + Conv_B (B'Address).all, Ld_B, + Info); + end; + + elsif Is_Double then + declare + subtype A_Type is + Double_Complex_Matrix (A'Range (1), A'Range (2)); + type B_Ptr is + access all Double_Complex_Matrix (B'Range (1), B'Range (2)); + function Conv_A is + new Unchecked_Conversion (Complex_Matrix, A_Type); + function Conv_B is new Unchecked_Conversion (Address, B_Ptr); + begin + zgetrs (Trans, N, N_Rhs, + Conv_A (A), Ld_A, + LAPACK.Integer_Vector (I_Piv), + Conv_B (B'Address).all, Ld_B, + Info); + end; + + else + declare + DP_A : Double_Complex_Matrix (A'Range (1), A'Range (2)); + DP_B : Double_Complex_Matrix (B'Range (1), B'Range (2)); + begin + DP_A := To_Double_Complex (A); + DP_B := To_Double_Complex (B); + zgetrs (Trans, N, N_Rhs, + DP_A, Ld_A, + LAPACK.Integer_Vector (I_Piv), + DP_B, Ld_B, + Info); + B := To_Complex (DP_B); + end; + end if; + end getrs; + + procedure heevr + (Job_Z : access constant Character; + Rng : access constant Character; + Uplo : access constant Character; + N : Natural; + A : in out Complex_Matrix; + Ld_A : Positive; + Vl, Vu : Real := 0.0; + Il, Iu : Integer := 1; + Abs_Tol : Real := 0.0; + M : out Integer; + W : out Real_Vector; + Z : out Complex_Matrix; + Ld_Z : Positive; + I_Supp_Z : out Integer_Vector; + Work : out Complex_Vector; + L_Work : Integer; + R_Work : out Real_Vector; + LR_Work : Integer; + I_Work : out Integer_Vector; + LI_Work : Integer; + Info : access Integer) + is + begin + if Is_Single then + declare + type A_Ptr is + access all BLAS.Complex_Matrix (A'Range (1), A'Range (2)); + type W_Ptr is + access all BLAS.Real_Vector (W'Range); + type Z_Ptr is + access all BLAS.Complex_Matrix (Z'Range (1), Z'Range (2)); + type Work_Ptr is access all BLAS.Complex_Vector (Work'Range); + type R_Work_Ptr is access all BLAS.Real_Vector (R_Work'Range); + + function Conv_A is new Unchecked_Conversion (Address, A_Ptr); + function Conv_W is new Unchecked_Conversion (Address, W_Ptr); + function Conv_Z is new Unchecked_Conversion (Address, Z_Ptr); + function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr); + function Conv_R_Work is + new Unchecked_Conversion (Address, R_Work_Ptr); + begin + cheevr (Job_Z, Rng, Uplo, N, + Conv_A (A'Address).all, Ld_A, + Fortran.Real (Vl), Fortran.Real (Vu), + Il, Iu, Fortran.Real (Abs_Tol), M, + Conv_W (W'Address).all, + Conv_Z (Z'Address).all, Ld_Z, + LAPACK.Integer_Vector (I_Supp_Z), + Conv_Work (Work'Address).all, L_Work, + Conv_R_Work (R_Work'Address).all, LR_Work, + LAPACK.Integer_Vector (I_Work), LI_Work, Info); + end; + + elsif Is_Double then + declare + type A_Ptr is + access all BLAS.Double_Complex_Matrix (A'Range (1), A'Range (2)); + type W_Ptr is + access all BLAS.Double_Precision_Vector (W'Range); + type Z_Ptr is + access all BLAS.Double_Complex_Matrix (Z'Range (1), Z'Range (2)); + type Work_Ptr is + access all BLAS.Double_Complex_Vector (Work'Range); + type R_Work_Ptr is + access all BLAS.Double_Precision_Vector (R_Work'Range); + + function Conv_A is new Unchecked_Conversion (Address, A_Ptr); + function Conv_W is new Unchecked_Conversion (Address, W_Ptr); + function Conv_Z is new Unchecked_Conversion (Address, Z_Ptr); + function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr); + function Conv_R_Work is + new Unchecked_Conversion (Address, R_Work_Ptr); + begin + zheevr (Job_Z, Rng, Uplo, N, + Conv_A (A'Address).all, Ld_A, + Double_Precision (Vl), Double_Precision (Vu), + Il, Iu, Double_Precision (Abs_Tol), M, + Conv_W (W'Address).all, + Conv_Z (Z'Address).all, Ld_Z, + LAPACK.Integer_Vector (I_Supp_Z), + Conv_Work (Work'Address).all, L_Work, + Conv_R_Work (R_Work'Address).all, LR_Work, + LAPACK.Integer_Vector (I_Work), LI_Work, Info); + end; + + else + declare + DP_A : Double_Complex_Matrix (A'Range (1), A'Range (2)); + DP_W : Double_Precision_Vector (W'Range); + DP_Z : Double_Complex_Matrix (Z'Range (1), Z'Range (2)); + DP_Work : Double_Complex_Vector (Work'Range); + DP_R_Work : Double_Precision_Vector (R_Work'Range); + + begin + DP_A := To_Double_Complex (A); + + zheevr (Job_Z, Rng, Uplo, N, + DP_A, Ld_A, + Double_Precision (Vl), Double_Precision (Vu), + Il, Iu, Double_Precision (Abs_Tol), M, + DP_W, DP_Z, Ld_Z, + LAPACK.Integer_Vector (I_Supp_Z), + DP_Work, L_Work, + DP_R_Work, LR_Work, + LAPACK.Integer_Vector (I_Work), LI_Work, Info); + + A := To_Complex (DP_A); + W := To_Real (DP_W); + Z := To_Complex (DP_Z); + + Work (1) := To_Complex (DP_Work (1)); + R_Work (1) := To_Real (DP_R_Work (1)); + end; + end if; + end heevr; + + ----------- + -- steqr -- + ----------- + + procedure steqr + (Comp_Z : access constant Character; + N : Natural; + D : in out Real_Vector; + E : in out Real_Vector; + Z : in out Complex_Matrix; + Ld_Z : Positive; + Work : out Real_Vector; + Info : access Integer) + is + begin + if Is_Single then + declare + type D_Ptr is access all BLAS.Real_Vector (D'Range); + type E_Ptr is access all BLAS.Real_Vector (E'Range); + type Z_Ptr is + access all BLAS.Complex_Matrix (Z'Range (1), Z'Range (2)); + type Work_Ptr is + access all BLAS.Real_Vector (Work'Range); + function Conv_D is new Unchecked_Conversion (Address, D_Ptr); + function Conv_E is new Unchecked_Conversion (Address, E_Ptr); + function Conv_Z is new Unchecked_Conversion (Address, Z_Ptr); + function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr); + begin + csteqr (Comp_Z, N, + Conv_D (D'Address).all, + Conv_E (E'Address).all, + Conv_Z (Z'Address).all, + Ld_Z, + Conv_Work (Work'Address).all, + Info); + end; + + elsif Is_Double then + declare + type D_Ptr is access all Double_Precision_Vector (D'Range); + type E_Ptr is access all Double_Precision_Vector (E'Range); + type Z_Ptr is + access all Double_Complex_Matrix (Z'Range (1), Z'Range (2)); + type Work_Ptr is + access all Double_Precision_Vector (Work'Range); + function Conv_D is new Unchecked_Conversion (Address, D_Ptr); + function Conv_E is new Unchecked_Conversion (Address, E_Ptr); + function Conv_Z is new Unchecked_Conversion (Address, Z_Ptr); + function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr); + begin + zsteqr (Comp_Z, N, + Conv_D (D'Address).all, + Conv_E (E'Address).all, + Conv_Z (Z'Address).all, + Ld_Z, + Conv_Work (Work'Address).all, + Info); + end; + + else + declare + DP_D : Double_Precision_Vector (D'Range); + DP_E : Double_Precision_Vector (E'Range); + DP_Z : Double_Complex_Matrix (Z'Range (1), Z'Range (2)); + DP_Work : Double_Precision_Vector (Work'Range); + begin + DP_D := To_Double_Precision (D); + DP_E := To_Double_Precision (E); + + if Comp_Z.all = 'V' then + DP_Z := To_Double_Complex (Z); + end if; + + zsteqr (Comp_Z, N, DP_D, DP_E, DP_Z, Ld_Z, DP_Work, Info); + + D := To_Real (DP_D); + E := To_Real (DP_E); + + if Comp_Z.all /= 'N' then + Z := To_Complex (DP_Z); + end if; + end; + end if; + end steqr; + +end System.Generic_Complex_LAPACK; diff --git a/gcc/ada/s-gecola.ads b/gcc/ada/s-gecola.ads new file mode 100644 index 0000000..2597a90 --- /dev/null +++ b/gcc/ada/s-gecola.ads @@ -0,0 +1,133 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- SYSTEM.GENERIC_COMPLEX_LAPACK -- +-- -- +-- S p e c -- +-- -- +-- Copyright (C) 2006, 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 2, 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. See the GNU General Public License -- +-- for more details. You should have received a copy of the GNU General -- +-- Public License distributed with GNAT; see file COPYING. If not, write -- +-- to the Free Software Foundation, 51 Franklin Street, Fifth Floor, -- +-- Boston, MA 02110-1301, USA. -- +-- -- +-- As a special exception, if other files instantiate generics from this -- +-- unit, or you link this unit with other files to produce an executable, -- +-- this unit does not by itself cause the resulting executable to be -- +-- covered by the GNU General Public License. This exception does not -- +-- however invalidate any other reasons why the executable file might be -- +-- covered by the GNU Public License. -- +-- -- +-- GNAT was originally developed by the GNAT team at New York University. -- +-- Extensive contributions were provided by Ada Core Technologies Inc. -- +-- -- +------------------------------------------------------------------------------ + +-- Package comment required ??? + +with Ada.Numerics.Generic_Complex_Types; +generic + type Real is digits <>; + type Real_Vector is array (Integer range <>) of Real; + + with package Complex_Types is new Ada.Numerics.Generic_Complex_Types (Real); + use Complex_Types; + + type Complex_Vector is array (Integer range <>) of Complex; + type Complex_Matrix is array (Integer range <>, Integer range <>) + of Complex; +package System.Generic_Complex_LAPACK is + pragma Pure; + + type Integer_Vector is array (Integer range <>) of Integer; + + Upper : aliased constant Character := 'U'; + Lower : aliased constant Character := 'L'; + + -- LAPACK Computational Routines + + -- getrf computes LU factorization of a general m-by-n matrix + + procedure getrf + (M : Natural; + N : Natural; + A : in out Complex_Matrix; + Ld_A : Positive; + I_Piv : out Integer_Vector; + Info : access Integer); + + -- getri computes inverse of an LU-factored square matrix, + -- with multiple right-hand sides + + procedure getri + (N : Natural; + A : in out Complex_Matrix; + Ld_A : Positive; + I_Piv : Integer_Vector; + Work : in out Complex_Vector; + L_Work : Integer; + Info : access Integer); + + -- getrs solves a system of linear equations with an LU-factored + -- square matrix, with multiple right-hand sides + + procedure getrs + (Trans : access constant Character; + N : Natural; + N_Rhs : Natural; + A : Complex_Matrix; + Ld_A : Positive; + I_Piv : Integer_Vector; + B : in out Complex_Matrix; + Ld_B : Positive; + Info : access Integer); + + -- heevr computes selected eigenvalues and, optionally, + -- eigenvectors of a Hermitian matrix using the Relatively + -- Robust Representations + + procedure heevr + (Job_Z : access constant Character; + Rng : access constant Character; + Uplo : access constant Character; + N : Natural; + A : in out Complex_Matrix; + Ld_A : Positive; + Vl, Vu : Real := 0.0; + Il, Iu : Integer := 1; + Abs_Tol : Real := 0.0; + M : out Integer; + W : out Real_Vector; + Z : out Complex_Matrix; + Ld_Z : Positive; + I_Supp_Z : out Integer_Vector; + Work : out Complex_Vector; + L_Work : Integer; + R_Work : out Real_Vector; + LR_Work : Integer; + I_Work : out Integer_Vector; + LI_Work : Integer; + Info : access Integer); + + -- steqr computes all eigenvalues and eigenvectors of a symmetric or + -- Hermitian matrix reduced to tridiagonal form (QR algorithm) + + procedure steqr + (Comp_Z : access constant Character; + N : Natural; + D : in out Real_Vector; + E : in out Real_Vector; + Z : in out Complex_Matrix; + Ld_Z : Positive; + Work : out Real_Vector; + Info : access Integer); + +end System.Generic_Complex_LAPACK; diff --git a/gcc/ada/s-gerebl.adb b/gcc/ada/s-gerebl.adb new file mode 100644 index 0000000..810a167 --- /dev/null +++ b/gcc/ada/s-gerebl.adb @@ -0,0 +1,313 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- SYSTEM.GENERIC_REAL_BLAS -- +-- -- +-- B o d y -- +-- -- +-- Copyright (C) 2006, 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 2, 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. See the GNU General Public License -- +-- for more details. You should have received a copy of the GNU General -- +-- Public License distributed with GNAT; see file COPYING. If not, write -- +-- to the Free Software Foundation, 51 Franklin Street, Fifth Floor, -- +-- Boston, MA 02110-1301, USA. -- +-- -- +-- As a special exception, if other files instantiate generics from this -- +-- unit, or you link this unit with other files to produce an executable, -- +-- this unit does not by itself cause the resulting executable to be -- +-- covered by the GNU General Public License. This exception does not -- +-- however invalidate any other reasons why the executable file might be -- +-- covered by the GNU Public License. -- +-- -- +-- GNAT was originally developed by the GNAT team at New York University. -- +-- Extensive contributions were provided by Ada Core Technologies Inc. -- +-- -- +------------------------------------------------------------------------------ + +with Ada.Unchecked_Conversion; use Ada; +with Interfaces; use Interfaces; +with Interfaces.Fortran; use Interfaces.Fortran; +with Interfaces.Fortran.BLAS; use Interfaces.Fortran.BLAS; +with System.Generic_Array_Operations; use System.Generic_Array_Operations; + +package body System.Generic_Real_BLAS is + + Is_Single : constant Boolean := + Real'Machine_Mantissa = Fortran.Real'Machine_Mantissa + and then Fortran.Real (Real'First) = Fortran.Real'First + and then Fortran.Real (Real'Last) = Fortran.Real'Last; + + Is_Double : constant Boolean := + Real'Machine_Mantissa = Double_Precision'Machine_Mantissa + and then + Double_Precision (Real'First) = Double_Precision'First + and then + Double_Precision (Real'Last) = Double_Precision'Last; + + -- Local subprograms + + function To_Double_Precision (X : Real) return Double_Precision; + pragma Inline_Always (To_Double_Precision); + + function To_Real (X : Double_Precision) return Real; + pragma Inline_Always (To_Real); + + -- Instantiations + + function To_Double_Precision is new + Vector_Elementwise_Operation + (X_Scalar => Real, + Result_Scalar => Double_Precision, + X_Vector => Real_Vector, + Result_Vector => Double_Precision_Vector, + Operation => To_Double_Precision); + + function To_Real is new + Vector_Elementwise_Operation + (X_Scalar => Double_Precision, + Result_Scalar => Real, + X_Vector => Double_Precision_Vector, + Result_Vector => Real_Vector, + Operation => To_Real); + + function To_Double_Precision is new + Matrix_Elementwise_Operation + (X_Scalar => Real, + Result_Scalar => Double_Precision, + X_Matrix => Real_Matrix, + Result_Matrix => Double_Precision_Matrix, + Operation => To_Double_Precision); + + function To_Real is new + Matrix_Elementwise_Operation + (X_Scalar => Double_Precision, + Result_Scalar => Real, + X_Matrix => Double_Precision_Matrix, + Result_Matrix => Real_Matrix, + Operation => To_Real); + + function To_Double_Precision (X : Real) return Double_Precision is + begin + return Double_Precision (X); + end To_Double_Precision; + + function To_Real (X : Double_Precision) return Real is + begin + return Real (X); + end To_Real; + + --------- + -- dot -- + --------- + + function dot + (N : Positive; + X : Real_Vector; + Inc_X : Integer := 1; + Y : Real_Vector; + Inc_Y : Integer := 1) return Real + is + begin + if Is_Single then + declare + type X_Ptr is access all BLAS.Real_Vector (X'Range); + type Y_Ptr is access all BLAS.Real_Vector (Y'Range); + function Conv_X is new Unchecked_Conversion (Address, X_Ptr); + function Conv_Y is new Unchecked_Conversion (Address, Y_Ptr); + begin + return Real (sdot (N, Conv_X (X'Address).all, Inc_X, + Conv_Y (Y'Address).all, Inc_Y)); + end; + + elsif Is_Double then + declare + type X_Ptr is access all BLAS.Double_Precision_Vector (X'Range); + type Y_Ptr is access all BLAS.Double_Precision_Vector (Y'Range); + function Conv_X is new Unchecked_Conversion (Address, X_Ptr); + function Conv_Y is new Unchecked_Conversion (Address, Y_Ptr); + begin + return Real (ddot (N, Conv_X (X'Address).all, Inc_X, + Conv_Y (Y'Address).all, Inc_Y)); + end; + + else + return Real (ddot (N, To_Double_Precision (X), Inc_X, + To_Double_Precision (Y), Inc_Y)); + end if; + end dot; + + ---------- + -- gemm -- + ---------- + + procedure gemm + (Trans_A : access constant Character; + Trans_B : access constant Character; + M : Positive; + N : Positive; + K : Positive; + Alpha : Real := 1.0; + A : Real_Matrix; + Ld_A : Integer; + B : Real_Matrix; + Ld_B : Integer; + Beta : Real := 0.0; + C : in out Real_Matrix; + Ld_C : Integer) + is + begin + if Is_Single then + declare + subtype A_Type is BLAS.Real_Matrix (A'Range (1), A'Range (2)); + subtype B_Type is BLAS.Real_Matrix (B'Range (1), B'Range (2)); + type C_Ptr is + access all BLAS.Real_Matrix (C'Range (1), C'Range (2)); + function Conv_A is new Unchecked_Conversion (Real_Matrix, A_Type); + function Conv_B is new Unchecked_Conversion (Real_Matrix, B_Type); + function Conv_C is new Unchecked_Conversion (Address, C_Ptr); + begin + sgemm (Trans_A, Trans_B, M, N, K, Fortran.Real (Alpha), + Conv_A (A), Ld_A, Conv_B (B), Ld_B, Fortran.Real (Beta), + Conv_C (C'Address).all, Ld_C); + end; + + elsif Is_Double then + declare + subtype A_Type is + Double_Precision_Matrix (A'Range (1), A'Range (2)); + subtype B_Type is + Double_Precision_Matrix (B'Range (1), B'Range (2)); + type C_Ptr is + access all Double_Precision_Matrix (C'Range (1), C'Range (2)); + function Conv_A is new Unchecked_Conversion (Real_Matrix, A_Type); + function Conv_B is new Unchecked_Conversion (Real_Matrix, B_Type); + function Conv_C is new Unchecked_Conversion (Address, C_Ptr); + begin + dgemm (Trans_A, Trans_B, M, N, K, Double_Precision (Alpha), + Conv_A (A), Ld_A, Conv_B (B), Ld_B, Double_Precision (Beta), + Conv_C (C'Address).all, Ld_C); + end; + + else + declare + DP_C : Double_Precision_Matrix (C'Range (1), C'Range (2)); + begin + if Beta /= 0.0 then + DP_C := To_Double_Precision (C); + end if; + + dgemm (Trans_A, Trans_B, M, N, K, Double_Precision (Alpha), + To_Double_Precision (A), Ld_A, + To_Double_Precision (B), Ld_B, Double_Precision (Beta), + DP_C, Ld_C); + + C := To_Real (DP_C); + end; + end if; + end gemm; + + ---------- + -- gemv -- + ---------- + + procedure gemv + (Trans : access constant Character; + M : Natural := 0; + N : Natural := 0; + Alpha : Real := 1.0; + A : Real_Matrix; + Ld_A : Positive; + X : Real_Vector; + Inc_X : Integer := 1; + Beta : Real := 0.0; + Y : in out Real_Vector; + Inc_Y : Integer := 1) + is + begin + if Is_Single then + declare + subtype A_Type is BLAS.Real_Matrix (A'Range (1), A'Range (2)); + subtype X_Type is BLAS.Real_Vector (X'Range); + type Y_Ptr is access all BLAS.Real_Vector (Y'Range); + function Conv_A is new Unchecked_Conversion (Real_Matrix, A_Type); + function Conv_X is new Unchecked_Conversion (Real_Vector, X_Type); + function Conv_Y is new Unchecked_Conversion (Address, Y_Ptr); + begin + sgemv (Trans, M, N, Fortran.Real (Alpha), + Conv_A (A), Ld_A, Conv_X (X), Inc_X, Fortran.Real (Beta), + Conv_Y (Y'Address).all, Inc_Y); + end; + + elsif Is_Double then + declare + subtype A_Type is + Double_Precision_Matrix (A'Range (1), A'Range (2)); + subtype X_Type is Double_Precision_Vector (X'Range); + type Y_Ptr is access all Double_Precision_Vector (Y'Range); + function Conv_A is new Unchecked_Conversion (Real_Matrix, A_Type); + function Conv_X is new Unchecked_Conversion (Real_Vector, X_Type); + function Conv_Y is new Unchecked_Conversion (Address, Y_Ptr); + begin + dgemv (Trans, M, N, Double_Precision (Alpha), + Conv_A (A), Ld_A, Conv_X (X), Inc_X, + Double_Precision (Beta), + Conv_Y (Y'Address).all, Inc_Y); + end; + + else + declare + DP_Y : Double_Precision_Vector (Y'Range); + begin + if Beta /= 0.0 then + DP_Y := To_Double_Precision (Y); + end if; + + dgemv (Trans, M, N, Double_Precision (Alpha), + To_Double_Precision (A), Ld_A, + To_Double_Precision (X), Inc_X, Double_Precision (Beta), + DP_Y, Inc_Y); + + Y := To_Real (DP_Y); + end; + end if; + end gemv; + + ---------- + -- nrm2 -- + ---------- + + function nrm2 + (N : Natural; + X : Real_Vector; + Inc_X : Integer := 1) return Real + is + begin + if Is_Single then + declare + subtype X_Type is BLAS.Real_Vector (X'Range); + function Conv_X is new Unchecked_Conversion (Real_Vector, X_Type); + begin + return Real (snrm2 (N, Conv_X (X), Inc_X)); + end; + + elsif Is_Double then + declare + subtype X_Type is Double_Precision_Vector (X'Range); + function Conv_X is new Unchecked_Conversion (Real_Vector, X_Type); + begin + return Real (dnrm2 (N, Conv_X (X), Inc_X)); + end; + + else + return Real (dnrm2 (N, To_Double_Precision (X), Inc_X)); + end if; + end nrm2; + +end System.Generic_Real_BLAS; diff --git a/gcc/ada/s-gerebl.ads b/gcc/ada/s-gerebl.ads new file mode 100644 index 0000000..adb62eb --- /dev/null +++ b/gcc/ada/s-gerebl.ads @@ -0,0 +1,98 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- SYSTEM.GENERIC_REAL_BLAS -- +-- -- +-- S p e c -- +-- -- +-- Copyright (C) 2006, 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 2, 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. See the GNU General Public License -- +-- for more details. You should have received a copy of the GNU General -- +-- Public License distributed with GNAT; see file COPYING. If not, write -- +-- to the Free Software Foundation, 51 Franklin Street, Fifth Floor, -- +-- Boston, MA 02110-1301, USA. -- +-- -- +-- As a special exception, if other files instantiate generics from this -- +-- unit, or you link this unit with other files to produce an executable, -- +-- this unit does not by itself cause the resulting executable to be -- +-- covered by the GNU General Public License. This exception does not -- +-- however invalidate any other reasons why the executable file might be -- +-- covered by the GNU Public License. -- +-- -- +-- GNAT was originally developed by the GNAT team at New York University. -- +-- Extensive contributions were provided by Ada Core Technologies Inc. -- +-- -- +------------------------------------------------------------------------------ + +-- Package comment required ??? + +generic + type Real is digits <>; + type Real_Vector is array (Integer range <>) of Real; + type Real_Matrix is array (Integer range <>, Integer range <>) of Real; +package System.Generic_Real_BLAS is + pragma Pure; + + -- Although BLAS support is only available for IEEE single and double + -- compatible floating-point types, this unit will accept any type + -- and apply conversions as necessary, with possible loss of + -- precision and range. + + No_Trans : aliased constant Character := 'N'; + Trans : aliased constant Character := 'T'; + Conj_Trans : aliased constant Character := 'C'; + + -- BLAS Level 1 Subprograms and Types + + function dot + (N : Positive; + X : Real_Vector; + Inc_X : Integer := 1; + Y : Real_Vector; + Inc_Y : Integer := 1) return Real; + + function nrm2 + (N : Natural; + X : Real_Vector; + Inc_X : Integer := 1) return Real; + + procedure gemv + (Trans : access constant Character; + M : Natural := 0; + N : Natural := 0; + Alpha : Real := 1.0; + A : Real_Matrix; + Ld_A : Positive; + X : Real_Vector; + Inc_X : Integer := 1; -- must be non-zero + Beta : Real := 0.0; + Y : in out Real_Vector; + Inc_Y : Integer := 1); -- must be non-zero + + -- BLAS Level 3 + + -- gemm s, d, c, z Matrix-matrix product of general matrices + + procedure gemm + (Trans_A : access constant Character; + Trans_B : access constant Character; + M : Positive; + N : Positive; + K : Positive; + Alpha : Real := 1.0; + A : Real_Matrix; + Ld_A : Integer; + B : Real_Matrix; + Ld_B : Integer; + Beta : Real := 0.0; + C : in out Real_Matrix; + Ld_C : Integer); + +end System.Generic_Real_BLAS; diff --git a/gcc/ada/s-gerela.adb b/gcc/ada/s-gerela.adb new file mode 100644 index 0000000..aa0e964 --- /dev/null +++ b/gcc/ada/s-gerela.adb @@ -0,0 +1,566 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- SYSTEM.GENERIC_REAL_LAPACK -- +-- -- +-- B o d y -- +-- -- +-- Copyright (C) 2006, 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 2, 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. See the GNU General Public License -- +-- for more details. You should have received a copy of the GNU General -- +-- Public License distributed with GNAT; see file COPYING. If not, write -- +-- to the Free Software Foundation, 51 Franklin Street, Fifth Floor, -- +-- Boston, MA 02110-1301, USA. -- +-- -- +-- As a special exception, if other files instantiate generics from this -- +-- unit, or you link this unit with other files to produce an executable, -- +-- this unit does not by itself cause the resulting executable to be -- +-- covered by the GNU General Public License. This exception does not -- +-- however invalidate any other reasons why the executable file might be -- +-- covered by the GNU Public License. -- +-- -- +-- GNAT was originally developed by the GNAT team at New York University. -- +-- Extensive contributions were provided by Ada Core Technologies Inc. -- +-- -- +------------------------------------------------------------------------------ + +with Ada.Unchecked_Conversion; use Ada; +with Interfaces; use Interfaces; +with Interfaces.Fortran; use Interfaces.Fortran; +with Interfaces.Fortran.BLAS; use Interfaces.Fortran.BLAS; +with Interfaces.Fortran.LAPACK; use Interfaces.Fortran.LAPACK; +with System.Generic_Array_Operations; use System.Generic_Array_Operations; + +package body System.Generic_Real_LAPACK is + + Is_Real : constant Boolean := + Real'Machine_Mantissa = Fortran.Real'Machine_Mantissa + and then Fortran.Real (Real'First) = Fortran.Real'First + and then Fortran.Real (Real'Last) = Fortran.Real'Last; + + Is_Double_Precision : constant Boolean := + Real'Machine_Mantissa = + Double_Precision'Machine_Mantissa + and then + Double_Precision (Real'First) = + Double_Precision'First + and then + Double_Precision (Real'Last) = + Double_Precision'Last; + + -- Local subprograms + + function To_Double_Precision (X : Real) return Double_Precision; + pragma Inline_Always (To_Double_Precision); + + function To_Real (X : Double_Precision) return Real; + pragma Inline_Always (To_Real); + + -- Instantiations + + function To_Double_Precision is new + Vector_Elementwise_Operation + (X_Scalar => Real, + Result_Scalar => Double_Precision, + X_Vector => Real_Vector, + Result_Vector => Double_Precision_Vector, + Operation => To_Double_Precision); + + function To_Real is new + Vector_Elementwise_Operation + (X_Scalar => Double_Precision, + Result_Scalar => Real, + X_Vector => Double_Precision_Vector, + Result_Vector => Real_Vector, + Operation => To_Real); + + function To_Double_Precision is new + Matrix_Elementwise_Operation + (X_Scalar => Real, + Result_Scalar => Double_Precision, + X_Matrix => Real_Matrix, + Result_Matrix => Double_Precision_Matrix, + Operation => To_Double_Precision); + + function To_Real is new + Matrix_Elementwise_Operation + (X_Scalar => Double_Precision, + Result_Scalar => Real, + X_Matrix => Double_Precision_Matrix, + Result_Matrix => Real_Matrix, + Operation => To_Real); + + function To_Double_Precision (X : Real) return Double_Precision is + begin + return Double_Precision (X); + end To_Double_Precision; + + function To_Real (X : Double_Precision) return Real is + begin + return Real (X); + end To_Real; + + ----------- + -- getrf -- + ----------- + + procedure getrf + (M : Natural; + N : Natural; + A : in out Real_Matrix; + Ld_A : Positive; + I_Piv : out Integer_Vector; + Info : access Integer) + is + begin + if Is_Real then + declare + type A_Ptr is + access all BLAS.Real_Matrix (A'Range (1), A'Range (2)); + function Conv_A is new Unchecked_Conversion (Address, A_Ptr); + begin + sgetrf (M, N, Conv_A (A'Address).all, Ld_A, + LAPACK.Integer_Vector (I_Piv), Info); + end; + + elsif Is_Double_Precision then + declare + type A_Ptr is + access all Double_Precision_Matrix (A'Range (1), A'Range (2)); + function Conv_A is new Unchecked_Conversion (Address, A_Ptr); + begin + dgetrf (M, N, Conv_A (A'Address).all, Ld_A, + LAPACK.Integer_Vector (I_Piv), Info); + end; + + else + declare + DP_A : Double_Precision_Matrix (A'Range (1), A'Range (2)); + begin + DP_A := To_Double_Precision (A); + dgetrf (M, N, DP_A, Ld_A, LAPACK.Integer_Vector (I_Piv), Info); + A := To_Real (DP_A); + end; + end if; + end getrf; + + ----------- + -- getri -- + ----------- + + procedure getri + (N : Natural; + A : in out Real_Matrix; + Ld_A : Positive; + I_Piv : Integer_Vector; + Work : in out Real_Vector; + L_Work : Integer; + Info : access Integer) + is + begin + if Is_Real then + declare + type A_Ptr is + access all BLAS.Real_Matrix (A'Range (1), A'Range (2)); + type Work_Ptr is + access all BLAS.Real_Vector (Work'Range); + function Conv_A is new Unchecked_Conversion (Address, A_Ptr); + function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr); + begin + sgetri (N, Conv_A (A'Address).all, Ld_A, + LAPACK.Integer_Vector (I_Piv), + Conv_Work (Work'Address).all, L_Work, + Info); + end; + + elsif Is_Double_Precision then + declare + type A_Ptr is + access all Double_Precision_Matrix (A'Range (1), A'Range (2)); + type Work_Ptr is + access all Double_Precision_Vector (Work'Range); + function Conv_A is new Unchecked_Conversion (Address, A_Ptr); + function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr); + begin + dgetri (N, Conv_A (A'Address).all, Ld_A, + LAPACK.Integer_Vector (I_Piv), + Conv_Work (Work'Address).all, L_Work, + Info); + end; + + else + declare + DP_A : Double_Precision_Matrix (A'Range (1), A'Range (2)); + DP_Work : Double_Precision_Vector (Work'Range); + begin + DP_A := To_Double_Precision (A); + dgetri (N, DP_A, Ld_A, LAPACK.Integer_Vector (I_Piv), + DP_Work, L_Work, Info); + A := To_Real (DP_A); + Work (1) := To_Real (DP_Work (1)); + end; + end if; + end getri; + + ----------- + -- getrs -- + ----------- + + procedure getrs + (Trans : access constant Character; + N : Natural; + N_Rhs : Natural; + A : Real_Matrix; + Ld_A : Positive; + I_Piv : Integer_Vector; + B : in out Real_Matrix; + Ld_B : Positive; + Info : access Integer) + is + begin + if Is_Real then + declare + subtype A_Type is BLAS.Real_Matrix (A'Range (1), A'Range (2)); + type B_Ptr is + access all BLAS.Real_Matrix (B'Range (1), B'Range (2)); + function Conv_A is new Unchecked_Conversion (Real_Matrix, A_Type); + function Conv_B is new Unchecked_Conversion (Address, B_Ptr); + begin + sgetrs (Trans, N, N_Rhs, + Conv_A (A), Ld_A, + LAPACK.Integer_Vector (I_Piv), + Conv_B (B'Address).all, Ld_B, + Info); + end; + + elsif Is_Double_Precision then + declare + subtype A_Type is + Double_Precision_Matrix (A'Range (1), A'Range (2)); + type B_Ptr is + access all Double_Precision_Matrix (B'Range (1), B'Range (2)); + function Conv_A is new Unchecked_Conversion (Real_Matrix, A_Type); + function Conv_B is new Unchecked_Conversion (Address, B_Ptr); + begin + dgetrs (Trans, N, N_Rhs, + Conv_A (A), Ld_A, + LAPACK.Integer_Vector (I_Piv), + Conv_B (B'Address).all, Ld_B, + Info); + end; + + else + declare + DP_A : Double_Precision_Matrix (A'Range (1), A'Range (2)); + DP_B : Double_Precision_Matrix (B'Range (1), B'Range (2)); + begin + DP_A := To_Double_Precision (A); + DP_B := To_Double_Precision (B); + dgetrs (Trans, N, N_Rhs, + DP_A, Ld_A, + LAPACK.Integer_Vector (I_Piv), + DP_B, Ld_B, + Info); + B := To_Real (DP_B); + end; + end if; + end getrs; + + ----------- + -- orgtr -- + ----------- + + procedure orgtr + (Uplo : access constant Character; + N : Natural; + A : in out Real_Matrix; + Ld_A : Positive; + Tau : in Real_Vector; + Work : out Real_Vector; + L_Work : Integer; + Info : access Integer) + is + begin + if Is_Real then + declare + type A_Ptr is + access all BLAS.Real_Matrix (A'Range (1), A'Range (2)); + subtype Tau_Type is BLAS.Real_Vector (Tau'Range); + type Work_Ptr is + access all BLAS.Real_Vector (Work'Range); + function Conv_A is new Unchecked_Conversion (Address, A_Ptr); + function Conv_Tau is + new Unchecked_Conversion (Real_Vector, Tau_Type); + function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr); + begin + sorgtr (Uplo, N, + Conv_A (A'Address).all, Ld_A, + Conv_Tau (Tau), + Conv_Work (Work'Address).all, L_Work, + Info); + end; + + elsif Is_Double_Precision then + declare + type A_Ptr is + access all Double_Precision_Matrix (A'Range (1), A'Range (2)); + subtype Tau_Type is Double_Precision_Vector (Tau'Range); + type Work_Ptr is + access all Double_Precision_Vector (Work'Range); + function Conv_A is new Unchecked_Conversion (Address, A_Ptr); + function Conv_Tau is + new Unchecked_Conversion (Real_Vector, Tau_Type); + function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr); + begin + dorgtr (Uplo, N, + Conv_A (A'Address).all, Ld_A, + Conv_Tau (Tau), + Conv_Work (Work'Address).all, L_Work, + Info); + end; + + else + declare + DP_A : Double_Precision_Matrix (A'Range (1), A'Range (2)); + DP_Work : Double_Precision_Vector (Work'Range); + DP_Tau : Double_Precision_Vector (Tau'Range); + begin + DP_A := To_Double_Precision (A); + DP_Tau := To_Double_Precision (Tau); + dorgtr (Uplo, N, DP_A, Ld_A, DP_Tau, DP_Work, L_Work, Info); + A := To_Real (DP_A); + Work (1) := To_Real (DP_Work (1)); + end; + end if; + end orgtr; + + ----------- + -- steqr -- + ----------- + + procedure steqr + (Comp_Z : access constant Character; + N : Natural; + D : in out Real_Vector; + E : in out Real_Vector; + Z : in out Real_Matrix; + Ld_Z : Positive; + Work : out Real_Vector; + Info : access Integer) + is + begin + if Is_Real then + declare + type D_Ptr is access all BLAS.Real_Vector (D'Range); + type E_Ptr is access all BLAS.Real_Vector (E'Range); + type Z_Ptr is + access all BLAS.Real_Matrix (Z'Range (1), Z'Range (2)); + type Work_Ptr is + access all BLAS.Real_Vector (Work'Range); + function Conv_D is new Unchecked_Conversion (Address, D_Ptr); + function Conv_E is new Unchecked_Conversion (Address, E_Ptr); + function Conv_Z is new Unchecked_Conversion (Address, Z_Ptr); + function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr); + begin + ssteqr (Comp_Z, N, + Conv_D (D'Address).all, + Conv_E (E'Address).all, + Conv_Z (Z'Address).all, + Ld_Z, + Conv_Work (Work'Address).all, + Info); + end; + + elsif Is_Double_Precision then + declare + type D_Ptr is access all Double_Precision_Vector (D'Range); + type E_Ptr is access all Double_Precision_Vector (E'Range); + type Z_Ptr is + access all Double_Precision_Matrix (Z'Range (1), Z'Range (2)); + type Work_Ptr is + access all Double_Precision_Vector (Work'Range); + function Conv_D is new Unchecked_Conversion (Address, D_Ptr); + function Conv_E is new Unchecked_Conversion (Address, E_Ptr); + function Conv_Z is new Unchecked_Conversion (Address, Z_Ptr); + function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr); + begin + dsteqr (Comp_Z, N, + Conv_D (D'Address).all, + Conv_E (E'Address).all, + Conv_Z (Z'Address).all, + Ld_Z, + Conv_Work (Work'Address).all, + Info); + end; + + else + declare + DP_D : Double_Precision_Vector (D'Range); + DP_E : Double_Precision_Vector (E'Range); + DP_Z : Double_Precision_Matrix (Z'Range (1), Z'Range (2)); + DP_Work : Double_Precision_Vector (Work'Range); + begin + DP_D := To_Double_Precision (D); + DP_E := To_Double_Precision (E); + + if Comp_Z.all = 'V' then + DP_Z := To_Double_Precision (Z); + end if; + + dsteqr (Comp_Z, N, DP_D, DP_E, DP_Z, Ld_Z, DP_Work, Info); + + D := To_Real (DP_D); + E := To_Real (DP_E); + Z := To_Real (DP_Z); + end; + end if; + end steqr; + + ----------- + -- sterf -- + ----------- + + procedure sterf + (N : Natural; + D : in out Real_Vector; + E : in out Real_Vector; + Info : access Integer) + is + begin + if Is_Real then + declare + type D_Ptr is access all BLAS.Real_Vector (D'Range); + type E_Ptr is access all BLAS.Real_Vector (E'Range); + function Conv_D is new Unchecked_Conversion (Address, D_Ptr); + function Conv_E is new Unchecked_Conversion (Address, E_Ptr); + begin + ssterf (N, Conv_D (D'Address).all, Conv_E (E'Address).all, Info); + end; + + elsif Is_Double_Precision then + declare + type D_Ptr is access all Double_Precision_Vector (D'Range); + type E_Ptr is access all Double_Precision_Vector (E'Range); + function Conv_D is new Unchecked_Conversion (Address, D_Ptr); + function Conv_E is new Unchecked_Conversion (Address, E_Ptr); + begin + dsterf (N, Conv_D (D'Address).all, Conv_E (E'Address).all, Info); + end; + + else + declare + DP_D : Double_Precision_Vector (D'Range); + DP_E : Double_Precision_Vector (E'Range); + + begin + DP_D := To_Double_Precision (D); + DP_E := To_Double_Precision (E); + + dsterf (N, DP_D, DP_E, Info); + + D := To_Real (DP_D); + E := To_Real (DP_E); + end; + end if; + end sterf; + + ----------- + -- sytrd -- + ----------- + + procedure sytrd + (Uplo : access constant Character; + N : Natural; + A : in out Real_Matrix; + Ld_A : Positive; + D : out Real_Vector; + E : out Real_Vector; + Tau : out Real_Vector; + Work : out Real_Vector; + L_Work : Integer; + Info : access Integer) + is + begin + if Is_Real then + declare + type A_Ptr is + access all BLAS.Real_Matrix (A'Range (1), A'Range (2)); + type D_Ptr is access all BLAS.Real_Vector (D'Range); + type E_Ptr is access all BLAS.Real_Vector (E'Range); + type Tau_Ptr is access all BLAS.Real_Vector (Tau'Range); + type Work_Ptr is + access all BLAS.Real_Vector (Work'Range); + function Conv_A is new Unchecked_Conversion (Address, A_Ptr); + function Conv_D is new Unchecked_Conversion (Address, D_Ptr); + function Conv_E is new Unchecked_Conversion (Address, E_Ptr); + function Conv_Tau is new Unchecked_Conversion (Address, Tau_Ptr); + function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr); + begin + ssytrd (Uplo, N, + Conv_A (A'Address).all, Ld_A, + Conv_D (D'Address).all, + Conv_E (E'Address).all, + Conv_Tau (Tau'Address).all, + Conv_Work (Work'Address).all, + L_Work, + Info); + end; + + elsif Is_Double_Precision then + declare + type A_Ptr is + access all Double_Precision_Matrix (A'Range (1), A'Range (2)); + type D_Ptr is access all Double_Precision_Vector (D'Range); + type E_Ptr is access all Double_Precision_Vector (E'Range); + type Tau_Ptr is access all Double_Precision_Vector (Tau'Range); + type Work_Ptr is + access all Double_Precision_Vector (Work'Range); + function Conv_A is new Unchecked_Conversion (Address, A_Ptr); + function Conv_D is new Unchecked_Conversion (Address, D_Ptr); + function Conv_E is new Unchecked_Conversion (Address, E_Ptr); + function Conv_Tau is new Unchecked_Conversion (Address, Tau_Ptr); + function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr); + begin + dsytrd (Uplo, N, + Conv_A (A'Address).all, Ld_A, + Conv_D (D'Address).all, + Conv_E (E'Address).all, + Conv_Tau (Tau'Address).all, + Conv_Work (Work'Address).all, + L_Work, + Info); + end; + + else + declare + DP_A : Double_Precision_Matrix (A'Range (1), A'Range (2)); + DP_D : Double_Precision_Vector (D'Range); + DP_E : Double_Precision_Vector (E'Range); + DP_Tau : Double_Precision_Vector (Tau'Range); + DP_Work : Double_Precision_Vector (Work'Range); + begin + DP_A := To_Double_Precision (A); + + dsytrd (Uplo, N, DP_A, Ld_A, DP_D, DP_E, DP_Tau, + DP_Work, L_Work, Info); + + if L_Work /= -1 then + A := To_Real (DP_A); + D := To_Real (DP_D); + E := To_Real (DP_E); + Tau := To_Real (DP_Tau); + end if; + + Work (1) := To_Real (DP_Work (1)); + end; + end if; + end sytrd; + +end System.Generic_Real_LAPACK; diff --git a/gcc/ada/s-gerela.ads b/gcc/ada/s-gerela.ads new file mode 100644 index 0000000..c9f29f2 --- /dev/null +++ b/gcc/ada/s-gerela.ads @@ -0,0 +1,130 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- SYSTEM.GENERIC_REAL_LAPACK -- +-- -- +-- S p e c -- +-- -- +-- Copyright (C) 2006, 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 2, 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. See the GNU General Public License -- +-- for more details. You should have received a copy of the GNU General -- +-- Public License distributed with GNAT; see file COPYING. If not, write -- +-- to the Free Software Foundation, 51 Franklin Street, Fifth Floor, -- +-- Boston, MA 02110-1301, USA. -- +-- -- +-- As a special exception, if other files instantiate generics from this -- +-- unit, or you link this unit with other files to produce an executable, -- +-- this unit does not by itself cause the resulting executable to be -- +-- covered by the GNU General Public License. This exception does not -- +-- however invalidate any other reasons why the executable file might be -- +-- covered by the GNU Public License. -- +-- -- +-- GNAT was originally developed by the GNAT team at New York University. -- +-- Extensive contributions were provided by Ada Core Technologies Inc. -- +-- -- +------------------------------------------------------------------------------ + +-- Package comment required ??? + +generic + type Real is digits <>; + type Real_Vector is array (Integer range <>) of Real; + type Real_Matrix is array (Integer range <>, Integer range <>) of Real; +package System.Generic_Real_LAPACK is + pragma Pure; + + type Integer_Vector is array (Integer range <>) of Integer; + + Upper : aliased constant Character := 'U'; + Lower : aliased constant Character := 'L'; + + -- LAPACK Computational Routines + + -- gerfs Refines the solution of a system of linear equations with + -- a general matrix and estimates its error + -- getrf Computes LU factorization of a general m-by-n matrix + -- getri Computes inverse of an LU-factored general matrix + -- square matrix, with multiple right-hand sides + -- getrs Solves a system of linear equations with an LU-factored + -- square matrix, with multiple right-hand sides + -- orgtr Generates the Float orthogonal matrix Q determined by sytrd + -- steqr Computes all eigenvalues and eigenvectors of a symmetric or + -- Hermitian matrix reduced to tridiagonal form (QR algorithm) + -- sterf Computes all eigenvalues of a Float symmetric + -- tridiagonal matrix using QR algorithm + -- sytrd Reduces a Float symmetric matrix to tridiagonal form + + procedure getrf + (M : Natural; + N : Natural; + A : in out Real_Matrix; + Ld_A : Positive; + I_Piv : out Integer_Vector; + Info : access Integer); + + procedure getri + (N : Natural; + A : in out Real_Matrix; + Ld_A : Positive; + I_Piv : Integer_Vector; + Work : in out Real_Vector; + L_Work : Integer; + Info : access Integer); + + procedure getrs + (Trans : access constant Character; + N : Natural; + N_Rhs : Natural; + A : Real_Matrix; + Ld_A : Positive; + I_Piv : Integer_Vector; + B : in out Real_Matrix; + Ld_B : Positive; + Info : access Integer); + + procedure orgtr + (Uplo : access constant Character; + N : Natural; + A : in out Real_Matrix; + Ld_A : Positive; + Tau : in Real_Vector; + Work : out Real_Vector; + L_Work : Integer; + Info : access Integer); + + procedure sterf + (N : Natural; + D : in out Real_Vector; + E : in out Real_Vector; + Info : access Integer); + + procedure steqr + (Comp_Z : access constant Character; + N : Natural; + D : in out Real_Vector; + E : in out Real_Vector; + Z : in out Real_Matrix; + Ld_Z : Positive; + Work : out Real_Vector; + Info : access Integer); + + procedure sytrd + (Uplo : access constant Character; + N : Natural; + A : in out Real_Matrix; + Ld_A : Positive; + D : out Real_Vector; + E : out Real_Vector; + Tau : out Real_Vector; + Work : out Real_Vector; + L_Work : Integer; + Info : access Integer); + +end System.Generic_Real_LAPACK; -- cgit v1.1