libgfortran/generated/product_r4.c

   1 /* Implementation of the PRODUCT intrinsic
   2    Copyright (C) 2002-2024 Free Software Foundation, Inc.
   3    Contributed by Paul Brook <paul@nowt.org>
   4
   5 This file is part of the GNU Fortran runtime library (libgfortran).
   6
   7 Libgfortran is free software; you can redistribute it and/or
   8 modify it under the terms of the GNU General Public
   9 License as published by the Free Software Foundation; either
  10 version 3 of the License, or (at your option) any later version.
  11
  12 Libgfortran is distributed in the hope that it will be useful,
  13 but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 GNU General Public License for more details.
  16
  17 Under Section 7 of GPL version 3, you are granted additional
  18 permissions described in the GCC Runtime Library Exception, version
  19 3.1, as published by the Free Software Foundation.
  20
  21 You should have received a copy of the GNU General Public License and
  22 a copy of the GCC Runtime Library Exception along with this program;
  23 see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
  24 <http://www.gnu.org/licenses/>.  */
  25
  26 #include "libgfortran.h"
  27
  28
  29 #if defined (HAVE_GFC_REAL_4) && defined (HAVE_GFC_REAL_4)
  30
  31
  32 extern void product_r4 (gfc_array_r4 * const restrict,
  33         gfc_array_r4 * const restrict, const index_type * const restrict);
  34 export_proto(product_r4);
  35
  36 void
  37 product_r4 (gfc_array_r4 * const restrict retarray,
  38         gfc_array_r4 * const restrict array,
  39         const index_type * const restrict pdim)
  40 {
  41   index_type count[GFC_MAX_DIMENSIONS];
  42   index_type extent[GFC_MAX_DIMENSIONS];
  43   index_type sstride[GFC_MAX_DIMENSIONS];
  44   index_type dstride[GFC_MAX_DIMENSIONS];
  45   const GFC_REAL_4 * restrict base;
  46   GFC_REAL_4 * restrict dest;
  47   index_type rank;
  48   index_type n;
  49   index_type len;
  50   index_type delta;
  51   index_type dim;
  52   int continue_loop;
  53
  54   /* Make dim zero based to avoid confusion.  */
  55   rank = GFC_DESCRIPTOR_RANK (array) - 1;
  56   dim = (*pdim) - 1;
  57
  58   if (unlikely (dim < 0 || dim > rank))
  59     {
  60       runtime_error ("Dim argument incorrect in PRODUCT intrinsic: "
  61                      "is %ld, should be between 1 and %ld",
  62                      (long int) dim + 1, (long int) rank + 1);
  63     }
  64
  65   len = GFC_DESCRIPTOR_EXTENT(array,dim);
  66   if (len < 0)
  67     len = 0;
  68   delta = GFC_DESCRIPTOR_STRIDE(array,dim);
  69
  70   for (n = 0; n < dim; n++)
  71     {
  72       sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n);
  73       extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);
  74
  75       if (extent[n] < 0)
  76         extent[n] = 0;
  77     }
  78   for (n = dim; n < rank; n++)
  79     {
  80       sstride[n] = GFC_DESCRIPTOR_STRIDE(array, n + 1);
  81       extent[n] = GFC_DESCRIPTOR_EXTENT(array, n + 1);
  82
  83       if (extent[n] < 0)
  84         extent[n] = 0;
  85     }
  86
  87   if (retarray->base_addr == NULL)
  88     {
  89       size_t alloc_size, str;
  90
  91       for (n = 0; n < rank; n++)
  92         {
  93           if (n == 0)
  94             str = 1;
  95           else
  96             str = GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1];
  97
  98           GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str);
  99
 100         }
 101
 102       retarray->offset = 0;
 103       retarray->dtype.rank = rank;
 104
 105       alloc_size = GFC_DESCRIPTOR_STRIDE(retarray,rank-1) * extent[rank-1];
 106
 107       retarray->base_addr = xmallocarray (alloc_size, sizeof (GFC_REAL_4));
 108       if (alloc_size == 0)
 109         return;
 110     }
 111   else
 112     {
 113       if (rank != GFC_DESCRIPTOR_RANK (retarray))
 114         runtime_error ("rank of return array incorrect in"
 115                        " PRODUCT intrinsic: is %ld, should be %ld",
 116                        (long int) (GFC_DESCRIPTOR_RANK (retarray)),
 117                        (long int) rank);
 118
 119       if (unlikely (compile_options.bounds_check))
 120         bounds_ifunction_return ((array_t *) retarray, extent,
 121                                  "return value", "PRODUCT");
 122     }
 123
 124   for (n = 0; n < rank; n++)
 125     {
 126       count[n] = 0;
 127       dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);
 128       if (extent[n] <= 0)
 129         return;
 130     }
 131
 132   base = array->base_addr;
 133   dest = retarray->base_addr;
 134
 135   continue_loop = 1;
 136   while (continue_loop)
 137     {
 138       const GFC_REAL_4 * restrict src;
 139       GFC_REAL_4 result;
 140       src = base;
 141       {
 142
 143   result = 1;
 144         if (len <= 0)
 145           *dest = 1;
 146         else
 147           {
 148 #if ! defined HAVE_BACK_ARG
 149             for (n = 0; n < len; n++, src += delta)
 150               {
 151 #endif
 152
 153   result *= *src;
 154               }
 155
 156             *dest = result;
 157           }
 158       }
 159       /* Advance to the next element.  */
 160       count[0]++;
 161       base += sstride[0];
 162       dest += dstride[0];
 163       n = 0;
 164       while (count[n] == extent[n])
 165         {
 166           /* When we get to the end of a dimension, reset it and increment
 167              the next dimension.  */
 168           count[n] = 0;
 169           /* We could precalculate these products, but this is a less
 170              frequently used path so probably not worth it.  */
 171           base -= sstride[n] * extent[n];
 172           dest -= dstride[n] * extent[n];
 173           n++;
 174           if (n >= rank)
 175             {
 176               /* Break out of the loop.  */
 177               continue_loop = 0;
 178               break;
 179             }
 180           else
 181             {
 182               count[n]++;
 183               base += sstride[n];
 184               dest += dstride[n];
 185             }
 186         }
 187     }
 188 }
 189
 190
 191 extern void mproduct_r4 (gfc_array_r4 * const restrict,
 192         gfc_array_r4 * const restrict, const index_type * const restrict,
 193         gfc_array_l1 * const restrict);
 194 export_proto(mproduct_r4);
 195
 196 void
 197 mproduct_r4 (gfc_array_r4 * const restrict retarray,
 198         gfc_array_r4 * const restrict array,
 199         const index_type * const restrict pdim,
 200         gfc_array_l1 * const restrict mask)
 201 {
 202   index_type count[GFC_MAX_DIMENSIONS];
 203   index_type extent[GFC_MAX_DIMENSIONS];
 204   index_type sstride[GFC_MAX_DIMENSIONS];
 205   index_type dstride[GFC_MAX_DIMENSIONS];
 206   index_type mstride[GFC_MAX_DIMENSIONS];
 207   GFC_REAL_4 * restrict dest;
 208   const GFC_REAL_4 * restrict base;
 209   const GFC_LOGICAL_1 * restrict mbase;
 210   index_type rank;
 211   index_type dim;
 212   index_type n;
 213   index_type len;
 214   index_type delta;
 215   index_type mdelta;
 216   int mask_kind;
 217
 218   if (mask == NULL)
 219     {
 220 #ifdef HAVE_BACK_ARG
 221       product_r4 (retarray, array, pdim, back);
 222 #else
 223       product_r4 (retarray, array, pdim);
 224 #endif
 225       return;
 226     }
 227
 228   dim = (*pdim) - 1;
 229   rank = GFC_DESCRIPTOR_RANK (array) - 1;
 230
 231
 232   if (unlikely (dim < 0 || dim > rank))
 233     {
 234       runtime_error ("Dim argument incorrect in PRODUCT intrinsic: "
 235                      "is %ld, should be between 1 and %ld",
 236                      (long int) dim + 1, (long int) rank + 1);
 237     }
 238
 239   len = GFC_DESCRIPTOR_EXTENT(array,dim);
 240   if (len < 0)
 241     len = 0;
 242
 243   mbase = mask->base_addr;
 244
 245   mask_kind = GFC_DESCRIPTOR_SIZE (mask);
 246
 247   if (mask_kind == 1 || mask_kind == 2 || mask_kind == 4 || mask_kind == 8
 248 #ifdef HAVE_GFC_LOGICAL_16
 249       || mask_kind == 16
 250 #endif
 251       )
 252     mbase = GFOR_POINTER_TO_L1 (mbase, mask_kind);
 253   else
 254     runtime_error ("Funny sized logical array");
 255
 256   delta = GFC_DESCRIPTOR_STRIDE(array,dim);
 257   mdelta = GFC_DESCRIPTOR_STRIDE_BYTES(mask,dim);
 258
 259   for (n = 0; n < dim; n++)
 260     {
 261       sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n);
 262       mstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(mask,n);
 263       extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);
 264
 265       if (extent[n] < 0)
 266         extent[n] = 0;
 267
 268     }
 269   for (n = dim; n < rank; n++)
 270     {
 271       sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n + 1);
 272       mstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(mask, n + 1);
 273       extent[n] = GFC_DESCRIPTOR_EXTENT(array, n + 1);
 274
 275       if (extent[n] < 0)
 276         extent[n] = 0;
 277     }
 278
 279   if (retarray->base_addr == NULL)
 280     {
 281       size_t alloc_size, str;
 282
 283       for (n = 0; n < rank; n++)
 284         {
 285           if (n == 0)
 286             str = 1;
 287           else
 288             str= GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1];
 289
 290           GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str);
 291
 292         }
 293
 294       alloc_size = GFC_DESCRIPTOR_STRIDE(retarray,rank-1) * extent[rank-1];
 295
 296       retarray->offset = 0;
 297       retarray->dtype.rank = rank;
 298
 299       retarray->base_addr = xmallocarray (alloc_size, sizeof (GFC_REAL_4));
 300       if (alloc_size == 0)
 301         return;
 302     }
 303   else
 304     {
 305       if (rank != GFC_DESCRIPTOR_RANK (retarray))
 306         runtime_error ("rank of return array incorrect in PRODUCT intrinsic");
 307
 308       if (unlikely (compile_options.bounds_check))
 309         {
 310           bounds_ifunction_return ((array_t *) retarray, extent,
 311                                    "return value", "PRODUCT");
 312           bounds_equal_extents ((array_t *) mask, (array_t *) array,
 313                                 "MASK argument", "PRODUCT");
 314         }
 315     }
 316
 317   for (n = 0; n < rank; n++)
 318     {
 319       count[n] = 0;
 320       dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);
 321       if (extent[n] <= 0)
 322         return;
 323     }
 324
 325   dest = retarray->base_addr;
 326   base = array->base_addr;
 327
 328   while (base)
 329     {
 330       const GFC_REAL_4 * restrict src;
 331       const GFC_LOGICAL_1 * restrict msrc;
 332       GFC_REAL_4 result;
 333       src = base;
 334       msrc = mbase;
 335       {
 336
 337   result = 1;
 338         for (n = 0; n < len; n++, src += delta, msrc += mdelta)
 339           {
 340
 341   if (*msrc)
 342     result *= *src;
 343           }
 344         *dest = result;
 345       }
 346       /* Advance to the next element.  */
 347       count[0]++;
 348       base += sstride[0];
 349       mbase += mstride[0];
 350       dest += dstride[0];
 351       n = 0;
 352       while (count[n] == extent[n])
 353         {
 354           /* When we get to the end of a dimension, reset it and increment
 355              the next dimension.  */
 356           count[n] = 0;
 357           /* We could precalculate these products, but this is a less
 358              frequently used path so probably not worth it.  */
 359           base -= sstride[n] * extent[n];
 360           mbase -= mstride[n] * extent[n];
 361           dest -= dstride[n] * extent[n];
 362           n++;
 363           if (n >= rank)
 364             {
 365               /* Break out of the loop.  */
 366               base = NULL;
 367               break;
 368             }
 369           else
 370             {
 371               count[n]++;
 372               base += sstride[n];
 373               mbase += mstride[n];
 374               dest += dstride[n];
 375             }
 376         }
 377     }
 378 }
 379
 380
 381 extern void sproduct_r4 (gfc_array_r4 * const restrict,
 382         gfc_array_r4 * const restrict, const index_type * const restrict,
 383         GFC_LOGICAL_4 *);
 384 export_proto(sproduct_r4);
 385
 386 void
 387 sproduct_r4 (gfc_array_r4 * const restrict retarray,
 388         gfc_array_r4 * const restrict array,
 389         const index_type * const restrict pdim,
 390         GFC_LOGICAL_4 * mask)
 391 {
 392   index_type count[GFC_MAX_DIMENSIONS];
 393   index_type extent[GFC_MAX_DIMENSIONS];
 394   index_type dstride[GFC_MAX_DIMENSIONS];
 395   GFC_REAL_4 * restrict dest;
 396   index_type rank;
 397   index_type n;
 398   index_type dim;
 399
 400
 401   if (mask == NULL || *mask)
 402     {
 403 #ifdef HAVE_BACK_ARG
 404       product_r4 (retarray, array, pdim, back);
 405 #else
 406       product_r4 (retarray, array, pdim);
 407 #endif
 408       return;
 409     }
 410   /* Make dim zero based to avoid confusion.  */
 411   dim = (*pdim) - 1;
 412   rank = GFC_DESCRIPTOR_RANK (array) - 1;
 413
 414   if (unlikely (dim < 0 || dim > rank))
 415     {
 416       runtime_error ("Dim argument incorrect in PRODUCT intrinsic: "
 417                      "is %ld, should be between 1 and %ld",
 418                      (long int) dim + 1, (long int) rank + 1);
 419     }
 420
 421   for (n = 0; n < dim; n++)
 422     {
 423       extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);
 424
 425       if (extent[n] <= 0)
 426         extent[n] = 0;
 427     }
 428
 429   for (n = dim; n < rank; n++)
 430     {
 431       extent[n] =
 432         GFC_DESCRIPTOR_EXTENT(array,n + 1);
 433
 434       if (extent[n] <= 0)
 435         extent[n] = 0;
 436     }
 437
 438   if (retarray->base_addr == NULL)
 439     {
 440       size_t alloc_size, str;
 441
 442       for (n = 0; n < rank; n++)
 443         {
 444           if (n == 0)
 445             str = 1;
 446           else
 447             str = GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1];
 448
 449           GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str);
 450
 451         }
 452
 453       retarray->offset = 0;
 454       retarray->dtype.rank = rank;
 455
 456       alloc_size = GFC_DESCRIPTOR_STRIDE(retarray,rank-1) * extent[rank-1];
 457
 458       retarray->base_addr = xmallocarray (alloc_size, sizeof (GFC_REAL_4));
 459       if (alloc_size == 0)
 460         return;
 461     }
 462   else
 463     {
 464       if (rank != GFC_DESCRIPTOR_RANK (retarray))
 465         runtime_error ("rank of return array incorrect in"
 466                        " PRODUCT intrinsic: is %ld, should be %ld",
 467                        (long int) (GFC_DESCRIPTOR_RANK (retarray)),
 468                        (long int) rank);
 469
 470       if (unlikely (compile_options.bounds_check))
 471         {
 472           for (n=0; n < rank; n++)
 473             {
 474               index_type ret_extent;
 475
 476               ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,n);
 477               if (extent[n] != ret_extent)
 478                 runtime_error ("Incorrect extent in return value of"
 479                                " PRODUCT intrinsic in dimension %ld:"
 480                                " is %ld, should be %ld", (long int) n + 1,
 481                                (long int) ret_extent, (long int) extent[n]);
 482             }
 483         }
 484     }
 485
 486   for (n = 0; n < rank; n++)
 487     {
 488       count[n] = 0;
 489       dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);
 490     }
 491
 492   dest = retarray->base_addr;
 493
 494   while(1)
 495     {
 496       *dest = 1;
 497       count[0]++;
 498       dest += dstride[0];
 499       n = 0;
 500       while (count[n] == extent[n])
 501         {
 502           /* When we get to the end of a dimension, reset it and increment
 503              the next dimension.  */
 504           count[n] = 0;
 505           /* We could precalculate these products, but this is a less
 506              frequently used path so probably not worth it.  */
 507           dest -= dstride[n] * extent[n];
 508           n++;
 509           if (n >= rank)
 510             return;
 511           else
 512             {
 513               count[n]++;
 514               dest += dstride[n];
 515             }
 516         }
 517     }
 518 }
 519
 520 #endif