------------------------------------------------------------------------------ -- -- -- GNAT RUN-TIME COMPONENTS -- -- -- -- S Y S T E M . G E N E R I C _ C O M P L E X _ L A P A C K -- -- -- -- B o d y -- -- -- -- Copyright (C) 2006-2007, 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;