libgfortran/generated/pack_r16.c

   1 /* Specific implementation of the PACK intrinsic
   2    Copyright (C) 2002, 2004, 2005, 2006, 2007, 2008 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 Ligbfortran 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 "libgfortran.h"
  32 #include <stdlib.h>
  33 #include <assert.h>
  34 #include <string.h>
  35
  36
  37 #if defined (HAVE_GFC_REAL_16)
  38
  39 /* PACK is specified as follows:
  40
  41    13.14.80 PACK (ARRAY, MASK, [VECTOR])
  42
  43    Description: Pack an array into an array of rank one under the
  44    control of a mask.
  45
  46    Class: Transformational function.
  47
  48    Arguments:
  49       ARRAY   may be of any type. It shall not be scalar.
  50       MASK    shall be of type LOGICAL. It shall be conformable with ARRAY.
  51       VECTOR  (optional) shall be of the same type and type parameters
  52               as ARRAY. VECTOR shall have at least as many elements as
  53               there are true elements in MASK. If MASK is a scalar
  54               with the value true, VECTOR shall have at least as many
  55               elements as there are in ARRAY.
  56
  57    Result Characteristics: The result is an array of rank one with the
  58    same type and type parameters as ARRAY. If VECTOR is present, the
  59    result size is that of VECTOR; otherwise, the result size is the
  60    number /t/ of true elements in MASK unless MASK is scalar with the
  61    value true, in which case the result size is the size of ARRAY.
  62
  63    Result Value: Element /i/ of the result is the element of ARRAY
  64    that corresponds to the /i/th true element of MASK, taking elements
  65    in array element order, for /i/ = 1, 2, ..., /t/. If VECTOR is
  66    present and has size /n/ > /t/, element /i/ of the result has the
  67    value VECTOR(/i/), for /i/ = /t/ + 1, ..., /n/.
  68
  69    Examples: The nonzero elements of an array M with the value
  70    | 0 0 0 |
  71    | 9 0 0 | may be "gathered" by the function PACK. The result of
  72    | 0 0 7 |
  73    PACK (M, MASK = M.NE.0) is [9,7] and the result of PACK (M, M.NE.0,
  74    VECTOR = (/ 2,4,6,8,10,12 /)) is [9,7,6,8,10,12].
  75
  76 There are two variants of the PACK intrinsic: one, where MASK is
  77 array valued, and the other one where MASK is scalar.  */
  78
  79 void
  80 pack_r16 (gfc_array_r16 *ret, const gfc_array_r16 *array,
  81                const gfc_array_l1 *mask, const gfc_array_r16 *vector)
  82 {
  83   /* r.* indicates the return array.  */
  84   index_type rstride0;
  85   GFC_REAL_16 * restrict rptr;
  86   /* s.* indicates the source array.  */
  87   index_type sstride[GFC_MAX_DIMENSIONS];
  88   index_type sstride0;
  89   const GFC_REAL_16 *sptr;
  90   /* m.* indicates the mask array.  */
  91   index_type mstride[GFC_MAX_DIMENSIONS];
  92   index_type mstride0;
  93   const GFC_LOGICAL_1 *mptr;
  94
  95   index_type count[GFC_MAX_DIMENSIONS];
  96   index_type extent[GFC_MAX_DIMENSIONS];
  97   int zero_sized;
  98   index_type n;
  99   index_type dim;
 100   index_type nelem;
 101   index_type total;
 102   int mask_kind;
 103
 104   dim = GFC_DESCRIPTOR_RANK (array);
 105
 106   mptr = mask->data;
 107
 108   /* Use the same loop for all logical types, by using GFC_LOGICAL_1
 109      and using shifting to address size and endian issues.  */
 110
 111   mask_kind = GFC_DESCRIPTOR_SIZE (mask);
 112
 113   if (mask_kind == 1 || mask_kind == 2 || mask_kind == 4 || mask_kind == 8
 114 #ifdef HAVE_GFC_LOGICAL_16
 115       || mask_kind == 16
 116 #endif
 117       )
 118     {
 119       /*  Do not convert a NULL pointer as we use test for NULL below.  */
 120       if (mptr)
 121         mptr = GFOR_POINTER_TO_L1 (mptr, mask_kind);
 122     }
 123   else
 124     runtime_error ("Funny sized logical array");
 125
 126   zero_sized = 0;
 127   for (n = 0; n < dim; n++)
 128     {
 129       count[n] = 0;
 130       extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;
 131       if (extent[n] <= 0)
 132        zero_sized = 1;
 133       sstride[n] = array->dim[n].stride;
 134       mstride[n] = mask->dim[n].stride * mask_kind;
 135     }
 136   if (sstride[0] == 0)
 137     sstride[0] = 1;
 138   if (mstride[0] == 0)
 139     mstride[0] = mask_kind;
 140
 141   if (zero_sized)
 142     sptr = NULL;
 143   else
 144     sptr = array->data;
 145
 146   if (ret->data == NULL || compile_options.bounds_check)
 147     {
 148       /* Count the elements, either for allocating memory or
 149          for bounds checking.  */
 150
 151       if (vector != NULL)
 152         {
 153           /* The return array will have as many
 154              elements as there are in VECTOR.  */
 155           total = vector->dim[0].ubound + 1 - vector->dim[0].lbound;
 156           if (total < 0)
 157             {
 158               total = 0;
 159               vector = NULL;
 160             }
 161         }
 162       else
 163         {
 164           /* We have to count the true elements in MASK.  */
 165
 166           /* TODO: We could speed up pack easily in the case of only
 167              few .TRUE. entries in MASK, by keeping track of where we
 168              would be in the source array during the initial traversal
 169              of MASK, and caching the pointers to those elements. Then,
 170              supposed the number of elements is small enough, we would
 171              only have to traverse the list, and copy those elements
 172              into the result array. In the case of datatypes which fit
 173              in one of the integer types we could also cache the
 174              value instead of a pointer to it.
 175              This approach might be bad from the point of view of
 176              cache behavior in the case where our cache is not big
 177              enough to hold all elements that have to be copied.  */
 178
 179           const GFC_LOGICAL_1 *m = mptr;
 180
 181           total = 0;
 182           if (zero_sized)
 183             m = NULL;
 184
 185           while (m)
 186             {
 187               /* Test this element.  */
 188               if (*m)
 189                 total++;
 190
 191               /* Advance to the next element.  */
 192               m += mstride[0];
 193               count[0]++;
 194               n = 0;
 195               while (count[n] == extent[n])
 196                 {
 197                   /* When we get to the end of a dimension, reset it
 198                      and increment the next dimension.  */
 199                   count[n] = 0;
 200                   /* We could precalculate this product, but this is a
 201                      less frequently used path so probably not worth
 202                      it.  */
 203                   m -= mstride[n] * extent[n];
 204                   n++;
 205                   if (n >= dim)
 206                     {
 207                       /* Break out of the loop.  */
 208                       m = NULL;
 209                       break;
 210                     }
 211                   else
 212                     {
 213                       count[n]++;
 214                       m += mstride[n];
 215                     }
 216                 }
 217             }
 218         }
 219
 220       if (ret->data == NULL)
 221         {
 222           /* Setup the array descriptor.  */
 223           ret->dim[0].lbound = 0;
 224           ret->dim[0].ubound = total - 1;
 225           ret->dim[0].stride = 1;
 226
 227           ret->offset = 0;
 228           if (total == 0)
 229             {
 230               /* In this case, nothing remains to be done.  */
 231               ret->data = internal_malloc_size (1);
 232               return;
 233             }
 234           else
 235             ret->data = internal_malloc_size (sizeof (GFC_REAL_16) * total);
 236         }
 237       else
 238         {
 239           /* We come here because of range checking.  */
 240           index_type ret_extent;
 241
 242           ret_extent = ret->dim[0].ubound + 1 - ret->dim[0].lbound;
 243           if (total != ret_extent)
 244             runtime_error ("Incorrect extent in return value of PACK intrinsic;"
 245                            " is %ld, should be %ld", (long int) total,
 246                            (long int) ret_extent);
 247         }
 248     }
 249
 250   rstride0 = ret->dim[0].stride;
 251   if (rstride0 == 0)
 252     rstride0 = 1;
 253   sstride0 = sstride[0];
 254   mstride0 = mstride[0];
 255   rptr = ret->data;
 256
 257   while (sptr && mptr)
 258     {
 259       /* Test this element.  */
 260       if (*mptr)
 261         {
 262           /* Add it.  */
 263           *rptr = *sptr;
 264           rptr += rstride0;
 265         }
 266       /* Advance to the next element.  */
 267       sptr += sstride0;
 268       mptr += mstride0;
 269       count[0]++;
 270       n = 0;
 271       while (count[n] == extent[n])
 272         {
 273           /* When we get to the end of a dimension, reset it and increment
 274              the next dimension.  */
 275           count[n] = 0;
 276           /* We could precalculate these products, but this is a less
 277              frequently used path so probably not worth it.  */
 278           sptr -= sstride[n] * extent[n];
 279           mptr -= mstride[n] * extent[n];
 280           n++;
 281           if (n >= dim)
 282             {
 283               /* Break out of the loop.  */
 284               sptr = NULL;
 285               break;
 286             }
 287           else
 288             {
 289               count[n]++;
 290               sptr += sstride[n];
 291               mptr += mstride[n];
 292             }
 293         }
 294     }
 295
 296   /* Add any remaining elements from VECTOR.  */
 297   if (vector)
 298     {
 299       n = vector->dim[0].ubound + 1 - vector->dim[0].lbound;
 300       nelem = ((rptr - ret->data) / rstride0);
 301       if (n > nelem)
 302         {
 303           sstride0 = vector->dim[0].stride;
 304           if (sstride0 == 0)
 305             sstride0 = 1;
 306
 307           sptr = vector->data + sstride0 * nelem;
 308           n -= nelem;
 309           while (n--)
 310             {
 311               *rptr = *sptr;
 312               rptr += rstride0;
 313               sptr += sstride0;
 314             }
 315         }
 316     }
 317 }
 318
 319 #endif
 320