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