libgfortran/generated/maxloc1_16_r8.c

   1 /* Implementation of the MAXLOC intrinsic
   2    Copyright 2002 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 Libgfortran 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 <float.h>
  35 #include <limits.h>
  36 #include "libgfortran.h"
  37
  38
  39 #if defined (HAVE_GFC_REAL_8) && defined (HAVE_GFC_INTEGER_16)
  40
  41
  42 extern void maxloc1_16_r8 (gfc_array_i16 * const restrict,
  43         gfc_array_r8 * const restrict, const index_type * const restrict);
  44 export_proto(maxloc1_16_r8);
  45
  46 void
  47 maxloc1_16_r8 (gfc_array_i16 * const restrict retarray,
  48         gfc_array_r8 * const restrict array,
  49         const index_type * const restrict pdim)
  50 {
  51   index_type count[GFC_MAX_DIMENSIONS];
  52   index_type extent[GFC_MAX_DIMENSIONS];
  53   index_type sstride[GFC_MAX_DIMENSIONS];
  54   index_type dstride[GFC_MAX_DIMENSIONS];
  55   const GFC_REAL_8 * restrict base;
  56   GFC_INTEGER_16 * restrict dest;
  57   index_type rank;
  58   index_type n;
  59   index_type len;
  60   index_type delta;
  61   index_type dim;
  62
  63   /* Make dim zero based to avoid confusion.  */
  64   dim = (*pdim) - 1;
  65   rank = GFC_DESCRIPTOR_RANK (array) - 1;
  66
  67   /* TODO:  It should be a front end job to correctly set the strides.  */
  68
  69   if (array->dim[0].stride == 0)
  70     array->dim[0].stride = 1;
  71
  72   len = array->dim[dim].ubound + 1 - array->dim[dim].lbound;
  73   delta = array->dim[dim].stride;
  74
  75   for (n = 0; n < dim; n++)
  76     {
  77       sstride[n] = array->dim[n].stride;
  78       extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;
  79     }
  80   for (n = dim; n < rank; n++)
  81     {
  82       sstride[n] = array->dim[n + 1].stride;
  83       extent[n] =
  84         array->dim[n + 1].ubound + 1 - array->dim[n + 1].lbound;
  85     }
  86
  87   if (retarray->data == NULL)
  88     {
  89       for (n = 0; n < rank; n++)
  90         {
  91           retarray->dim[n].lbound = 0;
  92           retarray->dim[n].ubound = extent[n]-1;
  93           if (n == 0)
  94             retarray->dim[n].stride = 1;
  95           else
  96             retarray->dim[n].stride = retarray->dim[n-1].stride * extent[n-1];
  97         }
  98
  99       retarray->data
 100          = internal_malloc_size (sizeof (GFC_INTEGER_16)
 101                                  * retarray->dim[rank-1].stride
 102                                  * extent[rank-1]);
 103       retarray->offset = 0;
 104       retarray->dtype = (array->dtype & ~GFC_DTYPE_RANK_MASK) | rank;
 105     }
 106   else
 107     {
 108       if (retarray->dim[0].stride == 0)
 109         retarray->dim[0].stride = 1;
 110
 111       if (rank != GFC_DESCRIPTOR_RANK (retarray))
 112         runtime_error ("rank of return array incorrect");
 113     }
 114
 115   for (n = 0; n < rank; n++)
 116     {
 117       count[n] = 0;
 118       dstride[n] = retarray->dim[n].stride;
 119       if (extent[n] <= 0)
 120         len = 0;
 121     }
 122
 123   base = array->data;
 124   dest = retarray->data;
 125
 126   while (base)
 127     {
 128       const GFC_REAL_8 * restrict src;
 129       GFC_INTEGER_16 result;
 130       src = base;
 131       {
 132
 133   GFC_REAL_8 maxval;
 134   maxval = -GFC_REAL_8_HUGE;
 135   result = 1;
 136         if (len <= 0)
 137           *dest = 0;
 138         else
 139           {
 140             for (n = 0; n < len; n++, src += delta)
 141               {
 142
 143   if (*src > maxval)
 144     {
 145       maxval = *src;
 146       result = (GFC_INTEGER_16)n + 1;
 147     }
 148           }
 149             *dest = result;
 150           }
 151       }
 152       /* Advance to the next element.  */
 153       count[0]++;
 154       base += sstride[0];
 155       dest += dstride[0];
 156       n = 0;
 157       while (count[n] == extent[n])
 158         {
 159           /* When we get to the end of a dimension, reset it and increment
 160              the next dimension.  */
 161           count[n] = 0;
 162           /* We could precalculate these products, but this is a less
 163              frequently used path so proabably not worth it.  */
 164           base -= sstride[n] * extent[n];
 165           dest -= dstride[n] * extent[n];
 166           n++;
 167           if (n == rank)
 168             {
 169               /* Break out of the look.  */
 170               base = NULL;
 171               break;
 172             }
 173           else
 174             {
 175               count[n]++;
 176               base += sstride[n];
 177               dest += dstride[n];
 178             }
 179         }
 180     }
 181 }
 182
 183
 184 extern void mmaxloc1_16_r8 (gfc_array_i16 * const restrict,
 185         gfc_array_r8 * const restrict, const index_type * const restrict,
 186         gfc_array_l4 * const restrict);
 187 export_proto(mmaxloc1_16_r8);
 188
 189 void
 190 mmaxloc1_16_r8 (gfc_array_i16 * const restrict retarray,
 191         gfc_array_r8 * const restrict array,
 192         const index_type * const restrict pdim,
 193         gfc_array_l4 * const restrict mask)
 194 {
 195   index_type count[GFC_MAX_DIMENSIONS];
 196   index_type extent[GFC_MAX_DIMENSIONS];
 197   index_type sstride[GFC_MAX_DIMENSIONS];
 198   index_type dstride[GFC_MAX_DIMENSIONS];
 199   index_type mstride[GFC_MAX_DIMENSIONS];
 200   GFC_INTEGER_16 * restrict dest;
 201   const GFC_REAL_8 * restrict base;
 202   const GFC_LOGICAL_4 * restrict mbase;
 203   int rank;
 204   int dim;
 205   index_type n;
 206   index_type len;
 207   index_type delta;
 208   index_type mdelta;
 209
 210   dim = (*pdim) - 1;
 211   rank = GFC_DESCRIPTOR_RANK (array) - 1;
 212
 213   /* TODO:  It should be a front end job to correctly set the strides.  */
 214
 215   if (array->dim[0].stride == 0)
 216     array->dim[0].stride = 1;
 217
 218   if (mask->dim[0].stride == 0)
 219     mask->dim[0].stride = 1;
 220
 221   len = array->dim[dim].ubound + 1 - array->dim[dim].lbound;
 222   if (len <= 0)
 223     return;
 224   delta = array->dim[dim].stride;
 225   mdelta = mask->dim[dim].stride;
 226
 227   for (n = 0; n < dim; n++)
 228     {
 229       sstride[n] = array->dim[n].stride;
 230       mstride[n] = mask->dim[n].stride;
 231       extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;
 232     }
 233   for (n = dim; n < rank; n++)
 234     {
 235       sstride[n] = array->dim[n + 1].stride;
 236       mstride[n] = mask->dim[n + 1].stride;
 237       extent[n] =
 238         array->dim[n + 1].ubound + 1 - array->dim[n + 1].lbound;
 239     }
 240
 241   if (retarray->data == NULL)
 242     {
 243       for (n = 0; n < rank; n++)
 244         {
 245           retarray->dim[n].lbound = 0;
 246           retarray->dim[n].ubound = extent[n]-1;
 247           if (n == 0)
 248             retarray->dim[n].stride = 1;
 249           else
 250             retarray->dim[n].stride = retarray->dim[n-1].stride * extent[n-1];
 251         }
 252
 253       retarray->data
 254          = internal_malloc_size (sizeof (GFC_INTEGER_16)
 255                                  * retarray->dim[rank-1].stride
 256                                  * extent[rank-1]);
 257       retarray->offset = 0;
 258       retarray->dtype = (array->dtype & ~GFC_DTYPE_RANK_MASK) | rank;
 259     }
 260   else
 261     {
 262       if (retarray->dim[0].stride == 0)
 263         retarray->dim[0].stride = 1;
 264
 265       if (rank != GFC_DESCRIPTOR_RANK (retarray))
 266         runtime_error ("rank of return array incorrect");
 267     }
 268
 269   for (n = 0; n < rank; n++)
 270     {
 271       count[n] = 0;
 272       dstride[n] = retarray->dim[n].stride;
 273       if (extent[n] <= 0)
 274         return;
 275     }
 276
 277   dest = retarray->data;
 278   base = array->data;
 279   mbase = mask->data;
 280
 281   if (GFC_DESCRIPTOR_SIZE (mask) != 4)
 282     {
 283       /* This allows the same loop to be used for all logical types.  */
 284       assert (GFC_DESCRIPTOR_SIZE (mask) == 8);
 285       for (n = 0; n < rank; n++)
 286         mstride[n] <<= 1;
 287       mdelta <<= 1;
 288       mbase = (GFOR_POINTER_L8_TO_L4 (mbase));
 289     }
 290
 291   while (base)
 292     {
 293       const GFC_REAL_8 * restrict src;
 294       const GFC_LOGICAL_4 * restrict msrc;
 295       GFC_INTEGER_16 result;
 296       src = base;
 297       msrc = mbase;
 298       {
 299
 300   GFC_REAL_8 maxval;
 301   maxval = -GFC_REAL_8_HUGE;
 302   result = 1;
 303         if (len <= 0)
 304           *dest = 0;
 305         else
 306           {
 307             for (n = 0; n < len; n++, src += delta, msrc += mdelta)
 308               {
 309
 310   if (*msrc && *src > maxval)
 311     {
 312       maxval = *src;
 313       result = (GFC_INTEGER_16)n + 1;
 314     }
 315               }
 316             *dest = result;
 317           }
 318       }
 319       /* Advance to the next element.  */
 320       count[0]++;
 321       base += sstride[0];
 322       mbase += mstride[0];
 323       dest += dstride[0];
 324       n = 0;
 325       while (count[n] == extent[n])
 326         {
 327           /* When we get to the end of a dimension, reset it and increment
 328              the next dimension.  */
 329           count[n] = 0;
 330           /* We could precalculate these products, but this is a less
 331              frequently used path so proabably not worth it.  */
 332           base -= sstride[n] * extent[n];
 333           mbase -= mstride[n] * extent[n];
 334           dest -= dstride[n] * extent[n];
 335           n++;
 336           if (n == rank)
 337             {
 338               /* Break out of the look.  */
 339               base = NULL;
 340               break;
 341             }
 342           else
 343             {
 344               count[n]++;
 345               base += sstride[n];
 346               mbase += mstride[n];
 347               dest += dstride[n];
 348             }
 349         }
 350     }
 351 }
 352
 353 #endif