gcc/ada/libgnat/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-2018, Free Software Foundation, Inc.          --
  10 --                                                                          --
  11 -- GNAT is free software;  you can  redistribute it  and/or modify it under --
  12 -- terms of the  GNU General Public License as published  by the Free Soft- --
  13 -- ware  Foundation;  either version 3,  or (at your option) any later ver- --
  14 -- sion.  GNAT is distributed in the hope that it will be useful, but WITH- --
  15 -- OUT ANY WARRANTY;  without even the  implied warranty of MERCHANTABILITY --
  16 -- or FITNESS FOR A PARTICULAR PURPOSE.                                     --
  17 --                                                                          --
  18 -- As a special exception under Section 7 of GPL version 3, you are granted --
  19 -- additional permissions described in the GCC Runtime Library Exception,   --
  20 -- version 3.1, as published by the Free Software Foundation.               --
  21 --                                                                          --
  22 -- You should have received a copy of the GNU General Public License and    --
  23 -- a copy of the GCC Runtime Library Exception along with this program;     --
  24 -- see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see    --
  25 -- <http://www.gnu.org/licenses/>.                                          --
  26 --                                                                          --
  27 -- GNAT was originally developed  by the GNAT team at  New York University. --
  28 -- Extensive contributions were provided by Ada Core Technologies Inc.      --
  29 --                                                                          --
  30 ------------------------------------------------------------------------------
  31
  32 with Ada.Numerics; use Ada.Numerics;
  33 package body System.Generic_Array_Operations is
  34    function Check_Unit_Last
  35      (Index : Integer;
  36       Order : Positive;
  37       First : Integer) return Integer;
  38    pragma Inline_Always (Check_Unit_Last);
  39    --  Compute index of last element returned by Unit_Vector or Unit_Matrix.
  40    --  A separate function is needed to allow raising Constraint_Error before
  41    --  declaring the function result variable. The result variable needs to be
  42    --  declared first, to allow front-end inlining.
  43
  44    --------------
  45    -- Diagonal --
  46    --------------
  47
  48    function Diagonal (A : Matrix) return Vector is
  49       N : constant Natural := Natural'Min (A'Length (1), A'Length (2));
  50    begin
  51       return R : Vector (A'First (1) .. A'First (1) + N - 1) do
  52          for J in 0 .. N - 1 loop
  53             R (R'First + J) := A (A'First (1) + J, A'First (2) + J);
  54          end loop;
  55       end return;
  56    end Diagonal;
  57
  58    --------------------------
  59    -- Square_Matrix_Length --
  60    --------------------------
  61
  62    function Square_Matrix_Length (A : Matrix) return Natural is
  63    begin
  64       if A'Length (1) /= A'Length (2) then
  65          raise Constraint_Error with "matrix is not square";
  66       else
  67          return A'Length (1);
  68       end if;
  69    end Square_Matrix_Length;
  70
  71    ---------------------
  72    -- Check_Unit_Last --
  73    ---------------------
  74
  75    function Check_Unit_Last
  76       (Index : Integer;
  77        Order : Positive;
  78        First : Integer) return Integer
  79    is
  80    begin
  81       --  Order the tests carefully to avoid overflow
  82
  83       if Index < First
  84         or else First > Integer'Last - Order + 1
  85         or else Index > First + (Order - 1)
  86       then
  87          raise Constraint_Error;
  88       end if;
  89
  90       return First + (Order - 1);
  91    end Check_Unit_Last;
  92
  93    ---------------------
  94    -- Back_Substitute --
  95    ---------------------
  96
  97    procedure Back_Substitute (M, N : in out Matrix) is
  98       pragma Assert (M'First (1) = N'First (1)
  99                        and then
 100                      M'Last  (1) = N'Last (1));
 101
 102       procedure Sub_Row
 103         (M      : in out Matrix;
 104          Target : Integer;
 105          Source : Integer;
 106          Factor : Scalar);
 107       --  Elementary row operation that subtracts Factor * M (Source, <>) from
 108       --  M (Target, <>)
 109
 110       -------------
 111       -- Sub_Row --
 112       -------------
 113
 114       procedure Sub_Row
 115         (M      : in out Matrix;
 116          Target : Integer;
 117          Source : Integer;
 118          Factor : Scalar)
 119       is
 120       begin
 121          for J in M'Range (2) loop
 122             M (Target, J) := M (Target, J) - Factor * M (Source, J);
 123          end loop;
 124       end Sub_Row;
 125
 126       --  Local declarations
 127
 128       Max_Col : Integer := M'Last (2);
 129
 130    --  Start of processing for Back_Substitute
 131
 132    begin
 133       Do_Rows : for Row in reverse M'Range (1) loop
 134          Find_Non_Zero : for Col in reverse M'First (2) .. Max_Col loop
 135             if Is_Non_Zero (M (Row, Col)) then
 136
 137                --  Found first non-zero element, so subtract a multiple of this
 138                --  element  from all higher rows, to reduce all other elements
 139                --  in this column to zero.
 140
 141                declare
 142                   --  We can't use a for loop, as we'd need to iterate to
 143                   --  Row - 1, but that expression will overflow if M'First
 144                   --  equals Integer'First, which is true for aggregates
 145                   --  without explicit bounds..
 146
 147                   J : Integer := M'First (1);
 148
 149                begin
 150                   while J < Row loop
 151                      Sub_Row (N, J, Row, (M (J, Col) / M (Row, Col)));
 152                      Sub_Row (M, J, Row, (M (J, Col) / M (Row, Col)));
 153                      J := J + 1;
 154                   end loop;
 155                end;
 156
 157                --  Avoid potential overflow in the subtraction below
 158
 159                exit Do_Rows when Col = M'First (2);
 160
 161                Max_Col := Col - 1;
 162
 163                exit Find_Non_Zero;
 164             end if;
 165          end loop Find_Non_Zero;
 166       end loop Do_Rows;
 167    end Back_Substitute;
 168
 169    -----------------------
 170    -- Forward_Eliminate --
 171    -----------------------
 172
 173    procedure Forward_Eliminate
 174      (M   : in out Matrix;
 175       N   : in out Matrix;
 176       Det : out Scalar)
 177    is
 178       pragma Assert (M'First (1) = N'First (1)
 179                        and then
 180                      M'Last  (1) = N'Last (1));
 181
 182       --  The following are variations of the elementary matrix row operations:
 183       --  row switching, row multiplication and row addition. Because in this
 184       --  algorithm the addition factor is always a negated value, we chose to
 185       --  use  row subtraction instead. Similarly, instead of multiplying by
 186       --  a reciprocal, we divide.
 187
 188       procedure Sub_Row
 189         (M      : in out Matrix;
 190          Target : Integer;
 191          Source : Integer;
 192          Factor : Scalar);
 193       --  Subtrace Factor * M (Source, <>) from M (Target, <>)
 194
 195       procedure Divide_Row
 196         (M, N  : in out Matrix;
 197          Row   : Integer;
 198          Scale : Scalar);
 199       --  Divide M (Row) and N (Row) by Scale, and update Det
 200
 201       procedure Switch_Row
 202         (M, N  : in out Matrix;
 203          Row_1 : Integer;
 204          Row_2 : Integer);
 205       --  Exchange M (Row_1) and N (Row_1) with M (Row_2) and N (Row_2),
 206       --  negating Det in the process.
 207
 208       -------------
 209       -- Sub_Row --
 210       -------------
 211
 212       procedure Sub_Row
 213         (M      : in out Matrix;
 214          Target : Integer;
 215          Source : Integer;
 216          Factor : Scalar)
 217       is
 218       begin
 219          for J in M'Range (2) loop
 220             M (Target, J) := M (Target, J) - Factor * M (Source, J);
 221          end loop;
 222       end Sub_Row;
 223
 224       ----------------
 225       -- Divide_Row --
 226       ----------------
 227
 228       procedure Divide_Row
 229         (M, N  : in out Matrix;
 230          Row   : Integer;
 231          Scale : Scalar)
 232       is
 233       begin
 234          Det := Det * Scale;
 235
 236          for J in M'Range (2) loop
 237             M (Row, J) := M (Row, J) / Scale;
 238          end loop;
 239
 240          for J in N'Range (2) loop
 241             N (Row - M'First (1) + N'First (1), J) :=
 242               N (Row - M'First (1) + N'First (1), J) / Scale;
 243          end loop;
 244       end Divide_Row;
 245
 246       ----------------
 247       -- Switch_Row --
 248       ----------------
 249
 250       procedure Switch_Row
 251         (M, N  : in out Matrix;
 252          Row_1 : Integer;
 253          Row_2 : Integer)
 254       is
 255          procedure Swap (X, Y : in out Scalar);
 256          --  Exchange the values of X and Y
 257
 258          ----------
 259          -- Swap --
 260          ----------
 261
 262          procedure Swap (X, Y : in out Scalar) is
 263             T : constant Scalar := X;
 264          begin
 265             X := Y;
 266             Y := T;
 267          end Swap;
 268
 269       --  Start of processing for Switch_Row
 270
 271       begin
 272          if Row_1 /= Row_2 then
 273             Det := Zero - Det;
 274
 275             for J in M'Range (2) loop
 276                Swap (M (Row_1, J), M (Row_2, J));
 277             end loop;
 278
 279             for J in N'Range (2) loop
 280                Swap (N (Row_1 - M'First (1) + N'First (1), J),
 281                      N (Row_2 - M'First (1) + N'First (1), J));
 282             end loop;
 283          end if;
 284       end Switch_Row;
 285
 286       --  Local declarations
 287
 288       Row : Integer := M'First (1);
 289
 290    --  Start of processing for Forward_Eliminate
 291
 292    begin
 293       Det := One;
 294
 295       for J in M'Range (2) loop
 296          declare
 297             Max_Row : Integer := Row;
 298             Max_Abs : Real'Base := 0.0;
 299
 300          begin
 301             --  Find best pivot in column J, starting in row Row
 302
 303             for K in Row .. M'Last (1) loop
 304                declare
 305                   New_Abs : constant Real'Base := abs M (K, J);
 306                begin
 307                   if Max_Abs < New_Abs then
 308                      Max_Abs := New_Abs;
 309                      Max_Row := K;
 310                   end if;
 311                end;
 312             end loop;
 313
 314             if Max_Abs > 0.0 then
 315                Switch_Row (M, N, Row, Max_Row);
 316
 317                --  The temporaries below are necessary to force a copy of the
 318                --  value and avoid improper aliasing.
 319
 320                declare
 321                   Scale : constant Scalar := M (Row, J);
 322                begin
 323                   Divide_Row (M, N, Row, Scale);
 324                end;
 325
 326                for U in Row + 1 .. M'Last (1) loop
 327                   declare
 328                      Factor : constant Scalar := M (U, J);
 329                   begin
 330                      Sub_Row (N, U, Row, Factor);
 331                      Sub_Row (M, U, Row, Factor);
 332                   end;
 333                end loop;
 334
 335                exit when Row >= M'Last (1);
 336
 337                Row := Row + 1;
 338
 339             else
 340                --  Set zero (note that we do not have literals)
 341
 342                Det := Zero;
 343             end if;
 344          end;
 345       end loop;
 346    end Forward_Eliminate;
 347
 348    -------------------
 349    -- Inner_Product --
 350    -------------------
 351
 352    function Inner_Product
 353      (Left  : Left_Vector;
 354       Right : Right_Vector) return  Result_Scalar
 355    is
 356       R : Result_Scalar := Zero;
 357
 358    begin
 359       if Left'Length /= Right'Length then
 360          raise Constraint_Error with
 361             "vectors are of different length in inner product";
 362       end if;
 363
 364       for J in Left'Range loop
 365          R := R + Left (J) * Right (J - Left'First + Right'First);
 366       end loop;
 367
 368       return R;
 369    end Inner_Product;
 370
 371    -------------
 372    -- L2_Norm --
 373    -------------
 374
 375    function L2_Norm (X : X_Vector) return Result_Real'Base is
 376       Sum : Result_Real'Base := 0.0;
 377
 378    begin
 379       for J in X'Range loop
 380          Sum := Sum + Result_Real'Base (abs X (J))**2;
 381       end loop;
 382
 383       return Sqrt (Sum);
 384    end L2_Norm;
 385
 386    ----------------------------------
 387    -- Matrix_Elementwise_Operation --
 388    ----------------------------------
 389
 390    function Matrix_Elementwise_Operation (X : X_Matrix) return Result_Matrix is
 391    begin
 392       return R : Result_Matrix (X'Range (1), X'Range (2)) do
 393          for J in R'Range (1) loop
 394             for K in R'Range (2) loop
 395                R (J, K) := Operation (X (J, K));
 396             end loop;
 397          end loop;
 398       end return;
 399    end Matrix_Elementwise_Operation;
 400
 401    ----------------------------------
 402    -- Vector_Elementwise_Operation --
 403    ----------------------------------
 404
 405    function Vector_Elementwise_Operation (X : X_Vector) return Result_Vector is
 406    begin
 407       return R : Result_Vector (X'Range) do
 408          for J in R'Range loop
 409             R (J) := Operation (X (J));
 410          end loop;
 411       end return;
 412    end Vector_Elementwise_Operation;
 413
 414    -----------------------------------------
 415    -- Matrix_Matrix_Elementwise_Operation --
 416    -----------------------------------------
 417
 418    function Matrix_Matrix_Elementwise_Operation
 419      (Left  : Left_Matrix;
 420       Right : Right_Matrix) return Result_Matrix
 421    is
 422    begin
 423       return R : Result_Matrix (Left'Range (1), Left'Range (2)) do
 424          if Left'Length (1) /= Right'Length (1)
 425               or else
 426             Left'Length (2) /= Right'Length (2)
 427          then
 428             raise Constraint_Error with
 429               "matrices are of different dimension in elementwise operation";
 430          end if;
 431
 432          for J in R'Range (1) loop
 433             for K in R'Range (2) loop
 434                R (J, K) :=
 435                  Operation
 436                    (Left (J, K),
 437                     Right
 438                       (J - R'First (1) + Right'First (1),
 439                        K - R'First (2) + Right'First (2)));
 440             end loop;
 441          end loop;
 442       end return;
 443    end Matrix_Matrix_Elementwise_Operation;
 444
 445    ------------------------------------------------
 446    -- Matrix_Matrix_Scalar_Elementwise_Operation --
 447    ------------------------------------------------
 448
 449    function Matrix_Matrix_Scalar_Elementwise_Operation
 450      (X : X_Matrix;
 451       Y : Y_Matrix;
 452       Z : Z_Scalar) return Result_Matrix
 453    is
 454    begin
 455       return R : Result_Matrix (X'Range (1), X'Range (2)) do
 456          if X'Length (1) /= Y'Length (1)
 457               or else
 458             X'Length (2) /= Y'Length (2)
 459          then
 460             raise Constraint_Error with
 461               "matrices are of different dimension in elementwise operation";
 462          end if;
 463
 464          for J in R'Range (1) loop
 465             for K in R'Range (2) loop
 466                R (J, K) :=
 467                  Operation
 468                    (X (J, K),
 469                     Y (J - R'First (1) + Y'First (1),
 470                        K - R'First (2) + Y'First (2)),
 471                     Z);
 472             end loop;
 473          end loop;
 474       end return;
 475    end Matrix_Matrix_Scalar_Elementwise_Operation;
 476
 477    -----------------------------------------
 478    -- Vector_Vector_Elementwise_Operation --
 479    -----------------------------------------
 480
 481    function Vector_Vector_Elementwise_Operation
 482      (Left  : Left_Vector;
 483       Right : Right_Vector) return Result_Vector
 484    is
 485    begin
 486       return R : Result_Vector (Left'Range) do
 487          if Left'Length /= Right'Length then
 488             raise Constraint_Error with
 489               "vectors are of different length in elementwise operation";
 490          end if;
 491
 492          for J in R'Range loop
 493             R (J) := Operation (Left (J), Right (J - R'First + Right'First));
 494          end loop;
 495       end return;
 496    end Vector_Vector_Elementwise_Operation;
 497
 498    ------------------------------------------------
 499    -- Vector_Vector_Scalar_Elementwise_Operation --
 500    ------------------------------------------------
 501
 502    function Vector_Vector_Scalar_Elementwise_Operation
 503      (X : X_Vector;
 504       Y : Y_Vector;
 505       Z : Z_Scalar) return Result_Vector is
 506    begin
 507       return R : Result_Vector (X'Range) do
 508          if X'Length /= Y'Length then
 509             raise Constraint_Error with
 510               "vectors are of different length in elementwise operation";
 511          end if;
 512
 513          for J in R'Range loop
 514             R (J) := Operation (X (J), Y (J - X'First + Y'First), Z);
 515          end loop;
 516       end return;
 517    end Vector_Vector_Scalar_Elementwise_Operation;
 518
 519    -----------------------------------------
 520    -- Matrix_Scalar_Elementwise_Operation --
 521    -----------------------------------------
 522
 523    function Matrix_Scalar_Elementwise_Operation
 524      (Left  : Left_Matrix;
 525       Right : Right_Scalar) return Result_Matrix
 526    is
 527    begin
 528       return R : Result_Matrix (Left'Range (1), Left'Range (2)) do
 529          for J in R'Range (1) loop
 530             for K in R'Range (2) loop
 531                R (J, K) := Operation (Left (J, K), Right);
 532             end loop;
 533          end loop;
 534       end return;
 535    end Matrix_Scalar_Elementwise_Operation;
 536
 537    -----------------------------------------
 538    -- Vector_Scalar_Elementwise_Operation --
 539    -----------------------------------------
 540
 541    function Vector_Scalar_Elementwise_Operation
 542      (Left  : Left_Vector;
 543       Right : Right_Scalar) return Result_Vector
 544    is
 545    begin
 546       return R : Result_Vector (Left'Range) do
 547          for J in R'Range loop
 548             R (J) := Operation (Left (J), Right);
 549          end loop;
 550       end return;
 551    end Vector_Scalar_Elementwise_Operation;
 552
 553    -----------------------------------------
 554    -- Scalar_Matrix_Elementwise_Operation --
 555    -----------------------------------------
 556
 557    function Scalar_Matrix_Elementwise_Operation
 558      (Left  : Left_Scalar;
 559       Right : Right_Matrix) return Result_Matrix
 560    is
 561    begin
 562       return R : Result_Matrix (Right'Range (1), Right'Range (2)) do
 563          for J in R'Range (1) loop
 564             for K in R'Range (2) loop
 565                R (J, K) := Operation (Left, Right (J, K));
 566             end loop;
 567          end loop;
 568       end return;
 569    end Scalar_Matrix_Elementwise_Operation;
 570
 571    -----------------------------------------
 572    -- Scalar_Vector_Elementwise_Operation --
 573    -----------------------------------------
 574
 575    function Scalar_Vector_Elementwise_Operation
 576      (Left  : Left_Scalar;
 577       Right : Right_Vector) return Result_Vector
 578    is
 579    begin
 580       return R : Result_Vector (Right'Range) do
 581          for J in R'Range loop
 582             R (J) := Operation (Left, Right (J));
 583          end loop;
 584       end return;
 585    end Scalar_Vector_Elementwise_Operation;
 586
 587    ----------
 588    -- Sqrt --
 589    ----------
 590
 591    function Sqrt (X : Real'Base) return Real'Base is
 592       Root, Next : Real'Base;
 593
 594    begin
 595       --  Be defensive: any comparisons with NaN values will yield False.
 596
 597       if not (X > 0.0) then
 598          if X = 0.0 then
 599             return X;
 600          else
 601             raise Argument_Error;
 602          end if;
 603
 604       elsif X > Real'Base'Last then
 605
 606          --  X is infinity, which is its own square root
 607
 608          return X;
 609       end if;
 610
 611       --  Compute an initial estimate based on:
 612
 613       --     X = M * R**E and Sqrt (X) = Sqrt (M) * R**(E / 2.0),
 614
 615       --  where M is the mantissa, R is the radix and E the exponent.
 616
 617       --  By ignoring the mantissa and ignoring the case of an odd
 618       --  exponent, we get a final error that is at most R. In other words,
 619       --  the result has about a single bit precision.
 620
 621       Root := Real'Base (Real'Machine_Radix) ** (Real'Exponent (X) / 2);
 622
 623       --  Because of the poor initial estimate, use the Babylonian method of
 624       --  computing the square root, as it is stable for all inputs. Every step
 625       --  will roughly double the precision of the result. Just a few steps
 626       --  suffice in most cases. Eight iterations should give about 2**8 bits
 627       --  of precision.
 628
 629       for J in 1 .. 8 loop
 630          Next := (Root + X / Root) / 2.0;
 631          exit when Root = Next;
 632          Root := Next;
 633       end loop;
 634
 635       return Root;
 636    end Sqrt;
 637
 638    ---------------------------
 639    -- Matrix_Matrix_Product --
 640    ---------------------------
 641
 642    function Matrix_Matrix_Product
 643      (Left  : Left_Matrix;
 644       Right : Right_Matrix) return Result_Matrix
 645    is
 646    begin
 647       return R : Result_Matrix (Left'Range (1), Right'Range (2)) do
 648          if Left'Length (2) /= Right'Length (1) then
 649             raise Constraint_Error with
 650               "incompatible dimensions in matrix multiplication";
 651          end if;
 652
 653          for J in R'Range (1) loop
 654             for K in R'Range (2) loop
 655                declare
 656                   S : Result_Scalar := Zero;
 657
 658                begin
 659                   for M in Left'Range (2) loop
 660                      S := S + Left (J, M) *
 661                                 Right
 662                                   (M - Left'First (2) + Right'First (1), K);
 663                   end loop;
 664
 665                   R (J, K) := S;
 666                end;
 667             end loop;
 668          end loop;
 669       end return;
 670    end  Matrix_Matrix_Product;
 671
 672    ----------------------------
 673    -- Matrix_Vector_Solution --
 674    ----------------------------
 675
 676    function Matrix_Vector_Solution (A : Matrix; X : Vector) return Vector is
 677       N   : constant Natural := A'Length (1);
 678       MA  : Matrix := A;
 679       MX  : Matrix (A'Range (1), 1 .. 1);
 680       R   : Vector (A'Range (2));
 681       Det : Scalar;
 682
 683    begin
 684       if A'Length (2) /= N then
 685          raise Constraint_Error with "matrix is not square";
 686       end if;
 687
 688       if X'Length /= N then
 689          raise Constraint_Error with "incompatible vector length";
 690       end if;
 691
 692       for J in 0 .. MX'Length (1) - 1 loop
 693          MX (MX'First (1) + J, 1) := X (X'First + J);
 694       end loop;
 695
 696       Forward_Eliminate (MA, MX, Det);
 697
 698       if Det = Zero then
 699          raise Constraint_Error with "matrix is singular";
 700       end if;
 701
 702       Back_Substitute (MA, MX);
 703
 704       for J in 0 .. R'Length - 1 loop
 705          R (R'First + J) := MX (MX'First (1) + J, 1);
 706       end loop;
 707
 708       return R;
 709    end Matrix_Vector_Solution;
 710
 711    ----------------------------
 712    -- Matrix_Matrix_Solution --
 713    ----------------------------
 714
 715    function Matrix_Matrix_Solution (A, X : Matrix) return Matrix is
 716       N   : constant Natural := A'Length (1);
 717       MA  : Matrix (A'Range (2), A'Range (2));
 718       MB  : Matrix (A'Range (2), X'Range (2));
 719       Det : Scalar;
 720
 721    begin
 722       if A'Length (2) /= N then
 723          raise Constraint_Error with "matrix is not square";
 724       end if;
 725
 726       if X'Length (1) /= N then
 727          raise Constraint_Error with "matrices have unequal number of rows";
 728       end if;
 729
 730       for J in 0 .. A'Length (1) - 1 loop
 731          for K in MA'Range (2) loop
 732             MA (MA'First (1) + J, K) := A (A'First (1) + J, K);
 733          end loop;
 734
 735          for K in MB'Range (2) loop
 736             MB (MB'First (1) + J, K) := X (X'First (1) + J, K);
 737          end loop;
 738       end loop;
 739
 740       Forward_Eliminate (MA, MB, Det);
 741
 742       if Det = Zero then
 743          raise Constraint_Error with "matrix is singular";
 744       end if;
 745
 746       Back_Substitute (MA, MB);
 747
 748       return MB;
 749    end Matrix_Matrix_Solution;
 750
 751    ---------------------------
 752    -- Matrix_Vector_Product --
 753    ---------------------------
 754
 755    function Matrix_Vector_Product
 756      (Left  : Matrix;
 757       Right : Right_Vector) return Result_Vector
 758    is
 759    begin
 760       return R : Result_Vector (Left'Range (1)) do
 761          if Left'Length (2) /= Right'Length then
 762             raise Constraint_Error with
 763               "incompatible dimensions in matrix-vector multiplication";
 764          end if;
 765
 766          for J in Left'Range (1) loop
 767             declare
 768                S : Result_Scalar := Zero;
 769
 770             begin
 771                for K in Left'Range (2) loop
 772                   S := S + Left (J, K)
 773                          * Right (K - Left'First (2) + Right'First);
 774                end loop;
 775
 776                R (J) := S;
 777             end;
 778          end loop;
 779       end return;
 780    end Matrix_Vector_Product;
 781
 782    -------------------
 783    -- Outer_Product --
 784    -------------------
 785
 786    function Outer_Product
 787      (Left  : Left_Vector;
 788       Right : Right_Vector) return Matrix
 789    is
 790    begin
 791       return R : Matrix (Left'Range, Right'Range) do
 792          for J in R'Range (1) loop
 793             for K in R'Range (2) loop
 794                R (J, K) := Left (J) * Right (K);
 795             end loop;
 796          end loop;
 797       end return;
 798    end Outer_Product;
 799
 800    -----------------
 801    -- Swap_Column --
 802    -----------------
 803
 804    procedure Swap_Column (A : in out Matrix; Left, Right : Integer) is
 805       Temp : Scalar;
 806    begin
 807       for J in A'Range (1) loop
 808          Temp := A (J, Left);
 809          A (J, Left) := A (J, Right);
 810          A (J, Right) := Temp;
 811       end loop;
 812    end Swap_Column;
 813
 814    ---------------
 815    -- Transpose --
 816    ---------------
 817
 818    procedure Transpose (A : Matrix; R : out Matrix) is
 819    begin
 820       for J in R'Range (1) loop
 821          for K in R'Range (2) loop
 822             R (J, K) := A (K - R'First (2) + A'First (1),
 823                            J - R'First (1) + A'First (2));
 824          end loop;
 825       end loop;
 826    end Transpose;
 827
 828    -------------------------------
 829    -- Update_Matrix_With_Matrix --
 830    -------------------------------
 831
 832    procedure Update_Matrix_With_Matrix (X : in out X_Matrix; Y : Y_Matrix) is
 833    begin
 834       if X'Length (1) /= Y'Length (1)
 835            or else
 836          X'Length (2) /= Y'Length (2)
 837       then
 838          raise Constraint_Error with
 839            "matrices are of different dimension in update operation";
 840       end if;
 841
 842       for J in X'Range (1) loop
 843          for K in X'Range (2) loop
 844             Update (X (J, K), Y (J - X'First (1) + Y'First (1),
 845                                  K - X'First (2) + Y'First (2)));
 846          end loop;
 847       end loop;
 848    end Update_Matrix_With_Matrix;
 849
 850    -------------------------------
 851    -- Update_Vector_With_Vector --
 852    -------------------------------
 853
 854    procedure Update_Vector_With_Vector (X : in out X_Vector; Y : Y_Vector) is
 855    begin
 856       if X'Length /= Y'Length then
 857          raise Constraint_Error with
 858            "vectors are of different length in update operation";
 859       end if;
 860
 861       for J in X'Range loop
 862          Update (X (J), Y (J - X'First + Y'First));
 863       end loop;
 864    end Update_Vector_With_Vector;
 865
 866    -----------------
 867    -- Unit_Matrix --
 868    -----------------
 869
 870    function Unit_Matrix
 871      (Order   : Positive;
 872       First_1 : Integer := 1;
 873       First_2 : Integer := 1) return Matrix
 874    is
 875    begin
 876       return R : Matrix (First_1 .. Check_Unit_Last (First_1, Order, First_1),
 877                          First_2 .. Check_Unit_Last (First_2, Order, First_2))
 878       do
 879          R := (others => (others => Zero));
 880
 881          for J in 0 .. Order - 1 loop
 882             R (First_1 + J, First_2 + J) := One;
 883          end loop;
 884       end return;
 885    end Unit_Matrix;
 886
 887    -----------------
 888    -- Unit_Vector --
 889    -----------------
 890
 891    function Unit_Vector
 892      (Index : Integer;
 893       Order : Positive;
 894       First : Integer := 1) return Vector
 895    is
 896    begin
 897       return R : Vector (First .. Check_Unit_Last (Index, Order, First)) do
 898          R := (others => Zero);
 899          R (Index) := One;
 900       end return;
 901    end Unit_Vector;
 902
 903    ---------------------------
 904    -- Vector_Matrix_Product --
 905    ---------------------------
 906
 907    function Vector_Matrix_Product
 908      (Left  : Left_Vector;
 909       Right : Matrix) return Result_Vector
 910    is
 911    begin
 912       return R : Result_Vector (Right'Range (2)) do
 913          if Left'Length /= Right'Length (1) then
 914             raise Constraint_Error with
 915               "incompatible dimensions in vector-matrix multiplication";
 916          end if;
 917
 918          for J in Right'Range (2) loop
 919             declare
 920                S : Result_Scalar := Zero;
 921
 922             begin
 923                for K in Right'Range (1) loop
 924                   S := S + Left (K - Right'First (1)
 925                                    + Left'First) * Right (K, J);
 926                end loop;
 927
 928                R (J) := S;
 929             end;
 930          end loop;
 931       end return;
 932    end Vector_Matrix_Product;
 933
 934 end System.Generic_Array_Operations;