libgfortran/m4/spread.m4

   1 `/* Special implementation of the SPREAD intrinsic
   2    Copyright 2008, 2009 Free Software Foundation, Inc.
   3    Contributed by Thomas Koenig <tkoenig@gcc.gnu.org>, based on
   4    spread_generic.c written by Paul Brook <paul@nowt.org>
   5
   6 This file is part of the GNU Fortran 95 runtime library (libgfortran).
   7
   8 Libgfortran is free software; you can redistribute it and/or
   9 modify it under the terms of the GNU General Public
  10 License as published by the Free Software Foundation; either
  11 version 3 of the License, or (at your option) any later version.
  12
  13 Ligbfortran is distributed in the hope that it will be useful,
  14 but WITHOUT ANY WARRANTY; without even the implied warranty of
  15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16 GNU General Public License for more details.
  17
  18 Under Section 7 of GPL version 3, you are granted additional
  19 permissions described in the GCC Runtime Library Exception, version
  20 3.1, as published by the Free Software Foundation.
  21
  22 You should have received a copy of the GNU General Public License and
  23 a copy of the GCC Runtime Library Exception along with this program;
  24 see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
  25 <http://www.gnu.org/licenses/>.  */
  26
  27 #include "libgfortran.h"
  28 #include <stdlib.h>
  29 #include <assert.h>
  30 #include <string.h>'
  31
  32 include(iparm.m4)dnl
  33
  34 `#if defined (HAVE_'rtype_name`)
  35
  36 void
  37 spread_'rtype_code` ('rtype` *ret, const 'rtype` *source,
  38                  const index_type along, const index_type pncopies)
  39 {
  40   /* r.* indicates the return array.  */
  41   index_type rstride[GFC_MAX_DIMENSIONS];
  42   index_type rstride0;
  43   index_type rdelta = 0;
  44   index_type rrank;
  45   index_type rs;
  46   'rtype_name` *rptr;
  47   'rtype_name` * restrict dest;
  48   /* s.* indicates the source array.  */
  49   index_type sstride[GFC_MAX_DIMENSIONS];
  50   index_type sstride0;
  51   index_type srank;
  52   const 'rtype_name` *sptr;
  53
  54   index_type count[GFC_MAX_DIMENSIONS];
  55   index_type extent[GFC_MAX_DIMENSIONS];
  56   index_type n;
  57   index_type dim;
  58   index_type ncopies;
  59
  60   srank = GFC_DESCRIPTOR_RANK(source);
  61
  62   rrank = srank + 1;
  63   if (rrank > GFC_MAX_DIMENSIONS)
  64     runtime_error ("return rank too large in spread()");
  65
  66   if (along > rrank)
  67       runtime_error ("dim outside of rank in spread()");
  68
  69   ncopies = pncopies;
  70
  71   if (ret->data == NULL)
  72     {
  73       /* The front end has signalled that we need to populate the
  74          return array descriptor.  */
  75       ret->dtype = (source->dtype & ~GFC_DTYPE_RANK_MASK) | rrank;
  76       dim = 0;
  77       rs = 1;
  78       for (n = 0; n < rrank; n++)
  79         {
  80           ret->dim[n].stride = rs;
  81           ret->dim[n].lbound = 0;
  82           if (n == along - 1)
  83             {
  84               ret->dim[n].ubound = ncopies - 1;
  85               rdelta = rs;
  86               rs *= ncopies;
  87             }
  88           else
  89             {
  90               count[dim] = 0;
  91               extent[dim] = source->dim[dim].ubound + 1
  92                 - source->dim[dim].lbound;
  93               sstride[dim] = source->dim[dim].stride;
  94               rstride[dim] = rs;
  95
  96               ret->dim[n].ubound = extent[dim]-1;
  97               rs *= extent[dim];
  98               dim++;
  99             }
 100         }
 101       ret->offset = 0;
 102       if (rs > 0)
 103         ret->data = internal_malloc_size (rs * sizeof('rtype_name`));
 104       else
 105         {
 106           ret->data = internal_malloc_size (1);
 107           return;
 108         }
 109     }
 110   else
 111     {
 112       int zero_sized;
 113
 114       zero_sized = 0;
 115
 116       dim = 0;
 117       if (GFC_DESCRIPTOR_RANK(ret) != rrank)
 118         runtime_error ("rank mismatch in spread()");
 119
 120       if (unlikely (compile_options.bounds_check))
 121         {
 122           for (n = 0; n < rrank; n++)
 123             {
 124               index_type ret_extent;
 125
 126               ret_extent = ret->dim[n].ubound + 1 - ret->dim[n].lbound;
 127               if (n == along - 1)
 128                 {
 129                   rdelta = ret->dim[n].stride;
 130
 131                   if (ret_extent != ncopies)
 132                     runtime_error("Incorrect extent in return value of SPREAD"
 133                                   " intrinsic in dimension %ld: is %ld,"
 134                                   " should be %ld", (long int) n+1,
 135                                   (long int) ret_extent, (long int) ncopies);
 136                 }
 137               else
 138                 {
 139                   count[dim] = 0;
 140                   extent[dim] = source->dim[dim].ubound + 1
 141                     - source->dim[dim].lbound;
 142                   if (ret_extent != extent[dim])
 143                     runtime_error("Incorrect extent in return value of SPREAD"
 144                                   " intrinsic in dimension %ld: is %ld,"
 145                                   " should be %ld", (long int) n+1,
 146                                   (long int) ret_extent,
 147                                   (long int) extent[dim]);
 148
 149                   if (extent[dim] <= 0)
 150                     zero_sized = 1;
 151                   sstride[dim] = source->dim[dim].stride;
 152                   rstride[dim] = ret->dim[n].stride;
 153                   dim++;
 154                 }
 155             }
 156         }
 157       else
 158         {
 159           for (n = 0; n < rrank; n++)
 160             {
 161               if (n == along - 1)
 162                 {
 163                   rdelta = ret->dim[n].stride;
 164                 }
 165               else
 166                 {
 167                   count[dim] = 0;
 168                   extent[dim] = source->dim[dim].ubound + 1
 169                     - source->dim[dim].lbound;
 170                   if (extent[dim] <= 0)
 171                     zero_sized = 1;
 172                   sstride[dim] = source->dim[dim].stride;
 173                   rstride[dim] = ret->dim[n].stride;
 174                   dim++;
 175                 }
 176             }
 177         }
 178
 179       if (zero_sized)
 180         return;
 181
 182       if (sstride[0] == 0)
 183         sstride[0] = 1;
 184     }
 185   sstride0 = sstride[0];
 186   rstride0 = rstride[0];
 187   rptr = ret->data;
 188   sptr = source->data;
 189
 190   while (sptr)
 191     {
 192       /* Spread this element.  */
 193       dest = rptr;
 194       for (n = 0; n < ncopies; n++)
 195         {
 196           *dest = *sptr;
 197           dest += rdelta;
 198         }
 199       /* Advance to the next element.  */
 200       sptr += sstride0;
 201       rptr += rstride0;
 202       count[0]++;
 203       n = 0;
 204       while (count[n] == extent[n])
 205         {
 206           /* When we get to the end of a dimension, reset it and increment
 207              the next dimension.  */
 208           count[n] = 0;
 209           /* We could precalculate these products, but this is a less
 210              frequently used path so probably not worth it.  */
 211           sptr -= sstride[n] * extent[n];
 212           rptr -= rstride[n] * extent[n];
 213           n++;
 214           if (n >= srank)
 215             {
 216               /* Break out of the loop.  */
 217               sptr = NULL;
 218               break;
 219             }
 220           else
 221             {
 222               count[n]++;
 223               sptr += sstride[n];
 224               rptr += rstride[n];
 225             }
 226         }
 227     }
 228 }
 229
 230 /* This version of spread_internal treats the special case of a scalar
 231    source.  This is much simpler than the more general case above.  */
 232
 233 void
 234 spread_scalar_'rtype_code` ('rtype` *ret, const 'rtype_name` *source,
 235                         const index_type along, const index_type pncopies)
 236 {
 237   int n;
 238   int ncopies = pncopies;
 239   'rtype_name` * restrict dest;
 240   index_type stride;
 241
 242   if (GFC_DESCRIPTOR_RANK (ret) != 1)
 243     runtime_error ("incorrect destination rank in spread()");
 244
 245   if (along > 1)
 246     runtime_error ("dim outside of rank in spread()");
 247
 248   if (ret->data == NULL)
 249     {
 250       ret->data = internal_malloc_size (ncopies * sizeof ('rtype_name`));
 251       ret->offset = 0;
 252       ret->dim[0].stride = 1;
 253       ret->dim[0].lbound = 0;
 254       ret->dim[0].ubound = ncopies - 1;
 255     }
 256   else
 257     {
 258       if (ncopies - 1 > (ret->dim[0].ubound - ret->dim[0].lbound)
 259                            / ret->dim[0].stride)
 260         runtime_error ("dim too large in spread()");
 261     }
 262
 263   dest = ret->data;
 264   stride = ret->dim[0].stride;
 265
 266   for (n = 0; n < ncopies; n++)
 267     {
 268       *dest = *source;
 269       dest += stride;
 270     }
 271 }
 272
 273 #endif
 274 '