libgfortran/generated/pack_c4.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_COMPLEX_4)
  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_c4 (gfc_array_c4 *ret, const gfc_array_c4 *array,
  81                const gfc_array_l1 *mask, const gfc_array_c4 *vector)
  82 {
  83   /* r.* indicates the return array.  */
  84   index_type rstride0;
  85   GFC_COMPLEX_4 *rptr;
  86   /* s.* indicates the source array.  */
  87   index_type sstride[GFC_MAX_DIMENSIONS];
  88   index_type sstride0;
  89   const GFC_COMPLEX_4 *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   sptr = array->data;
 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 (ret->data == NULL || compile_options.bounds_check)
 143     {
 144       /* Count the elements, either for allocating memory or
 145          for bounds checking.  */
 146
 147       if (vector != NULL)
 148         {
 149           /* The return array will have as many
 150              elements as there are in VECTOR.  */
 151           total = vector->dim[0].ubound + 1 - vector->dim[0].lbound;
 152         }
 153       else
 154         {
 155           /* We have to count the true elements in MASK.  */
 156
 157           /* TODO: We could speed up pack easily in the case of only
 158              few .TRUE. entries in MASK, by keeping track of where we
 159              would be in the source array during the initial traversal
 160              of MASK, and caching the pointers to those elements. Then,
 161              supposed the number of elements is small enough, we would
 162              only have to traverse the list, and copy those elements
 163              into the result array. In the case of datatypes which fit
 164              in one of the integer types we could also cache the
 165              value instead of a pointer to it.
 166              This approach might be bad from the point of view of
 167              cache behavior in the case where our cache is not big
 168              enough to hold all elements that have to be copied.  */
 169
 170           const GFC_LOGICAL_1 *m = mptr;
 171
 172           total = 0;
 173           if (zero_sized)
 174             m = NULL;
 175
 176           while (m)
 177             {
 178               /* Test this element.  */
 179               if (*m)
 180                 total++;
 181
 182               /* Advance to the next element.  */
 183               m += mstride[0];
 184               count[0]++;
 185               n = 0;
 186               while (count[n] == extent[n])
 187                 {
 188                   /* When we get to the end of a dimension, reset it
 189                      and increment the next dimension.  */
 190                   count[n] = 0;
 191                   /* We could precalculate this product, but this is a
 192                      less frequently used path so probably not worth
 193                      it.  */
 194                   m -= mstride[n] * extent[n];
 195                   n++;
 196                   if (n >= dim)
 197                     {
 198                       /* Break out of the loop.  */
 199                       m = NULL;
 200                       break;
 201                     }
 202                   else
 203                     {
 204                       count[n]++;
 205                       m += mstride[n];
 206                     }
 207                 }
 208             }
 209         }
 210
 211       if (ret->data == NULL)
 212         {
 213           /* Setup the array descriptor.  */
 214           ret->dim[0].lbound = 0;
 215           ret->dim[0].ubound = total - 1;
 216           ret->dim[0].stride = 1;
 217
 218           ret->offset = 0;
 219           if (total == 0)
 220             {
 221               /* In this case, nothing remains to be done.  */
 222               ret->data = internal_malloc_size (1);
 223               return;
 224             }
 225           else
 226             ret->data = internal_malloc_size (sizeof (GFC_COMPLEX_4) * total);
 227         }
 228       else
 229         {
 230           /* We come here because of range checking.  */
 231           index_type ret_extent;
 232
 233           ret_extent = ret->dim[0].ubound + 1 - ret->dim[0].lbound;
 234           if (total != ret_extent)
 235             runtime_error ("Incorrect extent in return value of PACK intrinsic;"
 236                            " is %ld, should be %ld", (long int) total,
 237                            (long int) ret_extent);
 238         }
 239     }
 240
 241   rstride0 = ret->dim[0].stride;
 242   if (rstride0 == 0)
 243     rstride0 = 1;
 244   sstride0 = sstride[0];
 245   mstride0 = mstride[0];
 246   rptr = ret->data;
 247
 248   while (sptr && mptr)
 249     {
 250       /* Test this element.  */
 251       if (*mptr)
 252         {
 253           /* Add it.  */
 254           *rptr = *sptr;
 255           rptr += rstride0;
 256         }
 257       /* Advance to the next element.  */
 258       sptr += sstride0;
 259       mptr += mstride0;
 260       count[0]++;
 261       n = 0;
 262       while (count[n] == extent[n])
 263         {
 264           /* When we get to the end of a dimension, reset it and increment
 265              the next dimension.  */
 266           count[n] = 0;
 267           /* We could precalculate these products, but this is a less
 268              frequently used path so probably not worth it.  */
 269           sptr -= sstride[n] * extent[n];
 270           mptr -= mstride[n] * extent[n];
 271           n++;
 272           if (n >= dim)
 273             {
 274               /* Break out of the loop.  */
 275               sptr = NULL;
 276               break;
 277             }
 278           else
 279             {
 280               count[n]++;
 281               sptr += sstride[n];
 282               mptr += mstride[n];
 283             }
 284         }
 285     }
 286
 287   /* Add any remaining elements from VECTOR.  */
 288   if (vector)
 289     {
 290       n = vector->dim[0].ubound + 1 - vector->dim[0].lbound;
 291       nelem = ((rptr - ret->data) / rstride0);
 292       if (n > nelem)
 293         {
 294           sstride0 = vector->dim[0].stride;
 295           if (sstride0 == 0)
 296             sstride0 = 1;
 297
 298           sptr = vector->data + sstride0 * nelem;
 299           n -= nelem;
 300           while (n--)
 301             {
 302               *rptr = *sptr;
 303               rptr += rstride0;
 304               sptr += sstride0;
 305             }
 306         }
 307     }
 308 }
 309
 310 #endif