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