# Objects needed only for tasking
GNATRTL_TASKING_OBJS= \
- a-diroro$(objext) \
a-dispat$(objext) \
a-dynpri$(objext) \
a-interr$(objext) \
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) \
g-bubsor$(objext) \
g-busora$(objext) \
g-busorg$(objext) \
+ g-bytswa$(objext) \
g-calend$(objext) \
g-casuti$(objext) \
g-catiio$(objext) \
g-regexp$(objext) \
g-regpat$(objext) \
g-sestin$(objext) \
+ g-sha1$(objext) \
g-soccon$(objext) \
g-socket$(objext) \
g-socthi$(objext) \
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- 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);
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- 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);
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- 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;
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- 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;
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- 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;
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- 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;
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- 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);
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- 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);
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- 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);
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- 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);
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- 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);
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- 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);
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- 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;
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- 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;
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- 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;
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- 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;
-- --
-- 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. --
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;
-- --
------------------------------------------------------------------------------
-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
"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
"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
"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
"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
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;
-- --
-- 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- --
-- 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;
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- 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;
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- 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;
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- 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;
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- 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;
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- 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;
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- 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;
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- 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;
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- 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;
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- 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;
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- 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;
projects / platform / upstream / gcc.git / commitdiff