libgfortran/generated/product_c4.c

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