libgfortran/generated/product_c16.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_16) && defined (HAVE_GFC_COMPLEX_16)
  37
  38
  39 extern void product_c16 (gfc_array_c16 * const restrict,
  40         gfc_array_c16 * const restrict, const index_type * const restrict);
  41 export_proto(product_c16);
  42
  43 void
  44 product_c16 (gfc_array_c16 * const restrict retarray,
  45         gfc_array_c16 * 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_16 * restrict base;
  53   GFC_COMPLEX_16 * 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_16) * 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");
 119     }
 120
 121   for (n = 0; n < rank; n++)
 122     {
 123       count[n] = 0;
 124       dstride[n] = retarray->dim[n].stride;
 125       if (extent[n] <= 0)
 126         len = 0;
 127     }
 128
 129   base = array->data;
 130   dest = retarray->data;
 131
 132   while (base)
 133     {
 134       const GFC_COMPLEX_16 * restrict src;
 135       GFC_COMPLEX_16 result;
 136       src = base;
 137       {
 138
 139   result = 1;
 140         if (len <= 0)
 141           *dest = 1;
 142         else
 143           {
 144             for (n = 0; n < len; n++, src += delta)
 145               {
 146
 147   result *= *src;
 148           }
 149             *dest = result;
 150           }
 151       }
 152       /* Advance to the next element.  */
 153       count[0]++;
 154       base += sstride[0];
 155       dest += dstride[0];
 156       n = 0;
 157       while (count[n] == extent[n])
 158         {
 159           /* When we get to the end of a dimension, reset it and increment
 160              the next dimension.  */
 161           count[n] = 0;
 162           /* We could precalculate these products, but this is a less
 163              frequently used path so probably not worth it.  */
 164           base -= sstride[n] * extent[n];
 165           dest -= dstride[n] * extent[n];
 166           n++;
 167           if (n == rank)
 168             {
 169               /* Break out of the look.  */
 170               base = NULL;
 171               break;
 172             }
 173           else
 174             {
 175               count[n]++;
 176               base += sstride[n];
 177               dest += dstride[n];
 178             }
 179         }
 180     }
 181 }
 182
 183
 184 extern void mproduct_c16 (gfc_array_c16 * const restrict,
 185         gfc_array_c16 * const restrict, const index_type * const restrict,
 186         gfc_array_l1 * const restrict);
 187 export_proto(mproduct_c16);
 188
 189 void
 190 mproduct_c16 (gfc_array_c16 * const restrict retarray,
 191         gfc_array_c16 * const restrict array,
 192         const index_type * const restrict pdim,
 193         gfc_array_l1 * const restrict mask)
 194 {
 195   index_type count[GFC_MAX_DIMENSIONS];
 196   index_type extent[GFC_MAX_DIMENSIONS];
 197   index_type sstride[GFC_MAX_DIMENSIONS];
 198   index_type dstride[GFC_MAX_DIMENSIONS];
 199   index_type mstride[GFC_MAX_DIMENSIONS];
 200   GFC_COMPLEX_16 * restrict dest;
 201   const GFC_COMPLEX_16 * restrict base;
 202   const GFC_LOGICAL_1 * restrict mbase;
 203   int rank;
 204   int dim;
 205   index_type n;
 206   index_type len;
 207   index_type delta;
 208   index_type mdelta;
 209   int mask_kind;
 210
 211   dim = (*pdim) - 1;
 212   rank = GFC_DESCRIPTOR_RANK (array) - 1;
 213
 214   len = array->dim[dim].ubound + 1 - array->dim[dim].lbound;
 215   if (len <= 0)
 216     return;
 217
 218   mbase = mask->data;
 219
 220   mask_kind = GFC_DESCRIPTOR_SIZE (mask);
 221
 222   if (mask_kind == 1 || mask_kind == 2 || mask_kind == 4 || mask_kind == 8
 223 #ifdef HAVE_GFC_LOGICAL_16
 224       || mask_kind == 16
 225 #endif
 226       )
 227     mbase = GFOR_POINTER_TO_L1 (mbase, mask_kind);
 228   else
 229     runtime_error ("Funny sized logical array");
 230
 231   delta = array->dim[dim].stride;
 232   mdelta = mask->dim[dim].stride * mask_kind;
 233
 234   for (n = 0; n < dim; n++)
 235     {
 236       sstride[n] = array->dim[n].stride;
 237       mstride[n] = mask->dim[n].stride * mask_kind;
 238       extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;
 239
 240       if (extent[n] < 0)
 241         extent[n] = 0;
 242
 243     }
 244   for (n = dim; n < rank; n++)
 245     {
 246       sstride[n] = array->dim[n + 1].stride;
 247       mstride[n] = mask->dim[n + 1].stride * mask_kind;
 248       extent[n] =
 249         array->dim[n + 1].ubound + 1 - array->dim[n + 1].lbound;
 250
 251       if (extent[n] < 0)
 252         extent[n] = 0;
 253     }
 254
 255   if (retarray->data == NULL)
 256     {
 257       size_t alloc_size;
 258
 259       for (n = 0; n < rank; n++)
 260         {
 261           retarray->dim[n].lbound = 0;
 262           retarray->dim[n].ubound = extent[n]-1;
 263           if (n == 0)
 264             retarray->dim[n].stride = 1;
 265           else
 266             retarray->dim[n].stride = retarray->dim[n-1].stride * extent[n-1];
 267         }
 268
 269       alloc_size = sizeof (GFC_COMPLEX_16) * retarray->dim[rank-1].stride
 270                    * extent[rank-1];
 271
 272       retarray->offset = 0;
 273       retarray->dtype = (array->dtype & ~GFC_DTYPE_RANK_MASK) | rank;
 274
 275       if (alloc_size == 0)
 276         {
 277           /* Make sure we have a zero-sized array.  */
 278           retarray->dim[0].lbound = 0;
 279           retarray->dim[0].ubound = -1;
 280           return;
 281         }
 282       else
 283         retarray->data = internal_malloc_size (alloc_size);
 284
 285     }
 286   else
 287     {
 288       if (rank != GFC_DESCRIPTOR_RANK (retarray))
 289         runtime_error ("rank of return array incorrect");
 290     }
 291
 292   for (n = 0; n < rank; n++)
 293     {
 294       count[n] = 0;
 295       dstride[n] = retarray->dim[n].stride;
 296       if (extent[n] <= 0)
 297         return;
 298     }
 299
 300   dest = retarray->data;
 301   base = array->data;
 302
 303   while (base)
 304     {
 305       const GFC_COMPLEX_16 * restrict src;
 306       const GFC_LOGICAL_1 * restrict msrc;
 307       GFC_COMPLEX_16 result;
 308       src = base;
 309       msrc = mbase;
 310       {
 311
 312   result = 1;
 313         if (len <= 0)
 314           *dest = 1;
 315         else
 316           {
 317             for (n = 0; n < len; n++, src += delta, msrc += mdelta)
 318               {
 319
 320   if (*msrc)
 321     result *= *src;
 322               }
 323             *dest = result;
 324           }
 325       }
 326       /* Advance to the next element.  */
 327       count[0]++;
 328       base += sstride[0];
 329       mbase += mstride[0];
 330       dest += dstride[0];
 331       n = 0;
 332       while (count[n] == extent[n])
 333         {
 334           /* When we get to the end of a dimension, reset it and increment
 335              the next dimension.  */
 336           count[n] = 0;
 337           /* We could precalculate these products, but this is a less
 338              frequently used path so probably not worth it.  */
 339           base -= sstride[n] * extent[n];
 340           mbase -= mstride[n] * extent[n];
 341           dest -= dstride[n] * extent[n];
 342           n++;
 343           if (n == rank)
 344             {
 345               /* Break out of the look.  */
 346               base = NULL;
 347               break;
 348             }
 349           else
 350             {
 351               count[n]++;
 352               base += sstride[n];
 353               mbase += mstride[n];
 354               dest += dstride[n];
 355             }
 356         }
 357     }
 358 }
 359
 360
 361 extern void sproduct_c16 (gfc_array_c16 * const restrict,
 362         gfc_array_c16 * const restrict, const index_type * const restrict,
 363         GFC_LOGICAL_4 *);
 364 export_proto(sproduct_c16);
 365
 366 void
 367 sproduct_c16 (gfc_array_c16 * const restrict retarray,
 368         gfc_array_c16 * const restrict array,
 369         const index_type * const restrict pdim,
 370         GFC_LOGICAL_4 * mask)
 371 {
 372   index_type rank;
 373   index_type n;
 374   index_type dstride;
 375   GFC_COMPLEX_16 *dest;
 376
 377   if (*mask)
 378     {
 379       product_c16 (retarray, array, pdim);
 380       return;
 381     }
 382     rank = GFC_DESCRIPTOR_RANK (array);
 383   if (rank <= 0)
 384     runtime_error ("Rank of array needs to be > 0");
 385
 386   if (retarray->data == NULL)
 387     {
 388       retarray->dim[0].lbound = 0;
 389       retarray->dim[0].ubound = rank-1;
 390       retarray->dim[0].stride = 1;
 391       retarray->dtype = (retarray->dtype & ~GFC_DTYPE_RANK_MASK) | 1;
 392       retarray->offset = 0;
 393       retarray->data = internal_malloc_size (sizeof (GFC_COMPLEX_16) * rank);
 394     }
 395   else
 396     {
 397       if (GFC_DESCRIPTOR_RANK (retarray) != 1)
 398         runtime_error ("rank of return array does not equal 1");
 399
 400       if (retarray->dim[0].ubound + 1 - retarray->dim[0].lbound != rank)
 401         runtime_error ("dimension of return array incorrect");
 402     }
 403
 404     dstride = retarray->dim[0].stride;
 405     dest = retarray->data;
 406
 407     for (n = 0; n < rank; n++)
 408       dest[n * dstride] = 1 ;
 409 }
 410
 411 #endif