1 ------------------------------------------------------------------------------
3 -- GNAT RUN-TIME COMPONENTS --
5 -- SYSTEM.GENERIC_REAL_LAPACK --
9 -- Copyright (C) 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_Real_LAPACK is
41 Is_Real : 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_Precision : constant Boolean :=
47 Real'Machine_Mantissa =
48 Double_Precision'Machine_Mantissa
50 Double_Precision (Real'First) =
51 Double_Precision'First
53 Double_Precision (Real'Last) =
54 Double_Precision'Last;
58 function To_Double_Precision (X : Real) return Double_Precision;
59 pragma Inline_Always (To_Double_Precision);
61 function To_Real (X : Double_Precision) return Real;
62 pragma Inline_Always (To_Real);
66 function To_Double_Precision is new
67 Vector_Elementwise_Operation
69 Result_Scalar => Double_Precision,
70 X_Vector => Real_Vector,
71 Result_Vector => Double_Precision_Vector,
72 Operation => To_Double_Precision);
74 function To_Real is new
75 Vector_Elementwise_Operation
76 (X_Scalar => Double_Precision,
77 Result_Scalar => Real,
78 X_Vector => Double_Precision_Vector,
79 Result_Vector => Real_Vector,
80 Operation => To_Real);
82 function To_Double_Precision is new
83 Matrix_Elementwise_Operation
85 Result_Scalar => Double_Precision,
86 X_Matrix => Real_Matrix,
87 Result_Matrix => Double_Precision_Matrix,
88 Operation => To_Double_Precision);
90 function To_Real is new
91 Matrix_Elementwise_Operation
92 (X_Scalar => Double_Precision,
93 Result_Scalar => Real,
94 X_Matrix => Double_Precision_Matrix,
95 Result_Matrix => Real_Matrix,
96 Operation => To_Real);
98 function To_Double_Precision (X : Real) return Double_Precision is
100 return Double_Precision (X);
101 end To_Double_Precision;
103 function To_Real (X : Double_Precision) return Real is
115 A : in out Real_Matrix;
117 I_Piv : out Integer_Vector;
118 Info : access Integer)
124 access all BLAS.Real_Matrix (A'Range (1), A'Range (2));
125 function Conv_A is new Unchecked_Conversion (Address, A_Ptr);
127 sgetrf (M, N, Conv_A (A'Address).all, Ld_A,
128 LAPACK.Integer_Vector (I_Piv), Info);
131 elsif Is_Double_Precision then
134 access all Double_Precision_Matrix (A'Range (1), A'Range (2));
135 function Conv_A is new Unchecked_Conversion (Address, A_Ptr);
137 dgetrf (M, N, Conv_A (A'Address).all, Ld_A,
138 LAPACK.Integer_Vector (I_Piv), Info);
143 DP_A : Double_Precision_Matrix (A'Range (1), A'Range (2));
145 DP_A := To_Double_Precision (A);
146 dgetrf (M, N, DP_A, Ld_A, LAPACK.Integer_Vector (I_Piv), Info);
158 A : in out Real_Matrix;
160 I_Piv : Integer_Vector;
161 Work : in out Real_Vector;
163 Info : access Integer)
169 access all BLAS.Real_Matrix (A'Range (1), A'Range (2));
171 access all BLAS.Real_Vector (Work'Range);
172 function Conv_A is new Unchecked_Conversion (Address, A_Ptr);
173 function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr);
175 sgetri (N, Conv_A (A'Address).all, Ld_A,
176 LAPACK.Integer_Vector (I_Piv),
177 Conv_Work (Work'Address).all, L_Work,
181 elsif Is_Double_Precision then
184 access all Double_Precision_Matrix (A'Range (1), A'Range (2));
186 access all Double_Precision_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 dgetri (N, Conv_A (A'Address).all, Ld_A,
191 LAPACK.Integer_Vector (I_Piv),
192 Conv_Work (Work'Address).all, L_Work,
198 DP_A : Double_Precision_Matrix (A'Range (1), A'Range (2));
199 DP_Work : Double_Precision_Vector (Work'Range);
201 DP_A := To_Double_Precision (A);
202 dgetri (N, DP_A, Ld_A, LAPACK.Integer_Vector (I_Piv),
203 DP_Work, L_Work, Info);
205 Work (1) := To_Real (DP_Work (1));
215 (Trans : access constant Character;
220 I_Piv : Integer_Vector;
221 B : in out Real_Matrix;
223 Info : access Integer)
228 subtype A_Type is BLAS.Real_Matrix (A'Range (1), A'Range (2));
230 access all BLAS.Real_Matrix (B'Range (1), B'Range (2));
231 function Conv_A is new Unchecked_Conversion (Real_Matrix, A_Type);
232 function Conv_B is new Unchecked_Conversion (Address, B_Ptr);
234 sgetrs (Trans, N, N_Rhs,
236 LAPACK.Integer_Vector (I_Piv),
237 Conv_B (B'Address).all, Ld_B,
241 elsif Is_Double_Precision then
244 Double_Precision_Matrix (A'Range (1), A'Range (2));
246 access all Double_Precision_Matrix (B'Range (1), B'Range (2));
247 function Conv_A is new Unchecked_Conversion (Real_Matrix, A_Type);
248 function Conv_B is new Unchecked_Conversion (Address, B_Ptr);
250 dgetrs (Trans, N, N_Rhs,
252 LAPACK.Integer_Vector (I_Piv),
253 Conv_B (B'Address).all, Ld_B,
259 DP_A : Double_Precision_Matrix (A'Range (1), A'Range (2));
260 DP_B : Double_Precision_Matrix (B'Range (1), B'Range (2));
262 DP_A := To_Double_Precision (A);
263 DP_B := To_Double_Precision (B);
264 dgetrs (Trans, N, N_Rhs,
266 LAPACK.Integer_Vector (I_Piv),
279 (Uplo : access constant Character;
281 A : in out Real_Matrix;
284 Work : out Real_Vector;
286 Info : access Integer)
292 access all BLAS.Real_Matrix (A'Range (1), A'Range (2));
293 subtype Tau_Type is BLAS.Real_Vector (Tau'Range);
295 access all BLAS.Real_Vector (Work'Range);
296 function Conv_A is new Unchecked_Conversion (Address, A_Ptr);
298 new Unchecked_Conversion (Real_Vector, Tau_Type);
299 function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr);
302 Conv_A (A'Address).all, Ld_A,
304 Conv_Work (Work'Address).all, L_Work,
308 elsif Is_Double_Precision then
311 access all Double_Precision_Matrix (A'Range (1), A'Range (2));
312 subtype Tau_Type is Double_Precision_Vector (Tau'Range);
314 access all Double_Precision_Vector (Work'Range);
315 function Conv_A is new Unchecked_Conversion (Address, A_Ptr);
317 new Unchecked_Conversion (Real_Vector, Tau_Type);
318 function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr);
321 Conv_A (A'Address).all, Ld_A,
323 Conv_Work (Work'Address).all, L_Work,
329 DP_A : Double_Precision_Matrix (A'Range (1), A'Range (2));
330 DP_Work : Double_Precision_Vector (Work'Range);
331 DP_Tau : Double_Precision_Vector (Tau'Range);
333 DP_A := To_Double_Precision (A);
334 DP_Tau := To_Double_Precision (Tau);
335 dorgtr (Uplo, N, DP_A, Ld_A, DP_Tau, DP_Work, L_Work, Info);
337 Work (1) := To_Real (DP_Work (1));
347 (Comp_Z : access constant Character;
349 D : in out Real_Vector;
350 E : in out Real_Vector;
351 Z : in out Real_Matrix;
353 Work : out Real_Vector;
354 Info : access Integer)
359 type D_Ptr is access all BLAS.Real_Vector (D'Range);
360 type E_Ptr is access all BLAS.Real_Vector (E'Range);
362 access all BLAS.Real_Matrix (Z'Range (1), Z'Range (2));
364 access all BLAS.Real_Vector (Work'Range);
365 function Conv_D is new Unchecked_Conversion (Address, D_Ptr);
366 function Conv_E is new Unchecked_Conversion (Address, E_Ptr);
367 function Conv_Z is new Unchecked_Conversion (Address, Z_Ptr);
368 function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr);
371 Conv_D (D'Address).all,
372 Conv_E (E'Address).all,
373 Conv_Z (Z'Address).all,
375 Conv_Work (Work'Address).all,
379 elsif Is_Double_Precision then
381 type D_Ptr is access all Double_Precision_Vector (D'Range);
382 type E_Ptr is access all Double_Precision_Vector (E'Range);
384 access all Double_Precision_Matrix (Z'Range (1), Z'Range (2));
386 access all Double_Precision_Vector (Work'Range);
387 function Conv_D is new Unchecked_Conversion (Address, D_Ptr);
388 function Conv_E is new Unchecked_Conversion (Address, E_Ptr);
389 function Conv_Z is new Unchecked_Conversion (Address, Z_Ptr);
390 function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr);
393 Conv_D (D'Address).all,
394 Conv_E (E'Address).all,
395 Conv_Z (Z'Address).all,
397 Conv_Work (Work'Address).all,
403 DP_D : Double_Precision_Vector (D'Range);
404 DP_E : Double_Precision_Vector (E'Range);
405 DP_Z : Double_Precision_Matrix (Z'Range (1), Z'Range (2));
406 DP_Work : Double_Precision_Vector (Work'Range);
408 DP_D := To_Double_Precision (D);
409 DP_E := To_Double_Precision (E);
411 if Comp_Z.all = 'V' then
412 DP_Z := To_Double_Precision (Z);
415 dsteqr (Comp_Z, N, DP_D, DP_E, DP_Z, Ld_Z, DP_Work, Info);
430 D : in out Real_Vector;
431 E : in out Real_Vector;
432 Info : access Integer)
437 type D_Ptr is access all BLAS.Real_Vector (D'Range);
438 type E_Ptr is access all BLAS.Real_Vector (E'Range);
439 function Conv_D is new Unchecked_Conversion (Address, D_Ptr);
440 function Conv_E is new Unchecked_Conversion (Address, E_Ptr);
442 ssterf (N, Conv_D (D'Address).all, Conv_E (E'Address).all, Info);
445 elsif Is_Double_Precision then
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 function Conv_D is new Unchecked_Conversion (Address, D_Ptr);
450 function Conv_E is new Unchecked_Conversion (Address, E_Ptr);
452 dsterf (N, Conv_D (D'Address).all, Conv_E (E'Address).all, Info);
457 DP_D : Double_Precision_Vector (D'Range);
458 DP_E : Double_Precision_Vector (E'Range);
461 DP_D := To_Double_Precision (D);
462 DP_E := To_Double_Precision (E);
464 dsterf (N, DP_D, DP_E, Info);
477 (Uplo : access constant Character;
479 A : in out Real_Matrix;
483 Tau : out Real_Vector;
484 Work : out Real_Vector;
486 Info : access Integer)
492 access all BLAS.Real_Matrix (A'Range (1), A'Range (2));
493 type D_Ptr is access all BLAS.Real_Vector (D'Range);
494 type E_Ptr is access all BLAS.Real_Vector (E'Range);
495 type Tau_Ptr is access all BLAS.Real_Vector (Tau'Range);
497 access all BLAS.Real_Vector (Work'Range);
498 function Conv_A is new Unchecked_Conversion (Address, A_Ptr);
499 function Conv_D is new Unchecked_Conversion (Address, D_Ptr);
500 function Conv_E is new Unchecked_Conversion (Address, E_Ptr);
501 function Conv_Tau is new Unchecked_Conversion (Address, Tau_Ptr);
502 function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr);
505 Conv_A (A'Address).all, Ld_A,
506 Conv_D (D'Address).all,
507 Conv_E (E'Address).all,
508 Conv_Tau (Tau'Address).all,
509 Conv_Work (Work'Address).all,
514 elsif Is_Double_Precision then
517 access all Double_Precision_Matrix (A'Range (1), A'Range (2));
518 type D_Ptr is access all Double_Precision_Vector (D'Range);
519 type E_Ptr is access all Double_Precision_Vector (E'Range);
520 type Tau_Ptr is access all Double_Precision_Vector (Tau'Range);
522 access all Double_Precision_Vector (Work'Range);
523 function Conv_A is new Unchecked_Conversion (Address, A_Ptr);
524 function Conv_D is new Unchecked_Conversion (Address, D_Ptr);
525 function Conv_E is new Unchecked_Conversion (Address, E_Ptr);
526 function Conv_Tau is new Unchecked_Conversion (Address, Tau_Ptr);
527 function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr);
530 Conv_A (A'Address).all, Ld_A,
531 Conv_D (D'Address).all,
532 Conv_E (E'Address).all,
533 Conv_Tau (Tau'Address).all,
534 Conv_Work (Work'Address).all,
541 DP_A : Double_Precision_Matrix (A'Range (1), A'Range (2));
542 DP_D : Double_Precision_Vector (D'Range);
543 DP_E : Double_Precision_Vector (E'Range);
544 DP_Tau : Double_Precision_Vector (Tau'Range);
545 DP_Work : Double_Precision_Vector (Work'Range);
547 DP_A := To_Double_Precision (A);
549 dsytrd (Uplo, N, DP_A, Ld_A, DP_D, DP_E, DP_Tau,
550 DP_Work, L_Work, Info);
556 Tau := To_Real (DP_Tau);
559 Work (1) := To_Real (DP_Work (1));
564 end System.Generic_Real_LAPACK;