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