libgfortran/generated/product_r4.c

   1 /* Implementation of the PRODUCT intrinsic
   2    Copyright 2002 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 (libgfor).
   6
   7 Libgfortran is free software; you can redistribute it and/or
   8 modify it under the terms of the GNU Lesser General Public
   9 License as published by the Free Software Foundation; either
  10 version 2.1 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 Lesser General Public License for more details.
  16
  17 You should have received a copy of the GNU Lesser General Public
  18 License along with libgfor; see the file COPYING.LIB.  If not,
  19 write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  20 Boston, MA 02111-1307, USA.  */
  21
  22 #include "config.h"
  23 #include <stdlib.h>
  24 #include <assert.h>
  25 #include "libgfortran.h"
  26
  27
  28 extern void product_r4 (gfc_array_r4 *, gfc_array_r4 *, index_type *);
  29 export_proto(product_r4);
  30
  31 void
  32 product_r4 (gfc_array_r4 *retarray, gfc_array_r4 *array, index_type *pdim)
  33 {
  34   index_type count[GFC_MAX_DIMENSIONS - 1];
  35   index_type extent[GFC_MAX_DIMENSIONS - 1];
  36   index_type sstride[GFC_MAX_DIMENSIONS - 1];
  37   index_type dstride[GFC_MAX_DIMENSIONS - 1];
  38   GFC_REAL_4 *base;
  39   GFC_REAL_4 *dest;
  40   index_type rank;
  41   index_type n;
  42   index_type len;
  43   index_type delta;
  44   index_type dim;
  45
  46   /* Make dim zero based to avoid confusion.  */
  47   dim = (*pdim) - 1;
  48   rank = GFC_DESCRIPTOR_RANK (array) - 1;
  49   assert (rank == GFC_DESCRIPTOR_RANK (retarray));
  50   if (array->dim[0].stride == 0)
  51     array->dim[0].stride = 1;
  52   if (retarray->dim[0].stride == 0)
  53     retarray->dim[0].stride = 1;
  54
  55   len = array->dim[dim].ubound + 1 - array->dim[dim].lbound;
  56   delta = array->dim[dim].stride;
  57
  58   for (n = 0; n < dim; n++)
  59     {
  60       sstride[n] = array->dim[n].stride;
  61       extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;
  62     }
  63   for (n = dim; n < rank; n++)
  64     {
  65       sstride[n] = array->dim[n + 1].stride;
  66       extent[n] =
  67         array->dim[n + 1].ubound + 1 - array->dim[n + 1].lbound;
  68     }
  69
  70   if (retarray->data == NULL)
  71     {
  72       for (n = 0; n < rank; n++)
  73         {
  74           retarray->dim[n].lbound = 0;
  75           retarray->dim[n].ubound = extent[n]-1;
  76           if (n == 0)
  77             retarray->dim[n].stride = 1;
  78           else
  79             retarray->dim[n].stride = retarray->dim[n-1].stride * extent[n-1];
  80         }
  81
  82       retarray->data
  83          = internal_malloc_size (sizeof (GFC_REAL_4)
  84                                  * retarray->dim[rank-1].stride
  85                                  * extent[rank-1]);
  86       retarray->base = 0;
  87     }
  88
  89   for (n = 0; n < rank; n++)
  90     {
  91       count[n] = 0;
  92       dstride[n] = retarray->dim[n].stride;
  93       if (extent[n] <= 0)
  94         len = 0;
  95     }
  96
  97   base = array->data;
  98   dest = retarray->data;
  99
 100   while (base)
 101     {
 102       GFC_REAL_4 *src;
 103       GFC_REAL_4 result;
 104       src = base;
 105       {
 106
 107   result = 1;
 108         if (len <= 0)
 109           *dest = 1;
 110         else
 111           {
 112             for (n = 0; n < len; n++, src += delta)
 113               {
 114
 115   result *= *src;
 116           }
 117             *dest = result;
 118           }
 119       }
 120       /* Advance to the next element.  */
 121       count[0]++;
 122       base += sstride[0];
 123       dest += dstride[0];
 124       n = 0;
 125       while (count[n] == extent[n])
 126         {
 127           /* When we get to the end of a dimension, reset it and increment
 128              the next dimension.  */
 129           count[n] = 0;
 130           /* We could precalculate these products, but this is a less
 131              frequently used path so proabably not worth it.  */
 132           base -= sstride[n] * extent[n];
 133           dest -= dstride[n] * extent[n];
 134           n++;
 135           if (n == rank)
 136             {
 137               /* Break out of the look.  */
 138               base = NULL;
 139               break;
 140             }
 141           else
 142             {
 143               count[n]++;
 144               base += sstride[n];
 145               dest += dstride[n];
 146             }
 147         }
 148     }
 149 }
 150
 151
 152 extern void mproduct_r4 (gfc_array_r4 *, gfc_array_r4 *, index_type *,
 153                                                gfc_array_l4 *);
 154 export_proto(mproduct_r4);
 155
 156 void
 157 mproduct_r4 (gfc_array_r4 * retarray, gfc_array_r4 * array,
 158                                   index_type *pdim, gfc_array_l4 * mask)
 159 {
 160   index_type count[GFC_MAX_DIMENSIONS - 1];
 161   index_type extent[GFC_MAX_DIMENSIONS - 1];
 162   index_type sstride[GFC_MAX_DIMENSIONS - 1];
 163   index_type dstride[GFC_MAX_DIMENSIONS - 1];
 164   index_type mstride[GFC_MAX_DIMENSIONS - 1];
 165   GFC_REAL_4 *dest;
 166   GFC_REAL_4 *base;
 167   GFC_LOGICAL_4 *mbase;
 168   int rank;
 169   int dim;
 170   index_type n;
 171   index_type len;
 172   index_type delta;
 173   index_type mdelta;
 174
 175   dim = (*pdim) - 1;
 176   rank = GFC_DESCRIPTOR_RANK (array) - 1;
 177   assert (rank == GFC_DESCRIPTOR_RANK (retarray));
 178   if (array->dim[0].stride == 0)
 179     array->dim[0].stride = 1;
 180   if (retarray->dim[0].stride == 0)
 181     retarray->dim[0].stride = 1;
 182
 183   len = array->dim[dim].ubound + 1 - array->dim[dim].lbound;
 184   if (len <= 0)
 185     return;
 186   delta = array->dim[dim].stride;
 187   mdelta = mask->dim[dim].stride;
 188
 189   for (n = 0; n < dim; n++)
 190     {
 191       sstride[n] = array->dim[n].stride;
 192       mstride[n] = mask->dim[n].stride;
 193       extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;
 194     }
 195   for (n = dim; n < rank; n++)
 196     {
 197       sstride[n] = array->dim[n + 1].stride;
 198       mstride[n] = mask->dim[n + 1].stride;
 199       extent[n] =
 200         array->dim[n + 1].ubound + 1 - array->dim[n + 1].lbound;
 201     }
 202
 203   for (n = 0; n < rank; n++)
 204     {
 205       count[n] = 0;
 206       dstride[n] = retarray->dim[n].stride;
 207       if (extent[n] <= 0)
 208         return;
 209     }
 210
 211   dest = retarray->data;
 212   base = array->data;
 213   mbase = mask->data;
 214
 215   if (GFC_DESCRIPTOR_SIZE (mask) != 4)
 216     {
 217       /* This allows the same loop to be used for all logical types.  */
 218       assert (GFC_DESCRIPTOR_SIZE (mask) == 8);
 219       for (n = 0; n < rank; n++)
 220         mstride[n] <<= 1;
 221       mdelta <<= 1;
 222       mbase = (GFOR_POINTER_L8_TO_L4 (mbase));
 223     }
 224
 225   while (base)
 226     {
 227       GFC_REAL_4 *src;
 228       GFC_LOGICAL_4 *msrc;
 229       GFC_REAL_4 result;
 230       src = base;
 231       msrc = mbase;
 232       {
 233
 234   result = 1;
 235         if (len <= 0)
 236           *dest = 1;
 237         else
 238           {
 239             for (n = 0; n < len; n++, src += delta, msrc += mdelta)
 240               {
 241
 242   if (*msrc)
 243     result *= *src;
 244               }
 245             *dest = result;
 246           }
 247       }
 248       /* Advance to the next element.  */
 249       count[0]++;
 250       base += sstride[0];
 251       mbase += mstride[0];
 252       dest += dstride[0];
 253       n = 0;
 254       while (count[n] == extent[n])
 255         {
 256           /* When we get to the end of a dimension, reset it and increment
 257              the next dimension.  */
 258           count[n] = 0;
 259           /* We could precalculate these products, but this is a less
 260              frequently used path so proabably not worth it.  */
 261           base -= sstride[n] * extent[n];
 262           mbase -= mstride[n] * extent[n];
 263           dest -= dstride[n] * extent[n];
 264           n++;
 265           if (n == rank)
 266             {
 267               /* Break out of the look.  */
 268               base = NULL;
 269               break;
 270             }
 271           else
 272             {
 273               count[n]++;
 274               base += sstride[n];
 275               mbase += mstride[n];
 276               dest += dstride[n];
 277             }
 278         }
 279     }
 280 }
 281