libgfortran/intrinsics/pack_generic.c

   1 /* Generic implementation of the PACK intrinsic
   2    Copyright (C) 2002, 2004 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 (libgfor).
   6
   7 Libgfor is free software; you can redistribute it and/or
   8 modify it under the terms of the GNU Lesser General Public
   9 License as published by the Free Software Foundation; either
  10 version 2.1 of the License, or (at your option) any later version.
  11
  12 Ligbfor is distributed in the hope that it will be useful,
  13 but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 GNU Lesser General Public License for more details.
  16
  17 You should have received a copy of the GNU Lesser General Public
  18 License along with libgfor; see the file COPYING.LIB.  If not,
  19 write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  20 Boston, MA 02111-1307, USA.  */
  21
  22 #include "config.h"
  23 #include <stdlib.h>
  24 #include <assert.h>
  25 #include <string.h>
  26 #include "libgfortran.h"
  27
  28 /* PACK is specified as follows:
  29
  30    13.14.80 PACK (ARRAY, MASK, [VECTOR])
  31
  32    Description: Pack an array into an array of rank one under the
  33    control of a mask.
  34
  35    Class: Transformational fucntion.
  36
  37    Arguments:
  38       ARRAY   may be of any type. It shall not be scalar.
  39       MASK    shall be of type LOGICAL. It shall be conformable with ARRAY.
  40       VECTOR  (optional) shall be of the same type and type parameters
  41               as ARRAY. VECTOR shall have at least as many elements as
  42               there are true elements in MASK. If MASK is a scalar
  43               with the value true, VECTOR shall have at least as many
  44               elements as there are in ARRAY.
  45
  46    Result Characteristics: The result is an array of rank one with the
  47    same type and type parameters as ARRAY. If VECTOR is present, the
  48    result size is that of VECTOR; otherwise, the result size is the
  49    number /t/ of true elements in MASK unless MASK is scalar with the
  50    value true, in which case the result size is the size of ARRAY.
  51
  52    Result Value: Element /i/ of the result is the element of ARRAY
  53    that corresponds to the /i/th true element of MASK, taking elements
  54    in array element order, for /i/ = 1, 2, ..., /t/. If VECTOR is
  55    present and has size /n/ > /t/, element /i/ of the result has the
  56    value VECTOR(/i/), for /i/ = /t/ + 1, ..., /n/.
  57
  58    Examples: The nonzero elements of an array M with the value
  59    | 0 0 0 |
  60    | 9 0 0 | may be "gathered" by the function PACK. The result of
  61    | 0 0 7 |
  62    PACK (M, MASK = M.NE.0) is [9,7] and the result of PACK (M, M.NE.0,
  63    VECTOR = (/ 2,4,6,8,10,12 /)) is [9,7,6,8,10,12].
  64
  65 There are two variants of the PACK intrinsic: one, where MASK is
  66 array valued, and the other one where MASK is scalar.  */
  67
  68 extern void __pack (gfc_array_char *, const gfc_array_char *,
  69                     const gfc_array_l4 *, const gfc_array_char *);
  70 export_proto_np(__pack);
  71
  72 void
  73 __pack (gfc_array_char * ret, const gfc_array_char * array,
  74         const gfc_array_l4 * mask, const gfc_array_char * vector)
  75 {
  76   /* r.* indicates the return array.  */
  77   index_type rstride0;
  78   char *rptr;
  79   /* s.* indicates the source array.  */
  80   index_type sstride[GFC_MAX_DIMENSIONS];
  81   index_type sstride0;
  82   const char *sptr;
  83   /* m.* indicates the mask array.  */
  84   index_type mstride[GFC_MAX_DIMENSIONS];
  85   index_type mstride0;
  86   const GFC_LOGICAL_4 *mptr;
  87
  88   index_type count[GFC_MAX_DIMENSIONS];
  89   index_type extent[GFC_MAX_DIMENSIONS];
  90   index_type n;
  91   index_type dim;
  92   index_type size;
  93   index_type nelem;
  94
  95   size = GFC_DESCRIPTOR_SIZE (array);
  96   dim = GFC_DESCRIPTOR_RANK (array);
  97   for (n = 0; n < dim; n++)
  98     {
  99       count[n] = 0;
 100       extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;
 101       sstride[n] = array->dim[n].stride * size;
 102       mstride[n] = mask->dim[n].stride;
 103     }
 104   if (sstride[0] == 0)
 105     sstride[0] = size;
 106   if (mstride[0] == 0)
 107     mstride[0] = 1;
 108
 109   sptr = array->data;
 110   mptr = mask->data;
 111
 112   /* Use the same loop for both logical types. */
 113   if (GFC_DESCRIPTOR_SIZE (mask) != 4)
 114     {
 115       if (GFC_DESCRIPTOR_SIZE (mask) != 8)
 116         runtime_error ("Funny sized logical array");
 117       for (n = 0; n < dim; n++)
 118         mstride[n] <<= 1;
 119       mstride0 <<= 1;
 120       mptr = GFOR_POINTER_L8_TO_L4 (mptr);
 121     }
 122
 123   if (ret->data == NULL)
 124     {
 125       /* Allocate the memory for the result.  */
 126       int total;
 127
 128       if (vector != NULL)
 129         {
 130           /* The return array will have as many
 131              elements as there are in VECTOR.  */
 132           total = vector->dim[0].ubound + 1 - vector->dim[0].lbound;
 133         }
 134       else
 135         {
 136           /* We have to count the true elements in MASK.  */
 137
 138           /* TODO: We could speed up pack easily in the case of only
 139              few .TRUE. entries in MASK, by keeping track of where we
 140              would be in the source array during the initial traversal
 141              of MASK, and caching the pointers to those elements. Then,
 142              supposed the number of elements is small enough, we would
 143              only have to traverse the list, and copy those elements
 144              into the result array. In the case of datatypes which fit
 145              in one of the integer types we could also cache the
 146              value instead of a pointer to it.
 147              This approach might be bad from the point of view of
 148              cache behavior in the case where our cache is not big
 149              enough to hold all elements that have to be copied.  */
 150
 151           const GFC_LOGICAL_4 *m = mptr;
 152
 153           total = 0;
 154
 155           while (m)
 156             {
 157               /* Test this element.  */
 158               if (*m)
 159                 total++;
 160
 161               /* Advance to the next element.  */
 162               m += mstride[0];
 163               count[0]++;
 164               n = 0;
 165               while (count[n] == extent[n])
 166                 {
 167                   /* When we get to the end of a dimension, reset it
 168                      and increment the next dimension.  */
 169                   count[n] = 0;
 170                   /* We could precalculate this product, but this is a
 171                      less frequently used path so proabably not worth
 172                      it.  */
 173                   m -= mstride[n] * extent[n];
 174                   n++;
 175                   if (n >= dim)
 176                     {
 177                       /* Break out of the loop.  */
 178                       m = NULL;
 179                       break;
 180                     }
 181                   else
 182                     {
 183                       count[n]++;
 184                       mptr += mstride[n];
 185                     }
 186                 }
 187             }
 188         }
 189
 190       /* Setup the array descriptor.  */
 191       ret->dim[0].lbound = 0;
 192       ret->dim[0].ubound = total - 1;
 193       ret->dim[0].stride = 1;
 194
 195       ret->data = internal_malloc_size (size * total);
 196       ret->base = 0;
 197
 198       if (total == 0)
 199         /* In this case, nothing remains to be done.  */
 200         return;
 201     }
 202
 203   rstride0 = ret->dim[0].stride * size;
 204   if (rstride0 == 0)
 205     rstride0 = size;
 206   sstride0 = sstride[0];
 207   mstride0 = mstride[0];
 208   rptr = ret->data;
 209
 210   while (sptr)
 211     {
 212       /* Test this element.  */
 213       if (*mptr)
 214         {
 215           /* Add it.  */
 216           memcpy (rptr, sptr, size);
 217           rptr += rstride0;
 218         }
 219       /* Advance to the next element.  */
 220       sptr += sstride0;
 221       mptr += mstride0;
 222       count[0]++;
 223       n = 0;
 224       while (count[n] == extent[n])
 225         {
 226           /* When we get to the end of a dimension, reset it and increment
 227              the next dimension.  */
 228           count[n] = 0;
 229           /* We could precalculate these products, but this is a less
 230              frequently used path so proabably not worth it.  */
 231           sptr -= sstride[n] * extent[n];
 232           mptr -= mstride[n] * extent[n];
 233           n++;
 234           if (n >= dim)
 235             {
 236               /* Break out of the loop.  */
 237               sptr = NULL;
 238               break;
 239             }
 240           else
 241             {
 242               count[n]++;
 243               sptr += sstride[n];
 244               mptr += mstride[n];
 245             }
 246         }
 247     }
 248
 249   /* Add any remaining elements from VECTOR.  */
 250   if (vector)
 251     {
 252       n = vector->dim[0].ubound + 1 - vector->dim[0].lbound;
 253       nelem = ((rptr - ret->data) / rstride0);
 254       if (n > nelem)
 255         {
 256           sstride0 = vector->dim[0].stride * size;
 257           if (sstride0 == 0)
 258             sstride0 = size;
 259
 260           sptr = vector->data + sstride0 * nelem;
 261           n -= nelem;
 262           while (n--)
 263             {
 264               memcpy (rptr, sptr, size);
 265               rptr += rstride0;
 266               sptr += sstride0;
 267             }
 268         }
 269     }
 270 }
 271
 272 extern void __pack_s (gfc_array_char *ret, const gfc_array_char *array,
 273                       const GFC_LOGICAL_4 *, const gfc_array_char *);
 274 export_proto_np(__pack_s);
 275
 276 void
 277 __pack_s (gfc_array_char * ret, const gfc_array_char * array,
 278           const GFC_LOGICAL_4 * mask, const gfc_array_char * vector)
 279 {
 280   /* r.* indicates the return array.  */
 281   index_type rstride0;
 282   char *rptr;
 283   /* s.* indicates the source array.  */
 284   index_type sstride[GFC_MAX_DIMENSIONS];
 285   index_type sstride0;
 286   const char *sptr;
 287
 288   index_type count[GFC_MAX_DIMENSIONS];
 289   index_type extent[GFC_MAX_DIMENSIONS];
 290   index_type n;
 291   index_type dim;
 292   index_type size;
 293   index_type nelem;
 294
 295   size = GFC_DESCRIPTOR_SIZE (array);
 296   dim = GFC_DESCRIPTOR_RANK (array);
 297   for (n = 0; n < dim; n++)
 298     {
 299       count[n] = 0;
 300       extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;
 301       sstride[n] = array->dim[n].stride * size;
 302     }
 303   if (sstride[0] == 0)
 304     sstride[0] = size;
 305
 306   sstride0 = sstride[0];
 307   sptr = array->data;
 308
 309   if (ret->data == NULL)
 310     {
 311       /* Allocate the memory for the result.  */
 312       int total;
 313
 314       if (vector != NULL)
 315         {
 316           /* The return array will have as many elements as there are
 317              in vector.  */
 318           total = vector->dim[0].ubound + 1 - vector->dim[0].lbound;
 319         }
 320       else
 321         {
 322           if (*mask)
 323             {
 324               /* The result array will have as many elements as the input
 325                  array.  */
 326               total = extent[0];
 327               for (n = 1; n < dim; n++)
 328                 total *= extent[n];
 329             }
 330           else
 331             {
 332               /* The result array will be empty.  */
 333               ret->dim[0].lbound = 0;
 334               ret->dim[0].ubound = -1;
 335               ret->dim[0].stride = 1;
 336               ret->data = internal_malloc_size (0);
 337               ret->base = 0;
 338
 339               return;
 340             }
 341         }
 342
 343       /* Setup the array descriptor.  */
 344       ret->dim[0].lbound = 0;
 345       ret->dim[0].ubound = total - 1;
 346       ret->dim[0].stride = 1;
 347
 348       ret->data = internal_malloc_size (size * total);
 349       ret->base = 0;
 350     }
 351
 352   rstride0 = ret->dim[0].stride * size;
 353   if (rstride0 == 0)
 354     rstride0 = size;
 355   rptr = ret->data;
 356
 357   /* The remaining possibilities are now:
 358        If MASK is .TRUE., we have to copy the source array into the
 359      result array. We then have to fill it up with elements from VECTOR.
 360        If MASK is .FALSE., we have to copy VECTOR into the result
 361      array. If VECTOR were not present we would have already returned.  */
 362
 363   if (*mask)
 364     {
 365       while (sptr)
 366         {
 367           /* Add this element.  */
 368           memcpy (rptr, sptr, size);
 369           rptr += rstride0;
 370
 371           /* Advance to the next element.  */
 372           sptr += sstride0;
 373           count[0]++;
 374           n = 0;
 375           while (count[n] == extent[n])
 376             {
 377               /* When we get to the end of a dimension, reset it and
 378                  increment the next dimension.  */
 379               count[n] = 0;
 380               /* We could precalculate these products, but this is a
 381                  less frequently used path so proabably not worth it.  */
 382               sptr -= sstride[n] * extent[n];
 383               n++;
 384               if (n >= dim)
 385                 {
 386                   /* Break out of the loop.  */
 387                   sptr = NULL;
 388                   break;
 389                 }
 390               else
 391                 {
 392                   count[n]++;
 393                   sptr += sstride[n];
 394                 }
 395             }
 396         }
 397     }
 398
 399   /* Add any remaining elements from VECTOR.  */
 400   if (vector)
 401     {
 402       n = vector->dim[0].ubound + 1 - vector->dim[0].lbound;
 403       nelem = ((rptr - ret->data) / rstride0);
 404       if (n > nelem)
 405         {
 406           sstride0 = vector->dim[0].stride * size;
 407           if (sstride0 == 0)
 408             sstride0 = size;
 409
 410           sptr = vector->data + sstride0 * nelem;
 411           n -= nelem;
 412           while (n--)
 413             {
 414               memcpy (rptr, sptr, size);
 415               rptr += rstride0;
 416               sptr += sstride0;
 417             }
 418         }
 419     }
 420 }