1 ------------------------------------------------------------------------------
3 -- GNAT RUN-TIME COMPONENTS --
5 -- ADA.NUMERICS.GENERIC_REAL_ARRAYS --
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 System; use System;
33 with System.Generic_Real_BLAS;
34 with System.Generic_Real_LAPACK;
35 with System.Generic_Array_Operations; use System.Generic_Array_Operations;
37 package body Ada.Numerics.Generic_Real_Arrays is
39 -- Operations involving inner products use BLAS library implementations.
40 -- This allows larger matrices and vectors to be computed efficiently,
41 -- taking into account memory hierarchy issues and vector instructions
42 -- that vary widely between machines.
44 -- Operations that are defined in terms of operations on the type Real,
45 -- such as addition, subtraction and scaling, are computed in the canonical
46 -- way looping over all elements.
48 -- Operations for solving linear systems and computing determinant,
49 -- eigenvalues, eigensystem and inverse, are implemented using the
53 new Generic_Real_BLAS (Real'Base, Real_Vector, Real_Matrix);
56 new Generic_Real_LAPACK (Real'Base, Real_Vector, Real_Matrix);
60 -- Procedure versions of functions returning unconstrained values.
61 -- This allows for inlining the function wrapper.
63 procedure Eigenvalues (A : Real_Matrix; Values : out Real_Vector);
64 procedure Inverse (A : Real_Matrix; R : out Real_Matrix);
65 procedure Solve (A : Real_Matrix; X : Real_Vector; B : out Real_Vector);
66 procedure Solve (A : Real_Matrix; X : Real_Matrix; B : out Real_Matrix);
68 procedure Transpose is new
69 Generic_Array_Operations.Transpose
71 Matrix => Real_Matrix);
73 -- Helper function that raises a Constraint_Error is the argument is
74 -- not a square matrix, and otherwise returns its length.
76 function Length is new Square_Matrix_Length (Real'Base, Real_Matrix);
78 -- Instantiating the following subprograms directly would lead to
79 -- name clashes, so use a local package.
81 package Instantiations is
84 Vector_Elementwise_Operation
85 (X_Scalar => Real'Base,
86 Result_Scalar => Real'Base,
87 X_Vector => Real_Vector,
88 Result_Vector => Real_Vector,
92 Matrix_Elementwise_Operation
93 (X_Scalar => Real'Base,
94 Result_Scalar => Real'Base,
95 X_Matrix => Real_Matrix,
96 Result_Matrix => Real_Matrix,
100 Vector_Vector_Elementwise_Operation
101 (Left_Scalar => Real'Base,
102 Right_Scalar => Real'Base,
103 Result_Scalar => Real'Base,
104 Left_Vector => Real_Vector,
105 Right_Vector => Real_Vector,
106 Result_Vector => Real_Vector,
110 Matrix_Matrix_Elementwise_Operation
111 (Left_Scalar => Real'Base,
112 Right_Scalar => Real'Base,
113 Result_Scalar => Real'Base,
114 Left_Matrix => Real_Matrix,
115 Right_Matrix => Real_Matrix,
116 Result_Matrix => Real_Matrix,
120 Vector_Elementwise_Operation
121 (X_Scalar => Real'Base,
122 Result_Scalar => Real'Base,
123 X_Vector => Real_Vector,
124 Result_Vector => Real_Vector,
128 Matrix_Elementwise_Operation
129 (X_Scalar => Real'Base,
130 Result_Scalar => Real'Base,
131 X_Matrix => Real_Matrix,
132 Result_Matrix => Real_Matrix,
136 Vector_Vector_Elementwise_Operation
137 (Left_Scalar => Real'Base,
138 Right_Scalar => Real'Base,
139 Result_Scalar => Real'Base,
140 Left_Vector => Real_Vector,
141 Right_Vector => Real_Vector,
142 Result_Vector => Real_Vector,
146 Matrix_Matrix_Elementwise_Operation
147 (Left_Scalar => Real'Base,
148 Right_Scalar => Real'Base,
149 Result_Scalar => Real'Base,
150 Left_Matrix => Real_Matrix,
151 Right_Matrix => Real_Matrix,
152 Result_Matrix => Real_Matrix,
156 Scalar_Vector_Elementwise_Operation
157 (Left_Scalar => Real'Base,
158 Right_Scalar => Real'Base,
159 Result_Scalar => Real'Base,
160 Right_Vector => Real_Vector,
161 Result_Vector => Real_Vector,
165 Scalar_Matrix_Elementwise_Operation
166 (Left_Scalar => Real'Base,
167 Right_Scalar => Real'Base,
168 Result_Scalar => Real'Base,
169 Right_Matrix => Real_Matrix,
170 Result_Matrix => Real_Matrix,
174 Vector_Scalar_Elementwise_Operation
175 (Left_Scalar => Real'Base,
176 Right_Scalar => Real'Base,
177 Result_Scalar => Real'Base,
178 Left_Vector => Real_Vector,
179 Result_Vector => Real_Vector,
183 Matrix_Scalar_Elementwise_Operation
184 (Left_Scalar => Real'Base,
185 Right_Scalar => Real'Base,
186 Result_Scalar => Real'Base,
187 Left_Matrix => Real_Matrix,
188 Result_Matrix => Real_Matrix,
193 (Left_Scalar => Real'Base,
194 Right_Scalar => Real'Base,
195 Result_Scalar => Real'Base,
196 Left_Vector => Real_Vector,
197 Right_Vector => Real_Vector,
198 Matrix => Real_Matrix);
201 Vector_Scalar_Elementwise_Operation
202 (Left_Scalar => Real'Base,
203 Right_Scalar => Real'Base,
204 Result_Scalar => Real'Base,
205 Left_Vector => Real_Vector,
206 Result_Vector => Real_Vector,
210 Matrix_Scalar_Elementwise_Operation
211 (Left_Scalar => Real'Base,
212 Right_Scalar => Real'Base,
213 Result_Scalar => Real'Base,
214 Left_Matrix => Real_Matrix,
215 Result_Matrix => Real_Matrix,
218 function "abs" is new
219 Vector_Elementwise_Operation
220 (X_Scalar => Real'Base,
221 Result_Scalar => Real'Base,
222 X_Vector => Real_Vector,
223 Result_Vector => Real_Vector,
226 function "abs" is new
227 Matrix_Elementwise_Operation
228 (X_Scalar => Real'Base,
229 Result_Scalar => Real'Base,
230 X_Matrix => Real_Matrix,
231 Result_Matrix => Real_Matrix,
234 function Unit_Matrix is new
235 Generic_Array_Operations.Unit_Matrix
236 (Scalar => Real'Base,
237 Matrix => Real_Matrix,
241 function Unit_Vector is new
242 Generic_Array_Operations.Unit_Vector
243 (Scalar => Real'Base,
244 Vector => Real_Vector,
254 function "+" (Right : Real_Vector) return Real_Vector
255 renames Instantiations."+";
257 function "+" (Right : Real_Matrix) return Real_Matrix
258 renames Instantiations."+";
260 function "+" (Left, Right : Real_Vector) return Real_Vector
261 renames Instantiations."+";
263 function "+" (Left, Right : Real_Matrix) return Real_Matrix
264 renames Instantiations."+";
270 function "-" (Right : Real_Vector) return Real_Vector
271 renames Instantiations."-";
273 function "-" (Right : Real_Matrix) return Real_Matrix
274 renames Instantiations."-";
276 function "-" (Left, Right : Real_Vector) return Real_Vector
277 renames Instantiations."-";
279 function "-" (Left, Right : Real_Matrix) return Real_Matrix
280 renames Instantiations."-";
286 -- Scalar multiplication
288 function "*" (Left : Real'Base; Right : Real_Vector) return Real_Vector
289 renames Instantiations."*";
291 function "*" (Left : Real_Vector; Right : Real'Base) return Real_Vector
292 renames Instantiations."*";
294 function "*" (Left : Real'Base; Right : Real_Matrix) return Real_Matrix
295 renames Instantiations."*";
297 function "*" (Left : Real_Matrix; Right : Real'Base) return Real_Matrix
298 renames Instantiations."*";
300 -- Vector multiplication
302 function "*" (Left, Right : Real_Vector) return Real'Base is
304 if Left'Length /= Right'Length then
305 raise Constraint_Error with
306 "vectors are of different length in inner product";
309 return dot (Left'Length, X => Left, Y => Right);
312 function "*" (Left, Right : Real_Vector) return Real_Matrix
313 renames Instantiations."*";
317 Right : Real_Matrix) return Real_Vector
319 R : Real_Vector (Right'Range (2));
322 if Left'Length /= Right'Length (1) then
323 raise Constraint_Error with
324 "incompatible dimensions in vector-matrix multiplication";
327 gemv (Trans => No_Trans'Access,
328 M => Right'Length (2),
329 N => Right'Length (1),
331 Ld_A => Right'Length (2),
340 Right : Real_Vector) return Real_Vector
342 R : Real_Vector (Left'Range (1));
345 if Left'Length (2) /= Right'Length then
346 raise Constraint_Error with
347 "incompatible dimensions in matrix-vector multiplication";
350 gemv (Trans => Trans'Access,
351 M => Left'Length (2),
352 N => Left'Length (1),
354 Ld_A => Left'Length (2),
361 -- Matrix Multiplication
363 function "*" (Left, Right : Real_Matrix) return Real_Matrix is
364 R : Real_Matrix (Left'Range (1), Right'Range (2));
367 if Left'Length (2) /= Right'Length (1) then
368 raise Constraint_Error with
369 "incompatible dimensions in matrix-matrix multiplication";
372 gemm (Trans_A => No_Trans'Access,
373 Trans_B => No_Trans'Access,
374 M => Right'Length (2),
375 N => Left'Length (1),
376 K => Right'Length (1),
378 Ld_A => Right'Length (2),
380 Ld_B => Left'Length (2),
382 Ld_C => R'Length (2));
391 function "/" (Left : Real_Vector; Right : Real'Base) return Real_Vector
392 renames Instantiations."/";
394 function "/" (Left : Real_Matrix; Right : Real'Base) return Real_Matrix
395 renames Instantiations."/";
401 function "abs" (Right : Real_Vector) return Real'Base is
403 return nrm2 (Right'Length, Right);
406 function "abs" (Right : Real_Vector) return Real_Vector
407 renames Instantiations."abs";
409 function "abs" (Right : Real_Matrix) return Real_Matrix
410 renames Instantiations."abs";
416 function Determinant (A : Real_Matrix) return Real'Base is
417 N : constant Integer := Length (A);
418 LU : Real_Matrix (1 .. N, 1 .. N) := A;
419 Piv : Integer_Vector (1 .. N);
420 Info : aliased Integer := -1;
429 Info => Info'Access);
432 raise Constraint_Error with "ill-conditioned matrix";
436 Det := (if Piv (J) /= J then -Det * LU (J, J) else Det * LU (J, J));
446 procedure Eigensystem
448 Values : out Real_Vector;
449 Vectors : out Real_Matrix)
451 N : constant Natural := Length (A);
452 Tau : Real_Vector (1 .. N);
453 L_Work : Real_Vector (1 .. 1);
454 Info : aliased Integer;
456 E : Real_Vector (1 .. N);
457 pragma Warnings (Off, E);
460 if Values'Length /= N then
461 raise Constraint_Error with "wrong length for output vector";
468 -- Initialize working matrix and check for symmetric input matrix
470 Transpose (A, Vectors);
473 raise Argument_Error with "matrix not symmetric";
476 -- Compute size of additional working space
478 sytrd (Uplo => Lower'Access,
487 Info => Info'Access);
490 Work : Real_Vector (1 .. Integer'Max (Integer (L_Work (1)), 2 * N));
491 pragma Warnings (Off, Work);
493 Comp_Z : aliased constant Character := 'V';
496 -- Reduce matrix to tridiagonal form
498 sytrd (Uplo => Lower'Access,
501 Ld_A => A'Length (1),
506 L_Work => Work'Length,
507 Info => Info'Access);
513 -- Generate the real orthogonal matrix determined by sytrd
515 orgtr (Uplo => Lower'Access,
521 L_Work => Work'Length,
522 Info => Info'Access);
528 -- Compute all eigenvalues and eigenvectors using QR algorithm
530 steqr (Comp_Z => Comp_Z'Access,
537 Info => Info'Access);
540 raise Constraint_Error with
541 "eigensystem computation failed to converge";
550 procedure Eigenvalues
552 Values : out Real_Vector)
554 N : constant Natural := Length (A);
555 L_Work : Real_Vector (1 .. 1);
556 Info : aliased Integer;
558 B : Real_Matrix (1 .. N, 1 .. N);
559 Tau : Real_Vector (1 .. N);
560 E : Real_Vector (1 .. N);
561 pragma Warnings (Off, B);
562 pragma Warnings (Off, Tau);
563 pragma Warnings (Off, E);
566 if Values'Length /= N then
567 raise Constraint_Error with "wrong length for output vector";
574 -- Initialize working matrix and check for symmetric input matrix
579 raise Argument_Error with "matrix not symmetric";
582 -- Find size of work area
584 sytrd (Uplo => Lower'Access,
593 Info => Info'Access);
596 Work : Real_Vector (1 .. Integer'Min (Integer (L_Work (1)), 4 * N));
597 pragma Warnings (Off, Work);
600 -- Reduce matrix to tridiagonal form
602 sytrd (Uplo => Lower'Access,
605 Ld_A => A'Length (1),
610 L_Work => Work'Length,
611 Info => Info'Access);
614 raise Constraint_Error;
617 -- Compute all eigenvalues using QR algorithm
622 Info => Info'Access);
625 raise Constraint_Error with
626 "eigenvalues computation failed to converge";
631 function Eigenvalues (A : Real_Matrix) return Real_Vector is
632 R : Real_Vector (A'Range (1));
642 procedure Inverse (A : Real_Matrix; R : out Real_Matrix) is
643 N : constant Integer := Length (A);
644 Piv : Integer_Vector (1 .. N);
645 L_Work : Real_Vector (1 .. 1);
646 Info : aliased Integer := -1;
649 -- All computations are done using column-major order, but this works
650 -- out fine, because Transpose (Inverse (Transpose (A))) = Inverse (A).
654 -- Compute LU decomposition
661 Info => Info'Access);
664 raise Constraint_Error with "inverting singular matrix";
667 -- Determine size of work area
675 Info => Info'Access);
678 raise Constraint_Error;
682 Work : Real_Vector (1 .. Integer (L_Work (1)));
683 pragma Warnings (Off, Work);
686 -- Compute inverse from LU decomposition
693 L_Work => Work'Length,
694 Info => Info'Access);
697 raise Constraint_Error with "inverting singular matrix";
700 -- ??? Should iterate with gerfs, based on implementation advice
704 function Inverse (A : Real_Matrix) return Real_Matrix is
705 R : Real_Matrix (A'Range (2), A'Range (1));
715 procedure Solve (A : Real_Matrix; X : Real_Vector; B : out Real_Vector) is
717 if Length (A) /= X'Length then
718 raise Constraint_Error with
719 "incompatible matrix and vector dimensions";
722 -- ??? Should solve directly, is faster and more accurate
724 B := Inverse (A) * X;
727 procedure Solve (A : Real_Matrix; X : Real_Matrix; B : out Real_Matrix) is
729 if Length (A) /= X'Length (1) then
730 raise Constraint_Error with "incompatible matrix dimensions";
733 -- ??? Should solve directly, is faster and more accurate
735 B := Inverse (A) * X;
738 function Solve (A : Real_Matrix; X : Real_Vector) return Real_Vector is
739 B : Real_Vector (A'Range (2));
745 function Solve (A, X : Real_Matrix) return Real_Matrix is
746 B : Real_Matrix (A'Range (2), X'Range (2));
756 function Transpose (X : Real_Matrix) return Real_Matrix is
757 R : Real_Matrix (X'Range (2), X'Range (1));
770 First_1 : Integer := 1;
771 First_2 : Integer := 1) return Real_Matrix
772 renames Instantiations.Unit_Matrix;
781 First : Integer := 1) return Real_Vector
782 renames Instantiations.Unit_Vector;
784 end Ada.Numerics.Generic_Real_Arrays;