| ------------------------------------------------------------------------------ |
| -- -- |
| -- GNAT RUN-TIME COMPONENTS -- |
| -- -- |
| -- S Y S T E M . G E N E R I C _ C O M P L E X _ B L A S -- |
| -- -- |
| -- 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 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.cdotu (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.zdotu (N, Conv_X (X'Address).all, Inc_X, |
| Conv_Y (Y'Address).all, Inc_Y)); |
| end; |
| |
| else |
| return To_Complex (BLAS.zdotu (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, 0.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, 0.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; |