libgfortran/generated/sum_c4.c

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