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-2018, 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 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 934end System.Generic_Array_Operations; 935