1 ------------------------------------------------------------------------------
3 -- GNAT RUN-TIME COMPONENTS --
5 -- SYSTEM.GENERIC_COMPLEX_LAPACK --
9 -- Copyright (C) 2006, 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;