1 ------------------------------------------------------------------------------
3 -- GNAT RUN-TIME COMPONENTS --
5 -- S Y S T E M . G E N E R I C _ C O M P L E X _ B L A S --
9 -- Copyright (C) 2006-2007, Free Software Foundation, Inc. --
11 -- GNAT is free software; you can redistribute it and/or modify it under --
12 -- terms of the GNU General Public License as published by the Free Soft- --
13 -- ware Foundation; either version 2, or (at your option) any later ver- --
14 -- sion. GNAT is distributed in the hope that it will be useful, but WITH- --
15 -- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY --
16 -- or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License --
17 -- for more details. You should have received a copy of the GNU General --
18 -- Public License distributed with GNAT; see file COPYING. If not, write --
19 -- to the Free Software Foundation, 51 Franklin Street, Fifth Floor, --
20 -- Boston, MA 02110-1301, USA. --
22 -- As a special exception, if other files instantiate generics from this --
23 -- unit, or you link this unit with other files to produce an executable, --
24 -- this unit does not by itself cause the resulting executable to be --
25 -- covered by the GNU General Public License. This exception does not --
26 -- however invalidate any other reasons why the executable file might be --
27 -- covered by the GNU Public License. --
29 -- GNAT was originally developed by the GNAT team at New York University. --
30 -- Extensive contributions were provided by Ada Core Technologies Inc. --
32 ------------------------------------------------------------------------------
34 with Ada.Unchecked_Conversion; use Ada;
35 with Interfaces; use Interfaces;
36 with Interfaces.Fortran; use Interfaces.Fortran;
37 with Interfaces.Fortran.BLAS; use Interfaces.Fortran.BLAS;
38 with System.Generic_Array_Operations; use System.Generic_Array_Operations;
40 package body System.Generic_Complex_BLAS is
42 Is_Single : constant Boolean :=
43 Real'Machine_Mantissa = Fortran.Real'Machine_Mantissa
44 and then Fortran.Real (Real'First) = Fortran.Real'First
45 and then Fortran.Real (Real'Last) = Fortran.Real'Last;
47 Is_Double : constant Boolean :=
48 Real'Machine_Mantissa = Double_Precision'Machine_Mantissa
50 Double_Precision (Real'First) = Double_Precision'First
52 Double_Precision (Real'Last) = Double_Precision'Last;
54 subtype Complex is Complex_Types.Complex;
58 function To_Double_Precision (X : Real) return Double_Precision;
59 pragma Inline (To_Double_Precision);
61 function To_Double_Complex (X : Complex) return Double_Complex;
62 pragma Inline (To_Double_Complex);
64 function To_Complex (X : Double_Complex) return Complex;
65 function To_Complex (X : Fortran.Complex) return Complex;
66 pragma Inline (To_Complex);
68 function To_Fortran (X : Complex) return Fortran.Complex;
69 pragma Inline (To_Fortran);
73 function To_Double_Complex is new
74 Vector_Elementwise_Operation
75 (X_Scalar => Complex_Types.Complex,
76 Result_Scalar => Fortran.Double_Complex,
77 X_Vector => Complex_Vector,
78 Result_Vector => BLAS.Double_Complex_Vector,
79 Operation => To_Double_Complex);
81 function To_Complex is new
82 Vector_Elementwise_Operation
83 (X_Scalar => Fortran.Double_Complex,
84 Result_Scalar => Complex,
85 X_Vector => BLAS.Double_Complex_Vector,
86 Result_Vector => Complex_Vector,
87 Operation => To_Complex);
89 function To_Double_Complex is new
90 Matrix_Elementwise_Operation
92 Result_Scalar => Double_Complex,
93 X_Matrix => Complex_Matrix,
94 Result_Matrix => BLAS.Double_Complex_Matrix,
95 Operation => To_Double_Complex);
97 function To_Complex is new
98 Matrix_Elementwise_Operation
99 (X_Scalar => Double_Complex,
100 Result_Scalar => Complex,
101 X_Matrix => BLAS.Double_Complex_Matrix,
102 Result_Matrix => Complex_Matrix,
103 Operation => To_Complex);
105 function To_Double_Precision (X : Real) return Double_Precision is
107 return Double_Precision (X);
108 end To_Double_Precision;
110 function To_Double_Complex (X : Complex) return Double_Complex is
112 return (To_Double_Precision (X.Re), To_Double_Precision (X.Im));
113 end To_Double_Complex;
115 function To_Complex (X : Double_Complex) return Complex is
117 return (Real (X.Re), Real (X.Im));
120 function To_Complex (X : Fortran.Complex) return Complex is
122 return (Real (X.Re), Real (X.Im));
125 function To_Fortran (X : Complex) return Fortran.Complex is
127 return (Fortran.Real (X.Re), Fortran.Real (X.Im));
137 Inc_X : Integer := 1;
139 Inc_Y : Integer := 1) return Complex
144 type X_Ptr is access all BLAS.Complex_Vector (X'Range);
145 type Y_Ptr is access all BLAS.Complex_Vector (Y'Range);
146 function Conv_X is new Unchecked_Conversion (Address, X_Ptr);
147 function Conv_Y is new Unchecked_Conversion (Address, Y_Ptr);
149 return To_Complex (BLAS.cdotu (N, Conv_X (X'Address).all, Inc_X,
150 Conv_Y (Y'Address).all, Inc_Y));
155 type X_Ptr is access all BLAS.Double_Complex_Vector (X'Range);
156 type Y_Ptr is access all BLAS.Double_Complex_Vector (Y'Range);
157 function Conv_X is new Unchecked_Conversion (Address, X_Ptr);
158 function Conv_Y is new Unchecked_Conversion (Address, Y_Ptr);
160 return To_Complex (BLAS.zdotu (N, Conv_X (X'Address).all, Inc_X,
161 Conv_Y (Y'Address).all, Inc_Y));
165 return To_Complex (BLAS.zdotu (N, To_Double_Complex (X), Inc_X,
166 To_Double_Complex (Y), Inc_Y));
175 (Trans_A : access constant Character;
176 Trans_B : access constant Character;
180 Alpha : Complex := (1.0, 0.0);
185 Beta : Complex := (0.0, 0.0);
186 C : in out Complex_Matrix;
192 subtype A_Type is BLAS.Complex_Matrix (A'Range (1), A'Range (2));
193 subtype B_Type is BLAS.Complex_Matrix (B'Range (1), B'Range (2));
195 access all BLAS.Complex_Matrix (C'Range (1), C'Range (2));
197 new Unchecked_Conversion (Complex_Matrix, A_Type);
199 new Unchecked_Conversion (Complex_Matrix, B_Type);
201 new Unchecked_Conversion (Address, C_Ptr);
203 BLAS.cgemm (Trans_A, Trans_B, M, N, K, To_Fortran (Alpha),
204 Conv_A (A), Ld_A, Conv_B (B), Ld_B, To_Fortran (Beta),
205 Conv_C (C'Address).all, Ld_C);
211 BLAS.Double_Complex_Matrix (A'Range (1), A'Range (2));
213 BLAS.Double_Complex_Matrix (B'Range (1), B'Range (2));
214 type C_Ptr is access all
215 BLAS.Double_Complex_Matrix (C'Range (1), C'Range (2));
217 new Unchecked_Conversion (Complex_Matrix, A_Type);
219 new Unchecked_Conversion (Complex_Matrix, B_Type);
220 function Conv_C is new Unchecked_Conversion (Address, C_Ptr);
222 BLAS.zgemm (Trans_A, Trans_B, M, N, K, To_Double_Complex (Alpha),
223 Conv_A (A), Ld_A, Conv_B (B), Ld_B,
224 To_Double_Complex (Beta),
225 Conv_C (C'Address).all, Ld_C);
230 DP_C : BLAS.Double_Complex_Matrix (C'Range (1), C'Range (2));
232 if Beta.Re /= 0.0 or else Beta.Im /= 0.0 then
233 DP_C := To_Double_Complex (C);
236 BLAS.zgemm (Trans_A, Trans_B, M, N, K, To_Double_Complex (Alpha),
237 To_Double_Complex (A), Ld_A,
238 To_Double_Complex (B), Ld_B, To_Double_Complex (Beta),
241 C := To_Complex (DP_C);
251 (Trans : access constant Character;
254 Alpha : Complex := (1.0, 0.0);
258 Inc_X : Integer := 1;
259 Beta : Complex := (0.0, 0.0);
260 Y : in out Complex_Vector;
261 Inc_Y : Integer := 1)
266 subtype A_Type is BLAS.Complex_Matrix (A'Range (1), A'Range (2));
267 subtype X_Type is BLAS.Complex_Vector (X'Range);
268 type Y_Ptr is access all BLAS.Complex_Vector (Y'Range);
270 new Unchecked_Conversion (Complex_Matrix, A_Type);
272 new Unchecked_Conversion (Complex_Vector, X_Type);
274 new Unchecked_Conversion (Address, Y_Ptr);
276 BLAS.cgemv (Trans, M, N, To_Fortran (Alpha),
277 Conv_A (A), Ld_A, Conv_X (X), Inc_X, To_Fortran (Beta),
278 Conv_Y (Y'Address).all, Inc_Y);
284 BLAS.Double_Complex_Matrix (A'Range (1), A'Range (2));
286 BLAS.Double_Complex_Vector (X'Range);
287 type Y_Ptr is access all BLAS.Double_Complex_Vector (Y'Range);
289 new Unchecked_Conversion (Complex_Matrix, A_Type);
291 new Unchecked_Conversion (Complex_Vector, X_Type);
293 new Unchecked_Conversion (Address, Y_Ptr);
295 BLAS.zgemv (Trans, M, N, To_Double_Complex (Alpha),
296 Conv_A (A), Ld_A, Conv_X (X), Inc_X,
297 To_Double_Complex (Beta),
298 Conv_Y (Y'Address).all, Inc_Y);
303 DP_Y : BLAS.Double_Complex_Vector (Y'Range);
305 if Beta.Re /= 0.0 or else Beta.Im /= 0.0 then
306 DP_Y := To_Double_Complex (Y);
309 BLAS.zgemv (Trans, M, N, To_Double_Complex (Alpha),
310 To_Double_Complex (A), Ld_A,
311 To_Double_Complex (X), Inc_X, To_Double_Complex (Beta),
314 Y := To_Complex (DP_Y);
326 Inc_X : Integer := 1) return Real
331 subtype X_Type is BLAS.Complex_Vector (X'Range);
333 new Unchecked_Conversion (Complex_Vector, X_Type);
335 return Real (BLAS.scnrm2 (N, Conv_X (X), Inc_X));
340 subtype X_Type is BLAS.Double_Complex_Vector (X'Range);
342 new Unchecked_Conversion (Complex_Vector, X_Type);
344 return Real (BLAS.dznrm2 (N, Conv_X (X), Inc_X));
348 return Real (BLAS.dznrm2 (N, To_Double_Complex (X), Inc_X));
352 end System.Generic_Complex_BLAS;