1 ------------------------------------------------------------------------------
3 -- GNAT RUN-TIME COMPONENTS --
5 -- ADA.NUMERICS.GENERIC_REAL_ARRAYS --
9 -- Copyright (C) 2006-2007, 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 System; use System;
35 with System.Generic_Real_BLAS;
36 with System.Generic_Real_LAPACK;
37 with System.Generic_Array_Operations; use System.Generic_Array_Operations;
39 package body Ada.Numerics.Generic_Real_Arrays is
41 -- Operations involving inner products use BLAS library implementations.
42 -- This allows larger matrices and vectors to be computed efficiently,
43 -- taking into account memory hierarchy issues and vector instructions
44 -- that vary widely between machines.
46 -- Operations that are defined in terms of operations on the type Real,
47 -- such as addition, subtraction and scaling, are computed in the canonical
48 -- way looping over all elements.
50 -- Operations for solving linear systems and computing determinant,
51 -- eigenvalues, eigensystem and inverse, are implemented using the
55 new Generic_Real_BLAS (Real'Base, Real_Vector, Real_Matrix);
58 new Generic_Real_LAPACK (Real'Base, Real_Vector, Real_Matrix);
62 -- Procedure versions of functions returning unconstrained values.
63 -- This allows for inlining the function wrapper.
65 procedure Eigenvalues (A : Real_Matrix; Values : out Real_Vector);
66 procedure Inverse (A : Real_Matrix; R : out Real_Matrix);
67 procedure Solve (A : Real_Matrix; X : Real_Vector; B : out Real_Vector);
68 procedure Solve (A : Real_Matrix; X : Real_Matrix; B : out Real_Matrix);
70 procedure Transpose is new
71 Generic_Array_Operations.Transpose
73 Matrix => Real_Matrix);
75 -- Helper function that raises a Constraint_Error is the argument is
76 -- not a square matrix, and otherwise returns its length.
78 function Length is new Square_Matrix_Length (Real'Base, Real_Matrix);
80 -- Instantiating the following subprograms directly would lead to
81 -- name clashes, so use a local package.
83 package Instantiations is
86 Vector_Elementwise_Operation
87 (X_Scalar => Real'Base,
88 Result_Scalar => Real'Base,
89 X_Vector => Real_Vector,
90 Result_Vector => Real_Vector,
94 Matrix_Elementwise_Operation
95 (X_Scalar => Real'Base,
96 Result_Scalar => Real'Base,
97 X_Matrix => Real_Matrix,
98 Result_Matrix => Real_Matrix,
102 Vector_Vector_Elementwise_Operation
103 (Left_Scalar => Real'Base,
104 Right_Scalar => Real'Base,
105 Result_Scalar => Real'Base,
106 Left_Vector => Real_Vector,
107 Right_Vector => Real_Vector,
108 Result_Vector => Real_Vector,
112 Matrix_Matrix_Elementwise_Operation
113 (Left_Scalar => Real'Base,
114 Right_Scalar => Real'Base,
115 Result_Scalar => Real'Base,
116 Left_Matrix => Real_Matrix,
117 Right_Matrix => Real_Matrix,
118 Result_Matrix => Real_Matrix,
122 Vector_Elementwise_Operation
123 (X_Scalar => Real'Base,
124 Result_Scalar => Real'Base,
125 X_Vector => Real_Vector,
126 Result_Vector => Real_Vector,
130 Matrix_Elementwise_Operation
131 (X_Scalar => Real'Base,
132 Result_Scalar => Real'Base,
133 X_Matrix => Real_Matrix,
134 Result_Matrix => Real_Matrix,
138 Vector_Vector_Elementwise_Operation
139 (Left_Scalar => Real'Base,
140 Right_Scalar => Real'Base,
141 Result_Scalar => Real'Base,
142 Left_Vector => Real_Vector,
143 Right_Vector => Real_Vector,
144 Result_Vector => Real_Vector,
148 Matrix_Matrix_Elementwise_Operation
149 (Left_Scalar => Real'Base,
150 Right_Scalar => Real'Base,
151 Result_Scalar => Real'Base,
152 Left_Matrix => Real_Matrix,
153 Right_Matrix => Real_Matrix,
154 Result_Matrix => Real_Matrix,
158 Scalar_Vector_Elementwise_Operation
159 (Left_Scalar => Real'Base,
160 Right_Scalar => Real'Base,
161 Result_Scalar => Real'Base,
162 Right_Vector => Real_Vector,
163 Result_Vector => Real_Vector,
167 Scalar_Matrix_Elementwise_Operation
168 (Left_Scalar => Real'Base,
169 Right_Scalar => Real'Base,
170 Result_Scalar => Real'Base,
171 Right_Matrix => Real_Matrix,
172 Result_Matrix => Real_Matrix,
176 Vector_Scalar_Elementwise_Operation
177 (Left_Scalar => Real'Base,
178 Right_Scalar => Real'Base,
179 Result_Scalar => Real'Base,
180 Left_Vector => Real_Vector,
181 Result_Vector => Real_Vector,
185 Matrix_Scalar_Elementwise_Operation
186 (Left_Scalar => Real'Base,
187 Right_Scalar => Real'Base,
188 Result_Scalar => Real'Base,
189 Left_Matrix => Real_Matrix,
190 Result_Matrix => Real_Matrix,
195 (Left_Scalar => Real'Base,
196 Right_Scalar => Real'Base,
197 Result_Scalar => Real'Base,
198 Left_Vector => Real_Vector,
199 Right_Vector => Real_Vector,
200 Matrix => Real_Matrix);
203 Vector_Scalar_Elementwise_Operation
204 (Left_Scalar => Real'Base,
205 Right_Scalar => Real'Base,
206 Result_Scalar => Real'Base,
207 Left_Vector => Real_Vector,
208 Result_Vector => Real_Vector,
212 Matrix_Scalar_Elementwise_Operation
213 (Left_Scalar => Real'Base,
214 Right_Scalar => Real'Base,
215 Result_Scalar => Real'Base,
216 Left_Matrix => Real_Matrix,
217 Result_Matrix => Real_Matrix,
220 function "abs" is new
221 Vector_Elementwise_Operation
222 (X_Scalar => Real'Base,
223 Result_Scalar => Real'Base,
224 X_Vector => Real_Vector,
225 Result_Vector => Real_Vector,
228 function "abs" is new
229 Matrix_Elementwise_Operation
230 (X_Scalar => Real'Base,
231 Result_Scalar => Real'Base,
232 X_Matrix => Real_Matrix,
233 Result_Matrix => Real_Matrix,
236 function Unit_Matrix is new
237 Generic_Array_Operations.Unit_Matrix
238 (Scalar => Real'Base,
239 Matrix => Real_Matrix,
243 function Unit_Vector is new
244 Generic_Array_Operations.Unit_Vector
245 (Scalar => Real'Base,
246 Vector => Real_Vector,
256 function "+" (Right : Real_Vector) return Real_Vector
257 renames Instantiations."+";
259 function "+" (Right : Real_Matrix) return Real_Matrix
260 renames Instantiations."+";
262 function "+" (Left, Right : Real_Vector) return Real_Vector
263 renames Instantiations."+";
265 function "+" (Left, Right : Real_Matrix) return Real_Matrix
266 renames Instantiations."+";
272 function "-" (Right : Real_Vector) return Real_Vector
273 renames Instantiations."-";
275 function "-" (Right : Real_Matrix) return Real_Matrix
276 renames Instantiations."-";
278 function "-" (Left, Right : Real_Vector) return Real_Vector
279 renames Instantiations."-";
281 function "-" (Left, Right : Real_Matrix) return Real_Matrix
282 renames Instantiations."-";
288 -- Scalar multiplication
290 function "*" (Left : Real'Base; Right : Real_Vector) return Real_Vector
291 renames Instantiations."*";
293 function "*" (Left : Real_Vector; Right : Real'Base) return Real_Vector
294 renames Instantiations."*";
296 function "*" (Left : Real'Base; Right : Real_Matrix) return Real_Matrix
297 renames Instantiations."*";
299 function "*" (Left : Real_Matrix; Right : Real'Base) return Real_Matrix
300 renames Instantiations."*";
302 -- Vector multiplication
304 function "*" (Left, Right : Real_Vector) return Real'Base is
306 if Left'Length /= Right'Length then
307 raise Constraint_Error with
308 "vectors are of different length in inner product";
311 return dot (Left'Length, X => Left, Y => Right);
314 function "*" (Left, Right : Real_Vector) return Real_Matrix
315 renames Instantiations."*";
319 Right : Real_Matrix) return Real_Vector
321 R : Real_Vector (Right'Range (2));
324 if Left'Length /= Right'Length (1) then
325 raise Constraint_Error with
326 "incompatible dimensions in vector-matrix multiplication";
329 gemv (Trans => No_Trans'Access,
330 M => Right'Length (2),
331 N => Right'Length (1),
333 Ld_A => Right'Length (2),
342 Right : Real_Vector) return Real_Vector
344 R : Real_Vector (Left'Range (1));
347 if Left'Length (2) /= Right'Length then
348 raise Constraint_Error with
349 "incompatible dimensions in matrix-vector multiplication";
352 gemv (Trans => Trans'Access,
353 M => Left'Length (2),
354 N => Left'Length (1),
356 Ld_A => Left'Length (2),
363 -- Matrix Multiplication
365 function "*" (Left, Right : Real_Matrix) return Real_Matrix is
366 R : Real_Matrix (Left'Range (1), Right'Range (2));
369 if Left'Length (2) /= Right'Length (1) then
370 raise Constraint_Error with
371 "incompatible dimensions in matrix-matrix multipication";
374 gemm (Trans_A => No_Trans'Access,
375 Trans_B => No_Trans'Access,
376 M => Right'Length (2),
377 N => Left'Length (1),
378 K => Right'Length (1),
380 Ld_A => Right'Length (2),
382 Ld_B => Left'Length (2),
384 Ld_C => R'Length (2));
393 function "/" (Left : Real_Vector; Right : Real'Base) return Real_Vector
394 renames Instantiations."/";
396 function "/" (Left : Real_Matrix; Right : Real'Base) return Real_Matrix
397 renames Instantiations."/";
403 function "abs" (Right : Real_Vector) return Real'Base is
405 return nrm2 (Right'Length, Right);
408 function "abs" (Right : Real_Vector) return Real_Vector
409 renames Instantiations."abs";
411 function "abs" (Right : Real_Matrix) return Real_Matrix
412 renames Instantiations."abs";
418 function Determinant (A : Real_Matrix) return Real'Base is
419 N : constant Integer := Length (A);
420 LU : Real_Matrix (1 .. N, 1 .. N) := A;
421 Piv : Integer_Vector (1 .. N);
422 Info : aliased Integer := -1;
431 Info => Info'Access);
434 raise Constraint_Error with "ill-conditioned matrix";
439 Det := -Det * LU (J, J);
441 Det := Det * LU (J, J);
452 procedure Eigensystem
454 Values : out Real_Vector;
455 Vectors : out Real_Matrix)
457 N : constant Natural := Length (A);
458 Tau : Real_Vector (1 .. N);
459 L_Work : Real_Vector (1 .. 1);
460 Info : aliased Integer;
462 E : Real_Vector (1 .. N);
463 pragma Warnings (Off, E);
466 if Values'Length /= N then
467 raise Constraint_Error with "wrong length for output vector";
474 -- Initialize working matrix and check for symmetric input matrix
476 Transpose (A, Vectors);
479 raise Argument_Error with "matrix not symmetric";
482 -- Compute size of additional working space
484 sytrd (Uplo => Lower'Access,
493 Info => Info'Access);
496 Work : Real_Vector (1 .. Integer'Max (Integer (L_Work (1)), 2 * N));
497 pragma Warnings (Off, Work);
499 Comp_Z : aliased constant Character := 'V';
502 -- Reduce matrix to tridiagonal form
504 sytrd (Uplo => Lower'Access,
507 Ld_A => A'Length (1),
512 L_Work => Work'Length,
513 Info => Info'Access);
519 -- Generate the real orthogonal matrix determined by sytrd
521 orgtr (Uplo => Lower'Access,
527 L_Work => Work'Length,
528 Info => Info'Access);
534 -- Compute all eigenvalues and eigenvectors using QR algorithm
536 steqr (Comp_Z => Comp_Z'Access,
543 Info => Info'Access);
546 raise Constraint_Error with
547 "eigensystem computation failed to converge";
556 procedure Eigenvalues
558 Values : out Real_Vector)
560 N : constant Natural := Length (A);
561 L_Work : Real_Vector (1 .. 1);
562 Info : aliased Integer;
564 B : Real_Matrix (1 .. N, 1 .. N);
565 Tau : Real_Vector (1 .. N);
566 E : Real_Vector (1 .. N);
567 pragma Warnings (Off, B);
568 pragma Warnings (Off, Tau);
569 pragma Warnings (Off, E);
572 if Values'Length /= N then
573 raise Constraint_Error with "wrong length for output vector";
580 -- Initialize working matrix and check for symmetric input matrix
585 raise Argument_Error with "matrix not symmetric";
588 -- Find size of work area
590 sytrd (Uplo => Lower'Access,
599 Info => Info'Access);
602 Work : Real_Vector (1 .. Integer'Min (Integer (L_Work (1)), 4 * N));
603 pragma Warnings (Off, Work);
606 -- Reduce matrix to tridiagonal form
608 sytrd (Uplo => Lower'Access,
611 Ld_A => A'Length (1),
616 L_Work => Work'Length,
617 Info => Info'Access);
620 raise Constraint_Error;
623 -- Compute all eigenvalues using QR algorithm
628 Info => Info'Access);
631 raise Constraint_Error with
632 "eigenvalues computation failed to converge";
637 function Eigenvalues (A : Real_Matrix) return Real_Vector is
638 R : Real_Vector (A'Range (1));
648 procedure Inverse (A : Real_Matrix; R : out Real_Matrix) is
649 N : constant Integer := Length (A);
650 Piv : Integer_Vector (1 .. N);
651 L_Work : Real_Vector (1 .. 1);
652 Info : aliased Integer := -1;
655 -- All computations are done using column-major order, but this works
656 -- out fine, because Transpose (Inverse (Transpose (A))) = Inverse (A).
660 -- Compute LU decomposition
667 Info => Info'Access);
670 raise Constraint_Error with "inverting singular matrix";
673 -- Determine size of work area
681 Info => Info'Access);
684 raise Constraint_Error;
688 Work : Real_Vector (1 .. Integer (L_Work (1)));
689 pragma Warnings (Off, Work);
692 -- Compute inverse from LU decomposition
699 L_Work => Work'Length,
700 Info => Info'Access);
703 raise Constraint_Error with "inverting singular matrix";
706 -- ??? Should iterate with gerfs, based on implementation advice
710 function Inverse (A : Real_Matrix) return Real_Matrix is
711 R : Real_Matrix (A'Range (2), A'Range (1));
721 procedure Solve (A : Real_Matrix; X : Real_Vector; B : out Real_Vector) is
723 if Length (A) /= X'Length then
724 raise Constraint_Error with
725 "incompatible matrix and vector dimensions";
728 -- ??? Should solve directly, is faster and more accurate
730 B := Inverse (A) * X;
733 procedure Solve (A : Real_Matrix; X : Real_Matrix; B : out Real_Matrix) is
735 if Length (A) /= X'Length (1) then
736 raise Constraint_Error with "incompatible matrix dimensions";
739 -- ??? Should solve directly, is faster and more accurate
741 B := Inverse (A) * X;
744 function Solve (A : Real_Matrix; X : Real_Vector) return Real_Vector is
745 B : Real_Vector (A'Range (2));
751 function Solve (A, X : Real_Matrix) return Real_Matrix is
752 B : Real_Matrix (A'Range (2), X'Range (2));
762 function Transpose (X : Real_Matrix) return Real_Matrix is
763 R : Real_Matrix (X'Range (2), X'Range (1));
776 First_1 : Integer := 1;
777 First_2 : Integer := 1) return Real_Matrix
778 renames Instantiations.Unit_Matrix;
787 First : Integer := 1) return Real_Vector
788 renames Instantiations.Unit_Vector;
790 end Ada.Numerics.Generic_Real_Arrays;