gomp-20050608-branch/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 (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 "config.h"
  32 #include <stdlib.h>
  33 #include <assert.h>
  34 #include "libgfortran.h"
  35
  36
  37 #if defined (HAVE_GFC_REAL_8) && defined (HAVE_GFC_REAL_8)
  38
  39
  40 extern void sum_r8 (gfc_array_r8 * const restrict,
  41         gfc_array_r8 * const restrict, const index_type * const restrict);
  42 export_proto(sum_r8);
  43
  44 void
  45 sum_r8 (gfc_array_r8 * const restrict retarray,
  46         gfc_array_r8 * const restrict array,
  47         const index_type * const restrict pdim)
  48 {
  49   index_type count[GFC_MAX_DIMENSIONS];
  50   index_type extent[GFC_MAX_DIMENSIONS];
  51   index_type sstride[GFC_MAX_DIMENSIONS];
  52   index_type dstride[GFC_MAX_DIMENSIONS];
  53   const GFC_REAL_8 * restrict base;
  54   GFC_REAL_8 * restrict dest;
  55   index_type rank;
  56   index_type n;
  57   index_type len;
  58   index_type delta;
  59   index_type dim;
  60
  61   /* Make dim zero based to avoid confusion.  */
  62   dim = (*pdim) - 1;
  63   rank = GFC_DESCRIPTOR_RANK (array) - 1;
  64
  65   /* TODO:  It should be a front end job to correctly set the strides.  */
  66
  67   if (array->dim[0].stride == 0)
  68     array->dim[0].stride = 1;
  69
  70   len = array->dim[dim].ubound + 1 - array->dim[dim].lbound;
  71   delta = array->dim[dim].stride;
  72
  73   for (n = 0; n < dim; n++)
  74     {
  75       sstride[n] = array->dim[n].stride;
  76       extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;
  77     }
  78   for (n = dim; n < rank; n++)
  79     {
  80       sstride[n] = array->dim[n + 1].stride;
  81       extent[n] =
  82         array->dim[n + 1].ubound + 1 - array->dim[n + 1].lbound;
  83     }
  84
  85   if (retarray->data == NULL)
  86     {
  87       for (n = 0; n < rank; n++)
  88         {
  89           retarray->dim[n].lbound = 0;
  90           retarray->dim[n].ubound = extent[n]-1;
  91           if (n == 0)
  92             retarray->dim[n].stride = 1;
  93           else
  94             retarray->dim[n].stride = retarray->dim[n-1].stride * extent[n-1];
  95         }
  96
  97       retarray->data
  98          = internal_malloc_size (sizeof (GFC_REAL_8)
  99                                  * retarray->dim[rank-1].stride
 100                                  * extent[rank-1]);
 101       retarray->offset = 0;
 102       retarray->dtype = (array->dtype & ~GFC_DTYPE_RANK_MASK) | rank;
 103     }
 104   else
 105     {
 106       if (retarray->dim[0].stride == 0)
 107         retarray->dim[0].stride = 1;
 108
 109       if (rank != GFC_DESCRIPTOR_RANK (retarray))
 110         runtime_error ("rank of return array incorrect");
 111     }
 112
 113   for (n = 0; n < rank; n++)
 114     {
 115       count[n] = 0;
 116       dstride[n] = retarray->dim[n].stride;
 117       if (extent[n] <= 0)
 118         len = 0;
 119     }
 120
 121   base = array->data;
 122   dest = retarray->data;
 123
 124   while (base)
 125     {
 126       const GFC_REAL_8 * restrict src;
 127       GFC_REAL_8 result;
 128       src = base;
 129       {
 130
 131   result = 0;
 132         if (len <= 0)
 133           *dest = 0;
 134         else
 135           {
 136             for (n = 0; n < len; n++, src += delta)
 137               {
 138
 139   result += *src;
 140           }
 141             *dest = result;
 142           }
 143       }
 144       /* Advance to the next element.  */
 145       count[0]++;
 146       base += sstride[0];
 147       dest += dstride[0];
 148       n = 0;
 149       while (count[n] == extent[n])
 150         {
 151           /* When we get to the end of a dimension, reset it and increment
 152              the next dimension.  */
 153           count[n] = 0;
 154           /* We could precalculate these products, but this is a less
 155              frequently used path so proabably not worth it.  */
 156           base -= sstride[n] * extent[n];
 157           dest -= dstride[n] * extent[n];
 158           n++;
 159           if (n == rank)
 160             {
 161               /* Break out of the look.  */
 162               base = NULL;
 163               break;
 164             }
 165           else
 166             {
 167               count[n]++;
 168               base += sstride[n];
 169               dest += dstride[n];
 170             }
 171         }
 172     }
 173 }
 174
 175
 176 extern void msum_r8 (gfc_array_r8 * const restrict,
 177         gfc_array_r8 * const restrict, const index_type * const restrict,
 178         gfc_array_l4 * const restrict);
 179 export_proto(msum_r8);
 180
 181 void
 182 msum_r8 (gfc_array_r8 * const restrict retarray,
 183         gfc_array_r8 * const restrict array,
 184         const index_type * const restrict pdim,
 185         gfc_array_l4 * const restrict mask)
 186 {
 187   index_type count[GFC_MAX_DIMENSIONS];
 188   index_type extent[GFC_MAX_DIMENSIONS];
 189   index_type sstride[GFC_MAX_DIMENSIONS];
 190   index_type dstride[GFC_MAX_DIMENSIONS];
 191   index_type mstride[GFC_MAX_DIMENSIONS];
 192   GFC_REAL_8 * restrict dest;
 193   const GFC_REAL_8 * restrict base;
 194   const GFC_LOGICAL_4 * restrict mbase;
 195   int rank;
 196   int dim;
 197   index_type n;
 198   index_type len;
 199   index_type delta;
 200   index_type mdelta;
 201
 202   dim = (*pdim) - 1;
 203   rank = GFC_DESCRIPTOR_RANK (array) - 1;
 204
 205   /* TODO:  It should be a front end job to correctly set the strides.  */
 206
 207   if (array->dim[0].stride == 0)
 208     array->dim[0].stride = 1;
 209
 210   if (mask->dim[0].stride == 0)
 211     mask->dim[0].stride = 1;
 212
 213   len = array->dim[dim].ubound + 1 - array->dim[dim].lbound;
 214   if (len <= 0)
 215     return;
 216   delta = array->dim[dim].stride;
 217   mdelta = mask->dim[dim].stride;
 218
 219   for (n = 0; n < dim; n++)
 220     {
 221       sstride[n] = array->dim[n].stride;
 222       mstride[n] = mask->dim[n].stride;
 223       extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;
 224     }
 225   for (n = dim; n < rank; n++)
 226     {
 227       sstride[n] = array->dim[n + 1].stride;
 228       mstride[n] = mask->dim[n + 1].stride;
 229       extent[n] =
 230         array->dim[n + 1].ubound + 1 - array->dim[n + 1].lbound;
 231     }
 232
 233   if (retarray->data == NULL)
 234     {
 235       for (n = 0; n < rank; n++)
 236         {
 237           retarray->dim[n].lbound = 0;
 238           retarray->dim[n].ubound = extent[n]-1;
 239           if (n == 0)
 240             retarray->dim[n].stride = 1;
 241           else
 242             retarray->dim[n].stride = retarray->dim[n-1].stride * extent[n-1];
 243         }
 244
 245       retarray->data
 246          = internal_malloc_size (sizeof (GFC_REAL_8)
 247                                  * retarray->dim[rank-1].stride
 248                                  * extent[rank-1]);
 249       retarray->offset = 0;
 250       retarray->dtype = (array->dtype & ~GFC_DTYPE_RANK_MASK) | rank;
 251     }
 252   else
 253     {
 254       if (retarray->dim[0].stride == 0)
 255         retarray->dim[0].stride = 1;
 256
 257       if (rank != GFC_DESCRIPTOR_RANK (retarray))
 258         runtime_error ("rank of return array incorrect");
 259     }
 260
 261   for (n = 0; n < rank; n++)
 262     {
 263       count[n] = 0;
 264       dstride[n] = retarray->dim[n].stride;
 265       if (extent[n] <= 0)
 266         return;
 267     }
 268
 269   dest = retarray->data;
 270   base = array->data;
 271   mbase = mask->data;
 272
 273   if (GFC_DESCRIPTOR_SIZE (mask) != 4)
 274     {
 275       /* This allows the same loop to be used for all logical types.  */
 276       assert (GFC_DESCRIPTOR_SIZE (mask) == 8);
 277       for (n = 0; n < rank; n++)
 278         mstride[n] <<= 1;
 279       mdelta <<= 1;
 280       mbase = (GFOR_POINTER_L8_TO_L4 (mbase));
 281     }
 282
 283   while (base)
 284     {
 285       const GFC_REAL_8 * restrict src;
 286       const GFC_LOGICAL_4 * restrict msrc;
 287       GFC_REAL_8 result;
 288       src = base;
 289       msrc = mbase;
 290       {
 291
 292   result = 0;
 293         if (len <= 0)
 294           *dest = 0;
 295         else
 296           {
 297             for (n = 0; n < len; n++, src += delta, msrc += mdelta)
 298               {
 299
 300   if (*msrc)
 301     result += *src;
 302               }
 303             *dest = result;
 304           }
 305       }
 306       /* Advance to the next element.  */
 307       count[0]++;
 308       base += sstride[0];
 309       mbase += mstride[0];
 310       dest += dstride[0];
 311       n = 0;
 312       while (count[n] == extent[n])
 313         {
 314           /* When we get to the end of a dimension, reset it and increment
 315              the next dimension.  */
 316           count[n] = 0;
 317           /* We could precalculate these products, but this is a less
 318              frequently used path so proabably not worth it.  */
 319           base -= sstride[n] * extent[n];
 320           mbase -= mstride[n] * extent[n];
 321           dest -= dstride[n] * extent[n];
 322           n++;
 323           if (n == rank)
 324             {
 325               /* Break out of the look.  */
 326               base = NULL;
 327               break;
 328             }
 329           else
 330             {
 331               count[n]++;
 332               base += sstride[n];
 333               mbase += mstride[n];
 334               dest += dstride[n];
 335             }
 336         }
 337     }
 338 }
 339
 340 #endif