libgfortran/m4/pack.m4

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