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-- -- 9-- Copyright (C) 2006-2021, Free Software Foundation, Inc. -- 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 32with Ada.Numerics; use Ada.Numerics; 33package 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 pragma Annotate 244 (CodePeer, False_Positive, "divide by zero", "Scale /= 0"); 245 end loop; 246 end Divide_Row; 247 248 ---------------- 249 -- Switch_Row -- 250 ---------------- 251 252 procedure Switch_Row 253 (M, N : in out Matrix; 254 Row_1 : Integer; 255 Row_2 : Integer) 256 is 257 procedure Swap (X, Y : in out Scalar); 258 -- Exchange the values of X and Y 259 260 ---------- 261 -- Swap -- 262 ---------- 263 264 procedure Swap (X, Y : in out Scalar) is 265 T : constant Scalar := X; 266 begin 267 X := Y; 268 Y := T; 269 end Swap; 270 271 -- Start of processing for Switch_Row 272 273 begin 274 if Row_1 /= Row_2 then 275 Det := Zero - Det; 276 277 for J in M'Range (2) loop 278 Swap (M (Row_1, J), M (Row_2, J)); 279 end loop; 280 281 for J in N'Range (2) loop 282 Swap (N (Row_1 - M'First (1) + N'First (1), J), 283 N (Row_2 - M'First (1) + N'First (1), J)); 284 end loop; 285 end if; 286 end Switch_Row; 287 288 -- Local declarations 289 290 Row : Integer := M'First (1); 291 292 -- Start of processing for Forward_Eliminate 293 294 begin 295 Det := One; 296 297 for J in M'Range (2) loop 298 declare 299 Max_Row : Integer := Row; 300 Max_Abs : Real'Base := 0.0; 301 302 begin 303 -- Find best pivot in column J, starting in row Row 304 305 for K in Row .. M'Last (1) loop 306 declare 307 New_Abs : constant Real'Base := abs M (K, J); 308 begin 309 if Max_Abs < New_Abs then 310 Max_Abs := New_Abs; 311 Max_Row := K; 312 end if; 313 end; 314 end loop; 315 316 if Max_Abs > 0.0 then 317 Switch_Row (M, N, Row, Max_Row); 318 319 -- The temporaries below are necessary to force a copy of the 320 -- value and avoid improper aliasing. 321 322 declare 323 Scale : constant Scalar := M (Row, J); 324 begin 325 Divide_Row (M, N, Row, Scale); 326 end; 327 328 for U in Row + 1 .. M'Last (1) loop 329 declare 330 Factor : constant Scalar := M (U, J); 331 begin 332 Sub_Row (N, U, Row, Factor); 333 Sub_Row (M, U, Row, Factor); 334 end; 335 end loop; 336 337 exit when Row >= M'Last (1); 338 339 Row := Row + 1; 340 341 else 342 -- Set zero (note that we do not have literals) 343 344 Det := Zero; 345 end if; 346 end; 347 end loop; 348 end Forward_Eliminate; 349 350 ------------------- 351 -- Inner_Product -- 352 ------------------- 353 354 function Inner_Product 355 (Left : Left_Vector; 356 Right : Right_Vector) return Result_Scalar 357 is 358 R : Result_Scalar := Zero; 359 360 begin 361 if Left'Length /= Right'Length then 362 raise Constraint_Error with 363 "vectors are of different length in inner product"; 364 end if; 365 366 for J in Left'Range loop 367 R := R + Left (J) * Right (J - Left'First + Right'First); 368 end loop; 369 370 return R; 371 end Inner_Product; 372 373 ------------- 374 -- L2_Norm -- 375 ------------- 376 377 function L2_Norm (X : X_Vector) return Result_Real'Base is 378 Sum : Result_Real'Base := 0.0; 379 380 begin 381 for J in X'Range loop 382 Sum := Sum + Result_Real'Base (abs X (J))**2; 383 end loop; 384 385 return Sqrt (Sum); 386 end L2_Norm; 387 388 ---------------------------------- 389 -- Matrix_Elementwise_Operation -- 390 ---------------------------------- 391 392 function Matrix_Elementwise_Operation (X : X_Matrix) return Result_Matrix is 393 begin 394 return R : Result_Matrix (X'Range (1), X'Range (2)) do 395 for J in R'Range (1) loop 396 for K in R'Range (2) loop 397 R (J, K) := Operation (X (J, K)); 398 end loop; 399 end loop; 400 end return; 401 end Matrix_Elementwise_Operation; 402 403 ---------------------------------- 404 -- Vector_Elementwise_Operation -- 405 ---------------------------------- 406 407 function Vector_Elementwise_Operation (X : X_Vector) return Result_Vector is 408 begin 409 return R : Result_Vector (X'Range) do 410 for J in R'Range loop 411 R (J) := Operation (X (J)); 412 end loop; 413 end return; 414 end Vector_Elementwise_Operation; 415 416 ----------------------------------------- 417 -- Matrix_Matrix_Elementwise_Operation -- 418 ----------------------------------------- 419 420 function Matrix_Matrix_Elementwise_Operation 421 (Left : Left_Matrix; 422 Right : Right_Matrix) return Result_Matrix 423 is 424 begin 425 return R : Result_Matrix (Left'Range (1), Left'Range (2)) do 426 if Left'Length (1) /= Right'Length (1) 427 or else 428 Left'Length (2) /= Right'Length (2) 429 then 430 raise Constraint_Error with 431 "matrices are of different dimension in elementwise operation"; 432 end if; 433 434 for J in R'Range (1) loop 435 for K in R'Range (2) loop 436 R (J, K) := 437 Operation 438 (Left (J, K), 439 Right 440 (J - R'First (1) + Right'First (1), 441 K - R'First (2) + Right'First (2))); 442 end loop; 443 end loop; 444 end return; 445 end Matrix_Matrix_Elementwise_Operation; 446 447 ------------------------------------------------ 448 -- Matrix_Matrix_Scalar_Elementwise_Operation -- 449 ------------------------------------------------ 450 451 function Matrix_Matrix_Scalar_Elementwise_Operation 452 (X : X_Matrix; 453 Y : Y_Matrix; 454 Z : Z_Scalar) return Result_Matrix 455 is 456 begin 457 return R : Result_Matrix (X'Range (1), X'Range (2)) do 458 if X'Length (1) /= Y'Length (1) 459 or else 460 X'Length (2) /= Y'Length (2) 461 then 462 raise Constraint_Error with 463 "matrices are of different dimension in elementwise operation"; 464 end if; 465 466 for J in R'Range (1) loop 467 for K in R'Range (2) loop 468 R (J, K) := 469 Operation 470 (X (J, K), 471 Y (J - R'First (1) + Y'First (1), 472 K - R'First (2) + Y'First (2)), 473 Z); 474 end loop; 475 end loop; 476 end return; 477 end Matrix_Matrix_Scalar_Elementwise_Operation; 478 479 ----------------------------------------- 480 -- Vector_Vector_Elementwise_Operation -- 481 ----------------------------------------- 482 483 function Vector_Vector_Elementwise_Operation 484 (Left : Left_Vector; 485 Right : Right_Vector) return Result_Vector 486 is 487 begin 488 return R : Result_Vector (Left'Range) do 489 if Left'Length /= Right'Length then 490 raise Constraint_Error with 491 "vectors are of different length in elementwise operation"; 492 end if; 493 494 for J in R'Range loop 495 R (J) := Operation (Left (J), Right (J - R'First + Right'First)); 496 end loop; 497 end return; 498 end Vector_Vector_Elementwise_Operation; 499 500 ------------------------------------------------ 501 -- Vector_Vector_Scalar_Elementwise_Operation -- 502 ------------------------------------------------ 503 504 function Vector_Vector_Scalar_Elementwise_Operation 505 (X : X_Vector; 506 Y : Y_Vector; 507 Z : Z_Scalar) return Result_Vector is 508 begin 509 return R : Result_Vector (X'Range) do 510 if X'Length /= Y'Length then 511 raise Constraint_Error with 512 "vectors are of different length in elementwise operation"; 513 end if; 514 515 for J in R'Range loop 516 R (J) := Operation (X (J), Y (J - X'First + Y'First), Z); 517 end loop; 518 end return; 519 end Vector_Vector_Scalar_Elementwise_Operation; 520 521 ----------------------------------------- 522 -- Matrix_Scalar_Elementwise_Operation -- 523 ----------------------------------------- 524 525 function Matrix_Scalar_Elementwise_Operation 526 (Left : Left_Matrix; 527 Right : Right_Scalar) return Result_Matrix 528 is 529 begin 530 return R : Result_Matrix (Left'Range (1), Left'Range (2)) do 531 for J in R'Range (1) loop 532 for K in R'Range (2) loop 533 R (J, K) := Operation (Left (J, K), Right); 534 end loop; 535 end loop; 536 end return; 537 end Matrix_Scalar_Elementwise_Operation; 538 539 ----------------------------------------- 540 -- Vector_Scalar_Elementwise_Operation -- 541 ----------------------------------------- 542 543 function Vector_Scalar_Elementwise_Operation 544 (Left : Left_Vector; 545 Right : Right_Scalar) return Result_Vector 546 is 547 begin 548 return R : Result_Vector (Left'Range) do 549 for J in R'Range loop 550 R (J) := Operation (Left (J), Right); 551 end loop; 552 end return; 553 end Vector_Scalar_Elementwise_Operation; 554 555 ----------------------------------------- 556 -- Scalar_Matrix_Elementwise_Operation -- 557 ----------------------------------------- 558 559 function Scalar_Matrix_Elementwise_Operation 560 (Left : Left_Scalar; 561 Right : Right_Matrix) return Result_Matrix 562 is 563 begin 564 return R : Result_Matrix (Right'Range (1), Right'Range (2)) do 565 for J in R'Range (1) loop 566 for K in R'Range (2) loop 567 R (J, K) := Operation (Left, Right (J, K)); 568 end loop; 569 end loop; 570 end return; 571 end Scalar_Matrix_Elementwise_Operation; 572 573 ----------------------------------------- 574 -- Scalar_Vector_Elementwise_Operation -- 575 ----------------------------------------- 576 577 function Scalar_Vector_Elementwise_Operation 578 (Left : Left_Scalar; 579 Right : Right_Vector) return Result_Vector 580 is 581 begin 582 return R : Result_Vector (Right'Range) do 583 for J in R'Range loop 584 R (J) := Operation (Left, Right (J)); 585 end loop; 586 end return; 587 end Scalar_Vector_Elementwise_Operation; 588 589 ---------- 590 -- Sqrt -- 591 ---------- 592 593 function Sqrt (X : Real'Base) return Real'Base is 594 Root, Next : Real'Base; 595 596 begin 597 -- Be defensive: any comparisons with NaN values will yield False. 598 599 if not (X > 0.0) then 600 if X = 0.0 then 601 return X; 602 else 603 raise Argument_Error; 604 end if; 605 606 elsif X > Real'Base'Last then 607 pragma Annotate 608 (CodePeer, Intentional, 609 "test always false", "test for infinity"); 610 611 -- X is infinity, which is its own square root 612 613 return X; 614 end if; 615 616 -- Compute an initial estimate based on: 617 618 -- X = M * R**E and Sqrt (X) = Sqrt (M) * R**(E / 2.0), 619 620 -- where M is the mantissa, R is the radix and E the exponent. 621 622 -- By ignoring the mantissa and ignoring the case of an odd 623 -- exponent, we get a final error that is at most R. In other words, 624 -- the result has about a single bit precision. 625 626 Root := Real'Base (Real'Machine_Radix) ** (Real'Exponent (X) / 2); 627 628 -- Because of the poor initial estimate, use the Babylonian method of 629 -- computing the square root, as it is stable for all inputs. Every step 630 -- will roughly double the precision of the result. Just a few steps 631 -- suffice in most cases. Eight iterations should give about 2**8 bits 632 -- of precision. 633 634 for J in 1 .. 8 loop 635 pragma Assert (Root /= 0.0); 636 637 Next := (Root + X / Root) / 2.0; 638 exit when Root = Next; 639 Root := Next; 640 end loop; 641 642 return Root; 643 end Sqrt; 644 645 --------------------------- 646 -- Matrix_Matrix_Product -- 647 --------------------------- 648 649 function Matrix_Matrix_Product 650 (Left : Left_Matrix; 651 Right : Right_Matrix) return Result_Matrix 652 is 653 begin 654 return R : Result_Matrix (Left'Range (1), Right'Range (2)) do 655 if Left'Length (2) /= Right'Length (1) then 656 raise Constraint_Error with 657 "incompatible dimensions in matrix multiplication"; 658 end if; 659 660 for J in R'Range (1) loop 661 for K in R'Range (2) loop 662 declare 663 S : Result_Scalar := Zero; 664 665 begin 666 for M in Left'Range (2) loop 667 S := S + Left (J, M) * 668 Right 669 (M - Left'First (2) + Right'First (1), K); 670 end loop; 671 672 R (J, K) := S; 673 end; 674 end loop; 675 end loop; 676 end return; 677 end Matrix_Matrix_Product; 678 679 ---------------------------- 680 -- Matrix_Vector_Solution -- 681 ---------------------------- 682 683 function Matrix_Vector_Solution (A : Matrix; X : Vector) return Vector is 684 N : constant Natural := A'Length (1); 685 MA : Matrix := A; 686 MX : Matrix (A'Range (1), 1 .. 1); 687 R : Vector (A'Range (2)); 688 Det : Scalar; 689 690 begin 691 if A'Length (2) /= N then 692 raise Constraint_Error with "matrix is not square"; 693 end if; 694 695 if X'Length /= N then 696 raise Constraint_Error with "incompatible vector length"; 697 end if; 698 699 for J in 0 .. MX'Length (1) - 1 loop 700 MX (MX'First (1) + J, 1) := X (X'First + J); 701 end loop; 702 703 Forward_Eliminate (MA, MX, Det); 704 705 if Det = Zero then 706 raise Constraint_Error with "matrix is singular"; 707 end if; 708 709 Back_Substitute (MA, MX); 710 711 for J in 0 .. R'Length - 1 loop 712 R (R'First + J) := MX (MX'First (1) + J, 1); 713 end loop; 714 715 return R; 716 end Matrix_Vector_Solution; 717 718 ---------------------------- 719 -- Matrix_Matrix_Solution -- 720 ---------------------------- 721 722 function Matrix_Matrix_Solution (A, X : Matrix) return Matrix is 723 N : constant Natural := A'Length (1); 724 MA : Matrix (A'Range (2), A'Range (2)); 725 MB : Matrix (A'Range (2), X'Range (2)); 726 Det : Scalar; 727 728 begin 729 if A'Length (2) /= N then 730 raise Constraint_Error with "matrix is not square"; 731 end if; 732 733 if X'Length (1) /= N then 734 raise Constraint_Error with "matrices have unequal number of rows"; 735 end if; 736 737 for J in 0 .. A'Length (1) - 1 loop 738 for K in MA'Range (2) loop 739 MA (MA'First (1) + J, K) := A (A'First (1) + J, K); 740 end loop; 741 742 for K in MB'Range (2) loop 743 MB (MB'First (1) + J, K) := X (X'First (1) + J, K); 744 end loop; 745 end loop; 746 747 Forward_Eliminate (MA, MB, Det); 748 749 if Det = Zero then 750 raise Constraint_Error with "matrix is singular"; 751 end if; 752 753 Back_Substitute (MA, MB); 754 755 return MB; 756 end Matrix_Matrix_Solution; 757 758 --------------------------- 759 -- Matrix_Vector_Product -- 760 --------------------------- 761 762 function Matrix_Vector_Product 763 (Left : Matrix; 764 Right : Right_Vector) return Result_Vector 765 is 766 begin 767 return R : Result_Vector (Left'Range (1)) do 768 if Left'Length (2) /= Right'Length then 769 raise Constraint_Error with 770 "incompatible dimensions in matrix-vector multiplication"; 771 end if; 772 773 for J in Left'Range (1) loop 774 declare 775 S : Result_Scalar := Zero; 776 777 begin 778 for K in Left'Range (2) loop 779 S := S + Left (J, K) 780 * Right (K - Left'First (2) + Right'First); 781 end loop; 782 783 R (J) := S; 784 end; 785 end loop; 786 end return; 787 end Matrix_Vector_Product; 788 789 ------------------- 790 -- Outer_Product -- 791 ------------------- 792 793 function Outer_Product 794 (Left : Left_Vector; 795 Right : Right_Vector) return Matrix 796 is 797 begin 798 return R : Matrix (Left'Range, Right'Range) do 799 for J in R'Range (1) loop 800 for K in R'Range (2) loop 801 R (J, K) := Left (J) * Right (K); 802 end loop; 803 end loop; 804 end return; 805 end Outer_Product; 806 807 ----------------- 808 -- Swap_Column -- 809 ----------------- 810 811 procedure Swap_Column (A : in out Matrix; Left, Right : Integer) is 812 Temp : Scalar; 813 begin 814 for J in A'Range (1) loop 815 Temp := A (J, Left); 816 A (J, Left) := A (J, Right); 817 A (J, Right) := Temp; 818 end loop; 819 end Swap_Column; 820 821 --------------- 822 -- Transpose -- 823 --------------- 824 825 procedure Transpose (A : Matrix; R : out Matrix) is 826 begin 827 for J in R'Range (1) loop 828 for K in R'Range (2) loop 829 R (J, K) := A (K - R'First (2) + A'First (1), 830 J - R'First (1) + A'First (2)); 831 end loop; 832 end loop; 833 end Transpose; 834 835 ------------------------------- 836 -- Update_Matrix_With_Matrix -- 837 ------------------------------- 838 839 procedure Update_Matrix_With_Matrix (X : in out X_Matrix; Y : Y_Matrix) is 840 begin 841 if X'Length (1) /= Y'Length (1) 842 or else 843 X'Length (2) /= Y'Length (2) 844 then 845 raise Constraint_Error with 846 "matrices are of different dimension in update operation"; 847 end if; 848 849 for J in X'Range (1) loop 850 for K in X'Range (2) loop 851 Update (X (J, K), Y (J - X'First (1) + Y'First (1), 852 K - X'First (2) + Y'First (2))); 853 end loop; 854 end loop; 855 end Update_Matrix_With_Matrix; 856 857 ------------------------------- 858 -- Update_Vector_With_Vector -- 859 ------------------------------- 860 861 procedure Update_Vector_With_Vector (X : in out X_Vector; Y : Y_Vector) is 862 begin 863 if X'Length /= Y'Length then 864 raise Constraint_Error with 865 "vectors are of different length in update operation"; 866 end if; 867 868 for J in X'Range loop 869 Update (X (J), Y (J - X'First + Y'First)); 870 end loop; 871 end Update_Vector_With_Vector; 872 873 ----------------- 874 -- Unit_Matrix -- 875 ----------------- 876 877 function Unit_Matrix 878 (Order : Positive; 879 First_1 : Integer := 1; 880 First_2 : Integer := 1) return Matrix 881 is 882 begin 883 return R : Matrix (First_1 .. Check_Unit_Last (First_1, Order, First_1), 884 First_2 .. Check_Unit_Last (First_2, Order, First_2)) 885 do 886 R := [others => [others => Zero]]; 887 888 for J in 0 .. Order - 1 loop 889 R (First_1 + J, First_2 + J) := One; 890 end loop; 891 end return; 892 end Unit_Matrix; 893 894 ----------------- 895 -- Unit_Vector -- 896 ----------------- 897 898 function Unit_Vector 899 (Index : Integer; 900 Order : Positive; 901 First : Integer := 1) return Vector 902 is 903 begin 904 return R : Vector (First .. Check_Unit_Last (Index, Order, First)) do 905 R := [others => Zero]; 906 R (Index) := One; 907 end return; 908 end Unit_Vector; 909 910 --------------------------- 911 -- Vector_Matrix_Product -- 912 --------------------------- 913 914 function Vector_Matrix_Product 915 (Left : Left_Vector; 916 Right : Matrix) return Result_Vector 917 is 918 begin 919 return R : Result_Vector (Right'Range (2)) do 920 if Left'Length /= Right'Length (1) then 921 raise Constraint_Error with 922 "incompatible dimensions in vector-matrix multiplication"; 923 end if; 924 925 for J in Right'Range (2) loop 926 declare 927 S : Result_Scalar := Zero; 928 929 begin 930 for K in Right'Range (1) loop 931 S := S + Left (K - Right'First (1) 932 + Left'First) * Right (K, J); 933 end loop; 934 935 R (J) := S; 936 end; 937 end loop; 938 end return; 939 end Vector_Matrix_Product; 940 941end System.Generic_Array_Operations; 942