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    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
 926 end System.Generic_Array_Operations;