fixing pr42337
[official-gcc.git] / gcc / ada / s-gecola.adb
blobad69fee9bc5ff2fcb49c85f0107ad171bc0d241f
1 ------------------------------------------------------------------------------
2 -- --
3 -- GNAT RUN-TIME COMPONENTS --
4 -- --
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 --
6 -- --
7 -- B o d y --
8 -- --
9 -- Copyright (C) 2006-2009, Free Software Foundation, Inc. --
10 -- --
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. --
17 -- --
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. --
21 -- --
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/>. --
26 -- --
27 -- GNAT was originally developed by the GNAT team at New York University. --
28 -- Extensive contributions were provided by Ada Core Technologies Inc. --
29 -- --
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
48 and then
49 Double_Precision (Real'First) = Double_Precision'First
50 and then
51 Double_Precision (Real'Last) = Double_Precision'Last;
53 subtype Complex is Complex_Types.Complex;
55 -- Local subprograms
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);
69 -- Instantiations
71 function To_Double_Precision is new
72 Vector_Elementwise_Operation
73 (X_Scalar => Real,
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
89 (X_Scalar => Complex,
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
104 begin
105 return Double_Precision (X);
106 end To_Double_Precision;
108 function To_Real (X : Double_Precision) return Real is
109 begin
110 return Real (X);
111 end To_Real;
113 function To_Double_Complex (X : Complex) return Double_Complex is
114 begin
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
119 begin
120 return (Real (X.Re), Real (X.Im));
121 end To_Complex;
123 -----------
124 -- getrf --
125 -----------
127 procedure getrf
128 (M : Natural;
129 N : Natural;
130 A : in out Complex_Matrix;
131 Ld_A : Positive;
132 I_Piv : out Integer_Vector;
133 Info : access Integer)
135 begin
136 if Is_Single then
137 declare
138 type A_Ptr is
139 access all BLAS.Complex_Matrix (A'Range (1), A'Range (2));
140 function Conv_A is new Unchecked_Conversion (Address, A_Ptr);
141 begin
142 cgetrf (M, N, Conv_A (A'Address).all, Ld_A,
143 LAPACK.Integer_Vector (I_Piv), Info);
144 end;
146 elsif Is_Double then
147 declare
148 type A_Ptr is
149 access all Double_Complex_Matrix (A'Range (1), A'Range (2));
150 function Conv_A is new Unchecked_Conversion (Address, A_Ptr);
151 begin
152 zgetrf (M, N, Conv_A (A'Address).all, Ld_A,
153 LAPACK.Integer_Vector (I_Piv), Info);
154 end;
156 else
157 declare
158 DP_A : Double_Complex_Matrix (A'Range (1), A'Range (2));
159 begin
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);
163 end;
164 end if;
165 end getrf;
167 -----------
168 -- getri --
169 -----------
171 procedure getri
172 (N : Natural;
173 A : in out Complex_Matrix;
174 Ld_A : Positive;
175 I_Piv : Integer_Vector;
176 Work : in out Complex_Vector;
177 L_Work : Integer;
178 Info : access Integer)
180 begin
181 if Is_Single then
182 declare
183 type A_Ptr is
184 access all BLAS.Complex_Matrix (A'Range (1), A'Range (2));
185 type Work_Ptr is
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);
189 begin
190 cgetri (N, Conv_A (A'Address).all, Ld_A,
191 LAPACK.Integer_Vector (I_Piv),
192 Conv_Work (Work'Address).all, L_Work,
193 Info);
194 end;
196 elsif Is_Double then
197 declare
198 type A_Ptr is
199 access all Double_Complex_Matrix (A'Range (1), A'Range (2));
200 type Work_Ptr is
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);
204 begin
205 zgetri (N, Conv_A (A'Address).all, Ld_A,
206 LAPACK.Integer_Vector (I_Piv),
207 Conv_Work (Work'Address).all, L_Work,
208 Info);
209 end;
211 else
212 declare
213 DP_A : Double_Complex_Matrix (A'Range (1), A'Range (2));
214 DP_Work : Double_Complex_Vector (Work'Range);
215 begin
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));
221 end;
222 end if;
223 end getri;
225 -----------
226 -- getrs --
227 -----------
229 procedure getrs
230 (Trans : access constant Character;
231 N : Natural;
232 N_Rhs : Natural;
233 A : Complex_Matrix;
234 Ld_A : Positive;
235 I_Piv : Integer_Vector;
236 B : in out Complex_Matrix;
237 Ld_B : Positive;
238 Info : access Integer)
240 begin
241 if Is_Single then
242 declare
243 subtype A_Type is BLAS.Complex_Matrix (A'Range (1), A'Range (2));
244 type B_Ptr is
245 access all BLAS.Complex_Matrix (B'Range (1), B'Range (2));
246 function Conv_A is
247 new Unchecked_Conversion (Complex_Matrix, A_Type);
248 function Conv_B is new Unchecked_Conversion (Address, B_Ptr);
249 begin
250 cgetrs (Trans, N, N_Rhs,
251 Conv_A (A), Ld_A,
252 LAPACK.Integer_Vector (I_Piv),
253 Conv_B (B'Address).all, Ld_B,
254 Info);
255 end;
257 elsif Is_Double then
258 declare
259 subtype A_Type is
260 Double_Complex_Matrix (A'Range (1), A'Range (2));
261 type B_Ptr is
262 access all Double_Complex_Matrix (B'Range (1), B'Range (2));
263 function Conv_A is
264 new Unchecked_Conversion (Complex_Matrix, A_Type);
265 function Conv_B is new Unchecked_Conversion (Address, B_Ptr);
266 begin
267 zgetrs (Trans, N, N_Rhs,
268 Conv_A (A), Ld_A,
269 LAPACK.Integer_Vector (I_Piv),
270 Conv_B (B'Address).all, Ld_B,
271 Info);
272 end;
274 else
275 declare
276 DP_A : Double_Complex_Matrix (A'Range (1), A'Range (2));
277 DP_B : Double_Complex_Matrix (B'Range (1), B'Range (2));
278 begin
279 DP_A := To_Double_Complex (A);
280 DP_B := To_Double_Complex (B);
281 zgetrs (Trans, N, N_Rhs,
282 DP_A, Ld_A,
283 LAPACK.Integer_Vector (I_Piv),
284 DP_B, Ld_B,
285 Info);
286 B := To_Complex (DP_B);
287 end;
288 end if;
289 end getrs;
291 procedure heevr
292 (Job_Z : access constant Character;
293 Rng : access constant Character;
294 Uplo : access constant Character;
295 N : Natural;
296 A : in out Complex_Matrix;
297 Ld_A : Positive;
298 Vl, Vu : Real := 0.0;
299 Il, Iu : Integer := 1;
300 Abs_Tol : Real := 0.0;
301 M : out Integer;
302 W : out Real_Vector;
303 Z : out Complex_Matrix;
304 Ld_Z : Positive;
305 I_Supp_Z : out Integer_Vector;
306 Work : out Complex_Vector;
307 L_Work : Integer;
308 R_Work : out Real_Vector;
309 LR_Work : Integer;
310 I_Work : out Integer_Vector;
311 LI_Work : Integer;
312 Info : access Integer)
314 begin
315 if Is_Single then
316 declare
317 type A_Ptr is
318 access all BLAS.Complex_Matrix (A'Range (1), A'Range (2));
319 type W_Ptr is
320 access all BLAS.Real_Vector (W'Range);
321 type Z_Ptr is
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);
332 begin
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);
343 end;
345 elsif Is_Double then
346 declare
347 type A_Ptr is
348 access all BLAS.Double_Complex_Matrix (A'Range (1), A'Range (2));
349 type W_Ptr is
350 access all BLAS.Double_Precision_Vector (W'Range);
351 type Z_Ptr is
352 access all BLAS.Double_Complex_Matrix (Z'Range (1), Z'Range (2));
353 type Work_Ptr is
354 access all BLAS.Double_Complex_Vector (Work'Range);
355 type R_Work_Ptr is
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);
364 begin
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);
375 end;
377 else
378 declare
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);
385 begin
386 DP_A := To_Double_Complex (A);
388 zheevr (Job_Z, Rng, Uplo, N,
389 DP_A, Ld_A,
390 Double_Precision (Vl), Double_Precision (Vu),
391 Il, Iu, Double_Precision (Abs_Tol), M,
392 DP_W, DP_Z, Ld_Z,
393 LAPACK.Integer_Vector (I_Supp_Z),
394 DP_Work, L_Work,
395 DP_R_Work, LR_Work,
396 LAPACK.Integer_Vector (I_Work), LI_Work, Info);
398 A := To_Complex (DP_A);
399 W := To_Real (DP_W);
400 Z := To_Complex (DP_Z);
402 Work (1) := To_Complex (DP_Work (1));
403 R_Work (1) := To_Real (DP_R_Work (1));
404 end;
405 end if;
406 end heevr;
408 -----------
409 -- steqr --
410 -----------
412 procedure steqr
413 (Comp_Z : access constant Character;
414 N : Natural;
415 D : in out Real_Vector;
416 E : in out Real_Vector;
417 Z : in out Complex_Matrix;
418 Ld_Z : Positive;
419 Work : out Real_Vector;
420 Info : access Integer)
422 begin
423 if Is_Single then
424 declare
425 type D_Ptr is access all BLAS.Real_Vector (D'Range);
426 type E_Ptr is access all BLAS.Real_Vector (E'Range);
427 type Z_Ptr is
428 access all BLAS.Complex_Matrix (Z'Range (1), Z'Range (2));
429 type Work_Ptr is
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);
435 begin
436 csteqr (Comp_Z, N,
437 Conv_D (D'Address).all,
438 Conv_E (E'Address).all,
439 Conv_Z (Z'Address).all,
440 Ld_Z,
441 Conv_Work (Work'Address).all,
442 Info);
443 end;
445 elsif Is_Double then
446 declare
447 type D_Ptr is access all Double_Precision_Vector (D'Range);
448 type E_Ptr is access all Double_Precision_Vector (E'Range);
449 type Z_Ptr is
450 access all Double_Complex_Matrix (Z'Range (1), Z'Range (2));
451 type Work_Ptr is
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);
457 begin
458 zsteqr (Comp_Z, N,
459 Conv_D (D'Address).all,
460 Conv_E (E'Address).all,
461 Conv_Z (Z'Address).all,
462 Ld_Z,
463 Conv_Work (Work'Address).all,
464 Info);
465 end;
467 else
468 declare
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);
473 begin
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);
479 end if;
481 zsteqr (Comp_Z, N, DP_D, DP_E, DP_Z, Ld_Z, DP_Work, Info);
483 D := To_Real (DP_D);
484 E := To_Real (DP_E);
486 if Comp_Z.all /= 'N' then
487 Z := To_Complex (DP_Z);
488 end if;
489 end;
490 end if;
491 end steqr;
493 end System.Generic_Complex_LAPACK;