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 _ L A P A C K --
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 Interfaces
.Fortran
.LAPACK
; use Interfaces
.Fortran
.LAPACK
;
37 with System
.Generic_Array_Operations
; use System
.Generic_Array_Operations
;
39 package body System
.Generic_Complex_LAPACK
is
41 Is_Single
: constant Boolean :=
42 Real
'Machine_Mantissa = Fortran
.Real
'Machine_Mantissa
43 and then Fortran
.Real
(Real
'First) = Fortran
.Real
'First
44 and then Fortran
.Real
(Real
'Last) = Fortran
.Real
'Last;
46 Is_Double
: constant Boolean :=
47 Real
'Machine_Mantissa = Double_Precision
'Machine_Mantissa
49 Double_Precision
(Real
'First) = Double_Precision
'First
51 Double_Precision
(Real
'Last) = Double_Precision
'Last;
53 subtype Complex
is Complex_Types
.Complex
;
57 function To_Double_Precision
(X
: Real
) return Double_Precision
;
58 pragma Inline
(To_Double_Precision
);
60 function To_Real
(X
: Double_Precision
) return Real
;
61 pragma Inline
(To_Real
);
63 function To_Double_Complex
(X
: Complex
) return Double_Complex
;
64 pragma Inline
(To_Double_Complex
);
66 function To_Complex
(X
: Double_Complex
) return Complex
;
67 pragma Inline
(To_Complex
);
71 function To_Double_Precision
is new
72 Vector_Elementwise_Operation
74 Result_Scalar
=> Double_Precision
,
75 X_Vector
=> Real_Vector
,
76 Result_Vector
=> Double_Precision_Vector
,
77 Operation
=> To_Double_Precision
);
79 function To_Real
is new
80 Vector_Elementwise_Operation
81 (X_Scalar
=> Double_Precision
,
82 Result_Scalar
=> Real
,
83 X_Vector
=> Double_Precision_Vector
,
84 Result_Vector
=> Real_Vector
,
85 Operation
=> To_Real
);
87 function To_Double_Complex
is new
88 Matrix_Elementwise_Operation
90 Result_Scalar
=> Double_Complex
,
91 X_Matrix
=> Complex_Matrix
,
92 Result_Matrix
=> 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
=> 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_Real
(X
: Double_Precision
) return Real
is
113 function To_Double_Complex
(X
: Complex
) return Double_Complex
is
115 return (To_Double_Precision
(X
.Re
), To_Double_Precision
(X
.Im
));
116 end To_Double_Complex
;
118 function To_Complex
(X
: Double_Complex
) return Complex
is
120 return (Real
(X
.Re
), Real
(X
.Im
));
130 A
: in out Complex_Matrix
;
132 I_Piv
: out Integer_Vector
;
133 Info
: access Integer)
139 access all BLAS
.Complex_Matrix
(A
'Range (1), A
'Range (2));
140 function Conv_A
is new Unchecked_Conversion
(Address
, A_Ptr
);
142 cgetrf
(M
, N
, Conv_A
(A
'Address).all, Ld_A
,
143 LAPACK
.Integer_Vector
(I_Piv
), Info
);
149 access all Double_Complex_Matrix
(A
'Range (1), A
'Range (2));
150 function Conv_A
is new Unchecked_Conversion
(Address
, A_Ptr
);
152 zgetrf
(M
, N
, Conv_A
(A
'Address).all, Ld_A
,
153 LAPACK
.Integer_Vector
(I_Piv
), Info
);
158 DP_A
: Double_Complex_Matrix
(A
'Range (1), A
'Range (2));
160 DP_A
:= To_Double_Complex
(A
);
161 zgetrf
(M
, N
, DP_A
, Ld_A
, LAPACK
.Integer_Vector
(I_Piv
), Info
);
162 A
:= To_Complex
(DP_A
);
173 A
: in out Complex_Matrix
;
175 I_Piv
: Integer_Vector
;
176 Work
: in out Complex_Vector
;
178 Info
: access Integer)
184 access all BLAS
.Complex_Matrix
(A
'Range (1), A
'Range (2));
186 access all BLAS
.Complex_Vector
(Work
'Range);
187 function Conv_A
is new Unchecked_Conversion
(Address
, A_Ptr
);
188 function Conv_Work
is new Unchecked_Conversion
(Address
, Work_Ptr
);
190 cgetri
(N
, Conv_A
(A
'Address).all, Ld_A
,
191 LAPACK
.Integer_Vector
(I_Piv
),
192 Conv_Work
(Work
'Address).all, L_Work
,
199 access all Double_Complex_Matrix
(A
'Range (1), A
'Range (2));
201 access all Double_Complex_Vector
(Work
'Range);
202 function Conv_A
is new Unchecked_Conversion
(Address
, A_Ptr
);
203 function Conv_Work
is new Unchecked_Conversion
(Address
, Work_Ptr
);
205 zgetri
(N
, Conv_A
(A
'Address).all, Ld_A
,
206 LAPACK
.Integer_Vector
(I_Piv
),
207 Conv_Work
(Work
'Address).all, L_Work
,
213 DP_A
: Double_Complex_Matrix
(A
'Range (1), A
'Range (2));
214 DP_Work
: Double_Complex_Vector
(Work
'Range);
216 DP_A
:= To_Double_Complex
(A
);
217 zgetri
(N
, DP_A
, Ld_A
, LAPACK
.Integer_Vector
(I_Piv
),
218 DP_Work
, L_Work
, Info
);
219 A
:= To_Complex
(DP_A
);
220 Work
(1) := To_Complex
(DP_Work
(1));
230 (Trans
: access constant Character;
235 I_Piv
: Integer_Vector
;
236 B
: in out Complex_Matrix
;
238 Info
: access Integer)
243 subtype A_Type
is BLAS
.Complex_Matrix
(A
'Range (1), A
'Range (2));
245 access all BLAS
.Complex_Matrix
(B
'Range (1), B
'Range (2));
247 new Unchecked_Conversion
(Complex_Matrix
, A_Type
);
248 function Conv_B
is new Unchecked_Conversion
(Address
, B_Ptr
);
250 cgetrs
(Trans
, N
, N_Rhs
,
252 LAPACK
.Integer_Vector
(I_Piv
),
253 Conv_B
(B
'Address).all, Ld_B
,
260 Double_Complex_Matrix
(A
'Range (1), A
'Range (2));
262 access all Double_Complex_Matrix
(B
'Range (1), B
'Range (2));
264 new Unchecked_Conversion
(Complex_Matrix
, A_Type
);
265 function Conv_B
is new Unchecked_Conversion
(Address
, B_Ptr
);
267 zgetrs
(Trans
, N
, N_Rhs
,
269 LAPACK
.Integer_Vector
(I_Piv
),
270 Conv_B
(B
'Address).all, Ld_B
,
276 DP_A
: Double_Complex_Matrix
(A
'Range (1), A
'Range (2));
277 DP_B
: Double_Complex_Matrix
(B
'Range (1), B
'Range (2));
279 DP_A
:= To_Double_Complex
(A
);
280 DP_B
:= To_Double_Complex
(B
);
281 zgetrs
(Trans
, N
, N_Rhs
,
283 LAPACK
.Integer_Vector
(I_Piv
),
286 B
:= To_Complex
(DP_B
);
292 (Job_Z
: access constant Character;
293 Rng
: access constant Character;
294 Uplo
: access constant Character;
296 A
: in out Complex_Matrix
;
298 Vl
, Vu
: Real
:= 0.0;
299 Il
, Iu
: Integer := 1;
300 Abs_Tol
: Real
:= 0.0;
303 Z
: out Complex_Matrix
;
305 I_Supp_Z
: out Integer_Vector
;
306 Work
: out Complex_Vector
;
308 R_Work
: out Real_Vector
;
310 I_Work
: out Integer_Vector
;
312 Info
: access Integer)
318 access all BLAS
.Complex_Matrix
(A
'Range (1), A
'Range (2));
320 access all BLAS
.Real_Vector
(W
'Range);
322 access all BLAS
.Complex_Matrix
(Z
'Range (1), Z
'Range (2));
323 type Work_Ptr
is access all BLAS
.Complex_Vector
(Work
'Range);
324 type R_Work_Ptr
is access all BLAS
.Real_Vector
(R_Work
'Range);
326 function Conv_A
is new Unchecked_Conversion
(Address
, A_Ptr
);
327 function Conv_W
is new Unchecked_Conversion
(Address
, W_Ptr
);
328 function Conv_Z
is new Unchecked_Conversion
(Address
, Z_Ptr
);
329 function Conv_Work
is new Unchecked_Conversion
(Address
, Work_Ptr
);
330 function Conv_R_Work
is
331 new Unchecked_Conversion
(Address
, R_Work_Ptr
);
333 cheevr
(Job_Z
, Rng
, Uplo
, N
,
334 Conv_A
(A
'Address).all, Ld_A
,
335 Fortran
.Real
(Vl
), Fortran
.Real
(Vu
),
336 Il
, Iu
, Fortran
.Real
(Abs_Tol
), M
,
337 Conv_W
(W
'Address).all,
338 Conv_Z
(Z
'Address).all, Ld_Z
,
339 LAPACK
.Integer_Vector
(I_Supp_Z
),
340 Conv_Work
(Work
'Address).all, L_Work
,
341 Conv_R_Work
(R_Work
'Address).all, LR_Work
,
342 LAPACK
.Integer_Vector
(I_Work
), LI_Work
, Info
);
348 access all BLAS
.Double_Complex_Matrix
(A
'Range (1), A
'Range (2));
350 access all BLAS
.Double_Precision_Vector
(W
'Range);
352 access all BLAS
.Double_Complex_Matrix
(Z
'Range (1), Z
'Range (2));
354 access all BLAS
.Double_Complex_Vector
(Work
'Range);
356 access all BLAS
.Double_Precision_Vector
(R_Work
'Range);
358 function Conv_A
is new Unchecked_Conversion
(Address
, A_Ptr
);
359 function Conv_W
is new Unchecked_Conversion
(Address
, W_Ptr
);
360 function Conv_Z
is new Unchecked_Conversion
(Address
, Z_Ptr
);
361 function Conv_Work
is new Unchecked_Conversion
(Address
, Work_Ptr
);
362 function Conv_R_Work
is
363 new Unchecked_Conversion
(Address
, R_Work_Ptr
);
365 zheevr
(Job_Z
, Rng
, Uplo
, N
,
366 Conv_A
(A
'Address).all, Ld_A
,
367 Double_Precision
(Vl
), Double_Precision
(Vu
),
368 Il
, Iu
, Double_Precision
(Abs_Tol
), M
,
369 Conv_W
(W
'Address).all,
370 Conv_Z
(Z
'Address).all, Ld_Z
,
371 LAPACK
.Integer_Vector
(I_Supp_Z
),
372 Conv_Work
(Work
'Address).all, L_Work
,
373 Conv_R_Work
(R_Work
'Address).all, LR_Work
,
374 LAPACK
.Integer_Vector
(I_Work
), LI_Work
, Info
);
379 DP_A
: Double_Complex_Matrix
(A
'Range (1), A
'Range (2));
380 DP_W
: Double_Precision_Vector
(W
'Range);
381 DP_Z
: Double_Complex_Matrix
(Z
'Range (1), Z
'Range (2));
382 DP_Work
: Double_Complex_Vector
(Work
'Range);
383 DP_R_Work
: Double_Precision_Vector
(R_Work
'Range);
386 DP_A
:= To_Double_Complex
(A
);
388 zheevr
(Job_Z
, Rng
, Uplo
, N
,
390 Double_Precision
(Vl
), Double_Precision
(Vu
),
391 Il
, Iu
, Double_Precision
(Abs_Tol
), M
,
393 LAPACK
.Integer_Vector
(I_Supp_Z
),
396 LAPACK
.Integer_Vector
(I_Work
), LI_Work
, Info
);
398 A
:= To_Complex
(DP_A
);
400 Z
:= To_Complex
(DP_Z
);
402 Work
(1) := To_Complex
(DP_Work
(1));
403 R_Work
(1) := To_Real
(DP_R_Work
(1));
413 (Comp_Z
: access constant Character;
415 D
: in out Real_Vector
;
416 E
: in out Real_Vector
;
417 Z
: in out Complex_Matrix
;
419 Work
: out Real_Vector
;
420 Info
: access Integer)
425 type D_Ptr
is access all BLAS
.Real_Vector
(D
'Range);
426 type E_Ptr
is access all BLAS
.Real_Vector
(E
'Range);
428 access all BLAS
.Complex_Matrix
(Z
'Range (1), Z
'Range (2));
430 access all BLAS
.Real_Vector
(Work
'Range);
431 function Conv_D
is new Unchecked_Conversion
(Address
, D_Ptr
);
432 function Conv_E
is new Unchecked_Conversion
(Address
, E_Ptr
);
433 function Conv_Z
is new Unchecked_Conversion
(Address
, Z_Ptr
);
434 function Conv_Work
is new Unchecked_Conversion
(Address
, Work_Ptr
);
437 Conv_D
(D
'Address).all,
438 Conv_E
(E
'Address).all,
439 Conv_Z
(Z
'Address).all,
441 Conv_Work
(Work
'Address).all,
447 type D_Ptr
is access all Double_Precision_Vector
(D
'Range);
448 type E_Ptr
is access all Double_Precision_Vector
(E
'Range);
450 access all Double_Complex_Matrix
(Z
'Range (1), Z
'Range (2));
452 access all Double_Precision_Vector
(Work
'Range);
453 function Conv_D
is new Unchecked_Conversion
(Address
, D_Ptr
);
454 function Conv_E
is new Unchecked_Conversion
(Address
, E_Ptr
);
455 function Conv_Z
is new Unchecked_Conversion
(Address
, Z_Ptr
);
456 function Conv_Work
is new Unchecked_Conversion
(Address
, Work_Ptr
);
459 Conv_D
(D
'Address).all,
460 Conv_E
(E
'Address).all,
461 Conv_Z
(Z
'Address).all,
463 Conv_Work
(Work
'Address).all,
469 DP_D
: Double_Precision_Vector
(D
'Range);
470 DP_E
: Double_Precision_Vector
(E
'Range);
471 DP_Z
: Double_Complex_Matrix
(Z
'Range (1), Z
'Range (2));
472 DP_Work
: Double_Precision_Vector
(Work
'Range);
474 DP_D
:= To_Double_Precision
(D
);
475 DP_E
:= To_Double_Precision
(E
);
477 if Comp_Z
.all = 'V' then
478 DP_Z
:= To_Double_Complex
(Z
);
481 zsteqr
(Comp_Z
, N
, DP_D
, DP_E
, DP_Z
, Ld_Z
, DP_Work
, Info
);
486 if Comp_Z
.all /= 'N' then
487 Z
:= To_Complex
(DP_Z
);
493 end System
.Generic_Complex_LAPACK
;