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-2009, 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 3, 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. --
18 -- As a special exception under Section 7 of GPL version 3, you are granted --
19 -- additional permissions described in the GCC Runtime Library Exception, --
20 -- version 3.1, as published by the Free Software Foundation. --
22 -- You should have received a copy of the GNU General Public License and --
23 -- a copy of the GCC Runtime Library Exception along with this program; --
24 -- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see --
25 -- <http://www.gnu.org/licenses/>. --
27 -- GNAT was originally developed by the GNAT team at New York University. --
28 -- Extensive contributions were provided by Ada Core Technologies Inc. --
30 ------------------------------------------------------------------------------
32 with Ada
.Unchecked_Conversion
; use Ada
;
33 with Interfaces
; use Interfaces
;
34 with Interfaces
.Fortran
; use Interfaces
.Fortran
;
35 with Interfaces
.Fortran
.BLAS
; use Interfaces
.Fortran
.BLAS
;
36 with System
.Generic_Array_Operations
; use System
.Generic_Array_Operations
;
38 package body System
.Generic_Complex_BLAS
is
40 Is_Single
: constant Boolean :=
41 Real
'Machine_Mantissa = Fortran
.Real
'Machine_Mantissa
42 and then Fortran
.Real
(Real
'First) = Fortran
.Real
'First
43 and then Fortran
.Real
(Real
'Last) = Fortran
.Real
'Last;
45 Is_Double
: constant Boolean :=
46 Real
'Machine_Mantissa = Double_Precision
'Machine_Mantissa
48 Double_Precision
(Real
'First) = Double_Precision
'First
50 Double_Precision
(Real
'Last) = Double_Precision
'Last;
52 subtype Complex
is Complex_Types
.Complex
;
56 function To_Double_Precision
(X
: Real
) return Double_Precision
;
57 pragma Inline
(To_Double_Precision
);
59 function To_Double_Complex
(X
: Complex
) return Double_Complex
;
60 pragma Inline
(To_Double_Complex
);
62 function To_Complex
(X
: Double_Complex
) return Complex
;
63 function To_Complex
(X
: Fortran
.Complex
) return Complex
;
64 pragma Inline
(To_Complex
);
66 function To_Fortran
(X
: Complex
) return Fortran
.Complex
;
67 pragma Inline
(To_Fortran
);
71 function To_Double_Complex
is new
72 Vector_Elementwise_Operation
73 (X_Scalar
=> Complex_Types
.Complex
,
74 Result_Scalar
=> Fortran
.Double_Complex
,
75 X_Vector
=> Complex_Vector
,
76 Result_Vector
=> BLAS
.Double_Complex_Vector
,
77 Operation
=> To_Double_Complex
);
79 function To_Complex
is new
80 Vector_Elementwise_Operation
81 (X_Scalar
=> Fortran
.Double_Complex
,
82 Result_Scalar
=> Complex
,
83 X_Vector
=> BLAS
.Double_Complex_Vector
,
84 Result_Vector
=> Complex_Vector
,
85 Operation
=> To_Complex
);
87 function To_Double_Complex
is new
88 Matrix_Elementwise_Operation
90 Result_Scalar
=> Double_Complex
,
91 X_Matrix
=> Complex_Matrix
,
92 Result_Matrix
=> BLAS
.Double_Complex_Matrix
,
93 Operation
=> To_Double_Complex
);
95 function To_Complex
is new
96 Matrix_Elementwise_Operation
97 (X_Scalar
=> Double_Complex
,
98 Result_Scalar
=> Complex
,
99 X_Matrix
=> BLAS
.Double_Complex_Matrix
,
100 Result_Matrix
=> Complex_Matrix
,
101 Operation
=> To_Complex
);
103 function To_Double_Precision
(X
: Real
) return Double_Precision
is
105 return Double_Precision
(X
);
106 end To_Double_Precision
;
108 function To_Double_Complex
(X
: Complex
) return Double_Complex
is
110 return (To_Double_Precision
(X
.Re
), To_Double_Precision
(X
.Im
));
111 end To_Double_Complex
;
113 function To_Complex
(X
: Double_Complex
) return Complex
is
115 return (Real
(X
.Re
), Real
(X
.Im
));
118 function To_Complex
(X
: Fortran
.Complex
) return Complex
is
120 return (Real
(X
.Re
), Real
(X
.Im
));
123 function To_Fortran
(X
: Complex
) return Fortran
.Complex
is
125 return (Fortran
.Real
(X
.Re
), Fortran
.Real
(X
.Im
));
135 Inc_X
: Integer := 1;
137 Inc_Y
: Integer := 1) return Complex
142 type X_Ptr
is access all BLAS
.Complex_Vector
(X
'Range);
143 type Y_Ptr
is access all BLAS
.Complex_Vector
(Y
'Range);
144 function Conv_X
is new Unchecked_Conversion
(Address
, X_Ptr
);
145 function Conv_Y
is new Unchecked_Conversion
(Address
, Y_Ptr
);
147 return To_Complex
(BLAS
.cdotu
(N
, Conv_X
(X
'Address).all, Inc_X
,
148 Conv_Y
(Y
'Address).all, Inc_Y
));
153 type X_Ptr
is access all BLAS
.Double_Complex_Vector
(X
'Range);
154 type Y_Ptr
is access all BLAS
.Double_Complex_Vector
(Y
'Range);
155 function Conv_X
is new Unchecked_Conversion
(Address
, X_Ptr
);
156 function Conv_Y
is new Unchecked_Conversion
(Address
, Y_Ptr
);
158 return To_Complex
(BLAS
.zdotu
(N
, Conv_X
(X
'Address).all, Inc_X
,
159 Conv_Y
(Y
'Address).all, Inc_Y
));
163 return To_Complex
(BLAS
.zdotu
(N
, To_Double_Complex
(X
), Inc_X
,
164 To_Double_Complex
(Y
), Inc_Y
));
173 (Trans_A
: access constant Character;
174 Trans_B
: access constant Character;
178 Alpha
: Complex
:= (1.0, 0.0);
183 Beta
: Complex
:= (0.0, 0.0);
184 C
: in out Complex_Matrix
;
190 subtype A_Type
is BLAS
.Complex_Matrix
(A
'Range (1), A
'Range (2));
191 subtype B_Type
is BLAS
.Complex_Matrix
(B
'Range (1), B
'Range (2));
193 access all BLAS
.Complex_Matrix
(C
'Range (1), C
'Range (2));
195 new Unchecked_Conversion
(Complex_Matrix
, A_Type
);
197 new Unchecked_Conversion
(Complex_Matrix
, B_Type
);
199 new Unchecked_Conversion
(Address
, C_Ptr
);
201 BLAS
.cgemm
(Trans_A
, Trans_B
, M
, N
, K
, To_Fortran
(Alpha
),
202 Conv_A
(A
), Ld_A
, Conv_B
(B
), Ld_B
, To_Fortran
(Beta
),
203 Conv_C
(C
'Address).all, Ld_C
);
209 BLAS
.Double_Complex_Matrix
(A
'Range (1), A
'Range (2));
211 BLAS
.Double_Complex_Matrix
(B
'Range (1), B
'Range (2));
212 type C_Ptr
is access all
213 BLAS
.Double_Complex_Matrix
(C
'Range (1), C
'Range (2));
215 new Unchecked_Conversion
(Complex_Matrix
, A_Type
);
217 new Unchecked_Conversion
(Complex_Matrix
, B_Type
);
218 function Conv_C
is new Unchecked_Conversion
(Address
, C_Ptr
);
220 BLAS
.zgemm
(Trans_A
, Trans_B
, M
, N
, K
, To_Double_Complex
(Alpha
),
221 Conv_A
(A
), Ld_A
, Conv_B
(B
), Ld_B
,
222 To_Double_Complex
(Beta
),
223 Conv_C
(C
'Address).all, Ld_C
);
228 DP_C
: BLAS
.Double_Complex_Matrix
(C
'Range (1), C
'Range (2));
230 if Beta
.Re
/= 0.0 or else Beta
.Im
/= 0.0 then
231 DP_C
:= To_Double_Complex
(C
);
234 BLAS
.zgemm
(Trans_A
, Trans_B
, M
, N
, K
, To_Double_Complex
(Alpha
),
235 To_Double_Complex
(A
), Ld_A
,
236 To_Double_Complex
(B
), Ld_B
, To_Double_Complex
(Beta
),
239 C
:= To_Complex
(DP_C
);
249 (Trans
: access constant Character;
252 Alpha
: Complex
:= (1.0, 0.0);
256 Inc_X
: Integer := 1;
257 Beta
: Complex
:= (0.0, 0.0);
258 Y
: in out Complex_Vector
;
259 Inc_Y
: Integer := 1)
264 subtype A_Type
is BLAS
.Complex_Matrix
(A
'Range (1), A
'Range (2));
265 subtype X_Type
is BLAS
.Complex_Vector
(X
'Range);
266 type Y_Ptr
is access all BLAS
.Complex_Vector
(Y
'Range);
268 new Unchecked_Conversion
(Complex_Matrix
, A_Type
);
270 new Unchecked_Conversion
(Complex_Vector
, X_Type
);
272 new Unchecked_Conversion
(Address
, Y_Ptr
);
274 BLAS
.cgemv
(Trans
, M
, N
, To_Fortran
(Alpha
),
275 Conv_A
(A
), Ld_A
, Conv_X
(X
), Inc_X
, To_Fortran
(Beta
),
276 Conv_Y
(Y
'Address).all, Inc_Y
);
282 BLAS
.Double_Complex_Matrix
(A
'Range (1), A
'Range (2));
284 BLAS
.Double_Complex_Vector
(X
'Range);
285 type Y_Ptr
is access all BLAS
.Double_Complex_Vector
(Y
'Range);
287 new Unchecked_Conversion
(Complex_Matrix
, A_Type
);
289 new Unchecked_Conversion
(Complex_Vector
, X_Type
);
291 new Unchecked_Conversion
(Address
, Y_Ptr
);
293 BLAS
.zgemv
(Trans
, M
, N
, To_Double_Complex
(Alpha
),
294 Conv_A
(A
), Ld_A
, Conv_X
(X
), Inc_X
,
295 To_Double_Complex
(Beta
),
296 Conv_Y
(Y
'Address).all, Inc_Y
);
301 DP_Y
: BLAS
.Double_Complex_Vector
(Y
'Range);
303 if Beta
.Re
/= 0.0 or else Beta
.Im
/= 0.0 then
304 DP_Y
:= To_Double_Complex
(Y
);
307 BLAS
.zgemv
(Trans
, M
, N
, To_Double_Complex
(Alpha
),
308 To_Double_Complex
(A
), Ld_A
,
309 To_Double_Complex
(X
), Inc_X
, To_Double_Complex
(Beta
),
312 Y
:= To_Complex
(DP_Y
);
324 Inc_X
: Integer := 1) return Real
329 subtype X_Type
is BLAS
.Complex_Vector
(X
'Range);
331 new Unchecked_Conversion
(Complex_Vector
, X_Type
);
333 return Real
(BLAS
.scnrm2
(N
, Conv_X
(X
), Inc_X
));
338 subtype X_Type
is BLAS
.Double_Complex_Vector
(X
'Range);
340 new Unchecked_Conversion
(Complex_Vector
, X_Type
);
342 return Real
(BLAS
.dznrm2
(N
, Conv_X
(X
), Inc_X
));
346 return Real
(BLAS
.dznrm2
(N
, To_Double_Complex
(X
), Inc_X
));
350 end System
.Generic_Complex_BLAS
;