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-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 Interfaces
.Fortran
.LAPACK
; use Interfaces
.Fortran
.LAPACK
;
39 with System
.Generic_Array_Operations
; use System
.Generic_Array_Operations
;
41 package body System
.Generic_Complex_LAPACK
is
43 Is_Single
: constant Boolean :=
44 Real
'Machine_Mantissa = Fortran
.Real
'Machine_Mantissa
45 and then Fortran
.Real
(Real
'First) = Fortran
.Real
'First
46 and then Fortran
.Real
(Real
'Last) = Fortran
.Real
'Last;
48 Is_Double
: constant Boolean :=
49 Real
'Machine_Mantissa = Double_Precision
'Machine_Mantissa
51 Double_Precision
(Real
'First) = Double_Precision
'First
53 Double_Precision
(Real
'Last) = Double_Precision
'Last;
55 subtype Complex
is Complex_Types
.Complex
;
59 function To_Double_Precision
(X
: Real
) return Double_Precision
;
60 pragma Inline
(To_Double_Precision
);
62 function To_Real
(X
: Double_Precision
) return Real
;
63 pragma Inline
(To_Real
);
65 function To_Double_Complex
(X
: Complex
) return Double_Complex
;
66 pragma Inline
(To_Double_Complex
);
68 function To_Complex
(X
: Double_Complex
) return Complex
;
69 pragma Inline
(To_Complex
);
73 function To_Double_Precision
is new
74 Vector_Elementwise_Operation
76 Result_Scalar
=> Double_Precision
,
77 X_Vector
=> Real_Vector
,
78 Result_Vector
=> Double_Precision_Vector
,
79 Operation
=> To_Double_Precision
);
81 function To_Real
is new
82 Vector_Elementwise_Operation
83 (X_Scalar
=> Double_Precision
,
84 Result_Scalar
=> Real
,
85 X_Vector
=> Double_Precision_Vector
,
86 Result_Vector
=> Real_Vector
,
87 Operation
=> To_Real
);
89 function To_Double_Complex
is new
90 Matrix_Elementwise_Operation
92 Result_Scalar
=> Double_Complex
,
93 X_Matrix
=> Complex_Matrix
,
94 Result_Matrix
=> 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
=> 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_Real
(X
: Double_Precision
) return Real
is
115 function To_Double_Complex
(X
: Complex
) return Double_Complex
is
117 return (To_Double_Precision
(X
.Re
), To_Double_Precision
(X
.Im
));
118 end To_Double_Complex
;
120 function To_Complex
(X
: Double_Complex
) return Complex
is
122 return (Real
(X
.Re
), Real
(X
.Im
));
132 A
: in out Complex_Matrix
;
134 I_Piv
: out Integer_Vector
;
135 Info
: access Integer)
141 access all BLAS
.Complex_Matrix
(A
'Range (1), A
'Range (2));
142 function Conv_A
is new Unchecked_Conversion
(Address
, A_Ptr
);
144 cgetrf
(M
, N
, Conv_A
(A
'Address).all, Ld_A
,
145 LAPACK
.Integer_Vector
(I_Piv
), Info
);
151 access all Double_Complex_Matrix
(A
'Range (1), A
'Range (2));
152 function Conv_A
is new Unchecked_Conversion
(Address
, A_Ptr
);
154 zgetrf
(M
, N
, Conv_A
(A
'Address).all, Ld_A
,
155 LAPACK
.Integer_Vector
(I_Piv
), Info
);
160 DP_A
: Double_Complex_Matrix
(A
'Range (1), A
'Range (2));
162 DP_A
:= To_Double_Complex
(A
);
163 zgetrf
(M
, N
, DP_A
, Ld_A
, LAPACK
.Integer_Vector
(I_Piv
), Info
);
164 A
:= To_Complex
(DP_A
);
175 A
: in out Complex_Matrix
;
177 I_Piv
: Integer_Vector
;
178 Work
: in out Complex_Vector
;
180 Info
: access Integer)
186 access all BLAS
.Complex_Matrix
(A
'Range (1), A
'Range (2));
188 access all BLAS
.Complex_Vector
(Work
'Range);
189 function Conv_A
is new Unchecked_Conversion
(Address
, A_Ptr
);
190 function Conv_Work
is new Unchecked_Conversion
(Address
, Work_Ptr
);
192 cgetri
(N
, Conv_A
(A
'Address).all, Ld_A
,
193 LAPACK
.Integer_Vector
(I_Piv
),
194 Conv_Work
(Work
'Address).all, L_Work
,
201 access all Double_Complex_Matrix
(A
'Range (1), A
'Range (2));
203 access all Double_Complex_Vector
(Work
'Range);
204 function Conv_A
is new Unchecked_Conversion
(Address
, A_Ptr
);
205 function Conv_Work
is new Unchecked_Conversion
(Address
, Work_Ptr
);
207 zgetri
(N
, Conv_A
(A
'Address).all, Ld_A
,
208 LAPACK
.Integer_Vector
(I_Piv
),
209 Conv_Work
(Work
'Address).all, L_Work
,
215 DP_A
: Double_Complex_Matrix
(A
'Range (1), A
'Range (2));
216 DP_Work
: Double_Complex_Vector
(Work
'Range);
218 DP_A
:= To_Double_Complex
(A
);
219 zgetri
(N
, DP_A
, Ld_A
, LAPACK
.Integer_Vector
(I_Piv
),
220 DP_Work
, L_Work
, Info
);
221 A
:= To_Complex
(DP_A
);
222 Work
(1) := To_Complex
(DP_Work
(1));
232 (Trans
: access constant Character;
237 I_Piv
: Integer_Vector
;
238 B
: in out Complex_Matrix
;
240 Info
: access Integer)
245 subtype A_Type
is BLAS
.Complex_Matrix
(A
'Range (1), A
'Range (2));
247 access all BLAS
.Complex_Matrix
(B
'Range (1), B
'Range (2));
249 new Unchecked_Conversion
(Complex_Matrix
, A_Type
);
250 function Conv_B
is new Unchecked_Conversion
(Address
, B_Ptr
);
252 cgetrs
(Trans
, N
, N_Rhs
,
254 LAPACK
.Integer_Vector
(I_Piv
),
255 Conv_B
(B
'Address).all, Ld_B
,
262 Double_Complex_Matrix
(A
'Range (1), A
'Range (2));
264 access all Double_Complex_Matrix
(B
'Range (1), B
'Range (2));
266 new Unchecked_Conversion
(Complex_Matrix
, A_Type
);
267 function Conv_B
is new Unchecked_Conversion
(Address
, B_Ptr
);
269 zgetrs
(Trans
, N
, N_Rhs
,
271 LAPACK
.Integer_Vector
(I_Piv
),
272 Conv_B
(B
'Address).all, Ld_B
,
278 DP_A
: Double_Complex_Matrix
(A
'Range (1), A
'Range (2));
279 DP_B
: Double_Complex_Matrix
(B
'Range (1), B
'Range (2));
281 DP_A
:= To_Double_Complex
(A
);
282 DP_B
:= To_Double_Complex
(B
);
283 zgetrs
(Trans
, N
, N_Rhs
,
285 LAPACK
.Integer_Vector
(I_Piv
),
288 B
:= To_Complex
(DP_B
);
294 (Job_Z
: access constant Character;
295 Rng
: access constant Character;
296 Uplo
: access constant Character;
298 A
: in out Complex_Matrix
;
300 Vl
, Vu
: Real
:= 0.0;
301 Il
, Iu
: Integer := 1;
302 Abs_Tol
: Real
:= 0.0;
305 Z
: out Complex_Matrix
;
307 I_Supp_Z
: out Integer_Vector
;
308 Work
: out Complex_Vector
;
310 R_Work
: out Real_Vector
;
312 I_Work
: out Integer_Vector
;
314 Info
: access Integer)
320 access all BLAS
.Complex_Matrix
(A
'Range (1), A
'Range (2));
322 access all BLAS
.Real_Vector
(W
'Range);
324 access all BLAS
.Complex_Matrix
(Z
'Range (1), Z
'Range (2));
325 type Work_Ptr
is access all BLAS
.Complex_Vector
(Work
'Range);
326 type R_Work_Ptr
is access all BLAS
.Real_Vector
(R_Work
'Range);
328 function Conv_A
is new Unchecked_Conversion
(Address
, A_Ptr
);
329 function Conv_W
is new Unchecked_Conversion
(Address
, W_Ptr
);
330 function Conv_Z
is new Unchecked_Conversion
(Address
, Z_Ptr
);
331 function Conv_Work
is new Unchecked_Conversion
(Address
, Work_Ptr
);
332 function Conv_R_Work
is
333 new Unchecked_Conversion
(Address
, R_Work_Ptr
);
335 cheevr
(Job_Z
, Rng
, Uplo
, N
,
336 Conv_A
(A
'Address).all, Ld_A
,
337 Fortran
.Real
(Vl
), Fortran
.Real
(Vu
),
338 Il
, Iu
, Fortran
.Real
(Abs_Tol
), M
,
339 Conv_W
(W
'Address).all,
340 Conv_Z
(Z
'Address).all, Ld_Z
,
341 LAPACK
.Integer_Vector
(I_Supp_Z
),
342 Conv_Work
(Work
'Address).all, L_Work
,
343 Conv_R_Work
(R_Work
'Address).all, LR_Work
,
344 LAPACK
.Integer_Vector
(I_Work
), LI_Work
, Info
);
350 access all BLAS
.Double_Complex_Matrix
(A
'Range (1), A
'Range (2));
352 access all BLAS
.Double_Precision_Vector
(W
'Range);
354 access all BLAS
.Double_Complex_Matrix
(Z
'Range (1), Z
'Range (2));
356 access all BLAS
.Double_Complex_Vector
(Work
'Range);
358 access all BLAS
.Double_Precision_Vector
(R_Work
'Range);
360 function Conv_A
is new Unchecked_Conversion
(Address
, A_Ptr
);
361 function Conv_W
is new Unchecked_Conversion
(Address
, W_Ptr
);
362 function Conv_Z
is new Unchecked_Conversion
(Address
, Z_Ptr
);
363 function Conv_Work
is new Unchecked_Conversion
(Address
, Work_Ptr
);
364 function Conv_R_Work
is
365 new Unchecked_Conversion
(Address
, R_Work_Ptr
);
367 zheevr
(Job_Z
, Rng
, Uplo
, N
,
368 Conv_A
(A
'Address).all, Ld_A
,
369 Double_Precision
(Vl
), Double_Precision
(Vu
),
370 Il
, Iu
, Double_Precision
(Abs_Tol
), M
,
371 Conv_W
(W
'Address).all,
372 Conv_Z
(Z
'Address).all, Ld_Z
,
373 LAPACK
.Integer_Vector
(I_Supp_Z
),
374 Conv_Work
(Work
'Address).all, L_Work
,
375 Conv_R_Work
(R_Work
'Address).all, LR_Work
,
376 LAPACK
.Integer_Vector
(I_Work
), LI_Work
, Info
);
381 DP_A
: Double_Complex_Matrix
(A
'Range (1), A
'Range (2));
382 DP_W
: Double_Precision_Vector
(W
'Range);
383 DP_Z
: Double_Complex_Matrix
(Z
'Range (1), Z
'Range (2));
384 DP_Work
: Double_Complex_Vector
(Work
'Range);
385 DP_R_Work
: Double_Precision_Vector
(R_Work
'Range);
388 DP_A
:= To_Double_Complex
(A
);
390 zheevr
(Job_Z
, Rng
, Uplo
, N
,
392 Double_Precision
(Vl
), Double_Precision
(Vu
),
393 Il
, Iu
, Double_Precision
(Abs_Tol
), M
,
395 LAPACK
.Integer_Vector
(I_Supp_Z
),
398 LAPACK
.Integer_Vector
(I_Work
), LI_Work
, Info
);
400 A
:= To_Complex
(DP_A
);
402 Z
:= To_Complex
(DP_Z
);
404 Work
(1) := To_Complex
(DP_Work
(1));
405 R_Work
(1) := To_Real
(DP_R_Work
(1));
415 (Comp_Z
: access constant Character;
417 D
: in out Real_Vector
;
418 E
: in out Real_Vector
;
419 Z
: in out Complex_Matrix
;
421 Work
: out Real_Vector
;
422 Info
: access Integer)
427 type D_Ptr
is access all BLAS
.Real_Vector
(D
'Range);
428 type E_Ptr
is access all BLAS
.Real_Vector
(E
'Range);
430 access all BLAS
.Complex_Matrix
(Z
'Range (1), Z
'Range (2));
432 access all BLAS
.Real_Vector
(Work
'Range);
433 function Conv_D
is new Unchecked_Conversion
(Address
, D_Ptr
);
434 function Conv_E
is new Unchecked_Conversion
(Address
, E_Ptr
);
435 function Conv_Z
is new Unchecked_Conversion
(Address
, Z_Ptr
);
436 function Conv_Work
is new Unchecked_Conversion
(Address
, Work_Ptr
);
439 Conv_D
(D
'Address).all,
440 Conv_E
(E
'Address).all,
441 Conv_Z
(Z
'Address).all,
443 Conv_Work
(Work
'Address).all,
449 type D_Ptr
is access all Double_Precision_Vector
(D
'Range);
450 type E_Ptr
is access all Double_Precision_Vector
(E
'Range);
452 access all Double_Complex_Matrix
(Z
'Range (1), Z
'Range (2));
454 access all Double_Precision_Vector
(Work
'Range);
455 function Conv_D
is new Unchecked_Conversion
(Address
, D_Ptr
);
456 function Conv_E
is new Unchecked_Conversion
(Address
, E_Ptr
);
457 function Conv_Z
is new Unchecked_Conversion
(Address
, Z_Ptr
);
458 function Conv_Work
is new Unchecked_Conversion
(Address
, Work_Ptr
);
461 Conv_D
(D
'Address).all,
462 Conv_E
(E
'Address).all,
463 Conv_Z
(Z
'Address).all,
465 Conv_Work
(Work
'Address).all,
471 DP_D
: Double_Precision_Vector
(D
'Range);
472 DP_E
: Double_Precision_Vector
(E
'Range);
473 DP_Z
: Double_Complex_Matrix
(Z
'Range (1), Z
'Range (2));
474 DP_Work
: Double_Precision_Vector
(Work
'Range);
476 DP_D
:= To_Double_Precision
(D
);
477 DP_E
:= To_Double_Precision
(E
);
479 if Comp_Z
.all = 'V' then
480 DP_Z
:= To_Double_Complex
(Z
);
483 zsteqr
(Comp_Z
, N
, DP_D
, DP_E
, DP_Z
, Ld_Z
, DP_Work
, Info
);
488 if Comp_Z
.all /= 'N' then
489 Z
:= To_Complex
(DP_Z
);
495 end System
.Generic_Complex_LAPACK
;