libgfortran/generated/product_c8.c

   1 /* Implementation of the PRODUCT intrinsic
   2    Copyright 2002, 2007 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 2 of the License, or (at your option) any later version.
  11
  12 In addition to the permissions in the GNU General Public License, the
  13 Free Software Foundation gives you unlimited permission to link the
  14 compiled version of this file into combinations with other programs,
  15 and to distribute those combinations without any restriction coming
  16 from the use of this file.  (The General Public License restrictions
  17 do apply in other respects; for example, they cover modification of
  18 the file, and distribution when not linked into a combine
  19 executable.)
  20
  21 Libgfortran is distributed in the hope that it will be useful,
  22 but WITHOUT ANY WARRANTY; without even the implied warranty of
  23 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  24 GNU General Public License for more details.
  25
  26 You should have received a copy of the GNU General Public
  27 License along with libgfortran; see the file COPYING.  If not,
  28 write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
  29 Boston, MA 02110-1301, USA.  */
  30
  31 #include "libgfortran.h"
  32 #include <stdlib.h>
  33 #include <assert.h>
  34
  35
  36 #if defined (HAVE_GFC_COMPLEX_8) && defined (HAVE_GFC_COMPLEX_8)
  37
  38
  39 extern void product_c8 (gfc_array_c8 * const restrict,
  40         gfc_array_c8 * const restrict, const index_type * const restrict);
  41 export_proto(product_c8);
  42
  43 void
  44 product_c8 (gfc_array_c8 * const restrict retarray,
  45         gfc_array_c8 * const restrict array,
  46         const index_type * const restrict pdim)
  47 {
  48   index_type count[GFC_MAX_DIMENSIONS];
  49   index_type extent[GFC_MAX_DIMENSIONS];
  50   index_type sstride[GFC_MAX_DIMENSIONS];
  51   index_type dstride[GFC_MAX_DIMENSIONS];
  52   const GFC_COMPLEX_8 * restrict base;
  53   GFC_COMPLEX_8 * restrict dest;
  54   index_type rank;
  55   index_type n;
  56   index_type len;
  57   index_type delta;
  58   index_type dim;
  59
  60   /* Make dim zero based to avoid confusion.  */
  61   dim = (*pdim) - 1;
  62   rank = GFC_DESCRIPTOR_RANK (array) - 1;
  63
  64   len = array->dim[dim].ubound + 1 - array->dim[dim].lbound;
  65   delta = array->dim[dim].stride;
  66
  67   for (n = 0; n < dim; n++)
  68     {
  69       sstride[n] = array->dim[n].stride;
  70       extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;
  71
  72       if (extent[n] < 0)
  73         extent[n] = 0;
  74     }
  75   for (n = dim; n < rank; n++)
  76     {
  77       sstride[n] = array->dim[n + 1].stride;
  78       extent[n] =
  79         array->dim[n + 1].ubound + 1 - array->dim[n + 1].lbound;
  80
  81       if (extent[n] < 0)
  82         extent[n] = 0;
  83     }
  84
  85   if (retarray->data == NULL)
  86     {
  87       size_t alloc_size;
  88
  89       for (n = 0; n < rank; n++)
  90         {
  91           retarray->dim[n].lbound = 0;
  92           retarray->dim[n].ubound = extent[n]-1;
  93           if (n == 0)
  94             retarray->dim[n].stride = 1;
  95           else
  96             retarray->dim[n].stride = retarray->dim[n-1].stride * extent[n-1];
  97         }
  98
  99       retarray->offset = 0;
 100       retarray->dtype = (array->dtype & ~GFC_DTYPE_RANK_MASK) | rank;
 101
 102       alloc_size = sizeof (GFC_COMPLEX_8) * retarray->dim[rank-1].stride
 103                    * extent[rank-1];
 104
 105       if (alloc_size == 0)
 106         {
 107           /* Make sure we have a zero-sized array.  */
 108           retarray->dim[0].lbound = 0;
 109           retarray->dim[0].ubound = -1;
 110           return;
 111         }
 112       else
 113         retarray->data = internal_malloc_size (alloc_size);
 114     }
 115   else
 116     {
 117       if (rank != GFC_DESCRIPTOR_RANK (retarray))
 118         runtime_error ("rank of return array incorrect in"
 119                        " PRODUCT intrinsic: is %d, should be %d",
 120                        GFC_DESCRIPTOR_RANK (retarray), rank);
 121
 122       if (compile_options.bounds_check)
 123         {
 124           for (n=0; n < rank; n++)
 125             {
 126               index_type ret_extent;
 127
 128               ret_extent = retarray->dim[n].ubound + 1
 129                 - retarray->dim[n].lbound;
 130               if (extent[n] != ret_extent)
 131                 runtime_error ("Incorrect extent in return value of"
 132                                " PRODUCT intrinsic in dimension %d:"
 133                                " is %ld, should be %ld", n + 1,
 134                                (long int) ret_extent, (long int) extent[n]);
 135             }
 136         }
 137     }
 138
 139   for (n = 0; n < rank; n++)
 140     {
 141       count[n] = 0;
 142       dstride[n] = retarray->dim[n].stride;
 143       if (extent[n] <= 0)
 144         len = 0;
 145     }
 146
 147   base = array->data;
 148   dest = retarray->data;
 149
 150   while (base)
 151     {
 152       const GFC_COMPLEX_8 * restrict src;
 153       GFC_COMPLEX_8 result;
 154       src = base;
 155       {
 156
 157   result = 1;
 158         if (len <= 0)
 159           *dest = 1;
 160         else
 161           {
 162             for (n = 0; n < len; n++, src += delta)
 163               {
 164
 165   result *= *src;
 166           }
 167             *dest = result;
 168           }
 169       }
 170       /* Advance to the next element.  */
 171       count[0]++;
 172       base += sstride[0];
 173       dest += dstride[0];
 174       n = 0;
 175       while (count[n] == extent[n])
 176         {
 177           /* When we get to the end of a dimension, reset it and increment
 178              the next dimension.  */
 179           count[n] = 0;
 180           /* We could precalculate these products, but this is a less
 181              frequently used path so probably not worth it.  */
 182           base -= sstride[n] * extent[n];
 183           dest -= dstride[n] * extent[n];
 184           n++;
 185           if (n == rank)
 186             {
 187               /* Break out of the look.  */
 188               base = NULL;
 189               break;
 190             }
 191           else
 192             {
 193               count[n]++;
 194               base += sstride[n];
 195               dest += dstride[n];
 196             }
 197         }
 198     }
 199 }
 200
 201
 202 extern void mproduct_c8 (gfc_array_c8 * const restrict,
 203         gfc_array_c8 * const restrict, const index_type * const restrict,
 204         gfc_array_l1 * const restrict);
 205 export_proto(mproduct_c8);
 206
 207 void
 208 mproduct_c8 (gfc_array_c8 * const restrict retarray,
 209         gfc_array_c8 * const restrict array,
 210         const index_type * const restrict pdim,
 211         gfc_array_l1 * const restrict mask)
 212 {
 213   index_type count[GFC_MAX_DIMENSIONS];
 214   index_type extent[GFC_MAX_DIMENSIONS];
 215   index_type sstride[GFC_MAX_DIMENSIONS];
 216   index_type dstride[GFC_MAX_DIMENSIONS];
 217   index_type mstride[GFC_MAX_DIMENSIONS];
 218   GFC_COMPLEX_8 * restrict dest;
 219   const GFC_COMPLEX_8 * restrict base;
 220   const GFC_LOGICAL_1 * restrict mbase;
 221   int rank;
 222   int dim;
 223   index_type n;
 224   index_type len;
 225   index_type delta;
 226   index_type mdelta;
 227   int mask_kind;
 228
 229   dim = (*pdim) - 1;
 230   rank = GFC_DESCRIPTOR_RANK (array) - 1;
 231
 232   len = array->dim[dim].ubound + 1 - array->dim[dim].lbound;
 233   if (len <= 0)
 234     return;
 235
 236   mbase = mask->data;
 237
 238   mask_kind = GFC_DESCRIPTOR_SIZE (mask);
 239
 240   if (mask_kind == 1 || mask_kind == 2 || mask_kind == 4 || mask_kind == 8
 241 #ifdef HAVE_GFC_LOGICAL_16
 242       || mask_kind == 16
 243 #endif
 244       )
 245     mbase = GFOR_POINTER_TO_L1 (mbase, mask_kind);
 246   else
 247     runtime_error ("Funny sized logical array");
 248
 249   delta = array->dim[dim].stride;
 250   mdelta = mask->dim[dim].stride * mask_kind;
 251
 252   for (n = 0; n < dim; n++)
 253     {
 254       sstride[n] = array->dim[n].stride;
 255       mstride[n] = mask->dim[n].stride * mask_kind;
 256       extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;
 257
 258       if (extent[n] < 0)
 259         extent[n] = 0;
 260
 261     }
 262   for (n = dim; n < rank; n++)
 263     {
 264       sstride[n] = array->dim[n + 1].stride;
 265       mstride[n] = mask->dim[n + 1].stride * mask_kind;
 266       extent[n] =
 267         array->dim[n + 1].ubound + 1 - array->dim[n + 1].lbound;
 268
 269       if (extent[n] < 0)
 270         extent[n] = 0;
 271     }
 272
 273   if (retarray->data == NULL)
 274     {
 275       size_t alloc_size;
 276
 277       for (n = 0; n < rank; n++)
 278         {
 279           retarray->dim[n].lbound = 0;
 280           retarray->dim[n].ubound = extent[n]-1;
 281           if (n == 0)
 282             retarray->dim[n].stride = 1;
 283           else
 284             retarray->dim[n].stride = retarray->dim[n-1].stride * extent[n-1];
 285         }
 286
 287       alloc_size = sizeof (GFC_COMPLEX_8) * retarray->dim[rank-1].stride
 288                    * extent[rank-1];
 289
 290       retarray->offset = 0;
 291       retarray->dtype = (array->dtype & ~GFC_DTYPE_RANK_MASK) | rank;
 292
 293       if (alloc_size == 0)
 294         {
 295           /* Make sure we have a zero-sized array.  */
 296           retarray->dim[0].lbound = 0;
 297           retarray->dim[0].ubound = -1;
 298           return;
 299         }
 300       else
 301         retarray->data = internal_malloc_size (alloc_size);
 302
 303     }
 304   else
 305     {
 306       if (rank != GFC_DESCRIPTOR_RANK (retarray))
 307         runtime_error ("rank of return array incorrect in PRODUCT intrinsic");
 308
 309       if (compile_options.bounds_check)
 310         {
 311           for (n=0; n < rank; n++)
 312             {
 313               index_type ret_extent;
 314
 315               ret_extent = retarray->dim[n].ubound + 1
 316                 - retarray->dim[n].lbound;
 317               if (extent[n] != ret_extent)
 318                 runtime_error ("Incorrect extent in return value of"
 319                                " PRODUCT intrinsic in dimension %d:"
 320                                " is %ld, should be %ld", n + 1,
 321                                (long int) ret_extent, (long int) extent[n]);
 322             }
 323           for (n=0; n<= rank; n++)
 324             {
 325               index_type mask_extent, array_extent;
 326
 327               array_extent = array->dim[n].ubound + 1 - array->dim[n].lbound;
 328               mask_extent = mask->dim[n].ubound + 1 - mask->dim[n].lbound;
 329               if (array_extent != mask_extent)
 330                 runtime_error ("Incorrect extent in MASK argument of"
 331                                " PRODUCT intrinsic in dimension %d:"
 332                                " is %ld, should be %ld", n + 1,
 333                                (long int) mask_extent, (long int) array_extent);
 334             }
 335         }
 336     }
 337
 338   for (n = 0; n < rank; n++)
 339     {
 340       count[n] = 0;
 341       dstride[n] = retarray->dim[n].stride;
 342       if (extent[n] <= 0)
 343         return;
 344     }
 345
 346   dest = retarray->data;
 347   base = array->data;
 348
 349   while (base)
 350     {
 351       const GFC_COMPLEX_8 * restrict src;
 352       const GFC_LOGICAL_1 * restrict msrc;
 353       GFC_COMPLEX_8 result;
 354       src = base;
 355       msrc = mbase;
 356       {
 357
 358   result = 1;
 359         if (len <= 0)
 360           *dest = 1;
 361         else
 362           {
 363             for (n = 0; n < len; n++, src += delta, msrc += mdelta)
 364               {
 365
 366   if (*msrc)
 367     result *= *src;
 368               }
 369             *dest = result;
 370           }
 371       }
 372       /* Advance to the next element.  */
 373       count[0]++;
 374       base += sstride[0];
 375       mbase += mstride[0];
 376       dest += dstride[0];
 377       n = 0;
 378       while (count[n] == extent[n])
 379         {
 380           /* When we get to the end of a dimension, reset it and increment
 381              the next dimension.  */
 382           count[n] = 0;
 383           /* We could precalculate these products, but this is a less
 384              frequently used path so probably not worth it.  */
 385           base -= sstride[n] * extent[n];
 386           mbase -= mstride[n] * extent[n];
 387           dest -= dstride[n] * extent[n];
 388           n++;
 389           if (n == rank)
 390             {
 391               /* Break out of the look.  */
 392               base = NULL;
 393               break;
 394             }
 395           else
 396             {
 397               count[n]++;
 398               base += sstride[n];
 399               mbase += mstride[n];
 400               dest += dstride[n];
 401             }
 402         }
 403     }
 404 }
 405
 406
 407 extern void sproduct_c8 (gfc_array_c8 * const restrict,
 408         gfc_array_c8 * const restrict, const index_type * const restrict,
 409         GFC_LOGICAL_4 *);
 410 export_proto(sproduct_c8);
 411
 412 void
 413 sproduct_c8 (gfc_array_c8 * const restrict retarray,
 414         gfc_array_c8 * const restrict array,
 415         const index_type * const restrict pdim,
 416         GFC_LOGICAL_4 * mask)
 417 {
 418   index_type rank;
 419   index_type n;
 420   index_type dstride;
 421   GFC_COMPLEX_8 *dest;
 422
 423   if (*mask)
 424     {
 425       product_c8 (retarray, array, pdim);
 426       return;
 427     }
 428     rank = GFC_DESCRIPTOR_RANK (array);
 429   if (rank <= 0)
 430     runtime_error ("Rank of array needs to be > 0");
 431
 432   if (retarray->data == NULL)
 433     {
 434       retarray->dim[0].lbound = 0;
 435       retarray->dim[0].ubound = rank-1;
 436       retarray->dim[0].stride = 1;
 437       retarray->dtype = (retarray->dtype & ~GFC_DTYPE_RANK_MASK) | 1;
 438       retarray->offset = 0;
 439       retarray->data = internal_malloc_size (sizeof (GFC_COMPLEX_8) * rank);
 440     }
 441   else
 442     {
 443       if (compile_options.bounds_check)
 444         {
 445           int ret_rank;
 446           index_type ret_extent;
 447
 448           ret_rank = GFC_DESCRIPTOR_RANK (retarray);
 449           if (ret_rank != 1)
 450             runtime_error ("rank of return array in PRODUCT intrinsic"
 451                            " should be 1, is %d", ret_rank);
 452
 453           ret_extent = retarray->dim[0].ubound + 1 - retarray->dim[0].lbound;
 454             if (ret_extent != rank)
 455               runtime_error ("dimension of return array incorrect");
 456         }
 457     }
 458     dstride = retarray->dim[0].stride;
 459     dest = retarray->data;
 460
 461     for (n = 0; n < rank; n++)
 462       dest[n * dstride] = 1 ;
 463 }
 464
 465 #endif