gcc/ada/s-gearop.adb

   1 ------------------------------------------------------------------------------
   2 --                                                                          --
   3 --                         GNAT RUN-TIME COMPONENTS                         --
   4 --                                                                          --
   5 --                     SYSTEM.GENERIC_ARRAY_OPERATIONS                      --
   6 --                                                                          --
   7 --                                 B o d y                                  --
   8 --                                                                          --
   9 --            Copyright (C) 2006, 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 2,  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.  See the GNU General Public License --
  17 -- for  more details.  You should have  received  a copy of the GNU General --
  18 -- Public License  distributed with GNAT;  see file COPYING.  If not, write --
  19 -- to  the  Free Software Foundation,  51  Franklin  Street,  Fifth  Floor, --
  20 -- Boston, MA 02110-1301, USA.                                              --
  21 --                                                                          --
  22 -- As a special exception,  if other files  instantiate  generics from this --
  23 -- unit, or you link  this unit with other files  to produce an executable, --
  24 -- this  unit  does not  by itself cause  the resulting  executable  to  be --
  25 -- covered  by the  GNU  General  Public  License.  This exception does not --
  26 -- however invalidate  any other reasons why  the executable file  might be --
  27 -- covered by the  GNU Public License.                                      --
  28 --                                                                          --
  29 -- GNAT was originally developed  by the GNAT team at  New York University. --
  30 -- Extensive contributions were provided by Ada Core Technologies Inc.      --
  31 --                                                                          --
  32 ------------------------------------------------------------------------------
  33
  34 package body System.Generic_Array_Operations is
  35
  36    --  The local function Check_Unit_Last computes the index
  37    --  of the last element returned by Unit_Vector or Unit_Matrix.
  38    --  A separate function is needed to allow raising Constraint_Error
  39    --  before declaring the function result variable. The result variable
  40    --  needs to be declared first, to 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    function Square_Matrix_Length (A : Matrix) return Natural is
  49    begin
  50       if A'Length (1) /= A'Length (2) then
  51          raise Constraint_Error with "matrix is not square";
  52       end if;
  53
  54       return A'Length (1);
  55    end Square_Matrix_Length;
  56
  57    ---------------------
  58    -- Check_Unit_Last --
  59    ---------------------
  60
  61    function Check_Unit_Last
  62       (Index : Integer;
  63        Order : Positive;
  64        First : Integer) return Integer is
  65    begin
  66       --  Order the tests carefully to avoid overflow
  67
  68       if Index < First
  69            or else First > Integer'Last - Order + 1
  70            or else Index > First + (Order - 1)
  71       then
  72          raise Constraint_Error;
  73       end if;
  74
  75       return First + (Order - 1);
  76    end Check_Unit_Last;
  77
  78    -------------------
  79    -- Inner_Product --
  80    -------------------
  81
  82    function Inner_Product
  83      (Left  : Left_Vector;
  84       Right : Right_Vector)
  85       return  Result_Scalar
  86    is
  87       R : Result_Scalar := Zero;
  88
  89    begin
  90       if Left'Length /= Right'Length then
  91          raise Constraint_Error with
  92             "vectors are of different length in inner product";
  93       end if;
  94
  95       for J in Left'Range loop
  96          R := R + Left (J) * Right (J - Left'First + Right'First);
  97       end loop;
  98
  99       return R;
 100    end Inner_Product;
 101
 102    ----------------------------------
 103    -- Matrix_Elementwise_Operation --
 104    ----------------------------------
 105
 106    function Matrix_Elementwise_Operation (X : X_Matrix) return Result_Matrix is
 107       R : Result_Matrix (X'Range (1), X'Range (2));
 108
 109    begin
 110       for J in R'Range (1) loop
 111          for K in R'Range (2) loop
 112             R (J, K) := Operation (X (J, K));
 113          end loop;
 114       end loop;
 115
 116       return R;
 117    end Matrix_Elementwise_Operation;
 118
 119    ----------------------------------
 120    -- Vector_Elementwise_Operation --
 121    ----------------------------------
 122
 123    function Vector_Elementwise_Operation (X : X_Vector) return Result_Vector is
 124       R : Result_Vector (X'Range);
 125
 126    begin
 127       for J in R'Range loop
 128          R (J) := Operation (X (J));
 129       end loop;
 130
 131       return R;
 132    end Vector_Elementwise_Operation;
 133
 134    -----------------------------------------
 135    -- Matrix_Matrix_Elementwise_Operation --
 136    -----------------------------------------
 137
 138    function Matrix_Matrix_Elementwise_Operation
 139      (Left  : Left_Matrix;
 140       Right : Right_Matrix)
 141       return Result_Matrix
 142    is
 143       R : Result_Matrix (Left'Range (1), Left'Range (2));
 144    begin
 145       if Left'Length (1) /= Right'Length (1)
 146         or else Left'Length (2) /= Right'Length (2)
 147       then
 148          raise Constraint_Error with
 149             "matrices are of different dimension in elementwise operation";
 150       end if;
 151
 152       for J in R'Range (1) loop
 153          for K in R'Range (2) loop
 154             R (J, K) := Operation (Left (J, K), Right (J, K));
 155          end loop;
 156       end loop;
 157
 158       return R;
 159    end Matrix_Matrix_Elementwise_Operation;
 160
 161    ------------------------------------------------
 162    -- Matrix_Matrix_Scalar_Elementwise_Operation --
 163    ------------------------------------------------
 164
 165    function Matrix_Matrix_Scalar_Elementwise_Operation
 166      (X    : X_Matrix;
 167       Y    : Y_Matrix;
 168       Z    : Z_Scalar) return Result_Matrix
 169    is
 170       R : Result_Matrix (X'Range (1), X'Range (2));
 171
 172    begin
 173       if X'Length (1) /= Y'Length (1)
 174         or else X'Length (2) /= Y'Length (2)
 175       then
 176          raise Constraint_Error with
 177             "matrices are of different dimension in elementwise operation";
 178       end if;
 179
 180       for J in R'Range (1) loop
 181          for K in R'Range (2) loop
 182             R (J, K) := Operation (X (J, K), Y (J, K), Z);
 183          end loop;
 184       end loop;
 185
 186       return R;
 187    end Matrix_Matrix_Scalar_Elementwise_Operation;
 188
 189    -----------------------------------------
 190    -- Vector_Vector_Elementwise_Operation --
 191    -----------------------------------------
 192
 193    function Vector_Vector_Elementwise_Operation
 194      (Left  : Left_Vector;
 195       Right : Right_Vector) return Result_Vector
 196    is
 197       R : Result_Vector (Left'Range);
 198
 199    begin
 200       if Left'Length /= Right'Length then
 201          raise Constraint_Error with
 202             "vectors are of different length in elementwise operation";
 203       end if;
 204
 205       for J in R'Range loop
 206          R (J) := Operation (Left (J), Right (J));
 207       end loop;
 208
 209       return R;
 210    end Vector_Vector_Elementwise_Operation;
 211
 212    ------------------------------------------------
 213    -- Vector_Vector_Scalar_Elementwise_Operation --
 214    ------------------------------------------------
 215
 216    function Vector_Vector_Scalar_Elementwise_Operation
 217      (X : X_Vector;
 218       Y : Y_Vector;
 219       Z : Z_Scalar) return Result_Vector
 220    is
 221       R : Result_Vector (X'Range);
 222
 223    begin
 224       if X'Length /= Y'Length then
 225          raise Constraint_Error with
 226             "vectors are of different length in elementwise operation";
 227       end if;
 228
 229       for J in R'Range loop
 230          R (J) := Operation (X (J), Y (J), Z);
 231       end loop;
 232
 233       return R;
 234    end Vector_Vector_Scalar_Elementwise_Operation;
 235
 236    -----------------------------------------
 237    -- Matrix_Scalar_Elementwise_Operation --
 238    -----------------------------------------
 239
 240    function Matrix_Scalar_Elementwise_Operation
 241      (Left  : Left_Matrix;
 242       Right : Right_Scalar) return Result_Matrix
 243    is
 244       R : Result_Matrix (Left'Range (1), Left'Range (2));
 245
 246    begin
 247       for J in R'Range (1) loop
 248          for K in R'Range (2) loop
 249             R (J, K) := Operation (Left (J, K), Right);
 250          end loop;
 251       end loop;
 252
 253       return R;
 254    end Matrix_Scalar_Elementwise_Operation;
 255
 256    -----------------------------------------
 257    -- Vector_Scalar_Elementwise_Operation --
 258    -----------------------------------------
 259
 260    function Vector_Scalar_Elementwise_Operation
 261      (Left  : Left_Vector;
 262       Right : Right_Scalar) return Result_Vector
 263    is
 264       R : Result_Vector (Left'Range);
 265
 266    begin
 267       for J in R'Range loop
 268          R (J) := Operation (Left (J), Right);
 269       end loop;
 270
 271       return R;
 272    end Vector_Scalar_Elementwise_Operation;
 273
 274    -----------------------------------------
 275    -- Scalar_Matrix_Elementwise_Operation --
 276    -----------------------------------------
 277
 278    function Scalar_Matrix_Elementwise_Operation
 279      (Left  : Left_Scalar;
 280       Right : Right_Matrix) return Result_Matrix
 281    is
 282       R : Result_Matrix (Right'Range (1), Right'Range (2));
 283
 284    begin
 285       for J in R'Range (1) loop
 286          for K in R'Range (2) loop
 287             R (J, K) := Operation (Left, Right (J, K));
 288          end loop;
 289       end loop;
 290
 291       return R;
 292    end Scalar_Matrix_Elementwise_Operation;
 293
 294    -----------------------------------------
 295    -- Scalar_Vector_Elementwise_Operation --
 296    -----------------------------------------
 297
 298    function Scalar_Vector_Elementwise_Operation
 299      (Left  : Left_Scalar;
 300       Right : Right_Vector) return Result_Vector
 301    is
 302       R : Result_Vector (Right'Range);
 303
 304    begin
 305       for J in R'Range loop
 306          R (J) := Operation (Left, Right (J));
 307       end loop;
 308
 309       return R;
 310    end Scalar_Vector_Elementwise_Operation;
 311
 312    ---------------------------
 313    -- Matrix_Matrix_Product --
 314    ---------------------------
 315
 316    function Matrix_Matrix_Product
 317      (Left  : Left_Matrix;
 318       Right : Right_Matrix) return Result_Matrix
 319    is
 320       R : Result_Matrix (Left'Range (1), Right'Range (2));
 321
 322    begin
 323       if Left'Length (2) /= Right'Length (1) then
 324          raise Constraint_Error with
 325             "incompatible dimensions in matrix multiplication";
 326       end if;
 327
 328       for J in R'Range (1) loop
 329          for K in R'Range (2) loop
 330             declare
 331                S : Result_Scalar := Zero;
 332             begin
 333                for M in Left'Range (2) loop
 334                   S := S + Left (J, M)
 335                             * Right (M - Left'First (2) + Right'First (1), K);
 336                end loop;
 337
 338                R (J, K) := S;
 339             end;
 340          end loop;
 341       end loop;
 342
 343       return R;
 344    end  Matrix_Matrix_Product;
 345
 346    ---------------------------
 347    -- Matrix_Vector_Product --
 348    ---------------------------
 349
 350    function Matrix_Vector_Product
 351      (Left  : Matrix;
 352       Right : Right_Vector) return Result_Vector
 353    is
 354       R : Result_Vector (Left'Range (1));
 355
 356    begin
 357       if Left'Length (2) /= Right'Length then
 358          raise Constraint_Error with
 359             "incompatible dimensions in matrix-vector multiplication";
 360       end if;
 361
 362       for J in Left'Range (1) loop
 363          declare
 364             S : Result_Scalar := Zero;
 365          begin
 366             for K in Left'Range (2) loop
 367                S := S + Left (J, K) * Right (K - Left'First (2) + Right'First);
 368             end loop;
 369
 370             R (J) := S;
 371          end;
 372       end loop;
 373
 374       return R;
 375    end Matrix_Vector_Product;
 376
 377    -------------------
 378    -- Outer_Product --
 379    -------------------
 380
 381    function Outer_Product
 382      (Left  : Left_Vector;
 383       Right : Right_Vector) return Matrix
 384    is
 385       R : Matrix (Left'Range, Right'Range);
 386
 387    begin
 388       for J in R'Range (1) loop
 389          for K in R'Range (2) loop
 390             R (J, K) := Left (J) * Right (K);
 391          end loop;
 392       end loop;
 393
 394       return R;
 395    end Outer_Product;
 396
 397    ---------------
 398    -- Transpose --
 399    ---------------
 400
 401    procedure Transpose (A : Matrix; R : out Matrix) is
 402    begin
 403       for J in R'Range (1) loop
 404          for K in R'Range (2) loop
 405             R (J, K) := A (J - R'First (1) + A'First (1),
 406                            K - R'First (2) + A'First (2));
 407          end loop;
 408       end loop;
 409    end Transpose;
 410
 411    -------------------------------
 412    -- Update_Matrix_With_Matrix --
 413    -------------------------------
 414
 415    procedure Update_Matrix_With_Matrix (X : in out X_Matrix; Y : Y_Matrix) is
 416    begin
 417       if X'Length (1) /= Y'Length (1)
 418         or else X'Length (2) /= Y'Length (2)
 419       then
 420          raise Constraint_Error with
 421             "matrices are of different dimension in update operation";
 422       end if;
 423
 424       for J in X'Range (1) loop
 425          for K in X'Range (2) loop
 426             Update (X (J, K), Y (J - X'First (1) + Y'First (1),
 427                                  K - X'First (2) + Y'First (2)));
 428          end loop;
 429       end loop;
 430    end Update_Matrix_With_Matrix;
 431
 432    -------------------------------
 433    -- Update_Vector_With_Vector --
 434    -------------------------------
 435
 436    procedure Update_Vector_With_Vector (X : in out X_Vector; Y : Y_Vector) is
 437    begin
 438       if X'Length /= Y'Length then
 439          raise Constraint_Error with
 440             "vectors are of different length in update operation";
 441       end if;
 442
 443       for J in X'Range loop
 444          Update (X (J), Y (J - X'First + Y'First));
 445       end loop;
 446    end Update_Vector_With_Vector;
 447
 448    -----------------
 449    -- Unit_Matrix --
 450    -----------------
 451
 452    function Unit_Matrix
 453      (Order   : Positive;
 454       First_1 : Integer := 1;
 455       First_2 : Integer := 1) return Matrix
 456    is
 457       R : Matrix (First_1 .. Check_Unit_Last (First_1, Order, First_1),
 458                   First_2 .. Check_Unit_Last (First_2, Order, First_2));
 459
 460    begin
 461       R := (others => (others => Zero));
 462
 463       for J in 0 .. Order - 1 loop
 464          R (First_1 + J, First_2 + J) := One;
 465       end loop;
 466
 467       return R;
 468    end Unit_Matrix;
 469
 470    -----------------
 471    -- Unit_Vector --
 472    -----------------
 473
 474    function Unit_Vector
 475      (Index : Integer;
 476       Order : Positive;
 477       First : Integer := 1) return Vector
 478    is
 479       R : Vector (First .. Check_Unit_Last (Index, Order, First));
 480    begin
 481       R := (others => Zero);
 482       R (Index) := One;
 483       return R;
 484    end Unit_Vector;
 485
 486    ---------------------------
 487    -- Vector_Matrix_Product --
 488    ---------------------------
 489
 490    function Vector_Matrix_Product
 491      (Left  : Left_Vector;
 492       Right : Matrix) return Result_Vector
 493    is
 494       R : Result_Vector (Right'Range (2));
 495
 496    begin
 497       if Left'Length /= Right'Length (2) then
 498          raise Constraint_Error with
 499             "incompatible dimensions in vector-matrix multiplication";
 500       end if;
 501
 502       for J in Right'Range (2) loop
 503          declare
 504             S : Result_Scalar := Zero;
 505
 506          begin
 507             for K in Right'Range (1) loop
 508                S := S + Left (J - Right'First (1) + Left'First) * Right (K, J);
 509             end loop;
 510
 511             R (J) := S;
 512          end;
 513       end loop;
 514
 515       return R;
 516    end Vector_Matrix_Product;
 517
 518 end System.Generic_Array_Operations;