libgfortran/generated/product_r8.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_r8 (gfc_array_r8 *, gfc_array_r8 *, index_type *);
  29 export_proto_np(__product_r8);
  30
  31 void
  32 __product_r8 (gfc_array_r8 *retarray, gfc_array_r8 *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_8 *base;
  39   GFC_REAL_8 *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_8)
  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_8 *src;
 103       GFC_REAL_8 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_r8 (gfc_array_r8 *, gfc_array_r8 *, index_type *,
 153                                                 gfc_array_l4 *);
 154 export_proto_np(__mproduct_r8);
 155
 156 void
 157 __mproduct_r8 (gfc_array_r8 * retarray, gfc_array_r8 * array, index_type *pdim, gfc_array_l4 * mask)
 158 {
 159   index_type count[GFC_MAX_DIMENSIONS - 1];
 160   index_type extent[GFC_MAX_DIMENSIONS - 1];
 161   index_type sstride[GFC_MAX_DIMENSIONS - 1];
 162   index_type dstride[GFC_MAX_DIMENSIONS - 1];
 163   index_type mstride[GFC_MAX_DIMENSIONS - 1];
 164   GFC_REAL_8 *dest;
 165   GFC_REAL_8 *base;
 166   GFC_LOGICAL_4 *mbase;
 167   int rank;
 168   int dim;
 169   index_type n;
 170   index_type len;
 171   index_type delta;
 172   index_type mdelta;
 173
 174   dim = (*pdim) - 1;
 175   rank = GFC_DESCRIPTOR_RANK (array) - 1;
 176   assert (rank == GFC_DESCRIPTOR_RANK (retarray));
 177   if (array->dim[0].stride == 0)
 178     array->dim[0].stride = 1;
 179   if (retarray->dim[0].stride == 0)
 180     retarray->dim[0].stride = 1;
 181
 182   len = array->dim[dim].ubound + 1 - array->dim[dim].lbound;
 183   if (len <= 0)
 184     return;
 185   delta = array->dim[dim].stride;
 186   mdelta = mask->dim[dim].stride;
 187
 188   for (n = 0; n < dim; n++)
 189     {
 190       sstride[n] = array->dim[n].stride;
 191       mstride[n] = mask->dim[n].stride;
 192       extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;
 193     }
 194   for (n = dim; n < rank; n++)
 195     {
 196       sstride[n] = array->dim[n + 1].stride;
 197       mstride[n] = mask->dim[n + 1].stride;
 198       extent[n] =
 199         array->dim[n + 1].ubound + 1 - array->dim[n + 1].lbound;
 200     }
 201
 202   for (n = 0; n < rank; n++)
 203     {
 204       count[n] = 0;
 205       dstride[n] = retarray->dim[n].stride;
 206       if (extent[n] <= 0)
 207         return;
 208     }
 209
 210   dest = retarray->data;
 211   base = array->data;
 212   mbase = mask->data;
 213
 214   if (GFC_DESCRIPTOR_SIZE (mask) != 4)
 215     {
 216       /* This allows the same loop to be used for all logical types.  */
 217       assert (GFC_DESCRIPTOR_SIZE (mask) == 8);
 218       for (n = 0; n < rank; n++)
 219         mstride[n] <<= 1;
 220       mdelta <<= 1;
 221       mbase = (GFOR_POINTER_L8_TO_L4 (mbase));
 222     }
 223
 224   while (base)
 225     {
 226       GFC_REAL_8 *src;
 227       GFC_LOGICAL_4 *msrc;
 228       GFC_REAL_8 result;
 229       src = base;
 230       msrc = mbase;
 231       {
 232
 233   result = 1;
 234         if (len <= 0)
 235           *dest = 1;
 236         else
 237           {
 238             for (n = 0; n < len; n++, src += delta, msrc += mdelta)
 239               {
 240
 241   if (*msrc)
 242     result *= *src;
 243               }
 244             *dest = result;
 245           }
 246       }
 247       /* Advance to the next element.  */
 248       count[0]++;
 249       base += sstride[0];
 250       mbase += mstride[0];
 251       dest += dstride[0];
 252       n = 0;
 253       while (count[n] == extent[n])
 254         {
 255           /* When we get to the end of a dimension, reset it and increment
 256              the next dimension.  */
 257           count[n] = 0;
 258           /* We could precalculate these products, but this is a less
 259              frequently used path so proabably not worth it.  */
 260           base -= sstride[n] * extent[n];
 261           mbase -= mstride[n] * extent[n];
 262           dest -= dstride[n] * extent[n];
 263           n++;
 264           if (n == rank)
 265             {
 266               /* Break out of the look.  */
 267               base = NULL;
 268               break;
 269             }
 270           else
 271             {
 272               count[n]++;
 273               base += sstride[n];
 274               mbase += mstride[n];
 275               dest += dstride[n];
 276             }
 277         }
 278     }
 279 }
 280