gcc/ada/s-gearop.adb

   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-2011, 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
  32 with Ada.Numerics; use Ada.Numerics;
  33
  34 package body System.Generic_Array_Operations is
  35
  36    --  The local function Check_Unit_Last computes the index of the last
  37    --  element returned by Unit_Vector or Unit_Matrix. A separate function is
  38    --  needed to allow raising Constraint_Error before declaring the function
  39    --  result variable. The result variable needs to be declared first, to
  40    --  allow front-end inlining.
  41
  42    function Check_Unit_Last
  43      (Index : Integer;
  44       Order : Positive;
  45       First : Integer) return Integer;
  46    pragma Inline_Always (Check_Unit_Last);
  47
  48    --------------
  49    -- Diagonal --
  50    --------------
  51
  52    function Diagonal (A : Matrix) return Vector is
  53       N : constant Natural := Natural'Min (A'Length (1), A'Length (2));
  54       R : Vector (A'First (1) .. A'First (1) + N - 1);
  55
  56    begin
  57       for J in 0 .. N - 1 loop
  58          R (R'First + J) := A (A'First (1) + J, A'First (2) + J);
  59       end loop;
  60
  61       return R;
  62    end Diagonal;
  63
  64    --------------------------
  65    -- Square_Matrix_Length --
  66    --------------------------
  67
  68    function Square_Matrix_Length (A : Matrix) return Natural is
  69    begin
  70       if A'Length (1) /= A'Length (2) then
  71          raise Constraint_Error with "matrix is not square";
  72       end if;
  73
  74       return A'Length (1);
  75    end Square_Matrix_Length;
  76
  77    ---------------------
  78    -- Check_Unit_Last --
  79    ---------------------
  80
  81    function Check_Unit_Last
  82       (Index : Integer;
  83        Order : Positive;
  84        First : Integer) return Integer
  85    is
  86    begin
  87       --  Order the tests carefully to avoid overflow
  88
  89       if Index < First
  90         or else First > Integer'Last - Order + 1
  91         or else Index > First + (Order - 1)
  92       then
  93          raise Constraint_Error;
  94       end if;
  95
  96       return First + (Order - 1);
  97    end Check_Unit_Last;
  98
  99    ---------------------
 100    -- Back_Substitute --
 101    ---------------------
 102
 103    procedure Back_Substitute (M, N : in out Matrix) is
 104       pragma Assert (M'First (1) = N'First (1)
 105                        and then
 106                      M'Last  (1) = N'Last (1));
 107
 108       procedure Sub_Row
 109         (M      : in out Matrix;
 110          Target : Integer;
 111          Source : Integer;
 112          Factor : Scalar);
 113       --  Elementary row operation that subtracts Factor * M (Source, <>) from
 114       --  M (Target, <>)
 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          procedure Swap (X, Y : in out Scalar) is
 261             T : constant Scalar := X;
 262          begin
 263             X := Y;
 264             Y := T;
 265          end Swap;
 266
 267       --  Start of processing for Switch_Row
 268
 269       begin
 270          if Row_1 /= Row_2 then
 271             Det := Zero - Det;
 272
 273             for J in M'Range (2) loop
 274                Swap (M (Row_1, J), M (Row_2, J));
 275             end loop;
 276
 277             for J in N'Range (2) loop
 278                Swap (N (Row_1 - M'First (1) + N'First (1), J),
 279                      N (Row_2 - M'First (1) + N'First (1), J));
 280             end loop;
 281          end if;
 282       end Switch_Row;
 283
 284       --  Local declarations
 285
 286       Row : Integer := M'First (1);
 287
 288    --  Start of processing for Forward_Eliminate
 289
 290    begin
 291       Det := One;
 292
 293       for J in M'Range (2) loop
 294          declare
 295             Max_Row : Integer := Row;
 296             Max_Abs : Real'Base := 0.0;
 297
 298          begin
 299             --  Find best pivot in column J, starting in row Row
 300
 301             for K in Row .. M'Last (1) loop
 302                declare
 303                   New_Abs : constant Real'Base := abs M (K, J);
 304                begin
 305                   if Max_Abs < New_Abs then
 306                      Max_Abs := New_Abs;
 307                      Max_Row := K;
 308                   end if;
 309                end;
 310             end loop;
 311
 312             if Max_Abs > 0.0 then
 313                Switch_Row (M, N, Row, Max_Row);
 314                Divide_Row (M, N, Row, M (Row, J));
 315
 316                for U in Row + 1 .. M'Last (1) loop
 317                   Sub_Row (N, U, Row, M (U, J));
 318                   Sub_Row (M, U, Row, M (U, J));
 319                end loop;
 320
 321                exit when Row >= M'Last (1);
 322
 323                Row := Row + 1;
 324
 325             else
 326                --  Set zero (note that we do not have literals)
 327
 328                Det := Zero;
 329             end if;
 330          end;
 331       end loop;
 332    end Forward_Eliminate;
 333
 334    -------------------
 335    -- Inner_Product --
 336    -------------------
 337
 338    function Inner_Product
 339      (Left  : Left_Vector;
 340       Right : Right_Vector) return  Result_Scalar
 341    is
 342       R : Result_Scalar := Zero;
 343
 344    begin
 345       if Left'Length /= Right'Length then
 346          raise Constraint_Error with
 347             "vectors are of different length in inner product";
 348       end if;
 349
 350       for J in Left'Range loop
 351          R := R + Left (J) * Right (J - Left'First + Right'First);
 352       end loop;
 353
 354       return R;
 355    end Inner_Product;
 356
 357    -------------
 358    -- L2_Norm --
 359    -------------
 360
 361    function L2_Norm (X : X_Vector) return Result_Real'Base is
 362       Sum : Result_Real'Base := 0.0;
 363
 364    begin
 365       for J in X'Range loop
 366          Sum := Sum + Result_Real'Base (abs X (J))**2;
 367       end loop;
 368
 369       return Sqrt (Sum);
 370    end L2_Norm;
 371
 372    ----------------------------------
 373    -- Matrix_Elementwise_Operation --
 374    ----------------------------------
 375
 376    function Matrix_Elementwise_Operation (X : X_Matrix) return Result_Matrix is
 377       R : Result_Matrix (X'Range (1), X'Range (2));
 378
 379    begin
 380       for J in R'Range (1) loop
 381          for K in R'Range (2) loop
 382             R (J, K) := Operation (X (J, K));
 383          end loop;
 384       end loop;
 385
 386       return R;
 387    end Matrix_Elementwise_Operation;
 388
 389    ----------------------------------
 390    -- Vector_Elementwise_Operation --
 391    ----------------------------------
 392
 393    function Vector_Elementwise_Operation (X : X_Vector) return Result_Vector is
 394       R : Result_Vector (X'Range);
 395
 396    begin
 397       for J in R'Range loop
 398          R (J) := Operation (X (J));
 399       end loop;
 400
 401       return R;
 402    end Vector_Elementwise_Operation;
 403
 404    -----------------------------------------
 405    -- Matrix_Matrix_Elementwise_Operation --
 406    -----------------------------------------
 407
 408    function Matrix_Matrix_Elementwise_Operation
 409      (Left  : Left_Matrix;
 410       Right : Right_Matrix) return Result_Matrix
 411    is
 412       R : Result_Matrix (Left'Range (1), Left'Range (2));
 413
 414    begin
 415       if Left'Length (1) /= Right'Length (1)
 416            or else
 417          Left'Length (2) /= Right'Length (2)
 418       then
 419          raise Constraint_Error with
 420            "matrices are of different dimension in elementwise operation";
 421       end if;
 422
 423       for J in R'Range (1) loop
 424          for K in R'Range (2) loop
 425             R (J, K) :=
 426               Operation
 427                 (Left (J, K),
 428                  Right
 429                    (J - R'First (1) + Right'First (1),
 430                     K - R'First (2) + Right'First (2)));
 431          end loop;
 432       end loop;
 433
 434       return R;
 435    end Matrix_Matrix_Elementwise_Operation;
 436
 437    ------------------------------------------------
 438    -- Matrix_Matrix_Scalar_Elementwise_Operation --
 439    ------------------------------------------------
 440
 441    function Matrix_Matrix_Scalar_Elementwise_Operation
 442      (X    : X_Matrix;
 443       Y    : Y_Matrix;
 444       Z    : Z_Scalar) return Result_Matrix
 445    is
 446       R : Result_Matrix (X'Range (1), X'Range (2));
 447
 448    begin
 449       if X'Length (1) /= Y'Length (1)
 450            or else
 451          X'Length (2) /= Y'Length (2)
 452       then
 453          raise Constraint_Error with
 454            "matrices are of different dimension in elementwise operation";
 455       end if;
 456
 457       for J in R'Range (1) loop
 458          for K in R'Range (2) loop
 459             R (J, K) :=
 460               Operation
 461                 (X (J, K),
 462                  Y (J - R'First (1) + Y'First (1),
 463                     K - R'First (2) + Y'First (2)),
 464                  Z);
 465          end loop;
 466       end loop;
 467
 468       return R;
 469    end Matrix_Matrix_Scalar_Elementwise_Operation;
 470
 471    -----------------------------------------
 472    -- Vector_Vector_Elementwise_Operation --
 473    -----------------------------------------
 474
 475    function Vector_Vector_Elementwise_Operation
 476      (Left  : Left_Vector;
 477       Right : Right_Vector) return Result_Vector
 478    is
 479       R : Result_Vector (Left'Range);
 480
 481    begin
 482       if Left'Length /= Right'Length then
 483          raise Constraint_Error with
 484            "vectors are of different length in elementwise operation";
 485       end if;
 486
 487       for J in R'Range loop
 488          R (J) := Operation (Left (J), Right (J - R'First + Right'First));
 489       end loop;
 490
 491       return R;
 492    end Vector_Vector_Elementwise_Operation;
 493
 494    ------------------------------------------------
 495    -- Vector_Vector_Scalar_Elementwise_Operation --
 496    ------------------------------------------------
 497
 498    function Vector_Vector_Scalar_Elementwise_Operation
 499      (X : X_Vector;
 500       Y : Y_Vector;
 501       Z : Z_Scalar) return Result_Vector
 502    is
 503       R : Result_Vector (X'Range);
 504
 505    begin
 506       if X'Length /= Y'Length then
 507          raise Constraint_Error with
 508            "vectors are of different length in elementwise operation";
 509       end if;
 510
 511       for J in R'Range loop
 512          R (J) := Operation (X (J), Y (J - X'First + Y'First), Z);
 513       end loop;
 514
 515       return R;
 516    end Vector_Vector_Scalar_Elementwise_Operation;
 517
 518    -----------------------------------------
 519    -- Matrix_Scalar_Elementwise_Operation --
 520    -----------------------------------------
 521
 522    function Matrix_Scalar_Elementwise_Operation
 523      (Left  : Left_Matrix;
 524       Right : Right_Scalar) return Result_Matrix
 525    is
 526       R : Result_Matrix (Left'Range (1), Left'Range (2));
 527
 528    begin
 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
 535       return R;
 536    end Matrix_Scalar_Elementwise_Operation;
 537
 538    -----------------------------------------
 539    -- Vector_Scalar_Elementwise_Operation --
 540    -----------------------------------------
 541
 542    function Vector_Scalar_Elementwise_Operation
 543      (Left  : Left_Vector;
 544       Right : Right_Scalar) return Result_Vector
 545    is
 546       R : Result_Vector (Left'Range);
 547
 548    begin
 549       for J in R'Range loop
 550          R (J) := Operation (Left (J), Right);
 551       end loop;
 552
 553       return R;
 554    end Vector_Scalar_Elementwise_Operation;
 555
 556    -----------------------------------------
 557    -- Scalar_Matrix_Elementwise_Operation --
 558    -----------------------------------------
 559
 560    function Scalar_Matrix_Elementwise_Operation
 561      (Left  : Left_Scalar;
 562       Right : Right_Matrix) return Result_Matrix
 563    is
 564       R : Result_Matrix (Right'Range (1), Right'Range (2));
 565
 566    begin
 567       for J in R'Range (1) loop
 568          for K in R'Range (2) loop
 569             R (J, K) := Operation (Left, Right (J, K));
 570          end loop;
 571       end loop;
 572
 573       return R;
 574    end Scalar_Matrix_Elementwise_Operation;
 575
 576    -----------------------------------------
 577    -- Scalar_Vector_Elementwise_Operation --
 578    -----------------------------------------
 579
 580    function Scalar_Vector_Elementwise_Operation
 581      (Left  : Left_Scalar;
 582       Right : Right_Vector) return Result_Vector
 583    is
 584       R : Result_Vector (Right'Range);
 585
 586    begin
 587       for J in R'Range loop
 588          R (J) := Operation (Left, Right (J));
 589       end loop;
 590
 591       return R;
 592    end Scalar_Vector_Elementwise_Operation;
 593
 594    ----------
 595    -- Sqrt --
 596    ----------
 597
 598    function Sqrt (X : Real'Base) return Real'Base is
 599       Root, Next : Real'Base;
 600
 601    begin
 602       --  Be defensive: any comparisons with NaN values will yield False.
 603
 604       if not (X > 0.0) then
 605          if X = 0.0 then
 606             return X;
 607          else
 608             raise Argument_Error;
 609          end if;
 610
 611       elsif X > Real'Base'Last then
 612
 613          --  X is infinity, which is its own square root
 614
 615          return X;
 616       end if;
 617
 618       --  Compute an initial estimate based on:
 619
 620       --     X = M * R**E and Sqrt (X) = Sqrt (M) * R**(E / 2.0),
 621
 622       --  where M is the mantissa, R is the radix and E the exponent.
 623
 624       --  By ignoring the mantissa and ignoring the case of an odd
 625       --  exponent, we get a final error that is at most R. In other words,
 626       --  the result has about a single bit precision.
 627
 628       Root := Real'Base (Real'Machine_Radix) ** (Real'Exponent (X) / 2);
 629
 630       --  Because of the poor initial estimate, use the Babylonian method of
 631       --  computing the square root, as it is stable for all inputs. Every step
 632       --  will roughly double the precision of the result. Just a few steps
 633       --  suffice in most cases. Eight iterations should give about 2**8 bits
 634       --  of precision.
 635
 636       for J in 1 .. 8 loop
 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       R : Result_Matrix (Left'Range (1), Right'Range (2));
 654
 655    begin
 656       if Left'Length (2) /= Right'Length (1) then
 657          raise Constraint_Error with
 658            "incompatible dimensions in matrix multiplication";
 659       end if;
 660
 661       for J in R'Range (1) loop
 662          for K in R'Range (2) loop
 663             declare
 664                S : Result_Scalar := Zero;
 665
 666             begin
 667                for M in Left'Range (2) loop
 668                   S := S + Left (J, M) *
 669                              Right (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
 677       return R;
 678    end  Matrix_Matrix_Product;
 679
 680    ----------------------------
 681    -- Matrix_Vector_Solution --
 682    ----------------------------
 683
 684    function Matrix_Vector_Solution (A : Matrix; X : Vector) return Vector is
 685       N   : constant Natural := A'Length (1);
 686       MA  : Matrix := A;
 687       MX  : Matrix (A'Range (1), 1 .. 1);
 688       R   : Vector (A'Range (2));
 689       Det : Scalar;
 690
 691    begin
 692       if A'Length (2) /= N then
 693          raise Constraint_Error with "matrix is not square";
 694       end if;
 695
 696       if X'Length /= N then
 697          raise Constraint_Error with "incompatible vector length";
 698       end if;
 699
 700       for J in 0 .. MX'Length (1) - 1 loop
 701          MX (MX'First (1) + J, 1) := X (X'First + J);
 702       end loop;
 703
 704       Forward_Eliminate (MA, MX, Det);
 705       Back_Substitute (MA, MX);
 706
 707       for J in 0 .. R'Length - 1 loop
 708          R (R'First + J) := MX (MX'First (1) + J, 1);
 709       end loop;
 710
 711       return R;
 712    end Matrix_Vector_Solution;
 713
 714    ----------------------------
 715    -- Matrix_Matrix_Solution --
 716    ----------------------------
 717
 718    function Matrix_Matrix_Solution (A, X : Matrix) return Matrix is
 719       N   : constant Natural := A'Length (1);
 720       MA  : Matrix (A'Range (2), A'Range (2));
 721       MB  : Matrix (A'Range (2), X'Range (2));
 722       Det : Scalar;
 723
 724    begin
 725       if A'Length (2) /= N then
 726          raise Constraint_Error with "matrix is not square";
 727       end if;
 728
 729       if X'Length (1) /= N then
 730          raise Constraint_Error with "matrices have unequal number of rows";
 731       end if;
 732
 733       for J in 0 .. A'Length (1) - 1 loop
 734          for K in MA'Range (2) loop
 735             MA (MA'First (1) + J, K) := A (A'First (1) + J, K);
 736          end loop;
 737
 738          for K in MB'Range (2) loop
 739             MB (MB'First (1) + J, K) := X (X'First (1) + J, K);
 740          end loop;
 741       end loop;
 742
 743       Forward_Eliminate (MA, MB, Det);
 744       Back_Substitute (MA, MB);
 745
 746       return MB;
 747    end Matrix_Matrix_Solution;
 748
 749    ---------------------------
 750    -- Matrix_Vector_Product --
 751    ---------------------------
 752
 753    function Matrix_Vector_Product
 754      (Left  : Matrix;
 755       Right : Right_Vector) return Result_Vector
 756    is
 757       R : Result_Vector (Left'Range (1));
 758
 759    begin
 760       if Left'Length (2) /= Right'Length then
 761          raise Constraint_Error with
 762             "incompatible dimensions in matrix-vector multiplication";
 763       end if;
 764
 765       for J in Left'Range (1) loop
 766          declare
 767             S : Result_Scalar := Zero;
 768
 769          begin
 770             for K in Left'Range (2) loop
 771                S := S + Left (J, K) * Right (K - Left'First (2) + Right'First);
 772             end loop;
 773
 774             R (J) := S;
 775          end;
 776       end loop;
 777
 778       return R;
 779    end Matrix_Vector_Product;
 780
 781    -------------------
 782    -- Outer_Product --
 783    -------------------
 784
 785    function Outer_Product
 786      (Left  : Left_Vector;
 787       Right : Right_Vector) return Matrix
 788    is
 789       R : Matrix (Left'Range, Right'Range);
 790
 791    begin
 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
 798       return R;
 799    end Outer_Product;
 800
 801    -----------------
 802    -- Swap_Column --
 803    -----------------
 804
 805    procedure Swap_Column (A : in out Matrix; Left, Right : Integer) is
 806       Temp : Scalar;
 807    begin
 808       for J in A'Range (1) loop
 809          Temp := A (J, Left);
 810          A (J, Left) := A (J, Right);
 811          A (J, Right) := Temp;
 812       end loop;
 813    end Swap_Column;
 814
 815    ---------------
 816    -- Transpose --
 817    ---------------
 818
 819    procedure Transpose (A : Matrix; R : out Matrix) is
 820    begin
 821       for J in R'Range (1) loop
 822          for K in R'Range (2) loop
 823             R (J, K) := A (K - R'First (2) + A'First (1),
 824                            J - R'First (1) + A'First (2));
 825          end loop;
 826       end loop;
 827    end Transpose;
 828
 829    -------------------------------
 830    -- Update_Matrix_With_Matrix --
 831    -------------------------------
 832
 833    procedure Update_Matrix_With_Matrix (X : in out X_Matrix; Y : Y_Matrix) is
 834    begin
 835       if X'Length (1) /= Y'Length (1)
 836         or else 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       R : Matrix (First_1 .. Check_Unit_Last (First_1, Order, First_1),
 876                   First_2 .. Check_Unit_Last (First_2, Order, First_2));
 877
 878    begin
 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
 885       return R;
 886    end Unit_Matrix;
 887
 888    -----------------
 889    -- Unit_Vector --
 890    -----------------
 891
 892    function Unit_Vector
 893      (Index : Integer;
 894       Order : Positive;
 895       First : Integer := 1) return Vector
 896    is
 897       R : Vector (First .. Check_Unit_Last (Index, Order, First));
 898    begin
 899       R := (others => Zero);
 900       R (Index) := One;
 901       return R;
 902    end Unit_Vector;
 903
 904    ---------------------------
 905    -- Vector_Matrix_Product --
 906    ---------------------------
 907
 908    function Vector_Matrix_Product
 909      (Left  : Left_Vector;
 910       Right : Matrix) return Result_Vector
 911    is
 912       R : Result_Vector (Right'Range (2));
 913
 914    begin
 915       if Left'Length /= Right'Length (2) then
 916          raise Constraint_Error with
 917            "incompatible dimensions in vector-matrix multiplication";
 918       end if;
 919
 920       for J in Right'Range (2) loop
 921          declare
 922             S : Result_Scalar := Zero;
 923
 924          begin
 925             for K in Right'Range (1) loop
 926                S := S + Left (J - Right'First (1) + Left'First) * Right (K, J);
 927             end loop;
 928
 929             R (J) := S;
 930          end;
 931       end loop;
 932
 933       return R;
 934    end Vector_Matrix_Product;
 935
 936 end System.Generic_Array_Operations;