i-fortra.ads: Add Double_Complex type.

2007-04-06  Geert Bosch  <bosch@adacore.com>
	    Robert Dewar  <dewar@adacore.com>

	* 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

@@ -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) \

gcc/ada/a-fzteio.ads

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

gcc/ada/a-izteio.ads

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

gcc/ada/a-ngcoar.adb

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

gcc/ada/a-ngcoar.ads

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

gcc/ada/a-ngrear.adb

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

gcc/ada/a-ngrear.ads

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

gcc/ada/a-nlcoar.ads

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

gcc/ada/a-nllcar.ads

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

gcc/ada/a-nllrar.ads

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

gcc/ada/a-nlrear.ads

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

gcc/ada/a-nucoar.ads

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

gcc/ada/a-nurear.ads

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

gcc/ada/g-sha1.adb

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

gcc/ada/g-sha1.ads

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

gcc/ada/i-forbla.ads

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

gcc/ada/i-forlap.ads

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

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;

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;

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;

gcc/ada/s-gearop.adb

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

gcc/ada/s-gearop.ads

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

gcc/ada/s-gecobl.adb

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

gcc/ada/s-gecobl.ads

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

gcc/ada/s-gecola.adb

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

gcc/ada/s-gecola.ads

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

gcc/ada/s-gerebl.adb

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

gcc/ada/s-gerebl.ads

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

gcc/ada/s-gerela.adb

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

gcc/ada/s-gerela.ads

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