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-2012, 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; 33 34package body System.Generic_Array_Operations is 35 36 function Check_Unit_Last 37 (Index : Integer; 38 Order : Positive; 39 First : Integer) return Integer; 40 pragma Inline_Always (Check_Unit_Last); 41 -- Compute index of last element returned by Unit_Vector or Unit_Matrix. 42 -- A separate function is needed to allow raising Constraint_Error before 43 -- declaring the function result variable. The result variable needs to be 44 -- declared first, to allow front-end inlining. 45 46 -------------- 47 -- Diagonal -- 48 -------------- 49 50 function Diagonal (A : Matrix) return Vector is 51 N : constant Natural := Natural'Min (A'Length (1), A'Length (2)); 52 begin 53 return R : Vector (A'First (1) .. A'First (1) + N - 1) do 54 for J in 0 .. N - 1 loop 55 R (R'First + J) := A (A'First (1) + J, A'First (2) + J); 56 end loop; 57 end return; 58 end Diagonal; 59 60 -------------------------- 61 -- Square_Matrix_Length -- 62 -------------------------- 63 64 function Square_Matrix_Length (A : Matrix) return Natural is 65 begin 66 if A'Length (1) /= A'Length (2) then 67 raise Constraint_Error with "matrix is not square"; 68 else 69 return A'Length (1); 70 end if; 71 end Square_Matrix_Length; 72 73 --------------------- 74 -- Check_Unit_Last -- 75 --------------------- 76 77 function Check_Unit_Last 78 (Index : Integer; 79 Order : Positive; 80 First : Integer) return Integer 81 is 82 begin 83 -- Order the tests carefully to avoid overflow 84 85 if Index < First 86 or else First > Integer'Last - Order + 1 87 or else Index > First + (Order - 1) 88 then 89 raise Constraint_Error; 90 end if; 91 92 return First + (Order - 1); 93 end Check_Unit_Last; 94 95 --------------------- 96 -- Back_Substitute -- 97 --------------------- 98 99 procedure Back_Substitute (M, N : in out Matrix) is 100 pragma Assert (M'First (1) = N'First (1) 101 and then 102 M'Last (1) = N'Last (1)); 103 104 procedure Sub_Row 105 (M : in out Matrix; 106 Target : Integer; 107 Source : Integer; 108 Factor : Scalar); 109 -- Elementary row operation that subtracts Factor * M (Source, <>) from 110 -- M (Target, <>) 111 112 ------------- 113 -- Sub_Row -- 114 ------------- 115 116 procedure Sub_Row 117 (M : in out Matrix; 118 Target : Integer; 119 Source : Integer; 120 Factor : Scalar) 121 is 122 begin 123 for J in M'Range (2) loop 124 M (Target, J) := M (Target, J) - Factor * M (Source, J); 125 end loop; 126 end Sub_Row; 127 128 -- Local declarations 129 130 Max_Col : Integer := M'Last (2); 131 132 -- Start of processing for Back_Substitute 133 134 begin 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 138 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. 142 143 declare 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.. 148 149 J : Integer := M'First (1); 150 151 begin 152 while J < Row loop 153 Sub_Row (N, J, Row, (M (J, Col) / M (Row, Col))); 154 Sub_Row (M, J, Row, (M (J, Col) / M (Row, Col))); 155 J := J + 1; 156 end loop; 157 end; 158 159 -- Avoid potential overflow in the subtraction below 160 161 exit Do_Rows when Col = M'First (2); 162 163 Max_Col := Col - 1; 164 165 exit Find_Non_Zero; 166 end if; 167 end loop Find_Non_Zero; 168 end loop Do_Rows; 169 end Back_Substitute; 170 171 ----------------------- 172 -- Forward_Eliminate -- 173 ----------------------- 174 175 procedure Forward_Eliminate 176 (M : in out Matrix; 177 N : in out Matrix; 178 Det : out Scalar) 179 is 180 pragma Assert (M'First (1) = N'First (1) 181 and then 182 M'Last (1) = N'Last (1)); 183 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. 189 190 procedure Sub_Row 191 (M : in out Matrix; 192 Target : Integer; 193 Source : Integer; 194 Factor : Scalar); 195 -- Subtrace Factor * M (Source, <>) from M (Target, <>) 196 197 procedure Divide_Row 198 (M, N : in out Matrix; 199 Row : Integer; 200 Scale : Scalar); 201 -- Divide M (Row) and N (Row) by Scale, and update Det 202 203 procedure Switch_Row 204 (M, N : in out Matrix; 205 Row_1 : Integer; 206 Row_2 : Integer); 207 -- Exchange M (Row_1) and N (Row_1) with M (Row_2) and N (Row_2), 208 -- negating Det in the process. 209 210 ------------- 211 -- Sub_Row -- 212 ------------- 213 214 procedure Sub_Row 215 (M : in out Matrix; 216 Target : Integer; 217 Source : Integer; 218 Factor : Scalar) 219 is 220 begin 221 for J in M'Range (2) loop 222 M (Target, J) := M (Target, J) - Factor * M (Source, J); 223 end loop; 224 end Sub_Row; 225 226 ---------------- 227 -- Divide_Row -- 228 ---------------- 229 230 procedure Divide_Row 231 (M, N : in out Matrix; 232 Row : Integer; 233 Scale : Scalar) 234 is 235 begin 236 Det := Det * Scale; 237 238 for J in M'Range (2) loop 239 M (Row, J) := M (Row, J) / Scale; 240 end loop; 241 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; 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 608 -- X is infinity, which is its own square root 609 610 return X; 611 end if; 612 613 -- Compute an initial estimate based on: 614 615 -- X = M * R**E and Sqrt (X) = Sqrt (M) * R**(E / 2.0), 616 617 -- where M is the mantissa, R is the radix and E the exponent. 618 619 -- By ignoring the mantissa and ignoring the case of an odd 620 -- exponent, we get a final error that is at most R. In other words, 621 -- the result has about a single bit precision. 622 623 Root := Real'Base (Real'Machine_Radix) ** (Real'Exponent (X) / 2); 624 625 -- Because of the poor initial estimate, use the Babylonian method of 626 -- computing the square root, as it is stable for all inputs. Every step 627 -- will roughly double the precision of the result. Just a few steps 628 -- suffice in most cases. Eight iterations should give about 2**8 bits 629 -- of precision. 630 631 for J in 1 .. 8 loop 632 Next := (Root + X / Root) / 2.0; 633 exit when Root = Next; 634 Root := Next; 635 end loop; 636 637 return Root; 638 end Sqrt; 639 640 --------------------------- 641 -- Matrix_Matrix_Product -- 642 --------------------------- 643 644 function Matrix_Matrix_Product 645 (Left : Left_Matrix; 646 Right : Right_Matrix) return Result_Matrix 647 is 648 begin 649 return R : Result_Matrix (Left'Range (1), Right'Range (2)) do 650 if Left'Length (2) /= Right'Length (1) then 651 raise Constraint_Error with 652 "incompatible dimensions in matrix multiplication"; 653 end if; 654 655 for J in R'Range (1) loop 656 for K in R'Range (2) loop 657 declare 658 S : Result_Scalar := Zero; 659 660 begin 661 for M in Left'Range (2) loop 662 S := S + Left (J, M) * 663 Right 664 (M - Left'First (2) + Right'First (1), K); 665 end loop; 666 667 R (J, K) := S; 668 end; 669 end loop; 670 end loop; 671 end return; 672 end Matrix_Matrix_Product; 673 674 ---------------------------- 675 -- Matrix_Vector_Solution -- 676 ---------------------------- 677 678 function Matrix_Vector_Solution (A : Matrix; X : Vector) return Vector is 679 N : constant Natural := A'Length (1); 680 MA : Matrix := A; 681 MX : Matrix (A'Range (1), 1 .. 1); 682 R : Vector (A'Range (2)); 683 Det : Scalar; 684 685 begin 686 if A'Length (2) /= N then 687 raise Constraint_Error with "matrix is not square"; 688 end if; 689 690 if X'Length /= N then 691 raise Constraint_Error with "incompatible vector length"; 692 end if; 693 694 for J in 0 .. MX'Length (1) - 1 loop 695 MX (MX'First (1) + J, 1) := X (X'First + J); 696 end loop; 697 698 Forward_Eliminate (MA, MX, Det); 699 Back_Substitute (MA, MX); 700 701 for J in 0 .. R'Length - 1 loop 702 R (R'First + J) := MX (MX'First (1) + J, 1); 703 end loop; 704 705 return R; 706 end Matrix_Vector_Solution; 707 708 ---------------------------- 709 -- Matrix_Matrix_Solution -- 710 ---------------------------- 711 712 function Matrix_Matrix_Solution (A, X : Matrix) return Matrix is 713 N : constant Natural := A'Length (1); 714 MA : Matrix (A'Range (2), A'Range (2)); 715 MB : Matrix (A'Range (2), X'Range (2)); 716 Det : Scalar; 717 718 begin 719 if A'Length (2) /= N then 720 raise Constraint_Error with "matrix is not square"; 721 end if; 722 723 if X'Length (1) /= N then 724 raise Constraint_Error with "matrices have unequal number of rows"; 725 end if; 726 727 for J in 0 .. A'Length (1) - 1 loop 728 for K in MA'Range (2) loop 729 MA (MA'First (1) + J, K) := A (A'First (1) + J, K); 730 end loop; 731 732 for K in MB'Range (2) loop 733 MB (MB'First (1) + J, K) := X (X'First (1) + J, K); 734 end loop; 735 end loop; 736 737 Forward_Eliminate (MA, MB, Det); 738 Back_Substitute (MA, MB); 739 740 return MB; 741 end Matrix_Matrix_Solution; 742 743 --------------------------- 744 -- Matrix_Vector_Product -- 745 --------------------------- 746 747 function Matrix_Vector_Product 748 (Left : Matrix; 749 Right : Right_Vector) return Result_Vector 750 is 751 begin 752 return R : Result_Vector (Left'Range (1)) do 753 if Left'Length (2) /= Right'Length then 754 raise Constraint_Error with 755 "incompatible dimensions in matrix-vector multiplication"; 756 end if; 757 758 for J in Left'Range (1) loop 759 declare 760 S : Result_Scalar := Zero; 761 762 begin 763 for K in Left'Range (2) loop 764 S := S + Left (J, K) 765 * Right (K - Left'First (2) + Right'First); 766 end loop; 767 768 R (J) := S; 769 end; 770 end loop; 771 end return; 772 end Matrix_Vector_Product; 773 774 ------------------- 775 -- Outer_Product -- 776 ------------------- 777 778 function Outer_Product 779 (Left : Left_Vector; 780 Right : Right_Vector) return Matrix 781 is 782 begin 783 return R : Matrix (Left'Range, Right'Range) do 784 for J in R'Range (1) loop 785 for K in R'Range (2) loop 786 R (J, K) := Left (J) * Right (K); 787 end loop; 788 end loop; 789 end return; 790 end Outer_Product; 791 792 ----------------- 793 -- Swap_Column -- 794 ----------------- 795 796 procedure Swap_Column (A : in out Matrix; Left, Right : Integer) is 797 Temp : Scalar; 798 begin 799 for J in A'Range (1) loop 800 Temp := A (J, Left); 801 A (J, Left) := A (J, Right); 802 A (J, Right) := Temp; 803 end loop; 804 end Swap_Column; 805 806 --------------- 807 -- Transpose -- 808 --------------- 809 810 procedure Transpose (A : Matrix; R : out Matrix) is 811 begin 812 for J in R'Range (1) loop 813 for K in R'Range (2) loop 814 R (J, K) := A (K - R'First (2) + A'First (1), 815 J - R'First (1) + A'First (2)); 816 end loop; 817 end loop; 818 end Transpose; 819 820 ------------------------------- 821 -- Update_Matrix_With_Matrix -- 822 ------------------------------- 823 824 procedure Update_Matrix_With_Matrix (X : in out X_Matrix; Y : Y_Matrix) is 825 begin 826 if X'Length (1) /= Y'Length (1) 827 or else 828 X'Length (2) /= Y'Length (2) 829 then 830 raise Constraint_Error with 831 "matrices are of different dimension in update operation"; 832 end if; 833 834 for J in X'Range (1) loop 835 for K in X'Range (2) loop 836 Update (X (J, K), Y (J - X'First (1) + Y'First (1), 837 K - X'First (2) + Y'First (2))); 838 end loop; 839 end loop; 840 end Update_Matrix_With_Matrix; 841 842 ------------------------------- 843 -- Update_Vector_With_Vector -- 844 ------------------------------- 845 846 procedure Update_Vector_With_Vector (X : in out X_Vector; Y : Y_Vector) is 847 begin 848 if X'Length /= Y'Length then 849 raise Constraint_Error with 850 "vectors are of different length in update operation"; 851 end if; 852 853 for J in X'Range loop 854 Update (X (J), Y (J - X'First + Y'First)); 855 end loop; 856 end Update_Vector_With_Vector; 857 858 ----------------- 859 -- Unit_Matrix -- 860 ----------------- 861 862 function Unit_Matrix 863 (Order : Positive; 864 First_1 : Integer := 1; 865 First_2 : Integer := 1) return Matrix 866 is 867 begin 868 return R : Matrix (First_1 .. Check_Unit_Last (First_1, Order, First_1), 869 First_2 .. Check_Unit_Last (First_2, Order, First_2)) 870 do 871 R := (others => (others => Zero)); 872 873 for J in 0 .. Order - 1 loop 874 R (First_1 + J, First_2 + J) := One; 875 end loop; 876 end return; 877 end Unit_Matrix; 878 879 ----------------- 880 -- Unit_Vector -- 881 ----------------- 882 883 function Unit_Vector 884 (Index : Integer; 885 Order : Positive; 886 First : Integer := 1) return Vector 887 is 888 begin 889 return R : Vector (First .. Check_Unit_Last (Index, Order, First)) do 890 R := (others => Zero); 891 R (Index) := One; 892 end return; 893 end Unit_Vector; 894 895 --------------------------- 896 -- Vector_Matrix_Product -- 897 --------------------------- 898 899 function Vector_Matrix_Product 900 (Left : Left_Vector; 901 Right : Matrix) return Result_Vector 902 is 903 begin 904 return R : Result_Vector (Right'Range (2)) do 905 if Left'Length /= Right'Length (1) then 906 raise Constraint_Error with 907 "incompatible dimensions in vector-matrix multiplication"; 908 end if; 909 910 for J in Right'Range (2) loop 911 declare 912 S : Result_Scalar := Zero; 913 914 begin 915 for K in Right'Range (1) loop 916 S := S + Left (K - Right'First (1) 917 + Left'First) * Right (K, J); 918 end loop; 919 920 R (J) := S; 921 end; 922 end loop; 923 end return; 924 end Vector_Matrix_Product; 925 926end System.Generic_Array_Operations; 927