123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556 |
- /* Implementation of the MINLOC intrinsic
- Copyright (C) 2002-2015 Free Software Foundation, Inc.
- Contributed by Paul Brook <paul@nowt.org>
- This file is part of the GNU Fortran runtime library (libgfortran).
- Libgfortran is free software; you can redistribute it and/or
- modify it under the terms of the GNU General Public
- License as published by the Free Software Foundation; either
- version 3 of the License, or (at your option) any later version.
- Libgfortran is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
- Under Section 7 of GPL version 3, you are granted additional
- permissions described in the GCC Runtime Library Exception, version
- 3.1, as published by the Free Software Foundation.
- You should have received a copy of the GNU General Public License and
- a copy of the GCC Runtime Library Exception along with this program;
- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
- <http://www.gnu.org/licenses/>. */
- #include "libgfortran.h"
- #include <stdlib.h>
- #include <assert.h>
- #include <limits.h>
- #if defined (HAVE_GFC_REAL_4) && defined (HAVE_GFC_INTEGER_16)
- extern void minloc1_16_r4 (gfc_array_i16 * const restrict,
- gfc_array_r4 * const restrict, const index_type * const restrict);
- export_proto(minloc1_16_r4);
- void
- minloc1_16_r4 (gfc_array_i16 * const restrict retarray,
- gfc_array_r4 * const restrict array,
- const index_type * const restrict pdim)
- {
- index_type count[GFC_MAX_DIMENSIONS];
- index_type extent[GFC_MAX_DIMENSIONS];
- index_type sstride[GFC_MAX_DIMENSIONS];
- index_type dstride[GFC_MAX_DIMENSIONS];
- const GFC_REAL_4 * restrict base;
- GFC_INTEGER_16 * restrict dest;
- index_type rank;
- index_type n;
- index_type len;
- index_type delta;
- index_type dim;
- int continue_loop;
- /* Make dim zero based to avoid confusion. */
- dim = (*pdim) - 1;
- rank = GFC_DESCRIPTOR_RANK (array) - 1;
- len = GFC_DESCRIPTOR_EXTENT(array,dim);
- if (len < 0)
- len = 0;
- delta = GFC_DESCRIPTOR_STRIDE(array,dim);
- for (n = 0; n < dim; n++)
- {
- sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n);
- extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);
- if (extent[n] < 0)
- extent[n] = 0;
- }
- for (n = dim; n < rank; n++)
- {
- sstride[n] = GFC_DESCRIPTOR_STRIDE(array, n + 1);
- extent[n] = GFC_DESCRIPTOR_EXTENT(array, n + 1);
- if (extent[n] < 0)
- extent[n] = 0;
- }
- if (retarray->base_addr == NULL)
- {
- size_t alloc_size, str;
- for (n = 0; n < rank; n++)
- {
- if (n == 0)
- str = 1;
- else
- str = GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1];
- GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str);
- }
- retarray->offset = 0;
- retarray->dtype = (array->dtype & ~GFC_DTYPE_RANK_MASK) | rank;
- alloc_size = GFC_DESCRIPTOR_STRIDE(retarray,rank-1) * extent[rank-1];
- retarray->base_addr = xmallocarray (alloc_size, sizeof (GFC_INTEGER_16));
- if (alloc_size == 0)
- {
- /* Make sure we have a zero-sized array. */
- GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 1);
- return;
- }
- }
- else
- {
- if (rank != GFC_DESCRIPTOR_RANK (retarray))
- runtime_error ("rank of return array incorrect in"
- " MINLOC intrinsic: is %ld, should be %ld",
- (long int) (GFC_DESCRIPTOR_RANK (retarray)),
- (long int) rank);
- if (unlikely (compile_options.bounds_check))
- bounds_ifunction_return ((array_t *) retarray, extent,
- "return value", "MINLOC");
- }
- for (n = 0; n < rank; n++)
- {
- count[n] = 0;
- dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);
- if (extent[n] <= 0)
- return;
- }
- base = array->base_addr;
- dest = retarray->base_addr;
- continue_loop = 1;
- while (continue_loop)
- {
- const GFC_REAL_4 * restrict src;
- GFC_INTEGER_16 result;
- src = base;
- {
- GFC_REAL_4 minval;
- #if defined (GFC_REAL_4_INFINITY)
- minval = GFC_REAL_4_INFINITY;
- #else
- minval = GFC_REAL_4_HUGE;
- #endif
- result = 1;
- if (len <= 0)
- *dest = 0;
- else
- {
- for (n = 0; n < len; n++, src += delta)
- {
- #if defined (GFC_REAL_4_QUIET_NAN)
- if (*src <= minval)
- {
- minval = *src;
- result = (GFC_INTEGER_16)n + 1;
- break;
- }
- }
- for (; n < len; n++, src += delta)
- {
- #endif
- if (*src < minval)
- {
- minval = *src;
- result = (GFC_INTEGER_16)n + 1;
- }
- }
-
- *dest = result;
- }
- }
- /* Advance to the next element. */
- count[0]++;
- base += sstride[0];
- dest += dstride[0];
- n = 0;
- while (count[n] == extent[n])
- {
- /* When we get to the end of a dimension, reset it and increment
- the next dimension. */
- count[n] = 0;
- /* We could precalculate these products, but this is a less
- frequently used path so probably not worth it. */
- base -= sstride[n] * extent[n];
- dest -= dstride[n] * extent[n];
- n++;
- if (n == rank)
- {
- /* Break out of the look. */
- continue_loop = 0;
- break;
- }
- else
- {
- count[n]++;
- base += sstride[n];
- dest += dstride[n];
- }
- }
- }
- }
- extern void mminloc1_16_r4 (gfc_array_i16 * const restrict,
- gfc_array_r4 * const restrict, const index_type * const restrict,
- gfc_array_l1 * const restrict);
- export_proto(mminloc1_16_r4);
- void
- mminloc1_16_r4 (gfc_array_i16 * const restrict retarray,
- gfc_array_r4 * const restrict array,
- const index_type * const restrict pdim,
- gfc_array_l1 * const restrict mask)
- {
- index_type count[GFC_MAX_DIMENSIONS];
- index_type extent[GFC_MAX_DIMENSIONS];
- index_type sstride[GFC_MAX_DIMENSIONS];
- index_type dstride[GFC_MAX_DIMENSIONS];
- index_type mstride[GFC_MAX_DIMENSIONS];
- GFC_INTEGER_16 * restrict dest;
- const GFC_REAL_4 * restrict base;
- const GFC_LOGICAL_1 * restrict mbase;
- int rank;
- int dim;
- index_type n;
- index_type len;
- index_type delta;
- index_type mdelta;
- int mask_kind;
- dim = (*pdim) - 1;
- rank = GFC_DESCRIPTOR_RANK (array) - 1;
- len = GFC_DESCRIPTOR_EXTENT(array,dim);
- if (len <= 0)
- return;
- mbase = mask->base_addr;
- mask_kind = GFC_DESCRIPTOR_SIZE (mask);
- if (mask_kind == 1 || mask_kind == 2 || mask_kind == 4 || mask_kind == 8
- #ifdef HAVE_GFC_LOGICAL_16
- || mask_kind == 16
- #endif
- )
- mbase = GFOR_POINTER_TO_L1 (mbase, mask_kind);
- else
- runtime_error ("Funny sized logical array");
- delta = GFC_DESCRIPTOR_STRIDE(array,dim);
- mdelta = GFC_DESCRIPTOR_STRIDE_BYTES(mask,dim);
- for (n = 0; n < dim; n++)
- {
- sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n);
- mstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(mask,n);
- extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);
- if (extent[n] < 0)
- extent[n] = 0;
- }
- for (n = dim; n < rank; n++)
- {
- sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n + 1);
- mstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(mask, n + 1);
- extent[n] = GFC_DESCRIPTOR_EXTENT(array, n + 1);
- if (extent[n] < 0)
- extent[n] = 0;
- }
- if (retarray->base_addr == NULL)
- {
- size_t alloc_size, str;
- for (n = 0; n < rank; n++)
- {
- if (n == 0)
- str = 1;
- else
- str= GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1];
- GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str);
- }
- alloc_size = GFC_DESCRIPTOR_STRIDE(retarray,rank-1) * extent[rank-1];
- retarray->offset = 0;
- retarray->dtype = (array->dtype & ~GFC_DTYPE_RANK_MASK) | rank;
- if (alloc_size == 0)
- {
- /* Make sure we have a zero-sized array. */
- GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 1);
- return;
- }
- else
- retarray->base_addr = xmallocarray (alloc_size, sizeof (GFC_INTEGER_16));
- }
- else
- {
- if (rank != GFC_DESCRIPTOR_RANK (retarray))
- runtime_error ("rank of return array incorrect in MINLOC intrinsic");
- if (unlikely (compile_options.bounds_check))
- {
- bounds_ifunction_return ((array_t *) retarray, extent,
- "return value", "MINLOC");
- bounds_equal_extents ((array_t *) mask, (array_t *) array,
- "MASK argument", "MINLOC");
- }
- }
- for (n = 0; n < rank; n++)
- {
- count[n] = 0;
- dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);
- if (extent[n] <= 0)
- return;
- }
- dest = retarray->base_addr;
- base = array->base_addr;
- while (base)
- {
- const GFC_REAL_4 * restrict src;
- const GFC_LOGICAL_1 * restrict msrc;
- GFC_INTEGER_16 result;
- src = base;
- msrc = mbase;
- {
- GFC_REAL_4 minval;
- #if defined (GFC_REAL_4_INFINITY)
- minval = GFC_REAL_4_INFINITY;
- #else
- minval = GFC_REAL_4_HUGE;
- #endif
- #if defined (GFC_REAL_4_QUIET_NAN)
- GFC_INTEGER_16 result2 = 0;
- #endif
- result = 0;
- for (n = 0; n < len; n++, src += delta, msrc += mdelta)
- {
- if (*msrc)
- {
- #if defined (GFC_REAL_4_QUIET_NAN)
- if (!result2)
- result2 = (GFC_INTEGER_16)n + 1;
- if (*src <= minval)
- #endif
- {
- minval = *src;
- result = (GFC_INTEGER_16)n + 1;
- break;
- }
- }
- }
- #if defined (GFC_REAL_4_QUIET_NAN)
- if (unlikely (n >= len))
- result = result2;
- else
- #endif
- for (; n < len; n++, src += delta, msrc += mdelta)
- {
- if (*msrc && *src < minval)
- {
- minval = *src;
- result = (GFC_INTEGER_16)n + 1;
- }
- }
- *dest = result;
- }
- /* Advance to the next element. */
- count[0]++;
- base += sstride[0];
- mbase += mstride[0];
- dest += dstride[0];
- n = 0;
- while (count[n] == extent[n])
- {
- /* When we get to the end of a dimension, reset it and increment
- the next dimension. */
- count[n] = 0;
- /* We could precalculate these products, but this is a less
- frequently used path so probably not worth it. */
- base -= sstride[n] * extent[n];
- mbase -= mstride[n] * extent[n];
- dest -= dstride[n] * extent[n];
- n++;
- if (n == rank)
- {
- /* Break out of the look. */
- base = NULL;
- break;
- }
- else
- {
- count[n]++;
- base += sstride[n];
- mbase += mstride[n];
- dest += dstride[n];
- }
- }
- }
- }
- extern void sminloc1_16_r4 (gfc_array_i16 * const restrict,
- gfc_array_r4 * const restrict, const index_type * const restrict,
- GFC_LOGICAL_4 *);
- export_proto(sminloc1_16_r4);
- void
- sminloc1_16_r4 (gfc_array_i16 * const restrict retarray,
- gfc_array_r4 * const restrict array,
- const index_type * const restrict pdim,
- GFC_LOGICAL_4 * mask)
- {
- index_type count[GFC_MAX_DIMENSIONS];
- index_type extent[GFC_MAX_DIMENSIONS];
- index_type dstride[GFC_MAX_DIMENSIONS];
- GFC_INTEGER_16 * restrict dest;
- index_type rank;
- index_type n;
- index_type dim;
- if (*mask)
- {
- minloc1_16_r4 (retarray, array, pdim);
- return;
- }
- /* Make dim zero based to avoid confusion. */
- dim = (*pdim) - 1;
- rank = GFC_DESCRIPTOR_RANK (array) - 1;
- for (n = 0; n < dim; n++)
- {
- extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);
- if (extent[n] <= 0)
- extent[n] = 0;
- }
- for (n = dim; n < rank; n++)
- {
- extent[n] =
- GFC_DESCRIPTOR_EXTENT(array,n + 1);
- if (extent[n] <= 0)
- extent[n] = 0;
- }
- if (retarray->base_addr == NULL)
- {
- size_t alloc_size, str;
- for (n = 0; n < rank; n++)
- {
- if (n == 0)
- str = 1;
- else
- str = GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1];
- GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str);
- }
- retarray->offset = 0;
- retarray->dtype = (array->dtype & ~GFC_DTYPE_RANK_MASK) | rank;
- alloc_size = GFC_DESCRIPTOR_STRIDE(retarray,rank-1) * extent[rank-1];
- if (alloc_size == 0)
- {
- /* Make sure we have a zero-sized array. */
- GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 1);
- return;
- }
- else
- retarray->base_addr = xmallocarray (alloc_size, sizeof (GFC_INTEGER_16));
- }
- else
- {
- if (rank != GFC_DESCRIPTOR_RANK (retarray))
- runtime_error ("rank of return array incorrect in"
- " MINLOC intrinsic: is %ld, should be %ld",
- (long int) (GFC_DESCRIPTOR_RANK (retarray)),
- (long int) rank);
- if (unlikely (compile_options.bounds_check))
- {
- for (n=0; n < rank; n++)
- {
- index_type ret_extent;
- ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,n);
- if (extent[n] != ret_extent)
- runtime_error ("Incorrect extent in return value of"
- " MINLOC intrinsic in dimension %ld:"
- " is %ld, should be %ld", (long int) n + 1,
- (long int) ret_extent, (long int) extent[n]);
- }
- }
- }
- for (n = 0; n < rank; n++)
- {
- count[n] = 0;
- dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);
- }
- dest = retarray->base_addr;
- while(1)
- {
- *dest = 0;
- count[0]++;
- dest += dstride[0];
- n = 0;
- while (count[n] == extent[n])
- {
- /* When we get to the end of a dimension, reset it and increment
- the next dimension. */
- count[n] = 0;
- /* We could precalculate these products, but this is a less
- frequently used path so probably not worth it. */
- dest -= dstride[n] * extent[n];
- n++;
- if (n == rank)
- return;
- else
- {
- count[n]++;
- dest += dstride[n];
- }
- }
- }
- }
- #endif
|