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;