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-2007, 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) :=
 155               Operation
 156                 (Left (J, K),
 157                  Right
 158                    (J - R'First (1) + Right'First (1),
 159                     K - R'First (2) + Right'First (2)));
 160          end loop;
 161       end loop;
 162
 163       return R;
 164    end Matrix_Matrix_Elementwise_Operation;
 165
 166    ------------------------------------------------
 167    -- Matrix_Matrix_Scalar_Elementwise_Operation --
 168    ------------------------------------------------
 169
 170    function Matrix_Matrix_Scalar_Elementwise_Operation
 171      (X    : X_Matrix;
 172       Y    : Y_Matrix;
 173       Z    : Z_Scalar) return Result_Matrix
 174    is
 175       R : Result_Matrix (X'Range (1), X'Range (2));
 176
 177    begin
 178       if X'Length (1) /= Y'Length (1)
 179         or else X'Length (2) /= Y'Length (2)
 180       then
 181          raise Constraint_Error with
 182             "matrices are of different dimension in elementwise operation";
 183       end if;
 184
 185       for J in R'Range (1) loop
 186          for K in R'Range (2) loop
 187             R (J, K) :=
 188               Operation
 189                 (X (J, K),
 190                  Y (J - R'First (1) + Y'First (1),
 191                     K - R'First (2) + Y'First (2)),
 192                  Z);
 193          end loop;
 194       end loop;
 195
 196       return R;
 197    end Matrix_Matrix_Scalar_Elementwise_Operation;
 198
 199    -----------------------------------------
 200    -- Vector_Vector_Elementwise_Operation --
 201    -----------------------------------------
 202
 203    function Vector_Vector_Elementwise_Operation
 204      (Left  : Left_Vector;
 205       Right : Right_Vector) return Result_Vector
 206    is
 207       R : Result_Vector (Left'Range);
 208
 209    begin
 210       if Left'Length /= Right'Length then
 211          raise Constraint_Error with
 212             "vectors are of different length in elementwise operation";
 213       end if;
 214
 215       for J in R'Range loop
 216          R (J) := Operation (Left (J), Right (J - R'First + Right'First));
 217       end loop;
 218
 219       return R;
 220    end Vector_Vector_Elementwise_Operation;
 221
 222    ------------------------------------------------
 223    -- Vector_Vector_Scalar_Elementwise_Operation --
 224    ------------------------------------------------
 225
 226    function Vector_Vector_Scalar_Elementwise_Operation
 227      (X : X_Vector;
 228       Y : Y_Vector;
 229       Z : Z_Scalar) return Result_Vector
 230    is
 231       R : Result_Vector (X'Range);
 232
 233    begin
 234       if X'Length /= Y'Length then
 235          raise Constraint_Error with
 236             "vectors are of different length in elementwise operation";
 237       end if;
 238
 239       for J in R'Range loop
 240          R (J) := Operation (X (J), Y (J - X'First + Y'First), Z);
 241       end loop;
 242
 243       return R;
 244    end Vector_Vector_Scalar_Elementwise_Operation;
 245
 246    -----------------------------------------
 247    -- Matrix_Scalar_Elementwise_Operation --
 248    -----------------------------------------
 249
 250    function Matrix_Scalar_Elementwise_Operation
 251      (Left  : Left_Matrix;
 252       Right : Right_Scalar) return Result_Matrix
 253    is
 254       R : Result_Matrix (Left'Range (1), Left'Range (2));
 255
 256    begin
 257       for J in R'Range (1) loop
 258          for K in R'Range (2) loop
 259             R (J, K) := Operation (Left (J, K), Right);
 260          end loop;
 261       end loop;
 262
 263       return R;
 264    end Matrix_Scalar_Elementwise_Operation;
 265
 266    -----------------------------------------
 267    -- Vector_Scalar_Elementwise_Operation --
 268    -----------------------------------------
 269
 270    function Vector_Scalar_Elementwise_Operation
 271      (Left  : Left_Vector;
 272       Right : Right_Scalar) return Result_Vector
 273    is
 274       R : Result_Vector (Left'Range);
 275
 276    begin
 277       for J in R'Range loop
 278          R (J) := Operation (Left (J), Right);
 279       end loop;
 280
 281       return R;
 282    end Vector_Scalar_Elementwise_Operation;
 283
 284    -----------------------------------------
 285    -- Scalar_Matrix_Elementwise_Operation --
 286    -----------------------------------------
 287
 288    function Scalar_Matrix_Elementwise_Operation
 289      (Left  : Left_Scalar;
 290       Right : Right_Matrix) return Result_Matrix
 291    is
 292       R : Result_Matrix (Right'Range (1), Right'Range (2));
 293
 294    begin
 295       for J in R'Range (1) loop
 296          for K in R'Range (2) loop
 297             R (J, K) := Operation (Left, Right (J, K));
 298          end loop;
 299       end loop;
 300
 301       return R;
 302    end Scalar_Matrix_Elementwise_Operation;
 303
 304    -----------------------------------------
 305    -- Scalar_Vector_Elementwise_Operation --
 306    -----------------------------------------
 307
 308    function Scalar_Vector_Elementwise_Operation
 309      (Left  : Left_Scalar;
 310       Right : Right_Vector) return Result_Vector
 311    is
 312       R : Result_Vector (Right'Range);
 313
 314    begin
 315       for J in R'Range loop
 316          R (J) := Operation (Left, Right (J));
 317       end loop;
 318
 319       return R;
 320    end Scalar_Vector_Elementwise_Operation;
 321
 322    ---------------------------
 323    -- Matrix_Matrix_Product --
 324    ---------------------------
 325
 326    function Matrix_Matrix_Product
 327      (Left  : Left_Matrix;
 328       Right : Right_Matrix) return Result_Matrix
 329    is
 330       R : Result_Matrix (Left'Range (1), Right'Range (2));
 331
 332    begin
 333       if Left'Length (2) /= Right'Length (1) then
 334          raise Constraint_Error with
 335             "incompatible dimensions in matrix multiplication";
 336       end if;
 337
 338       for J in R'Range (1) loop
 339          for K in R'Range (2) loop
 340             declare
 341                S : Result_Scalar := Zero;
 342             begin
 343                for M in Left'Range (2) loop
 344                   S := S + Left (J, M)
 345                             * Right (M - Left'First (2) + Right'First (1), K);
 346                end loop;
 347
 348                R (J, K) := S;
 349             end;
 350          end loop;
 351       end loop;
 352
 353       return R;
 354    end  Matrix_Matrix_Product;
 355
 356    ---------------------------
 357    -- Matrix_Vector_Product --
 358    ---------------------------
 359
 360    function Matrix_Vector_Product
 361      (Left  : Matrix;
 362       Right : Right_Vector) return Result_Vector
 363    is
 364       R : Result_Vector (Left'Range (1));
 365
 366    begin
 367       if Left'Length (2) /= Right'Length then
 368          raise Constraint_Error with
 369             "incompatible dimensions in matrix-vector multiplication";
 370       end if;
 371
 372       for J in Left'Range (1) loop
 373          declare
 374             S : Result_Scalar := Zero;
 375          begin
 376             for K in Left'Range (2) loop
 377                S := S + Left (J, K) * Right (K - Left'First (2) + Right'First);
 378             end loop;
 379
 380             R (J) := S;
 381          end;
 382       end loop;
 383
 384       return R;
 385    end Matrix_Vector_Product;
 386
 387    -------------------
 388    -- Outer_Product --
 389    -------------------
 390
 391    function Outer_Product
 392      (Left  : Left_Vector;
 393       Right : Right_Vector) return Matrix
 394    is
 395       R : Matrix (Left'Range, Right'Range);
 396
 397    begin
 398       for J in R'Range (1) loop
 399          for K in R'Range (2) loop
 400             R (J, K) := Left (J) * Right (K);
 401          end loop;
 402       end loop;
 403
 404       return R;
 405    end Outer_Product;
 406
 407    ---------------
 408    -- Transpose --
 409    ---------------
 410
 411    procedure Transpose (A : Matrix; R : out Matrix) is
 412    begin
 413       for J in R'Range (1) loop
 414          for K in R'Range (2) loop
 415             R (J, K) := A (K - R'First (2) + A'First (1),
 416                            J - R'First (1) + A'First (2));
 417          end loop;
 418       end loop;
 419    end Transpose;
 420
 421    -------------------------------
 422    -- Update_Matrix_With_Matrix --
 423    -------------------------------
 424
 425    procedure Update_Matrix_With_Matrix (X : in out X_Matrix; Y : Y_Matrix) is
 426    begin
 427       if X'Length (1) /= Y'Length (1)
 428         or else X'Length (2) /= Y'Length (2)
 429       then
 430          raise Constraint_Error with
 431             "matrices are of different dimension in update operation";
 432       end if;
 433
 434       for J in X'Range (1) loop
 435          for K in X'Range (2) loop
 436             Update (X (J, K), Y (J - X'First (1) + Y'First (1),
 437                                  K - X'First (2) + Y'First (2)));
 438          end loop;
 439       end loop;
 440    end Update_Matrix_With_Matrix;
 441
 442    -------------------------------
 443    -- Update_Vector_With_Vector --
 444    -------------------------------
 445
 446    procedure Update_Vector_With_Vector (X : in out X_Vector; Y : Y_Vector) is
 447    begin
 448       if X'Length /= Y'Length then
 449          raise Constraint_Error with
 450             "vectors are of different length in update operation";
 451       end if;
 452
 453       for J in X'Range loop
 454          Update (X (J), Y (J - X'First + Y'First));
 455       end loop;
 456    end Update_Vector_With_Vector;
 457
 458    -----------------
 459    -- Unit_Matrix --
 460    -----------------
 461
 462    function Unit_Matrix
 463      (Order   : Positive;
 464       First_1 : Integer := 1;
 465       First_2 : Integer := 1) return Matrix
 466    is
 467       R : Matrix (First_1 .. Check_Unit_Last (First_1, Order, First_1),
 468                   First_2 .. Check_Unit_Last (First_2, Order, First_2));
 469
 470    begin
 471       R := (others => (others => Zero));
 472
 473       for J in 0 .. Order - 1 loop
 474          R (First_1 + J, First_2 + J) := One;
 475       end loop;
 476
 477       return R;
 478    end Unit_Matrix;
 479
 480    -----------------
 481    -- Unit_Vector --
 482    -----------------
 483
 484    function Unit_Vector
 485      (Index : Integer;
 486       Order : Positive;
 487       First : Integer := 1) return Vector
 488    is
 489       R : Vector (First .. Check_Unit_Last (Index, Order, First));
 490    begin
 491       R := (others => Zero);
 492       R (Index) := One;
 493       return R;
 494    end Unit_Vector;
 495
 496    ---------------------------
 497    -- Vector_Matrix_Product --
 498    ---------------------------
 499
 500    function Vector_Matrix_Product
 501      (Left  : Left_Vector;
 502       Right : Matrix) return Result_Vector
 503    is
 504       R : Result_Vector (Right'Range (2));
 505
 506    begin
 507       if Left'Length /= Right'Length (2) then
 508          raise Constraint_Error with
 509             "incompatible dimensions in vector-matrix multiplication";
 510       end if;
 511
 512       for J in Right'Range (2) loop
 513          declare
 514             S : Result_Scalar := Zero;
 515
 516          begin
 517             for K in Right'Range (1) loop
 518                S := S + Left (J - Right'First (1) + Left'First) * Right (K, J);
 519             end loop;
 520
 521             R (J) := S;
 522          end;
 523       end loop;
 524
 525       return R;
 526    end Vector_Matrix_Product;
 527
 528 end System.Generic_Array_Operations;