diff options
Diffstat (limited to 'gcc-4.4.3/gcc/ada/a-ngcoar.adb')
| -rw-r--r-- | gcc-4.4.3/gcc/ada/a-ngcoar.adb | 1502 |
1 files changed, 0 insertions, 1502 deletions
diff --git a/gcc-4.4.3/gcc/ada/a-ngcoar.adb b/gcc-4.4.3/gcc/ada/a-ngcoar.adb deleted file mode 100644 index 9979d6bae..000000000 --- a/gcc-4.4.3/gcc/ada/a-ngcoar.adb +++ /dev/null @@ -1,1502 +0,0 @@ ------------------------------------------------------------------------------- --- -- --- GNAT COMPILER COMPONENTS -- --- -- --- ADA.NUMERICS.GENERIC_COMPLEX_ARRAYS -- --- -- --- B o d y -- --- -- --- Copyright (C) 2006-2009, Free Software Foundation, Inc. -- --- -- --- GNAT is free software; you can redistribute it and/or modify it under -- --- terms of the GNU General Public License as published by the Free Soft- -- --- ware Foundation; either version 3, or (at your option) any later ver- -- --- sion. GNAT is distributed in the hope that it will be useful, but WITH- -- --- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY -- --- or FITNESS FOR A PARTICULAR PURPOSE. -- --- -- --- As a special exception under Section 7 of GPL version 3, you are granted -- --- additional permissions described in the GCC Runtime Library Exception, -- --- version 3.1, as published by the Free Software Foundation. -- --- -- --- You should have received a copy of the GNU General Public License and -- --- a copy of the GCC Runtime Library Exception along with this program; -- --- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see -- --- <http://www.gnu.org/licenses/>. -- --- -- --- GNAT was originally developed by the GNAT team at New York University. -- --- Extensive contributions were provided by Ada Core Technologies Inc. -- --- -- ------------------------------------------------------------------------------- - -with 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); - - subtype Real is Real_Arrays.Real; - -- Work around visibility bug ??? - - 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 multiplication"; - end if; - - gemm (Trans_A => No_Trans'Access, - Trans_B => No_Trans'Access, - M => Right'Length (2), - N => Left'Length (1), - K => Right'Length (1), - A => Right, - Ld_A => Right'Length (2), - B => Left, - Ld_B => Left'Length (2), - C => R, - Ld_C => R'Length (2)); - - return R; - end "*"; - - function "*" - (Left : 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 : 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-Hermitian 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; |