libgfortran/generated/sum_r8.c

   1 /* Implementation of the SUM intrinsic
   2    Copyright (C) 2002-2019 Free Software Foundation, Inc.
   3    Contributed by Paul Brook <paul@nowt.org>
   4
   5 This file is part of the GNU Fortran 95 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_8) && defined (HAVE_GFC_REAL_8)
  30
  31
  32 extern void sum_r8 (gfc_array_r8 * const restrict,
  33         gfc_array_r8 * const restrict, const index_type * const restrict);
  34 export_proto(sum_r8);
  35
  36 void
  37 sum_r8 (gfc_array_r8 * const restrict retarray,
  38         gfc_array_r8 * 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_8 * restrict base;
  46   GFC_REAL_8 * 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 SUM 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_8));
 108       if (alloc_size == 0)
 109         {
 110           /* Make sure we have a zero-sized array.  */
 111           GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 1);
 112           return;
 113
 114         }
 115     }
 116   else
 117     {
 118       if (rank != GFC_DESCRIPTOR_RANK (retarray))
 119         runtime_error ("rank of return array incorrect in"
 120                        " SUM intrinsic: is %ld, should be %ld",
 121                        (long int) (GFC_DESCRIPTOR_RANK (retarray)),
 122                        (long int) rank);
 123
 124       if (unlikely (compile_options.bounds_check))
 125         bounds_ifunction_return ((array_t *) retarray, extent,
 126                                  "return value", "SUM");
 127     }
 128
 129   for (n = 0; n < rank; n++)
 130     {
 131       count[n] = 0;
 132       dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);
 133       if (extent[n] <= 0)
 134         return;
 135     }
 136
 137   base = array->base_addr;
 138   dest = retarray->base_addr;
 139
 140   continue_loop = 1;
 141   while (continue_loop)
 142     {
 143       const GFC_REAL_8 * restrict src;
 144       GFC_REAL_8 result;
 145       src = base;
 146       {
 147
 148   result = 0;
 149         if (len <= 0)
 150           *dest = 0;
 151         else
 152           {
 153 #if ! defined HAVE_BACK_ARG
 154             for (n = 0; n < len; n++, src += delta)
 155               {
 156 #endif
 157
 158   result += *src;
 159               }
 160
 161             *dest = result;
 162           }
 163       }
 164       /* Advance to the next element.  */
 165       count[0]++;
 166       base += sstride[0];
 167       dest += dstride[0];
 168       n = 0;
 169       while (count[n] == extent[n])
 170         {
 171           /* When we get to the end of a dimension, reset it and increment
 172              the next dimension.  */
 173           count[n] = 0;
 174           /* We could precalculate these products, but this is a less
 175              frequently used path so probably not worth it.  */
 176           base -= sstride[n] * extent[n];
 177           dest -= dstride[n] * extent[n];
 178           n++;
 179           if (n >= rank)
 180             {
 181               /* Break out of the loop.  */
 182               continue_loop = 0;
 183               break;
 184             }
 185           else
 186             {
 187               count[n]++;
 188               base += sstride[n];
 189               dest += dstride[n];
 190             }
 191         }
 192     }
 193 }
 194
 195
 196 extern void msum_r8 (gfc_array_r8 * const restrict,
 197         gfc_array_r8 * const restrict, const index_type * const restrict,
 198         gfc_array_l1 * const restrict);
 199 export_proto(msum_r8);
 200
 201 void
 202 msum_r8 (gfc_array_r8 * const restrict retarray,
 203         gfc_array_r8 * const restrict array,
 204         const index_type * const restrict pdim,
 205         gfc_array_l1 * const restrict mask)
 206 {
 207   index_type count[GFC_MAX_DIMENSIONS];
 208   index_type extent[GFC_MAX_DIMENSIONS];
 209   index_type sstride[GFC_MAX_DIMENSIONS];
 210   index_type dstride[GFC_MAX_DIMENSIONS];
 211   index_type mstride[GFC_MAX_DIMENSIONS];
 212   GFC_REAL_8 * restrict dest;
 213   const GFC_REAL_8 * restrict base;
 214   const GFC_LOGICAL_1 * restrict mbase;
 215   index_type rank;
 216   index_type dim;
 217   index_type n;
 218   index_type len;
 219   index_type delta;
 220   index_type mdelta;
 221   int mask_kind;
 222
 223   if (mask == NULL)
 224     {
 225 #ifdef HAVE_BACK_ARG
 226       sum_r8 (retarray, array, pdim, back);
 227 #else
 228       sum_r8 (retarray, array, pdim);
 229 #endif
 230       return;
 231     }
 232
 233   dim = (*pdim) - 1;
 234   rank = GFC_DESCRIPTOR_RANK (array) - 1;
 235
 236
 237   if (unlikely (dim < 0 || dim > rank))
 238     {
 239       runtime_error ("Dim argument incorrect in SUM intrinsic: "
 240                      "is %ld, should be between 1 and %ld",
 241                      (long int) dim + 1, (long int) rank + 1);
 242     }
 243
 244   len = GFC_DESCRIPTOR_EXTENT(array,dim);
 245   if (len <= 0)
 246     return;
 247
 248   mbase = mask->base_addr;
 249
 250   mask_kind = GFC_DESCRIPTOR_SIZE (mask);
 251
 252   if (mask_kind == 1 || mask_kind == 2 || mask_kind == 4 || mask_kind == 8
 253 #ifdef HAVE_GFC_LOGICAL_16
 254       || mask_kind == 16
 255 #endif
 256       )
 257     mbase = GFOR_POINTER_TO_L1 (mbase, mask_kind);
 258   else
 259     runtime_error ("Funny sized logical array");
 260
 261   delta = GFC_DESCRIPTOR_STRIDE(array,dim);
 262   mdelta = GFC_DESCRIPTOR_STRIDE_BYTES(mask,dim);
 263
 264   for (n = 0; n < dim; n++)
 265     {
 266       sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n);
 267       mstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(mask,n);
 268       extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);
 269
 270       if (extent[n] < 0)
 271         extent[n] = 0;
 272
 273     }
 274   for (n = dim; n < rank; n++)
 275     {
 276       sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n + 1);
 277       mstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(mask, n + 1);
 278       extent[n] = GFC_DESCRIPTOR_EXTENT(array, n + 1);
 279
 280       if (extent[n] < 0)
 281         extent[n] = 0;
 282     }
 283
 284   if (retarray->base_addr == NULL)
 285     {
 286       size_t alloc_size, str;
 287
 288       for (n = 0; n < rank; n++)
 289         {
 290           if (n == 0)
 291             str = 1;
 292           else
 293             str= GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1];
 294
 295           GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str);
 296
 297         }
 298
 299       alloc_size = GFC_DESCRIPTOR_STRIDE(retarray,rank-1) * extent[rank-1];
 300
 301       retarray->offset = 0;
 302       retarray->dtype.rank = rank;
 303
 304       if (alloc_size == 0)
 305         {
 306           /* Make sure we have a zero-sized array.  */
 307           GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 1);
 308           return;
 309         }
 310       else
 311         retarray->base_addr = xmallocarray (alloc_size, sizeof (GFC_REAL_8));
 312
 313     }
 314   else
 315     {
 316       if (rank != GFC_DESCRIPTOR_RANK (retarray))
 317         runtime_error ("rank of return array incorrect in SUM intrinsic");
 318
 319       if (unlikely (compile_options.bounds_check))
 320         {
 321           bounds_ifunction_return ((array_t *) retarray, extent,
 322                                    "return value", "SUM");
 323           bounds_equal_extents ((array_t *) mask, (array_t *) array,
 324                                 "MASK argument", "SUM");
 325         }
 326     }
 327
 328   for (n = 0; n < rank; n++)
 329     {
 330       count[n] = 0;
 331       dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);
 332       if (extent[n] <= 0)
 333         return;
 334     }
 335
 336   dest = retarray->base_addr;
 337   base = array->base_addr;
 338
 339   while (base)
 340     {
 341       const GFC_REAL_8 * restrict src;
 342       const GFC_LOGICAL_1 * restrict msrc;
 343       GFC_REAL_8 result;
 344       src = base;
 345       msrc = mbase;
 346       {
 347
 348   result = 0;
 349         for (n = 0; n < len; n++, src += delta, msrc += mdelta)
 350           {
 351
 352   if (*msrc)
 353     result += *src;
 354           }
 355         *dest = result;
 356       }
 357       /* Advance to the next element.  */
 358       count[0]++;
 359       base += sstride[0];
 360       mbase += mstride[0];
 361       dest += dstride[0];
 362       n = 0;
 363       while (count[n] == extent[n])
 364         {
 365           /* When we get to the end of a dimension, reset it and increment
 366              the next dimension.  */
 367           count[n] = 0;
 368           /* We could precalculate these products, but this is a less
 369              frequently used path so probably not worth it.  */
 370           base -= sstride[n] * extent[n];
 371           mbase -= mstride[n] * extent[n];
 372           dest -= dstride[n] * extent[n];
 373           n++;
 374           if (n >= rank)
 375             {
 376               /* Break out of the loop.  */
 377               base = NULL;
 378               break;
 379             }
 380           else
 381             {
 382               count[n]++;
 383               base += sstride[n];
 384               mbase += mstride[n];
 385               dest += dstride[n];
 386             }
 387         }
 388     }
 389 }
 390
 391
 392 extern void ssum_r8 (gfc_array_r8 * const restrict,
 393         gfc_array_r8 * const restrict, const index_type * const restrict,
 394         GFC_LOGICAL_4 *);
 395 export_proto(ssum_r8);
 396
 397 void
 398 ssum_r8 (gfc_array_r8 * const restrict retarray,
 399         gfc_array_r8 * const restrict array,
 400         const index_type * const restrict pdim,
 401         GFC_LOGICAL_4 * mask)
 402 {
 403   index_type count[GFC_MAX_DIMENSIONS];
 404   index_type extent[GFC_MAX_DIMENSIONS];
 405   index_type dstride[GFC_MAX_DIMENSIONS];
 406   GFC_REAL_8 * restrict dest;
 407   index_type rank;
 408   index_type n;
 409   index_type dim;
 410
 411
 412   if (mask == NULL || *mask)
 413     {
 414 #ifdef HAVE_BACK_ARG
 415       sum_r8 (retarray, array, pdim, back);
 416 #else
 417       sum_r8 (retarray, array, pdim);
 418 #endif
 419       return;
 420     }
 421   /* Make dim zero based to avoid confusion.  */
 422   dim = (*pdim) - 1;
 423   rank = GFC_DESCRIPTOR_RANK (array) - 1;
 424
 425   if (unlikely (dim < 0 || dim > rank))
 426     {
 427       runtime_error ("Dim argument incorrect in SUM intrinsic: "
 428                      "is %ld, should be between 1 and %ld",
 429                      (long int) dim + 1, (long int) rank + 1);
 430     }
 431
 432   for (n = 0; n < dim; n++)
 433     {
 434       extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);
 435
 436       if (extent[n] <= 0)
 437         extent[n] = 0;
 438     }
 439
 440   for (n = dim; n < rank; n++)
 441     {
 442       extent[n] =
 443         GFC_DESCRIPTOR_EXTENT(array,n + 1);
 444
 445       if (extent[n] <= 0)
 446         extent[n] = 0;
 447     }
 448
 449   if (retarray->base_addr == NULL)
 450     {
 451       size_t alloc_size, str;
 452
 453       for (n = 0; n < rank; n++)
 454         {
 455           if (n == 0)
 456             str = 1;
 457           else
 458             str = GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1];
 459
 460           GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str);
 461
 462         }
 463
 464       retarray->offset = 0;
 465       retarray->dtype.rank = rank;
 466
 467       alloc_size = GFC_DESCRIPTOR_STRIDE(retarray,rank-1) * extent[rank-1];
 468
 469       if (alloc_size == 0)
 470         {
 471           /* Make sure we have a zero-sized array.  */
 472           GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 1);
 473           return;
 474         }
 475       else
 476         retarray->base_addr = xmallocarray (alloc_size, sizeof (GFC_REAL_8));
 477     }
 478   else
 479     {
 480       if (rank != GFC_DESCRIPTOR_RANK (retarray))
 481         runtime_error ("rank of return array incorrect in"
 482                        " SUM intrinsic: is %ld, should be %ld",
 483                        (long int) (GFC_DESCRIPTOR_RANK (retarray)),
 484                        (long int) rank);
 485
 486       if (unlikely (compile_options.bounds_check))
 487         {
 488           for (n=0; n < rank; n++)
 489             {
 490               index_type ret_extent;
 491
 492               ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,n);
 493               if (extent[n] != ret_extent)
 494                 runtime_error ("Incorrect extent in return value of"
 495                                " SUM intrinsic in dimension %ld:"
 496                                " is %ld, should be %ld", (long int) n + 1,
 497                                (long int) ret_extent, (long int) extent[n]);
 498             }
 499         }
 500     }
 501
 502   for (n = 0; n < rank; n++)
 503     {
 504       count[n] = 0;
 505       dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);
 506     }
 507
 508   dest = retarray->base_addr;
 509
 510   while(1)
 511     {
 512       *dest = 0;
 513       count[0]++;
 514       dest += dstride[0];
 515       n = 0;
 516       while (count[n] == extent[n])
 517         {
 518           /* When we get to the end of a dimension, reset it and increment
 519              the next dimension.  */
 520           count[n] = 0;
 521           /* We could precalculate these products, but this is a less
 522              frequently used path so probably not worth it.  */
 523           dest -= dstride[n] * extent[n];
 524           n++;
 525           if (n >= rank)
 526             return;
 527           else
 528             {
 529               count[n]++;
 530               dest += dstride[n];
 531             }
 532         }
 533     }
 534 }
 535
 536 #endif