1 ------------------------------------------------------------------------------
3 -- GNAT RUN-TIME COMPONENTS --
5 -- S Y S T E M . G E N E R I C _ A R R A Y _ O P E R A T I O N S --
9 -- Copyright (C) 2006-2011, 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
.Numerics
; use Ada
.Numerics
;
34 package body System
.Generic_Array_Operations
is
36 -- The local function Check_Unit_Last computes the index of the last
37 -- element returned by Unit_Vector or Unit_Matrix. A separate function is
38 -- needed to allow raising Constraint_Error before declaring the function
39 -- result variable. The result variable needs to be declared first, to
40 -- allow front-end inlining.
42 function Check_Unit_Last
45 First
: Integer) return Integer;
46 pragma Inline_Always
(Check_Unit_Last
);
52 function Diagonal
(A
: Matrix
) return Vector
is
53 N
: constant Natural := Natural'Min (A
'Length (1), A
'Length (2));
54 R
: Vector
(A
'First (1) .. A
'First (1) + N
- 1);
57 for J
in 0 .. N
- 1 loop
58 R
(R
'First + J
) := A
(A
'First (1) + J
, A
'First (2) + J
);
64 --------------------------
65 -- Square_Matrix_Length --
66 --------------------------
68 function Square_Matrix_Length
(A
: Matrix
) return Natural is
70 if A
'Length (1) /= A
'Length (2) then
71 raise Constraint_Error
with "matrix is not square";
75 end Square_Matrix_Length
;
81 function Check_Unit_Last
84 First
: Integer) return Integer
87 -- Order the tests carefully to avoid overflow
90 or else First
> Integer'Last - Order
+ 1
91 or else Index
> First
+ (Order
- 1)
93 raise Constraint_Error
;
96 return First
+ (Order
- 1);
100 -- Back_Substitute --
101 ---------------------
103 procedure Back_Substitute
(M
, N
: in out Matrix
) is
104 pragma Assert
(M
'First (1) = N
'First (1)
106 M
'Last (1) = N
'Last (1));
113 -- Elementary row operation that subtracts Factor * M (Source, <>) from
123 for J
in M
'Range (2) loop
124 M
(Target
, J
) := M
(Target
, J
) - Factor
* M
(Source
, J
);
128 -- Local declarations
130 Max_Col
: Integer := M
'Last (2);
132 -- Start of processing for Back_Substitute
135 Do_Rows
: for Row
in reverse M
'Range (1) loop
136 Find_Non_Zero
: for Col
in reverse M
'First (2) .. Max_Col
loop
137 if Is_Non_Zero
(M
(Row
, Col
)) then
139 -- Found first non-zero element, so subtract a multiple of this
140 -- element from all higher rows, to reduce all other elements
141 -- in this column to zero.
144 -- We can't use a for loop, as we'd need to iterate to
145 -- Row - 1, but that expression will overflow if M'First
146 -- equals Integer'First, which is true for aggregates
147 -- without explicit bounds..
149 J
: Integer := M
'First (1);
153 Sub_Row
(N
, J
, Row
, (M
(J
, Col
) / M
(Row
, Col
)));
154 Sub_Row
(M
, J
, Row
, (M
(J
, Col
) / M
(Row
, Col
)));
159 -- Avoid potential overflow in the subtraction below
161 exit Do_Rows
when Col
= M
'First (2);
167 end loop Find_Non_Zero
;
171 -----------------------
172 -- Forward_Eliminate --
173 -----------------------
175 procedure Forward_Eliminate
180 pragma Assert
(M
'First (1) = N
'First (1)
182 M
'Last (1) = N
'Last (1));
184 -- The following are variations of the elementary matrix row operations:
185 -- row switching, row multiplication and row addition. Because in this
186 -- algorithm the addition factor is always a negated value, we chose to
187 -- use row subtraction instead. Similarly, instead of multiplying by
188 -- a reciprocal, we divide.
195 -- Subtrace Factor * M (Source, <>) from M (Target, <>)
198 (M
, N
: in out Matrix
;
201 -- Divide M (Row) and N (Row) by Scale, and update Det
204 (M
, N
: in out Matrix
;
207 -- Exchange M (Row_1) and N (Row_1) with M (Row_2) and N (Row_2),
208 -- negating Det in the process.
221 for J
in M
'Range (2) loop
222 M
(Target
, J
) := M
(Target
, J
) - Factor
* M
(Source
, J
);
231 (M
, N
: in out Matrix
;
238 for J
in M
'Range (2) loop
239 M
(Row
, J
) := M
(Row
, J
) / Scale
;
242 for J
in N
'Range (2) loop
243 N
(Row
- M
'First (1) + N
'First (1), J
) :=
244 N
(Row
- M
'First (1) + N
'First (1), J
) / Scale
;
253 (M
, N
: in out Matrix
;
257 procedure Swap
(X
, Y
: in out Scalar
);
258 -- Exchange the values of X and Y
260 procedure Swap
(X
, Y
: in out Scalar
) is
261 T
: constant Scalar
:= X
;
267 -- Start of processing for Switch_Row
270 if Row_1
/= Row_2
then
273 for J
in M
'Range (2) loop
274 Swap
(M
(Row_1
, J
), M
(Row_2
, J
));
277 for J
in N
'Range (2) loop
278 Swap
(N
(Row_1
- M
'First (1) + N
'First (1), J
),
279 N
(Row_2
- M
'First (1) + N
'First (1), J
));
284 -- Local declarations
286 Row
: Integer := M
'First (1);
288 -- Start of processing for Forward_Eliminate
293 for J
in M
'Range (2) loop
295 Max_Row
: Integer := Row
;
296 Max_Abs
: Real
'Base := 0.0;
299 -- Find best pivot in column J, starting in row Row
301 for K
in Row
.. M
'Last (1) loop
303 New_Abs
: constant Real
'Base := abs M
(K
, J
);
305 if Max_Abs
< New_Abs
then
312 if Max_Abs
> 0.0 then
313 Switch_Row
(M
, N
, Row
, Max_Row
);
314 Divide_Row
(M
, N
, Row
, M
(Row
, J
));
316 for U
in Row
+ 1 .. M
'Last (1) loop
317 Sub_Row
(N
, U
, Row
, M
(U
, J
));
318 Sub_Row
(M
, U
, Row
, M
(U
, J
));
321 exit when Row
>= M
'Last (1);
326 -- Set zero (note that we do not have literals)
332 end Forward_Eliminate
;
338 function Inner_Product
340 Right
: Right_Vector
) return Result_Scalar
342 R
: Result_Scalar
:= Zero
;
345 if Left
'Length /= Right
'Length then
346 raise Constraint_Error
with
347 "vectors are of different length in inner product";
350 for J
in Left
'Range loop
351 R
:= R
+ Left
(J
) * Right
(J
- Left
'First + Right
'First);
361 function L2_Norm
(X
: X_Vector
) return Result_Real
'Base is
362 Sum
: Result_Real
'Base := 0.0;
365 for J
in X
'Range loop
366 Sum
:= Sum
+ Result_Real
'Base (abs X
(J
))**2;
372 ----------------------------------
373 -- Matrix_Elementwise_Operation --
374 ----------------------------------
376 function Matrix_Elementwise_Operation
(X
: X_Matrix
) return Result_Matrix
is
377 R
: Result_Matrix
(X
'Range (1), X
'Range (2));
380 for J
in R
'Range (1) loop
381 for K
in R
'Range (2) loop
382 R
(J
, K
) := Operation
(X
(J
, K
));
387 end Matrix_Elementwise_Operation
;
389 ----------------------------------
390 -- Vector_Elementwise_Operation --
391 ----------------------------------
393 function Vector_Elementwise_Operation
(X
: X_Vector
) return Result_Vector
is
394 R
: Result_Vector
(X
'Range);
397 for J
in R
'Range loop
398 R
(J
) := Operation
(X
(J
));
402 end Vector_Elementwise_Operation
;
404 -----------------------------------------
405 -- Matrix_Matrix_Elementwise_Operation --
406 -----------------------------------------
408 function Matrix_Matrix_Elementwise_Operation
410 Right
: Right_Matrix
) return Result_Matrix
412 R
: Result_Matrix
(Left
'Range (1), Left
'Range (2));
415 if Left
'Length (1) /= Right
'Length (1)
417 Left
'Length (2) /= Right
'Length (2)
419 raise Constraint_Error
with
420 "matrices are of different dimension in elementwise operation";
423 for J
in R
'Range (1) loop
424 for K
in R
'Range (2) loop
429 (J
- R
'First (1) + Right
'First (1),
430 K
- R
'First (2) + Right
'First (2)));
435 end Matrix_Matrix_Elementwise_Operation
;
437 ------------------------------------------------
438 -- Matrix_Matrix_Scalar_Elementwise_Operation --
439 ------------------------------------------------
441 function Matrix_Matrix_Scalar_Elementwise_Operation
444 Z
: Z_Scalar
) return Result_Matrix
446 R
: Result_Matrix
(X
'Range (1), X
'Range (2));
449 if X
'Length (1) /= Y
'Length (1)
451 X
'Length (2) /= Y
'Length (2)
453 raise Constraint_Error
with
454 "matrices are of different dimension in elementwise operation";
457 for J
in R
'Range (1) loop
458 for K
in R
'Range (2) loop
462 Y
(J
- R
'First (1) + Y
'First (1),
463 K
- R
'First (2) + Y
'First (2)),
469 end Matrix_Matrix_Scalar_Elementwise_Operation
;
471 -----------------------------------------
472 -- Vector_Vector_Elementwise_Operation --
473 -----------------------------------------
475 function Vector_Vector_Elementwise_Operation
477 Right
: Right_Vector
) return Result_Vector
479 R
: Result_Vector
(Left
'Range);
482 if Left
'Length /= Right
'Length then
483 raise Constraint_Error
with
484 "vectors are of different length in elementwise operation";
487 for J
in R
'Range loop
488 R
(J
) := Operation
(Left
(J
), Right
(J
- R
'First + Right
'First));
492 end Vector_Vector_Elementwise_Operation
;
494 ------------------------------------------------
495 -- Vector_Vector_Scalar_Elementwise_Operation --
496 ------------------------------------------------
498 function Vector_Vector_Scalar_Elementwise_Operation
501 Z
: Z_Scalar
) return Result_Vector
503 R
: Result_Vector
(X
'Range);
506 if X
'Length /= Y
'Length then
507 raise Constraint_Error
with
508 "vectors are of different length in elementwise operation";
511 for J
in R
'Range loop
512 R
(J
) := Operation
(X
(J
), Y
(J
- X
'First + Y
'First), Z
);
516 end Vector_Vector_Scalar_Elementwise_Operation
;
518 -----------------------------------------
519 -- Matrix_Scalar_Elementwise_Operation --
520 -----------------------------------------
522 function Matrix_Scalar_Elementwise_Operation
524 Right
: Right_Scalar
) return Result_Matrix
526 R
: Result_Matrix
(Left
'Range (1), Left
'Range (2));
529 for J
in R
'Range (1) loop
530 for K
in R
'Range (2) loop
531 R
(J
, K
) := Operation
(Left
(J
, K
), Right
);
536 end Matrix_Scalar_Elementwise_Operation
;
538 -----------------------------------------
539 -- Vector_Scalar_Elementwise_Operation --
540 -----------------------------------------
542 function Vector_Scalar_Elementwise_Operation
544 Right
: Right_Scalar
) return Result_Vector
546 R
: Result_Vector
(Left
'Range);
549 for J
in R
'Range loop
550 R
(J
) := Operation
(Left
(J
), Right
);
554 end Vector_Scalar_Elementwise_Operation
;
556 -----------------------------------------
557 -- Scalar_Matrix_Elementwise_Operation --
558 -----------------------------------------
560 function Scalar_Matrix_Elementwise_Operation
562 Right
: Right_Matrix
) return Result_Matrix
564 R
: Result_Matrix
(Right
'Range (1), Right
'Range (2));
567 for J
in R
'Range (1) loop
568 for K
in R
'Range (2) loop
569 R
(J
, K
) := Operation
(Left
, Right
(J
, K
));
574 end Scalar_Matrix_Elementwise_Operation
;
576 -----------------------------------------
577 -- Scalar_Vector_Elementwise_Operation --
578 -----------------------------------------
580 function Scalar_Vector_Elementwise_Operation
582 Right
: Right_Vector
) return Result_Vector
584 R
: Result_Vector
(Right
'Range);
587 for J
in R
'Range loop
588 R
(J
) := Operation
(Left
, Right
(J
));
592 end Scalar_Vector_Elementwise_Operation
;
598 function Sqrt
(X
: Real
'Base) return Real
'Base is
599 Root
, Next
: Real
'Base;
602 -- Be defensive: any comparisons with NaN values will yield False.
604 if not (X
> 0.0) then
608 raise Argument_Error
;
611 elsif X
> Real
'Base'Last then
613 -- X is infinity, which is its own square root
618 -- Compute an initial estimate based on:
620 -- X = M * R**E and Sqrt (X) = Sqrt (M) * R**(E / 2.0),
622 -- where M is the mantissa, R is the radix and E the exponent.
624 -- By ignoring the mantissa and ignoring the case of an odd
625 -- exponent, we get a final error that is at most R. In other words,
626 -- the result has about a single bit precision.
628 Root := Real'Base (Real'Machine_Radix) ** (Real'Exponent (X) / 2);
630 -- Because of the poor initial estimate, use the Babylonian method of
631 -- computing the square root, as it is stable for all inputs. Every step
632 -- will roughly double the precision of the result. Just a few steps
633 -- suffice in most cases. Eight iterations should give about 2**8 bits
637 Next := (Root + X / Root) / 2.0;
638 exit when Root = Next;
645 ---------------------------
646 -- Matrix_Matrix_Product --
647 ---------------------------
649 function Matrix_Matrix_Product
651 Right : Right_Matrix) return Result_Matrix
653 R : Result_Matrix (Left'Range (1), Right'Range (2));
656 if Left'Length (2) /= Right'Length (1) then
657 raise Constraint_Error with
658 "incompatible dimensions in matrix multiplication";
661 for J in R'Range (1) loop
662 for K in R'Range (2) loop
664 S : Result_Scalar := Zero;
667 for M in Left'Range (2) loop
668 S := S + Left (J, M) *
669 Right (M - Left'First (2) + Right'First (1), K);
678 end Matrix_Matrix_Product;
680 ----------------------------
681 -- Matrix_Vector_Solution --
682 ----------------------------
684 function Matrix_Vector_Solution (A : Matrix; X : Vector) return Vector is
685 N : constant Natural := A'Length (1);
687 MX : Matrix (A'Range (1), 1 .. 1);
688 R : Vector (A'Range (2));
692 if A'Length (2) /= N then
693 raise Constraint_Error with "matrix is not square";
696 if X'Length /= N then
697 raise Constraint_Error with "incompatible vector length";
700 for J in 0 .. MX'Length (1) - 1 loop
701 MX (MX'First (1) + J, 1) := X (X'First + J);
704 Forward_Eliminate (MA, MX, Det);
705 Back_Substitute (MA, MX);
707 for J in 0 .. R'Length - 1 loop
708 R (R'First + J) := MX (MX'First (1) + J, 1);
712 end Matrix_Vector_Solution;
714 ----------------------------
715 -- Matrix_Matrix_Solution --
716 ----------------------------
718 function Matrix_Matrix_Solution (A, X : Matrix) return Matrix is
719 N : constant Natural := A'Length (1);
720 MA : Matrix (A'Range (2), A'Range (2));
721 MB : Matrix (A'Range (2), X'Range (2));
725 if A'Length (2) /= N then
726 raise Constraint_Error with "matrix is not square";
729 if X'Length (1) /= N then
730 raise Constraint_Error with "matrices have unequal number of rows";
733 for J in 0 .. A'Length (1) - 1 loop
734 for K in MA'Range (2) loop
735 MA (MA'First (1) + J, K) := A (A'First (1) + J, K);
738 for K in MB'Range (2) loop
739 MB (MB'First (1) + J, K) := X (X'First (1) + J, K);
743 Forward_Eliminate (MA, MB, Det);
744 Back_Substitute (MA, MB);
747 end Matrix_Matrix_Solution;
749 ---------------------------
750 -- Matrix_Vector_Product --
751 ---------------------------
753 function Matrix_Vector_Product
755 Right : Right_Vector) return Result_Vector
757 R : Result_Vector (Left'Range (1));
760 if Left'Length (2) /= Right'Length then
761 raise Constraint_Error with
762 "incompatible dimensions in matrix-vector multiplication";
765 for J in Left'Range (1) loop
767 S : Result_Scalar := Zero;
770 for K in Left'Range (2) loop
771 S := S + Left (J, K) * Right (K - Left'First (2) + Right'First);
779 end Matrix_Vector_Product;
785 function Outer_Product
787 Right : Right_Vector) return Matrix
789 R : Matrix (Left'Range, Right'Range);
792 for J in R'Range (1) loop
793 for K in R'Range (2) loop
794 R (J, K) := Left (J) * Right (K);
805 procedure Swap_Column (A : in out Matrix; Left, Right : Integer) is
808 for J in A'Range (1) loop
810 A (J, Left) := A (J, Right);
811 A (J, Right) := Temp;
819 procedure Transpose (A : Matrix; R : out Matrix) is
821 for J in R'Range (1) loop
822 for K in R'Range (2) loop
823 R (J, K) := A (K - R'First (2) + A'First (1),
824 J - R'First (1) + A'First (2));
829 -------------------------------
830 -- Update_Matrix_With_Matrix --
831 -------------------------------
833 procedure Update_Matrix_With_Matrix (X : in out X_Matrix; Y : Y_Matrix) is
835 if X'Length (1) /= Y'Length (1)
836 or else X'Length (2) /= Y'Length (2)
838 raise Constraint_Error with
839 "matrices are of different dimension in update operation";
842 for J in X'Range (1) loop
843 for K in X'Range (2) loop
844 Update (X (J, K), Y (J - X'First (1) + Y'First (1),
845 K - X'First (2) + Y'First (2)));
848 end Update_Matrix_With_Matrix;
850 -------------------------------
851 -- Update_Vector_With_Vector --
852 -------------------------------
854 procedure Update_Vector_With_Vector (X : in out X_Vector; Y : Y_Vector) is
856 if X'Length /= Y'Length then
857 raise Constraint_Error with
858 "vectors are of different length in update operation";
861 for J in X'Range loop
862 Update (X (J), Y (J - X'First + Y'First));
864 end Update_Vector_With_Vector;
872 First_1 : Integer := 1;
873 First_2 : Integer := 1) return Matrix
875 R : Matrix (First_1 .. Check_Unit_Last (First_1, Order, First_1),
876 First_2 .. Check_Unit_Last (First_2, Order, First_2));
879 R := (others => (others => Zero));
881 for J in 0 .. Order - 1 loop
882 R (First_1 + J, First_2 + J) := One;
895 First : Integer := 1) return Vector
897 R : Vector (First .. Check_Unit_Last (Index, Order, First));
899 R := (others => Zero);
904 ---------------------------
905 -- Vector_Matrix_Product --
906 ---------------------------
908 function Vector_Matrix_Product
910 Right : Matrix) return Result_Vector
912 R : Result_Vector (Right'Range (2));
915 if Left'Length /= Right'Length (2) then
916 raise Constraint_Error with
917 "incompatible dimensions in vector-matrix multiplication";
920 for J in Right'Range (2) loop
922 S : Result_Scalar := Zero;
925 for K in Right'Range (1) loop
926 S := S + Left (J - Right'First (1) + Left'First) * Right (K, J);
934 end Vector_Matrix_Product;
936 end System.Generic_Array_Operations;