libgfortran/generated/sum_r8.c

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