|
111
|
1 ------------------------------------------------------------------------------
|
|
|
2 -- --
|
|
|
3 -- GNAT RUN-TIME COMPONENTS --
|
|
|
4 -- --
|
|
|
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 --
|
|
|
6 -- --
|
|
|
7 -- B o d y --
|
|
|
8 -- --
|
|
145
|
9 -- Copyright (C) 2006-2019, Free Software Foundation, Inc. --
|
|
111
|
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 ------------------------------------------------------------------------------
|
|
|
31
|
|
|
32 with Ada.Numerics; use Ada.Numerics;
|
|
|
33 package body System.Generic_Array_Operations is
|
|
|
34 function Check_Unit_Last
|
|
|
35 (Index : Integer;
|
|
|
36 Order : Positive;
|
|
|
37 First : Integer) return Integer;
|
|
|
38 pragma Inline_Always (Check_Unit_Last);
|
|
|
39 -- Compute index of last element returned by Unit_Vector or Unit_Matrix.
|
|
|
40 -- A separate function is needed to allow raising Constraint_Error before
|
|
|
41 -- declaring the function result variable. The result variable needs to be
|
|
|
42 -- declared first, to allow front-end inlining.
|
|
|
43
|
|
|
44 --------------
|
|
|
45 -- Diagonal --
|
|
|
46 --------------
|
|
|
47
|
|
|
48 function Diagonal (A : Matrix) return Vector is
|
|
|
49 N : constant Natural := Natural'Min (A'Length (1), A'Length (2));
|
|
|
50 begin
|
|
|
51 return R : Vector (A'First (1) .. A'First (1) + N - 1) do
|
|
|
52 for J in 0 .. N - 1 loop
|
|
|
53 R (R'First + J) := A (A'First (1) + J, A'First (2) + J);
|
|
|
54 end loop;
|
|
|
55 end return;
|
|
|
56 end Diagonal;
|
|
|
57
|
|
|
58 --------------------------
|
|
|
59 -- Square_Matrix_Length --
|
|
|
60 --------------------------
|
|
|
61
|
|
|
62 function Square_Matrix_Length (A : Matrix) return Natural is
|
|
|
63 begin
|
|
|
64 if A'Length (1) /= A'Length (2) then
|
|
|
65 raise Constraint_Error with "matrix is not square";
|
|
|
66 else
|
|
|
67 return A'Length (1);
|
|
|
68 end if;
|
|
|
69 end Square_Matrix_Length;
|
|
|
70
|
|
|
71 ---------------------
|
|
|
72 -- Check_Unit_Last --
|
|
|
73 ---------------------
|
|
|
74
|
|
|
75 function Check_Unit_Last
|
|
|
76 (Index : Integer;
|
|
|
77 Order : Positive;
|
|
|
78 First : Integer) return Integer
|
|
|
79 is
|
|
|
80 begin
|
|
|
81 -- Order the tests carefully to avoid overflow
|
|
|
82
|
|
|
83 if Index < First
|
|
|
84 or else First > Integer'Last - Order + 1
|
|
|
85 or else Index > First + (Order - 1)
|
|
|
86 then
|
|
|
87 raise Constraint_Error;
|
|
|
88 end if;
|
|
|
89
|
|
|
90 return First + (Order - 1);
|
|
|
91 end Check_Unit_Last;
|
|
|
92
|
|
|
93 ---------------------
|
|
|
94 -- Back_Substitute --
|
|
|
95 ---------------------
|
|
|
96
|
|
|
97 procedure Back_Substitute (M, N : in out Matrix) is
|
|
|
98 pragma Assert (M'First (1) = N'First (1)
|
|
|
99 and then
|
|
|
100 M'Last (1) = N'Last (1));
|
|
|
101
|
|
|
102 procedure Sub_Row
|
|
|
103 (M : in out Matrix;
|
|
|
104 Target : Integer;
|
|
|
105 Source : Integer;
|
|
|
106 Factor : Scalar);
|
|
|
107 -- Elementary row operation that subtracts Factor * M (Source, <>) from
|
|
|
108 -- M (Target, <>)
|
|
|
109
|
|
|
110 -------------
|
|
|
111 -- Sub_Row --
|
|
|
112 -------------
|
|
|
113
|
|
|
114 procedure Sub_Row
|
|
|
115 (M : in out Matrix;
|
|
|
116 Target : Integer;
|
|
|
117 Source : Integer;
|
|
|
118 Factor : Scalar)
|
|
|
119 is
|
|
|
120 begin
|
|
|
121 for J in M'Range (2) loop
|
|
|
122 M (Target, J) := M (Target, J) - Factor * M (Source, J);
|
|
|
123 end loop;
|
|
|
124 end Sub_Row;
|
|
|
125
|
|
|
126 -- Local declarations
|
|
|
127
|
|
|
128 Max_Col : Integer := M'Last (2);
|
|
|
129
|
|
|
130 -- Start of processing for Back_Substitute
|
|
|
131
|
|
|
132 begin
|
|
|
133 Do_Rows : for Row in reverse M'Range (1) loop
|
|
|
134 Find_Non_Zero : for Col in reverse M'First (2) .. Max_Col loop
|
|
|
135 if Is_Non_Zero (M (Row, Col)) then
|
|
|
136
|
|
|
137 -- Found first non-zero element, so subtract a multiple of this
|
|
|
138 -- element from all higher rows, to reduce all other elements
|
|
|
139 -- in this column to zero.
|
|
|
140
|
|
|
141 declare
|
|
|
142 -- We can't use a for loop, as we'd need to iterate to
|
|
|
143 -- Row - 1, but that expression will overflow if M'First
|
|
|
144 -- equals Integer'First, which is true for aggregates
|
|
|
145 -- without explicit bounds..
|
|
|
146
|
|
|
147 J : Integer := M'First (1);
|
|
|
148
|
|
|
149 begin
|
|
|
150 while J < Row loop
|
|
|
151 Sub_Row (N, J, Row, (M (J, Col) / M (Row, Col)));
|
|
|
152 Sub_Row (M, J, Row, (M (J, Col) / M (Row, Col)));
|
|
|
153 J := J + 1;
|
|
|
154 end loop;
|
|
|
155 end;
|
|
|
156
|
|
|
157 -- Avoid potential overflow in the subtraction below
|
|
|
158
|
|
|
159 exit Do_Rows when Col = M'First (2);
|
|
|
160
|
|
|
161 Max_Col := Col - 1;
|
|
|
162
|
|
|
163 exit Find_Non_Zero;
|
|
|
164 end if;
|
|
|
165 end loop Find_Non_Zero;
|
|
|
166 end loop Do_Rows;
|
|
|
167 end Back_Substitute;
|
|
|
168
|
|
|
169 -----------------------
|
|
|
170 -- Forward_Eliminate --
|
|
|
171 -----------------------
|
|
|
172
|
|
|
173 procedure Forward_Eliminate
|
|
|
174 (M : in out Matrix;
|
|
|
175 N : in out Matrix;
|
|
|
176 Det : out Scalar)
|
|
|
177 is
|
|
|
178 pragma Assert (M'First (1) = N'First (1)
|
|
|
179 and then
|
|
|
180 M'Last (1) = N'Last (1));
|
|
|
181
|
|
|
182 -- The following are variations of the elementary matrix row operations:
|
|
|
183 -- row switching, row multiplication and row addition. Because in this
|
|
|
184 -- algorithm the addition factor is always a negated value, we chose to
|
|
|
185 -- use row subtraction instead. Similarly, instead of multiplying by
|
|
|
186 -- a reciprocal, we divide.
|
|
|
187
|
|
|
188 procedure Sub_Row
|
|
|
189 (M : in out Matrix;
|
|
|
190 Target : Integer;
|
|
|
191 Source : Integer;
|
|
|
192 Factor : Scalar);
|
|
|
193 -- Subtrace Factor * M (Source, <>) from M (Target, <>)
|
|
|
194
|
|
|
195 procedure Divide_Row
|
|
|
196 (M, N : in out Matrix;
|
|
|
197 Row : Integer;
|
|
|
198 Scale : Scalar);
|
|
|
199 -- Divide M (Row) and N (Row) by Scale, and update Det
|
|
|
200
|
|
|
201 procedure Switch_Row
|
|
|
202 (M, N : in out Matrix;
|
|
|
203 Row_1 : Integer;
|
|
|
204 Row_2 : Integer);
|
|
|
205 -- Exchange M (Row_1) and N (Row_1) with M (Row_2) and N (Row_2),
|
|
|
206 -- negating Det in the process.
|
|
|
207
|
|
|
208 -------------
|
|
|
209 -- Sub_Row --
|
|
|
210 -------------
|
|
|
211
|
|
|
212 procedure Sub_Row
|
|
|
213 (M : in out Matrix;
|
|
|
214 Target : Integer;
|
|
|
215 Source : Integer;
|
|
|
216 Factor : Scalar)
|
|
|
217 is
|
|
|
218 begin
|
|
|
219 for J in M'Range (2) loop
|
|
|
220 M (Target, J) := M (Target, J) - Factor * M (Source, J);
|
|
|
221 end loop;
|
|
|
222 end Sub_Row;
|
|
|
223
|
|
|
224 ----------------
|
|
|
225 -- Divide_Row --
|
|
|
226 ----------------
|
|
|
227
|
|
|
228 procedure Divide_Row
|
|
|
229 (M, N : in out Matrix;
|
|
|
230 Row : Integer;
|
|
|
231 Scale : Scalar)
|
|
|
232 is
|
|
|
233 begin
|
|
|
234 Det := Det * Scale;
|
|
|
235
|
|
|
236 for J in M'Range (2) loop
|
|
|
237 M (Row, J) := M (Row, J) / Scale;
|
|
|
238 end loop;
|
|
|
239
|
|
|
240 for J in N'Range (2) loop
|
|
|
241 N (Row - M'First (1) + N'First (1), J) :=
|
|
|
242 N (Row - M'First (1) + N'First (1), J) / Scale;
|
|
|
243 end loop;
|
|
|
244 end Divide_Row;
|
|
|
245
|
|
|
246 ----------------
|
|
|
247 -- Switch_Row --
|
|
|
248 ----------------
|
|
|
249
|
|
|
250 procedure Switch_Row
|
|
|
251 (M, N : in out Matrix;
|
|
|
252 Row_1 : Integer;
|
|
|
253 Row_2 : Integer)
|
|
|
254 is
|
|
|
255 procedure Swap (X, Y : in out Scalar);
|
|
|
256 -- Exchange the values of X and Y
|
|
|
257
|
|
|
258 ----------
|
|
|
259 -- Swap --
|
|
|
260 ----------
|
|
|
261
|
|
|
262 procedure Swap (X, Y : in out Scalar) is
|
|
|
263 T : constant Scalar := X;
|
|
|
264 begin
|
|
|
265 X := Y;
|
|
|
266 Y := T;
|
|
|
267 end Swap;
|
|
|
268
|
|
|
269 -- Start of processing for Switch_Row
|
|
|
270
|
|
|
271 begin
|
|
|
272 if Row_1 /= Row_2 then
|
|
|
273 Det := Zero - Det;
|
|
|
274
|
|
|
275 for J in M'Range (2) loop
|
|
|
276 Swap (M (Row_1, J), M (Row_2, J));
|
|
|
277 end loop;
|
|
|
278
|
|
|
279 for J in N'Range (2) loop
|
|
|
280 Swap (N (Row_1 - M'First (1) + N'First (1), J),
|
|
|
281 N (Row_2 - M'First (1) + N'First (1), J));
|
|
|
282 end loop;
|
|
|
283 end if;
|
|
|
284 end Switch_Row;
|
|
|
285
|
|
|
286 -- Local declarations
|
|
|
287
|
|
|
288 Row : Integer := M'First (1);
|
|
|
289
|
|
|
290 -- Start of processing for Forward_Eliminate
|
|
|
291
|
|
|
292 begin
|
|
|
293 Det := One;
|
|
|
294
|
|
|
295 for J in M'Range (2) loop
|
|
|
296 declare
|
|
|
297 Max_Row : Integer := Row;
|
|
|
298 Max_Abs : Real'Base := 0.0;
|
|
|
299
|
|
|
300 begin
|
|
|
301 -- Find best pivot in column J, starting in row Row
|
|
|
302
|
|
|
303 for K in Row .. M'Last (1) loop
|
|
|
304 declare
|
|
|
305 New_Abs : constant Real'Base := abs M (K, J);
|
|
|
306 begin
|
|
|
307 if Max_Abs < New_Abs then
|
|
|
308 Max_Abs := New_Abs;
|
|
|
309 Max_Row := K;
|
|
|
310 end if;
|
|
|
311 end;
|
|
|
312 end loop;
|
|
|
313
|
|
|
314 if Max_Abs > 0.0 then
|
|
|
315 Switch_Row (M, N, Row, Max_Row);
|
|
|
316
|
|
|
317 -- The temporaries below are necessary to force a copy of the
|
|
|
318 -- value and avoid improper aliasing.
|
|
|
319
|
|
|
320 declare
|
|
|
321 Scale : constant Scalar := M (Row, J);
|
|
|
322 begin
|
|
|
323 Divide_Row (M, N, Row, Scale);
|
|
|
324 end;
|
|
|
325
|
|
|
326 for U in Row + 1 .. M'Last (1) loop
|
|
|
327 declare
|
|
|
328 Factor : constant Scalar := M (U, J);
|
|
|
329 begin
|
|
|
330 Sub_Row (N, U, Row, Factor);
|
|
|
331 Sub_Row (M, U, Row, Factor);
|
|
|
332 end;
|
|
|
333 end loop;
|
|
|
334
|
|
|
335 exit when Row >= M'Last (1);
|
|
|
336
|
|
|
337 Row := Row + 1;
|
|
|
338
|
|
|
339 else
|
|
|
340 -- Set zero (note that we do not have literals)
|
|
|
341
|
|
|
342 Det := Zero;
|
|
|
343 end if;
|
|
|
344 end;
|
|
|
345 end loop;
|
|
|
346 end Forward_Eliminate;
|
|
|
347
|
|
|
348 -------------------
|
|
|
349 -- Inner_Product --
|
|
|
350 -------------------
|
|
|
351
|
|
|
352 function Inner_Product
|
|
|
353 (Left : Left_Vector;
|
|
|
354 Right : Right_Vector) return Result_Scalar
|
|
|
355 is
|
|
|
356 R : Result_Scalar := Zero;
|
|
|
357
|
|
|
358 begin
|
|
|
359 if Left'Length /= Right'Length then
|
|
|
360 raise Constraint_Error with
|
|
|
361 "vectors are of different length in inner product";
|
|
|
362 end if;
|
|
|
363
|
|
|
364 for J in Left'Range loop
|
|
|
365 R := R + Left (J) * Right (J - Left'First + Right'First);
|
|
|
366 end loop;
|
|
|
367
|
|
|
368 return R;
|
|
|
369 end Inner_Product;
|
|
|
370
|
|
|
371 -------------
|
|
|
372 -- L2_Norm --
|
|
|
373 -------------
|
|
|
374
|
|
|
375 function L2_Norm (X : X_Vector) return Result_Real'Base is
|
|
|
376 Sum : Result_Real'Base := 0.0;
|
|
|
377
|
|
|
378 begin
|
|
|
379 for J in X'Range loop
|
|
|
380 Sum := Sum + Result_Real'Base (abs X (J))**2;
|
|
|
381 end loop;
|
|
|
382
|
|
|
383 return Sqrt (Sum);
|
|
|
384 end L2_Norm;
|
|
|
385
|
|
|
386 ----------------------------------
|
|
|
387 -- Matrix_Elementwise_Operation --
|
|
|
388 ----------------------------------
|
|
|
389
|
|
|
390 function Matrix_Elementwise_Operation (X : X_Matrix) return Result_Matrix is
|
|
|
391 begin
|
|
|
392 return R : Result_Matrix (X'Range (1), X'Range (2)) do
|
|
|
393 for J in R'Range (1) loop
|
|
|
394 for K in R'Range (2) loop
|
|
|
395 R (J, K) := Operation (X (J, K));
|
|
|
396 end loop;
|
|
|
397 end loop;
|
|
|
398 end return;
|
|
|
399 end Matrix_Elementwise_Operation;
|
|
|
400
|
|
|
401 ----------------------------------
|
|
|
402 -- Vector_Elementwise_Operation --
|
|
|
403 ----------------------------------
|
|
|
404
|
|
|
405 function Vector_Elementwise_Operation (X : X_Vector) return Result_Vector is
|
|
|
406 begin
|
|
|
407 return R : Result_Vector (X'Range) do
|
|
|
408 for J in R'Range loop
|
|
|
409 R (J) := Operation (X (J));
|
|
|
410 end loop;
|
|
|
411 end return;
|
|
|
412 end Vector_Elementwise_Operation;
|
|
|
413
|
|
|
414 -----------------------------------------
|
|
|
415 -- Matrix_Matrix_Elementwise_Operation --
|
|
|
416 -----------------------------------------
|
|
|
417
|
|
|
418 function Matrix_Matrix_Elementwise_Operation
|
|
|
419 (Left : Left_Matrix;
|
|
|
420 Right : Right_Matrix) return Result_Matrix
|
|
|
421 is
|
|
|
422 begin
|
|
|
423 return R : Result_Matrix (Left'Range (1), Left'Range (2)) do
|
|
|
424 if Left'Length (1) /= Right'Length (1)
|
|
|
425 or else
|
|
|
426 Left'Length (2) /= Right'Length (2)
|
|
|
427 then
|
|
|
428 raise Constraint_Error with
|
|
|
429 "matrices are of different dimension in elementwise operation";
|
|
|
430 end if;
|
|
|
431
|
|
|
432 for J in R'Range (1) loop
|
|
|
433 for K in R'Range (2) loop
|
|
|
434 R (J, K) :=
|
|
|
435 Operation
|
|
|
436 (Left (J, K),
|
|
|
437 Right
|
|
|
438 (J - R'First (1) + Right'First (1),
|
|
|
439 K - R'First (2) + Right'First (2)));
|
|
|
440 end loop;
|
|
|
441 end loop;
|
|
|
442 end return;
|
|
|
443 end Matrix_Matrix_Elementwise_Operation;
|
|
|
444
|
|
|
445 ------------------------------------------------
|
|
|
446 -- Matrix_Matrix_Scalar_Elementwise_Operation --
|
|
|
447 ------------------------------------------------
|
|
|
448
|
|
|
449 function Matrix_Matrix_Scalar_Elementwise_Operation
|
|
|
450 (X : X_Matrix;
|
|
|
451 Y : Y_Matrix;
|
|
|
452 Z : Z_Scalar) return Result_Matrix
|
|
|
453 is
|
|
|
454 begin
|
|
|
455 return R : Result_Matrix (X'Range (1), X'Range (2)) do
|
|
|
456 if X'Length (1) /= Y'Length (1)
|
|
|
457 or else
|
|
|
458 X'Length (2) /= Y'Length (2)
|
|
|
459 then
|
|
|
460 raise Constraint_Error with
|
|
|
461 "matrices are of different dimension in elementwise operation";
|
|
|
462 end if;
|
|
|
463
|
|
|
464 for J in R'Range (1) loop
|
|
|
465 for K in R'Range (2) loop
|
|
|
466 R (J, K) :=
|
|
|
467 Operation
|
|
|
468 (X (J, K),
|
|
|
469 Y (J - R'First (1) + Y'First (1),
|
|
|
470 K - R'First (2) + Y'First (2)),
|
|
|
471 Z);
|
|
|
472 end loop;
|
|
|
473 end loop;
|
|
|
474 end return;
|
|
|
475 end Matrix_Matrix_Scalar_Elementwise_Operation;
|
|
|
476
|
|
|
477 -----------------------------------------
|
|
|
478 -- Vector_Vector_Elementwise_Operation --
|
|
|
479 -----------------------------------------
|
|
|
480
|
|
|
481 function Vector_Vector_Elementwise_Operation
|
|
|
482 (Left : Left_Vector;
|
|
|
483 Right : Right_Vector) return Result_Vector
|
|
|
484 is
|
|
|
485 begin
|
|
|
486 return R : Result_Vector (Left'Range) do
|
|
|
487 if Left'Length /= Right'Length then
|
|
|
488 raise Constraint_Error with
|
|
|
489 "vectors are of different length in elementwise operation";
|
|
|
490 end if;
|
|
|
491
|
|
|
492 for J in R'Range loop
|
|
|
493 R (J) := Operation (Left (J), Right (J - R'First + Right'First));
|
|
|
494 end loop;
|
|
|
495 end return;
|
|
|
496 end Vector_Vector_Elementwise_Operation;
|
|
|
497
|
|
|
498 ------------------------------------------------
|
|
|
499 -- Vector_Vector_Scalar_Elementwise_Operation --
|
|
|
500 ------------------------------------------------
|
|
|
501
|
|
|
502 function Vector_Vector_Scalar_Elementwise_Operation
|
|
|
503 (X : X_Vector;
|
|
|
504 Y : Y_Vector;
|
|
|
505 Z : Z_Scalar) return Result_Vector is
|
|
|
506 begin
|
|
|
507 return R : Result_Vector (X'Range) do
|
|
|
508 if X'Length /= Y'Length then
|
|
|
509 raise Constraint_Error with
|
|
|
510 "vectors are of different length in elementwise operation";
|
|
|
511 end if;
|
|
|
512
|
|
|
513 for J in R'Range loop
|
|
|
514 R (J) := Operation (X (J), Y (J - X'First + Y'First), Z);
|
|
|
515 end loop;
|
|
|
516 end return;
|
|
|
517 end Vector_Vector_Scalar_Elementwise_Operation;
|
|
|
518
|
|
|
519 -----------------------------------------
|
|
|
520 -- Matrix_Scalar_Elementwise_Operation --
|
|
|
521 -----------------------------------------
|
|
|
522
|
|
|
523 function Matrix_Scalar_Elementwise_Operation
|
|
|
524 (Left : Left_Matrix;
|
|
|
525 Right : Right_Scalar) return Result_Matrix
|
|
|
526 is
|
|
|
527 begin
|
|
|
528 return R : Result_Matrix (Left'Range (1), Left'Range (2)) do
|
|
|
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);
|
|
|
532 end loop;
|
|
|
533 end loop;
|
|
|
534 end return;
|
|
|
535 end Matrix_Scalar_Elementwise_Operation;
|
|
|
536
|
|
|
537 -----------------------------------------
|
|
|
538 -- Vector_Scalar_Elementwise_Operation --
|
|
|
539 -----------------------------------------
|
|
|
540
|
|
|
541 function Vector_Scalar_Elementwise_Operation
|
|
|
542 (Left : Left_Vector;
|
|
|
543 Right : Right_Scalar) return Result_Vector
|
|
|
544 is
|
|
|
545 begin
|
|
|
546 return R : Result_Vector (Left'Range) do
|
|
|
547 for J in R'Range loop
|
|
|
548 R (J) := Operation (Left (J), Right);
|
|
|
549 end loop;
|
|
|
550 end return;
|
|
|
551 end Vector_Scalar_Elementwise_Operation;
|
|
|
552
|
|
|
553 -----------------------------------------
|
|
|
554 -- Scalar_Matrix_Elementwise_Operation --
|
|
|
555 -----------------------------------------
|
|
|
556
|
|
|
557 function Scalar_Matrix_Elementwise_Operation
|
|
|
558 (Left : Left_Scalar;
|
|
|
559 Right : Right_Matrix) return Result_Matrix
|
|
|
560 is
|
|
|
561 begin
|
|
|
562 return R : Result_Matrix (Right'Range (1), Right'Range (2)) do
|
|
|
563 for J in R'Range (1) loop
|
|
|
564 for K in R'Range (2) loop
|
|
|
565 R (J, K) := Operation (Left, Right (J, K));
|
|
|
566 end loop;
|
|
|
567 end loop;
|
|
|
568 end return;
|
|
|
569 end Scalar_Matrix_Elementwise_Operation;
|
|
|
570
|
|
|
571 -----------------------------------------
|
|
|
572 -- Scalar_Vector_Elementwise_Operation --
|
|
|
573 -----------------------------------------
|
|
|
574
|
|
|
575 function Scalar_Vector_Elementwise_Operation
|
|
|
576 (Left : Left_Scalar;
|
|
|
577 Right : Right_Vector) return Result_Vector
|
|
|
578 is
|
|
|
579 begin
|
|
|
580 return R : Result_Vector (Right'Range) do
|
|
|
581 for J in R'Range loop
|
|
|
582 R (J) := Operation (Left, Right (J));
|
|
|
583 end loop;
|
|
|
584 end return;
|
|
|
585 end Scalar_Vector_Elementwise_Operation;
|
|
|
586
|
|
|
587 ----------
|
|
|
588 -- Sqrt --
|
|
|
589 ----------
|
|
|
590
|
|
|
591 function Sqrt (X : Real'Base) return Real'Base is
|
|
|
592 Root, Next : Real'Base;
|
|
|
593
|
|
|
594 begin
|
|
|
595 -- Be defensive: any comparisons with NaN values will yield False.
|
|
|
596
|
|
|
597 if not (X > 0.0) then
|
|
|
598 if X = 0.0 then
|
|
|
599 return X;
|
|
|
600 else
|
|
|
601 raise Argument_Error;
|
|
|
602 end if;
|
|
|
603
|
|
|
604 elsif X > Real'Base'Last then
|
|
|
605
|
|
|
606 -- X is infinity, which is its own square root
|
|
|
607
|
|
|
608 return X;
|
|
|
609 end if;
|
|
|
610
|
|
|
611 -- Compute an initial estimate based on:
|
|
|
612
|
|
|
613 -- X = M * R**E and Sqrt (X) = Sqrt (M) * R**(E / 2.0),
|
|
|
614
|
|
|
615 -- where M is the mantissa, R is the radix and E the exponent.
|
|
|
616
|
|
|
617 -- By ignoring the mantissa and ignoring the case of an odd
|
|
|
618 -- exponent, we get a final error that is at most R. In other words,
|
|
|
619 -- the result has about a single bit precision.
|
|
|
620
|
|
|
621 Root := Real'Base (Real'Machine_Radix) ** (Real'Exponent (X) / 2);
|
|
|
622
|
|
|
623 -- Because of the poor initial estimate, use the Babylonian method of
|
|
|
624 -- computing the square root, as it is stable for all inputs. Every step
|
|
|
625 -- will roughly double the precision of the result. Just a few steps
|
|
|
626 -- suffice in most cases. Eight iterations should give about 2**8 bits
|
|
|
627 -- of precision.
|
|
|
628
|
|
|
629 for J in 1 .. 8 loop
|
|
|
630 Next := (Root + X / Root) / 2.0;
|
|
|
631 exit when Root = Next;
|
|
|
632 Root := Next;
|
|
|
633 end loop;
|
|
|
634
|
|
|
635 return Root;
|
|
|
636 end Sqrt;
|
|
|
637
|
|
|
638 ---------------------------
|
|
|
639 -- Matrix_Matrix_Product --
|
|
|
640 ---------------------------
|
|
|
641
|
|
|
642 function Matrix_Matrix_Product
|
|
|
643 (Left : Left_Matrix;
|
|
|
644 Right : Right_Matrix) return Result_Matrix
|
|
|
645 is
|
|
|
646 begin
|
|
|
647 return R : Result_Matrix (Left'Range (1), Right'Range (2)) do
|
|
|
648 if Left'Length (2) /= Right'Length (1) then
|
|
|
649 raise Constraint_Error with
|
|
|
650 "incompatible dimensions in matrix multiplication";
|
|
|
651 end if;
|
|
|
652
|
|
|
653 for J in R'Range (1) loop
|
|
|
654 for K in R'Range (2) loop
|
|
|
655 declare
|
|
|
656 S : Result_Scalar := Zero;
|
|
|
657
|
|
|
658 begin
|
|
|
659 for M in Left'Range (2) loop
|
|
|
660 S := S + Left (J, M) *
|
|
|
661 Right
|
|
|
662 (M - Left'First (2) + Right'First (1), K);
|
|
|
663 end loop;
|
|
|
664
|
|
|
665 R (J, K) := S;
|
|
|
666 end;
|
|
|
667 end loop;
|
|
|
668 end loop;
|
|
|
669 end return;
|
|
|
670 end Matrix_Matrix_Product;
|
|
|
671
|
|
|
672 ----------------------------
|
|
|
673 -- Matrix_Vector_Solution --
|
|
|
674 ----------------------------
|
|
|
675
|
|
|
676 function Matrix_Vector_Solution (A : Matrix; X : Vector) return Vector is
|
|
|
677 N : constant Natural := A'Length (1);
|
|
|
678 MA : Matrix := A;
|
|
|
679 MX : Matrix (A'Range (1), 1 .. 1);
|
|
|
680 R : Vector (A'Range (2));
|
|
|
681 Det : Scalar;
|
|
|
682
|
|
|
683 begin
|
|
|
684 if A'Length (2) /= N then
|
|
|
685 raise Constraint_Error with "matrix is not square";
|
|
|
686 end if;
|
|
|
687
|
|
|
688 if X'Length /= N then
|
|
|
689 raise Constraint_Error with "incompatible vector length";
|
|
|
690 end if;
|
|
|
691
|
|
|
692 for J in 0 .. MX'Length (1) - 1 loop
|
|
|
693 MX (MX'First (1) + J, 1) := X (X'First + J);
|
|
|
694 end loop;
|
|
|
695
|
|
|
696 Forward_Eliminate (MA, MX, Det);
|
|
|
697
|
|
|
698 if Det = Zero then
|
|
|
699 raise Constraint_Error with "matrix is singular";
|
|
|
700 end if;
|
|
|
701
|
|
|
702 Back_Substitute (MA, MX);
|
|
|
703
|
|
|
704 for J in 0 .. R'Length - 1 loop
|
|
|
705 R (R'First + J) := MX (MX'First (1) + J, 1);
|
|
|
706 end loop;
|
|
|
707
|
|
|
708 return R;
|
|
|
709 end Matrix_Vector_Solution;
|
|
|
710
|
|
|
711 ----------------------------
|
|
|
712 -- Matrix_Matrix_Solution --
|
|
|
713 ----------------------------
|
|
|
714
|
|
|
715 function Matrix_Matrix_Solution (A, X : Matrix) return Matrix is
|
|
|
716 N : constant Natural := A'Length (1);
|
|
|
717 MA : Matrix (A'Range (2), A'Range (2));
|
|
|
718 MB : Matrix (A'Range (2), X'Range (2));
|
|
|
719 Det : Scalar;
|
|
|
720
|
|
|
721 begin
|
|
|
722 if A'Length (2) /= N then
|
|
|
723 raise Constraint_Error with "matrix is not square";
|
|
|
724 end if;
|
|
|
725
|
|
|
726 if X'Length (1) /= N then
|
|
|
727 raise Constraint_Error with "matrices have unequal number of rows";
|
|
|
728 end if;
|
|
|
729
|
|
|
730 for J in 0 .. A'Length (1) - 1 loop
|
|
|
731 for K in MA'Range (2) loop
|
|
|
732 MA (MA'First (1) + J, K) := A (A'First (1) + J, K);
|
|
|
733 end loop;
|
|
|
734
|
|
|
735 for K in MB'Range (2) loop
|
|
|
736 MB (MB'First (1) + J, K) := X (X'First (1) + J, K);
|
|
|
737 end loop;
|
|
|
738 end loop;
|
|
|
739
|
|
|
740 Forward_Eliminate (MA, MB, Det);
|
|
|
741
|
|
|
742 if Det = Zero then
|
|
|
743 raise Constraint_Error with "matrix is singular";
|
|
|
744 end if;
|
|
|
745
|
|
|
746 Back_Substitute (MA, MB);
|
|
|
747
|
|
|
748 return MB;
|
|
|
749 end Matrix_Matrix_Solution;
|
|
|
750
|
|
|
751 ---------------------------
|
|
|
752 -- Matrix_Vector_Product --
|
|
|
753 ---------------------------
|
|
|
754
|
|
|
755 function Matrix_Vector_Product
|
|
|
756 (Left : Matrix;
|
|
|
757 Right : Right_Vector) return Result_Vector
|
|
|
758 is
|
|
|
759 begin
|
|
|
760 return R : Result_Vector (Left'Range (1)) do
|
|
|
761 if Left'Length (2) /= Right'Length then
|
|
|
762 raise Constraint_Error with
|
|
|
763 "incompatible dimensions in matrix-vector multiplication";
|
|
|
764 end if;
|
|
|
765
|
|
|
766 for J in Left'Range (1) loop
|
|
|
767 declare
|
|
|
768 S : Result_Scalar := Zero;
|
|
|
769
|
|
|
770 begin
|
|
|
771 for K in Left'Range (2) loop
|
|
|
772 S := S + Left (J, K)
|
|
|
773 * Right (K - Left'First (2) + Right'First);
|
|
|
774 end loop;
|
|
|
775
|
|
|
776 R (J) := S;
|
|
|
777 end;
|
|
|
778 end loop;
|
|
|
779 end return;
|
|
|
780 end Matrix_Vector_Product;
|
|
|
781
|
|
|
782 -------------------
|
|
|
783 -- Outer_Product --
|
|
|
784 -------------------
|
|
|
785
|
|
|
786 function Outer_Product
|
|
|
787 (Left : Left_Vector;
|
|
|
788 Right : Right_Vector) return Matrix
|
|
|
789 is
|
|
|
790 begin
|
|
|
791 return R : Matrix (Left'Range, Right'Range) do
|
|
|
792 for J in R'Range (1) loop
|
|
|
793 for K in R'Range (2) loop
|
|
|
794 R (J, K) := Left (J) * Right (K);
|
|
|
795 end loop;
|
|
|
796 end loop;
|
|
|
797 end return;
|
|
|
798 end Outer_Product;
|
|
|
799
|
|
|
800 -----------------
|
|
|
801 -- Swap_Column --
|
|
|
802 -----------------
|
|
|
803
|
|
|
804 procedure Swap_Column (A : in out Matrix; Left, Right : Integer) is
|
|
|
805 Temp : Scalar;
|
|
|
806 begin
|
|
|
807 for J in A'Range (1) loop
|
|
|
808 Temp := A (J, Left);
|
|
|
809 A (J, Left) := A (J, Right);
|
|
|
810 A (J, Right) := Temp;
|
|
|
811 end loop;
|
|
|
812 end Swap_Column;
|
|
|
813
|
|
|
814 ---------------
|
|
|
815 -- Transpose --
|
|
|
816 ---------------
|
|
|
817
|
|
|
818 procedure Transpose (A : Matrix; R : out Matrix) is
|
|
|
819 begin
|
|
|
820 for J in R'Range (1) loop
|
|
|
821 for K in R'Range (2) loop
|
|
|
822 R (J, K) := A (K - R'First (2) + A'First (1),
|
|
|
823 J - R'First (1) + A'First (2));
|
|
|
824 end loop;
|
|
|
825 end loop;
|
|
|
826 end Transpose;
|
|
|
827
|
|
|
828 -------------------------------
|
|
|
829 -- Update_Matrix_With_Matrix --
|
|
|
830 -------------------------------
|
|
|
831
|
|
|
832 procedure Update_Matrix_With_Matrix (X : in out X_Matrix; Y : Y_Matrix) is
|
|
|
833 begin
|
|
|
834 if X'Length (1) /= Y'Length (1)
|
|
|
835 or else
|
|
|
836 X'Length (2) /= Y'Length (2)
|
|
|
837 then
|
|
|
838 raise Constraint_Error with
|
|
|
839 "matrices are of different dimension in update operation";
|
|
|
840 end if;
|
|
|
841
|
|
|
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)));
|
|
|
846 end loop;
|
|
|
847 end loop;
|
|
|
848 end Update_Matrix_With_Matrix;
|
|
|
849
|
|
|
850 -------------------------------
|
|
|
851 -- Update_Vector_With_Vector --
|
|
|
852 -------------------------------
|
|
|
853
|
|
|
854 procedure Update_Vector_With_Vector (X : in out X_Vector; Y : Y_Vector) is
|
|
|
855 begin
|
|
|
856 if X'Length /= Y'Length then
|
|
|
857 raise Constraint_Error with
|
|
|
858 "vectors are of different length in update operation";
|
|
|
859 end if;
|
|
|
860
|
|
|
861 for J in X'Range loop
|
|
|
862 Update (X (J), Y (J - X'First + Y'First));
|
|
|
863 end loop;
|
|
|
864 end Update_Vector_With_Vector;
|
|
|
865
|
|
|
866 -----------------
|
|
|
867 -- Unit_Matrix --
|
|
|
868 -----------------
|
|
|
869
|
|
|
870 function Unit_Matrix
|
|
|
871 (Order : Positive;
|
|
|
872 First_1 : Integer := 1;
|
|
|
873 First_2 : Integer := 1) return Matrix
|
|
|
874 is
|
|
|
875 begin
|
|
|
876 return R : Matrix (First_1 .. Check_Unit_Last (First_1, Order, First_1),
|
|
|
877 First_2 .. Check_Unit_Last (First_2, Order, First_2))
|
|
|
878 do
|
|
|
879 R := (others => (others => Zero));
|
|
|
880
|
|
|
881 for J in 0 .. Order - 1 loop
|
|
|
882 R (First_1 + J, First_2 + J) := One;
|
|
|
883 end loop;
|
|
|
884 end return;
|
|
|
885 end Unit_Matrix;
|
|
|
886
|
|
|
887 -----------------
|
|
|
888 -- Unit_Vector --
|
|
|
889 -----------------
|
|
|
890
|
|
|
891 function Unit_Vector
|
|
|
892 (Index : Integer;
|
|
|
893 Order : Positive;
|
|
|
894 First : Integer := 1) return Vector
|
|
|
895 is
|
|
|
896 begin
|
|
|
897 return R : Vector (First .. Check_Unit_Last (Index, Order, First)) do
|
|
|
898 R := (others => Zero);
|
|
|
899 R (Index) := One;
|
|
|
900 end return;
|
|
|
901 end Unit_Vector;
|
|
|
902
|
|
|
903 ---------------------------
|
|
|
904 -- Vector_Matrix_Product --
|
|
|
905 ---------------------------
|
|
|
906
|
|
|
907 function Vector_Matrix_Product
|
|
|
908 (Left : Left_Vector;
|
|
|
909 Right : Matrix) return Result_Vector
|
|
|
910 is
|
|
|
911 begin
|
|
|
912 return R : Result_Vector (Right'Range (2)) do
|
|
|
913 if Left'Length /= Right'Length (1) then
|
|
|
914 raise Constraint_Error with
|
|
|
915 "incompatible dimensions in vector-matrix multiplication";
|
|
|
916 end if;
|
|
|
917
|
|
|
918 for J in Right'Range (2) loop
|
|
|
919 declare
|
|
|
920 S : Result_Scalar := Zero;
|
|
|
921
|
|
|
922 begin
|
|
|
923 for K in Right'Range (1) loop
|
|
|
924 S := S + Left (K - Right'First (1)
|
|
|
925 + Left'First) * Right (K, J);
|
|
|
926 end loop;
|
|
|
927
|
|
|
928 R (J) := S;
|
|
|
929 end;
|
|
|
930 end loop;
|
|
|
931 end return;
|
|
|
932 end Vector_Matrix_Product;
|
|
|
933
|
|
|
934 end System.Generic_Array_Operations;
|
