libgfortran/intrinsics/pack_generic.c

   1 /* Generic implementation of the PACK intrinsic
   2    Copyright (C) 2002, 2004, 2005, 2006 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 "config.h"
  32 #include <stdlib.h>
  33 #include <assert.h>
  34 #include <string.h>
  35 #include "libgfortran.h"
  36
  37 /* PACK is specified as follows:
  38
  39    13.14.80 PACK (ARRAY, MASK, [VECTOR])
  40
  41    Description: Pack an array into an array of rank one under the
  42    control of a mask.
  43
  44    Class: Transformational function.
  45
  46    Arguments:
  47       ARRAY   may be of any type. It shall not be scalar.
  48       MASK    shall be of type LOGICAL. It shall be conformable with ARRAY.
  49       VECTOR  (optional) shall be of the same type and type parameters
  50               as ARRAY. VECTOR shall have at least as many elements as
  51               there are true elements in MASK. If MASK is a scalar
  52               with the value true, VECTOR shall have at least as many
  53               elements as there are in ARRAY.
  54
  55    Result Characteristics: The result is an array of rank one with the
  56    same type and type parameters as ARRAY. If VECTOR is present, the
  57    result size is that of VECTOR; otherwise, the result size is the
  58    number /t/ of true elements in MASK unless MASK is scalar with the
  59    value true, in which case the result size is the size of ARRAY.
  60
  61    Result Value: Element /i/ of the result is the element of ARRAY
  62    that corresponds to the /i/th true element of MASK, taking elements
  63    in array element order, for /i/ = 1, 2, ..., /t/. If VECTOR is
  64    present and has size /n/ > /t/, element /i/ of the result has the
  65    value VECTOR(/i/), for /i/ = /t/ + 1, ..., /n/.
  66
  67    Examples: The nonzero elements of an array M with the value
  68    | 0 0 0 |
  69    | 9 0 0 | may be "gathered" by the function PACK. The result of
  70    | 0 0 7 |
  71    PACK (M, MASK = M.NE.0) is [9,7] and the result of PACK (M, M.NE.0,
  72    VECTOR = (/ 2,4,6,8,10,12 /)) is [9,7,6,8,10,12].
  73
  74 There are two variants of the PACK intrinsic: one, where MASK is
  75 array valued, and the other one where MASK is scalar.  */
  76
  77 static void
  78 pack_internal (gfc_array_char *ret, const gfc_array_char *array,
  79                const gfc_array_l4 *mask, const gfc_array_char *vector,
  80                index_type size)
  81 {
  82   /* r.* indicates the return array.  */
  83   index_type rstride0;
  84   char *rptr;
  85   /* s.* indicates the source array.  */
  86   index_type sstride[GFC_MAX_DIMENSIONS];
  87   index_type sstride0;
  88   const char *sptr;
  89   /* m.* indicates the mask array.  */
  90   index_type mstride[GFC_MAX_DIMENSIONS];
  91   index_type mstride0;
  92   const GFC_LOGICAL_4 *mptr;
  93
  94   index_type count[GFC_MAX_DIMENSIONS];
  95   index_type extent[GFC_MAX_DIMENSIONS];
  96   int zero_sized;
  97   index_type n;
  98   index_type dim;
  99   index_type nelem;
 100
 101   dim = GFC_DESCRIPTOR_RANK (array);
 102   zero_sized = 0;
 103   for (n = 0; n < dim; n++)
 104     {
 105       count[n] = 0;
 106       extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;
 107       if (extent[n] <= 0)
 108        zero_sized = 1;
 109       sstride[n] = array->dim[n].stride * size;
 110       mstride[n] = mask->dim[n].stride;
 111     }
 112   if (sstride[0] == 0)
 113     sstride[0] = size;
 114   if (mstride[0] == 0)
 115     mstride[0] = 1;
 116
 117   sptr = array->data;
 118   mptr = mask->data;
 119
 120   /* Use the same loop for both logical types. */
 121   if (GFC_DESCRIPTOR_SIZE (mask) != 4)
 122     {
 123       if (GFC_DESCRIPTOR_SIZE (mask) != 8)
 124         runtime_error ("Funny sized logical array");
 125       for (n = 0; n < dim; n++)
 126         mstride[n] <<= 1;
 127       mptr = GFOR_POINTER_L8_TO_L4 (mptr);
 128     }
 129
 130   if (ret->data == NULL)
 131     {
 132       /* Allocate the memory for the result.  */
 133       int total;
 134
 135       if (vector != NULL)
 136         {
 137           /* The return array will have as many
 138              elements as there are in VECTOR.  */
 139           total = vector->dim[0].ubound + 1 - vector->dim[0].lbound;
 140         }
 141       else
 142         {
 143           /* We have to count the true elements in MASK.  */
 144
 145           /* TODO: We could speed up pack easily in the case of only
 146              few .TRUE. entries in MASK, by keeping track of where we
 147              would be in the source array during the initial traversal
 148              of MASK, and caching the pointers to those elements. Then,
 149              supposed the number of elements is small enough, we would
 150              only have to traverse the list, and copy those elements
 151              into the result array. In the case of datatypes which fit
 152              in one of the integer types we could also cache the
 153              value instead of a pointer to it.
 154              This approach might be bad from the point of view of
 155              cache behavior in the case where our cache is not big
 156              enough to hold all elements that have to be copied.  */
 157
 158           const GFC_LOGICAL_4 *m = mptr;
 159
 160           total = 0;
 161           if (zero_sized)
 162             m = NULL;
 163
 164           while (m)
 165             {
 166               /* Test this element.  */
 167               if (*m)
 168                 total++;
 169
 170               /* Advance to the next element.  */
 171               m += mstride[0];
 172               count[0]++;
 173               n = 0;
 174               while (count[n] == extent[n])
 175                 {
 176                   /* When we get to the end of a dimension, reset it
 177                      and increment the next dimension.  */
 178                   count[n] = 0;
 179                   /* We could precalculate this product, but this is a
 180                      less frequently used path so probably not worth
 181                      it.  */
 182                   m -= mstride[n] * extent[n];
 183                   n++;
 184                   if (n >= dim)
 185                     {
 186                       /* Break out of the loop.  */
 187                       m = NULL;
 188                       break;
 189                     }
 190                   else
 191                     {
 192                       count[n]++;
 193                       m += mstride[n];
 194                     }
 195                 }
 196             }
 197         }
 198
 199       /* Setup the array descriptor.  */
 200       ret->dim[0].lbound = 0;
 201       ret->dim[0].ubound = total - 1;
 202       ret->dim[0].stride = 1;
 203
 204       ret->offset = 0;
 205       if (total == 0)
 206         {
 207           /* In this case, nothing remains to be done.  */
 208           ret->data = internal_malloc_size (1);
 209           return;
 210         }
 211       else
 212         ret->data = internal_malloc_size (size * total);
 213     }
 214
 215   rstride0 = ret->dim[0].stride * size;
 216   if (rstride0 == 0)
 217     rstride0 = size;
 218   sstride0 = sstride[0];
 219   mstride0 = mstride[0];
 220   rptr = ret->data;
 221
 222   while (sptr && mptr)
 223     {
 224       /* Test this element.  */
 225       if (*mptr)
 226         {
 227           /* Add it.  */
 228           memcpy (rptr, sptr, size);
 229           rptr += rstride0;
 230         }
 231       /* Advance to the next element.  */
 232       sptr += sstride0;
 233       mptr += mstride0;
 234       count[0]++;
 235       n = 0;
 236       while (count[n] == extent[n])
 237         {
 238           /* When we get to the end of a dimension, reset it and increment
 239              the next dimension.  */
 240           count[n] = 0;
 241           /* We could precalculate these products, but this is a less
 242              frequently used path so probably not worth it.  */
 243           sptr -= sstride[n] * extent[n];
 244           mptr -= mstride[n] * extent[n];
 245           n++;
 246           if (n >= dim)
 247             {
 248               /* Break out of the loop.  */
 249               sptr = NULL;
 250               break;
 251             }
 252           else
 253             {
 254               count[n]++;
 255               sptr += sstride[n];
 256               mptr += mstride[n];
 257             }
 258         }
 259     }
 260
 261   /* Add any remaining elements from VECTOR.  */
 262   if (vector)
 263     {
 264       n = vector->dim[0].ubound + 1 - vector->dim[0].lbound;
 265       nelem = ((rptr - ret->data) / rstride0);
 266       if (n > nelem)
 267         {
 268           sstride0 = vector->dim[0].stride * size;
 269           if (sstride0 == 0)
 270             sstride0 = size;
 271
 272           sptr = vector->data + sstride0 * nelem;
 273           n -= nelem;
 274           while (n--)
 275             {
 276               memcpy (rptr, sptr, size);
 277               rptr += rstride0;
 278               sptr += sstride0;
 279             }
 280         }
 281     }
 282 }
 283
 284 extern void pack (gfc_array_char *, const gfc_array_char *,
 285                   const gfc_array_l4 *, const gfc_array_char *);
 286 export_proto(pack);
 287
 288 void
 289 pack (gfc_array_char *ret, const gfc_array_char *array,
 290       const gfc_array_l4 *mask, const gfc_array_char *vector)
 291 {
 292   pack_internal (ret, array, mask, vector, GFC_DESCRIPTOR_SIZE (array));
 293 }
 294
 295 extern void pack_char (gfc_array_char *, GFC_INTEGER_4, const gfc_array_char *,
 296                        const gfc_array_l4 *, const gfc_array_char *,
 297                        GFC_INTEGER_4, GFC_INTEGER_4);
 298 export_proto(pack_char);
 299
 300 void
 301 pack_char (gfc_array_char *ret,
 302            GFC_INTEGER_4 ret_length __attribute__((unused)),
 303            const gfc_array_char *array, const gfc_array_l4 *mask,
 304            const gfc_array_char *vector, GFC_INTEGER_4 array_length,
 305            GFC_INTEGER_4 vector_length __attribute__((unused)))
 306 {
 307   pack_internal (ret, array, mask, vector, array_length);
 308 }
 309
 310 static void
 311 pack_s_internal (gfc_array_char *ret, const gfc_array_char *array,
 312                  const GFC_LOGICAL_4 *mask, const gfc_array_char *vector,
 313                  index_type size)
 314 {
 315   /* r.* indicates the return array.  */
 316   index_type rstride0;
 317   char *rptr;
 318   /* s.* indicates the source array.  */
 319   index_type sstride[GFC_MAX_DIMENSIONS];
 320   index_type sstride0;
 321   const char *sptr;
 322
 323   index_type count[GFC_MAX_DIMENSIONS];
 324   index_type extent[GFC_MAX_DIMENSIONS];
 325   index_type n;
 326   index_type dim;
 327   index_type ssize;
 328   index_type nelem;
 329
 330   dim = GFC_DESCRIPTOR_RANK (array);
 331   ssize = 1;
 332   for (n = 0; n < dim; n++)
 333     {
 334       count[n] = 0;
 335       extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;
 336       sstride[n] = array->dim[n].stride * size;
 337       ssize *= extent[n];
 338     }
 339   if (sstride[0] == 0)
 340     sstride[0] = size;
 341
 342   sstride0 = sstride[0];
 343   sptr = array->data;
 344
 345   if (ret->data == NULL)
 346     {
 347       /* Allocate the memory for the result.  */
 348       int total;
 349
 350       if (vector != NULL)
 351         {
 352           /* The return array will have as many elements as there are
 353              in vector.  */
 354           total = vector->dim[0].ubound + 1 - vector->dim[0].lbound;
 355         }
 356       else
 357         {
 358           if (*mask)
 359             {
 360               /* The result array will have as many elements as the input
 361                  array.  */
 362               total = extent[0];
 363               for (n = 1; n < dim; n++)
 364                 total *= extent[n];
 365             }
 366           else
 367             /* The result array will be empty.  */
 368             total = 0;
 369         }
 370
 371       /* Setup the array descriptor.  */
 372       ret->dim[0].lbound = 0;
 373       ret->dim[0].ubound = total - 1;
 374       ret->dim[0].stride = 1;
 375       ret->offset = 0;
 376
 377       if (total == 0)
 378         {
 379           ret->data = internal_malloc_size (1);
 380           return;
 381         }
 382       else
 383         ret->data = internal_malloc_size (size * total);
 384     }
 385
 386   rstride0 = ret->dim[0].stride * size;
 387   if (rstride0 == 0)
 388     rstride0 = size;
 389   rptr = ret->data;
 390
 391   /* The remaining possibilities are now:
 392        If MASK is .TRUE., we have to copy the source array into the
 393      result array. We then have to fill it up with elements from VECTOR.
 394        If MASK is .FALSE., we have to copy VECTOR into the result
 395      array. If VECTOR were not present we would have already returned.  */
 396
 397   if (*mask && ssize != 0)
 398     {
 399       while (sptr)
 400         {
 401           /* Add this element.  */
 402           memcpy (rptr, sptr, size);
 403           rptr += rstride0;
 404
 405           /* Advance to the next element.  */
 406           sptr += sstride0;
 407           count[0]++;
 408           n = 0;
 409           while (count[n] == extent[n])
 410             {
 411               /* When we get to the end of a dimension, reset it and
 412                  increment the next dimension.  */
 413               count[n] = 0;
 414               /* We could precalculate these products, but this is a
 415                  less frequently used path so probably not worth it.  */
 416               sptr -= sstride[n] * extent[n];
 417               n++;
 418               if (n >= dim)
 419                 {
 420                   /* Break out of the loop.  */
 421                   sptr = NULL;
 422                   break;
 423                 }
 424               else
 425                 {
 426                   count[n]++;
 427                   sptr += sstride[n];
 428                 }
 429             }
 430         }
 431     }
 432
 433   /* Add any remaining elements from VECTOR.  */
 434   if (vector)
 435     {
 436       n = vector->dim[0].ubound + 1 - vector->dim[0].lbound;
 437       nelem = ((rptr - ret->data) / rstride0);
 438       if (n > nelem)
 439         {
 440           sstride0 = vector->dim[0].stride * size;
 441           if (sstride0 == 0)
 442             sstride0 = size;
 443
 444           sptr = vector->data + sstride0 * nelem;
 445           n -= nelem;
 446           while (n--)
 447             {
 448               memcpy (rptr, sptr, size);
 449               rptr += rstride0;
 450               sptr += sstride0;
 451             }
 452         }
 453     }
 454 }
 455
 456 extern void pack_s (gfc_array_char *ret, const gfc_array_char *array,
 457                     const GFC_LOGICAL_4 *, const gfc_array_char *);
 458 export_proto(pack_s);
 459
 460 void
 461 pack_s (gfc_array_char *ret, const gfc_array_char *array,
 462         const GFC_LOGICAL_4 *mask, const gfc_array_char *vector)
 463 {
 464   pack_s_internal (ret, array, mask, vector, GFC_DESCRIPTOR_SIZE (array));
 465 }
 466
 467 extern void pack_s_char (gfc_array_char *ret, GFC_INTEGER_4,
 468                          const gfc_array_char *array, const GFC_LOGICAL_4 *,
 469                          const gfc_array_char *, GFC_INTEGER_4,
 470                          GFC_INTEGER_4);
 471 export_proto(pack_s_char);
 472
 473 void
 474 pack_s_char (gfc_array_char *ret,
 475              GFC_INTEGER_4 ret_length __attribute__((unused)),
 476              const gfc_array_char *array, const GFC_LOGICAL_4 *mask,
 477              const gfc_array_char *vector, GFC_INTEGER_4 array_length,
 478              GFC_INTEGER_4 vector_length __attribute__((unused)))
 479 {
 480   pack_s_internal (ret, array, mask, vector, array_length);
 481 }