libgfortran/generated/spread_r8.c

   1 /* Special implementation of the SPREAD intrinsic
   2    Copyright 2008 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 2 of the License, or (at your option) any later version.
  12
  13 In addition to the permissions in the GNU General Public License, the
  14 Free Software Foundation gives you unlimited permission to link the
  15 compiled version of this file into combinations with other programs,
  16 and to distribute those combinations without any restriction coming
  17 from the use of this file.  (The General Public License restrictions
  18 do apply in other respects; for example, they cover modification of
  19 the file, and distribution when not linked into a combine
  20 executable.)
  21
  22 Ligbfortran is distributed in the hope that it will be useful,
  23 but WITHOUT ANY WARRANTY; without even the implied warranty of
  24 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  25 GNU General Public License for more details.
  26
  27 You should have received a copy of the GNU General Public
  28 License along with libgfortran; see the file COPYING.  If not,
  29 write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
  30 Boston, MA 02110-1301, USA.  */
  31
  32 #include "libgfortran.h"
  33 #include <stdlib.h>
  34 #include <assert.h>
  35 #include <string.h>
  36
  37
  38 #if defined (HAVE_GFC_REAL_8)
  39
  40 void
  41 spread_r8 (gfc_array_r8 *ret, const gfc_array_r8 *source,
  42                  const index_type along, const index_type pncopies)
  43 {
  44   /* r.* indicates the return array.  */
  45   index_type rstride[GFC_MAX_DIMENSIONS];
  46   index_type rstride0;
  47   index_type rdelta = 0;
  48   index_type rrank;
  49   index_type rs;
  50   GFC_REAL_8 *rptr;
  51   GFC_REAL_8 * restrict dest;
  52   /* s.* indicates the source array.  */
  53   index_type sstride[GFC_MAX_DIMENSIONS];
  54   index_type sstride0;
  55   index_type srank;
  56   const GFC_REAL_8 *sptr;
  57
  58   index_type count[GFC_MAX_DIMENSIONS];
  59   index_type extent[GFC_MAX_DIMENSIONS];
  60   index_type n;
  61   index_type dim;
  62   index_type ncopies;
  63
  64   srank = GFC_DESCRIPTOR_RANK(source);
  65
  66   rrank = srank + 1;
  67   if (rrank > GFC_MAX_DIMENSIONS)
  68     runtime_error ("return rank too large in spread()");
  69
  70   if (along > rrank)
  71       runtime_error ("dim outside of rank in spread()");
  72
  73   ncopies = pncopies;
  74
  75   if (ret->data == NULL)
  76     {
  77       /* The front end has signalled that we need to populate the
  78          return array descriptor.  */
  79       ret->dtype = (source->dtype & ~GFC_DTYPE_RANK_MASK) | rrank;
  80       dim = 0;
  81       rs = 1;
  82       for (n = 0; n < rrank; n++)
  83         {
  84           ret->dim[n].stride = rs;
  85           ret->dim[n].lbound = 0;
  86           if (n == along - 1)
  87             {
  88               ret->dim[n].ubound = ncopies - 1;
  89               rdelta = rs;
  90               rs *= ncopies;
  91             }
  92           else
  93             {
  94               count[dim] = 0;
  95               extent[dim] = source->dim[dim].ubound + 1
  96                 - source->dim[dim].lbound;
  97               sstride[dim] = source->dim[dim].stride;
  98               rstride[dim] = rs;
  99
 100               ret->dim[n].ubound = extent[dim]-1;
 101               rs *= extent[dim];
 102               dim++;
 103             }
 104         }
 105       ret->offset = 0;
 106       if (rs > 0)
 107         ret->data = internal_malloc_size (rs * sizeof(GFC_REAL_8));
 108       else
 109         {
 110           ret->data = internal_malloc_size (1);
 111           return;
 112         }
 113     }
 114   else
 115     {
 116       int zero_sized;
 117
 118       zero_sized = 0;
 119
 120       dim = 0;
 121       if (GFC_DESCRIPTOR_RANK(ret) != rrank)
 122         runtime_error ("rank mismatch in spread()");
 123
 124       if (compile_options.bounds_check)
 125         {
 126           for (n = 0; n < rrank; n++)
 127             {
 128               index_type ret_extent;
 129
 130               ret_extent = ret->dim[n].ubound + 1 - ret->dim[n].lbound;
 131               if (n == along - 1)
 132                 {
 133                   rdelta = ret->dim[n].stride;
 134
 135                   if (ret_extent != ncopies)
 136                     runtime_error("Incorrect extent in return value of SPREAD"
 137                                   " intrinsic in dimension %ld: is %ld,"
 138                                   " should be %ld", (long int) n+1,
 139                                   (long int) ret_extent, (long int) ncopies);
 140                 }
 141               else
 142                 {
 143                   count[dim] = 0;
 144                   extent[dim] = source->dim[dim].ubound + 1
 145                     - source->dim[dim].lbound;
 146                   if (ret_extent != extent[dim])
 147                     runtime_error("Incorrect extent in return value of SPREAD"
 148                                   " intrinsic in dimension %ld: is %ld,"
 149                                   " should be %ld", (long int) n+1,
 150                                   (long int) ret_extent,
 151                                   (long int) extent[dim]);
 152
 153                   if (extent[dim] <= 0)
 154                     zero_sized = 1;
 155                   sstride[dim] = source->dim[dim].stride;
 156                   rstride[dim] = ret->dim[n].stride;
 157                   dim++;
 158                 }
 159             }
 160         }
 161       else
 162         {
 163           for (n = 0; n < rrank; n++)
 164             {
 165               if (n == along - 1)
 166                 {
 167                   rdelta = ret->dim[n].stride;
 168                 }
 169               else
 170                 {
 171                   count[dim] = 0;
 172                   extent[dim] = source->dim[dim].ubound + 1
 173                     - source->dim[dim].lbound;
 174                   if (extent[dim] <= 0)
 175                     zero_sized = 1;
 176                   sstride[dim] = source->dim[dim].stride;
 177                   rstride[dim] = ret->dim[n].stride;
 178                   dim++;
 179                 }
 180             }
 181         }
 182
 183       if (zero_sized)
 184         return;
 185
 186       if (sstride[0] == 0)
 187         sstride[0] = 1;
 188     }
 189   sstride0 = sstride[0];
 190   rstride0 = rstride[0];
 191   rptr = ret->data;
 192   sptr = source->data;
 193
 194   while (sptr)
 195     {
 196       /* Spread this element.  */
 197       dest = rptr;
 198       for (n = 0; n < ncopies; n++)
 199         {
 200           *dest = *sptr;
 201           dest += rdelta;
 202         }
 203       /* Advance to the next element.  */
 204       sptr += sstride0;
 205       rptr += rstride0;
 206       count[0]++;
 207       n = 0;
 208       while (count[n] == extent[n])
 209         {
 210           /* When we get to the end of a dimension, reset it and increment
 211              the next dimension.  */
 212           count[n] = 0;
 213           /* We could precalculate these products, but this is a less
 214              frequently used path so probably not worth it.  */
 215           sptr -= sstride[n] * extent[n];
 216           rptr -= rstride[n] * extent[n];
 217           n++;
 218           if (n >= srank)
 219             {
 220               /* Break out of the loop.  */
 221               sptr = NULL;
 222               break;
 223             }
 224           else
 225             {
 226               count[n]++;
 227               sptr += sstride[n];
 228               rptr += rstride[n];
 229             }
 230         }
 231     }
 232 }
 233
 234 /* This version of spread_internal treats the special case of a scalar
 235    source.  This is much simpler than the more general case above.  */
 236
 237 void
 238 spread_scalar_r8 (gfc_array_r8 *ret, const GFC_REAL_8 *source,
 239                         const index_type along, const index_type pncopies)
 240 {
 241   int n;
 242   int ncopies = pncopies;
 243   GFC_REAL_8 * restrict dest;
 244   index_type stride;
 245
 246   if (GFC_DESCRIPTOR_RANK (ret) != 1)
 247     runtime_error ("incorrect destination rank in spread()");
 248
 249   if (along > 1)
 250     runtime_error ("dim outside of rank in spread()");
 251
 252   if (ret->data == NULL)
 253     {
 254       ret->data = internal_malloc_size (ncopies * sizeof (GFC_REAL_8));
 255       ret->offset = 0;
 256       ret->dim[0].stride = 1;
 257       ret->dim[0].lbound = 0;
 258       ret->dim[0].ubound = ncopies - 1;
 259     }
 260   else
 261     {
 262       if (ncopies - 1 > (ret->dim[0].ubound - ret->dim[0].lbound)
 263                            / ret->dim[0].stride)
 264         runtime_error ("dim too large in spread()");
 265     }
 266
 267   dest = ret->data;
 268   stride = ret->dim[0].stride;
 269
 270   for (n = 0; n < ncopies; n++)
 271     {
 272       *dest = *source;
 273       dest += stride;
 274     }
 275 }
 276
 277 #endif
 278