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