1 /* Implementation of the MINLOC intrinsic 2 Copyright (C) 2002-2024 Free Software Foundation, Inc. 3 Contributed by Paul Brook <paul (at) nowt.org> 4 5 This file is part of the GNU Fortran 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 3 of the License, or (at your option) any later version. 11 12 Libgfortran 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 General Public License for more details. 16 17 Under Section 7 of GPL version 3, you are granted additional 18 permissions described in the GCC Runtime Library Exception, version 19 3.1, as published by the Free Software Foundation. 20 21 You should have received a copy of the GNU General Public License and 22 a copy of the GCC Runtime Library Exception along with this program; 23 see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 24 <http://www.gnu.org/licenses/>. */ 25 26 #include "libgfortran.h" 27 #include <assert.h> 28 29 30 #if defined (HAVE_GFC_REAL_4) && defined (HAVE_GFC_INTEGER_16) 31 32 #define HAVE_BACK_ARG 1 33 34 35 extern void minloc1_16_r4 (gfc_array_i16 * const restrict, 36 gfc_array_r4 * const restrict, const index_type * const restrict, GFC_LOGICAL_4 back); 37 export_proto(minloc1_16_r4); 38 39 void 40 minloc1_16_r4 (gfc_array_i16 * const restrict retarray, 41 gfc_array_r4 * const restrict array, 42 const index_type * const restrict pdim, GFC_LOGICAL_4 back) 43 { 44 index_type count[GFC_MAX_DIMENSIONS]; 45 index_type extent[GFC_MAX_DIMENSIONS]; 46 index_type sstride[GFC_MAX_DIMENSIONS]; 47 index_type dstride[GFC_MAX_DIMENSIONS]; 48 const GFC_REAL_4 * restrict base; 49 GFC_INTEGER_16 * restrict dest; 50 index_type rank; 51 index_type n; 52 index_type len; 53 index_type delta; 54 index_type dim; 55 int continue_loop; 56 57 /* Make dim zero based to avoid confusion. */ 58 rank = GFC_DESCRIPTOR_RANK (array) - 1; 59 dim = (*pdim) - 1; 60 61 if (unlikely (dim < 0 || dim > rank)) 62 { 63 runtime_error ("Dim argument incorrect in MINLOC intrinsic: " 64 "is %ld, should be between 1 and %ld", 65 (long int) dim + 1, (long int) rank + 1); 66 } 67 68 len = GFC_DESCRIPTOR_EXTENT(array,dim); 69 if (len < 0) 70 len = 0; 71 delta = GFC_DESCRIPTOR_STRIDE(array,dim); 72 73 for (n = 0; n < dim; n++) 74 { 75 sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n); 76 extent[n] = GFC_DESCRIPTOR_EXTENT(array,n); 77 78 if (extent[n] < 0) 79 extent[n] = 0; 80 } 81 for (n = dim; n < rank; n++) 82 { 83 sstride[n] = GFC_DESCRIPTOR_STRIDE(array, n + 1); 84 extent[n] = GFC_DESCRIPTOR_EXTENT(array, n + 1); 85 86 if (extent[n] < 0) 87 extent[n] = 0; 88 } 89 90 if (retarray->base_addr == NULL) 91 { 92 size_t alloc_size, str; 93 94 for (n = 0; n < rank; n++) 95 { 96 if (n == 0) 97 str = 1; 98 else 99 str = GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1]; 100 101 GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str); 102 103 } 104 105 retarray->offset = 0; 106 retarray->dtype.rank = rank; 107 108 alloc_size = GFC_DESCRIPTOR_STRIDE(retarray,rank-1) * extent[rank-1]; 109 110 retarray->base_addr = xmallocarray (alloc_size, sizeof (GFC_INTEGER_16)); 111 if (alloc_size == 0) 112 return; 113 } 114 else 115 { 116 if (rank != GFC_DESCRIPTOR_RANK (retarray)) 117 runtime_error ("rank of return array incorrect in" 118 " MINLOC intrinsic: is %ld, should be %ld", 119 (long int) (GFC_DESCRIPTOR_RANK (retarray)), 120 (long int) rank); 121 122 if (unlikely (compile_options.bounds_check)) 123 bounds_ifunction_return ((array_t *) retarray, extent, 124 "return value", "MINLOC"); 125 } 126 127 for (n = 0; n < rank; n++) 128 { 129 count[n] = 0; 130 dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n); 131 if (extent[n] <= 0) 132 return; 133 } 134 135 base = array->base_addr; 136 dest = retarray->base_addr; 137 138 continue_loop = 1; 139 while (continue_loop) 140 { 141 const GFC_REAL_4 * restrict src; 142 GFC_INTEGER_16 result; 143 src = base; 144 { 145 146 GFC_REAL_4 minval; 147 #if defined (GFC_REAL_4_INFINITY) 148 minval = GFC_REAL_4_INFINITY; 149 #else 150 minval = GFC_REAL_4_HUGE; 151 #endif 152 result = 1; 153 if (len <= 0) 154 *dest = 0; 155 else 156 { 157 #if ! defined HAVE_BACK_ARG 158 for (n = 0; n < len; n++, src += delta) 159 { 160 #endif 161 162 #if defined (GFC_REAL_4_QUIET_NAN) 163 for (n = 0; n < len; n++, src += delta) 164 { 165 if (*src <= minval) 166 { 167 minval = *src; 168 result = (GFC_INTEGER_16)n + 1; 169 break; 170 } 171 } 172 #else 173 n = 0; 174 #endif 175 if (back) 176 for (; n < len; n++, src += delta) 177 { 178 if (unlikely (*src <= minval)) 179 { 180 minval = *src; 181 result = (GFC_INTEGER_16)n + 1; 182 } 183 } 184 else 185 for (; n < len; n++, src += delta) 186 { 187 if (unlikely (*src < minval)) 188 { 189 minval = *src; 190 result = (GFC_INTEGER_16) n + 1; 191 } 192 } 193 194 *dest = result; 195 } 196 } 197 /* Advance to the next element. */ 198 count[0]++; 199 base += sstride[0]; 200 dest += dstride[0]; 201 n = 0; 202 while (count[n] == extent[n]) 203 { 204 /* When we get to the end of a dimension, reset it and increment 205 the next dimension. */ 206 count[n] = 0; 207 /* We could precalculate these products, but this is a less 208 frequently used path so probably not worth it. */ 209 base -= sstride[n] * extent[n]; 210 dest -= dstride[n] * extent[n]; 211 n++; 212 if (n >= rank) 213 { 214 /* Break out of the loop. */ 215 continue_loop = 0; 216 break; 217 } 218 else 219 { 220 count[n]++; 221 base += sstride[n]; 222 dest += dstride[n]; 223 } 224 } 225 } 226 } 227 228 229 extern void mminloc1_16_r4 (gfc_array_i16 * const restrict, 230 gfc_array_r4 * const restrict, const index_type * const restrict, 231 gfc_array_l1 * const restrict, GFC_LOGICAL_4 back); 232 export_proto(mminloc1_16_r4); 233 234 void 235 mminloc1_16_r4 (gfc_array_i16 * const restrict retarray, 236 gfc_array_r4 * const restrict array, 237 const index_type * const restrict pdim, 238 gfc_array_l1 * const restrict mask, GFC_LOGICAL_4 back) 239 { 240 index_type count[GFC_MAX_DIMENSIONS]; 241 index_type extent[GFC_MAX_DIMENSIONS]; 242 index_type sstride[GFC_MAX_DIMENSIONS]; 243 index_type dstride[GFC_MAX_DIMENSIONS]; 244 index_type mstride[GFC_MAX_DIMENSIONS]; 245 GFC_INTEGER_16 * restrict dest; 246 const GFC_REAL_4 * restrict base; 247 const GFC_LOGICAL_1 * restrict mbase; 248 index_type rank; 249 index_type dim; 250 index_type n; 251 index_type len; 252 index_type delta; 253 index_type mdelta; 254 int mask_kind; 255 256 if (mask == NULL) 257 { 258 #ifdef HAVE_BACK_ARG 259 minloc1_16_r4 (retarray, array, pdim, back); 260 #else 261 minloc1_16_r4 (retarray, array, pdim); 262 #endif 263 return; 264 } 265 266 dim = (*pdim) - 1; 267 rank = GFC_DESCRIPTOR_RANK (array) - 1; 268 269 270 if (unlikely (dim < 0 || dim > rank)) 271 { 272 runtime_error ("Dim argument incorrect in MINLOC intrinsic: " 273 "is %ld, should be between 1 and %ld", 274 (long int) dim + 1, (long int) rank + 1); 275 } 276 277 len = GFC_DESCRIPTOR_EXTENT(array,dim); 278 if (len < 0) 279 len = 0; 280 281 mbase = mask->base_addr; 282 283 mask_kind = GFC_DESCRIPTOR_SIZE (mask); 284 285 if (mask_kind == 1 || mask_kind == 2 || mask_kind == 4 || mask_kind == 8 286 #ifdef HAVE_GFC_LOGICAL_16 287 || mask_kind == 16 288 #endif 289 ) 290 mbase = GFOR_POINTER_TO_L1 (mbase, mask_kind); 291 else 292 runtime_error ("Funny sized logical array"); 293 294 delta = GFC_DESCRIPTOR_STRIDE(array,dim); 295 mdelta = GFC_DESCRIPTOR_STRIDE_BYTES(mask,dim); 296 297 for (n = 0; n < dim; n++) 298 { 299 sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n); 300 mstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(mask,n); 301 extent[n] = GFC_DESCRIPTOR_EXTENT(array,n); 302 303 if (extent[n] < 0) 304 extent[n] = 0; 305 306 } 307 for (n = dim; n < rank; n++) 308 { 309 sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n + 1); 310 mstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(mask, n + 1); 311 extent[n] = GFC_DESCRIPTOR_EXTENT(array, n + 1); 312 313 if (extent[n] < 0) 314 extent[n] = 0; 315 } 316 317 if (retarray->base_addr == NULL) 318 { 319 size_t alloc_size, str; 320 321 for (n = 0; n < rank; n++) 322 { 323 if (n == 0) 324 str = 1; 325 else 326 str= GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1]; 327 328 GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str); 329 330 } 331 332 alloc_size = GFC_DESCRIPTOR_STRIDE(retarray,rank-1) * extent[rank-1]; 333 334 retarray->offset = 0; 335 retarray->dtype.rank = rank; 336 337 retarray->base_addr = xmallocarray (alloc_size, sizeof (GFC_INTEGER_16)); 338 if (alloc_size == 0) 339 return; 340 } 341 else 342 { 343 if (rank != GFC_DESCRIPTOR_RANK (retarray)) 344 runtime_error ("rank of return array incorrect in MINLOC intrinsic"); 345 346 if (unlikely (compile_options.bounds_check)) 347 { 348 bounds_ifunction_return ((array_t *) retarray, extent, 349 "return value", "MINLOC"); 350 bounds_equal_extents ((array_t *) mask, (array_t *) array, 351 "MASK argument", "MINLOC"); 352 } 353 } 354 355 for (n = 0; n < rank; n++) 356 { 357 count[n] = 0; 358 dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n); 359 if (extent[n] <= 0) 360 return; 361 } 362 363 dest = retarray->base_addr; 364 base = array->base_addr; 365 366 while (base) 367 { 368 const GFC_REAL_4 * restrict src; 369 const GFC_LOGICAL_1 * restrict msrc; 370 GFC_INTEGER_16 result; 371 src = base; 372 msrc = mbase; 373 { 374 375 GFC_REAL_4 minval; 376 #if defined (GFC_REAL_4_INFINITY) 377 minval = GFC_REAL_4_INFINITY; 378 #else 379 minval = GFC_REAL_4_HUGE; 380 #endif 381 #if defined (GFC_REAL_4_QUIET_NAN) 382 GFC_INTEGER_16 result2 = 0; 383 #endif 384 result = 0; 385 for (n = 0; n < len; n++, src += delta, msrc += mdelta) 386 { 387 388 if (*msrc) 389 { 390 #if defined (GFC_REAL_4_QUIET_NAN) 391 if (!result2) 392 result2 = (GFC_INTEGER_16)n + 1; 393 if (*src <= minval) 394 #endif 395 { 396 minval = *src; 397 result = (GFC_INTEGER_16)n + 1; 398 break; 399 } 400 } 401 } 402 #if defined (GFC_REAL_4_QUIET_NAN) 403 if (unlikely (n >= len)) 404 result = result2; 405 else 406 #endif 407 if (back) 408 for (; n < len; n++, src += delta, msrc += mdelta) 409 { 410 if (*msrc && unlikely (*src <= minval)) 411 { 412 minval = *src; 413 result = (GFC_INTEGER_16)n + 1; 414 } 415 } 416 else 417 for (; n < len; n++, src += delta, msrc += mdelta) 418 { 419 if (*msrc && unlikely (*src < minval)) 420 { 421 minval = *src; 422 result = (GFC_INTEGER_16) n + 1; 423 } 424 } 425 *dest = result; 426 } 427 /* Advance to the next element. */ 428 count[0]++; 429 base += sstride[0]; 430 mbase += mstride[0]; 431 dest += dstride[0]; 432 n = 0; 433 while (count[n] == extent[n]) 434 { 435 /* When we get to the end of a dimension, reset it and increment 436 the next dimension. */ 437 count[n] = 0; 438 /* We could precalculate these products, but this is a less 439 frequently used path so probably not worth it. */ 440 base -= sstride[n] * extent[n]; 441 mbase -= mstride[n] * extent[n]; 442 dest -= dstride[n] * extent[n]; 443 n++; 444 if (n >= rank) 445 { 446 /* Break out of the loop. */ 447 base = NULL; 448 break; 449 } 450 else 451 { 452 count[n]++; 453 base += sstride[n]; 454 mbase += mstride[n]; 455 dest += dstride[n]; 456 } 457 } 458 } 459 } 460 461 462 extern void sminloc1_16_r4 (gfc_array_i16 * const restrict, 463 gfc_array_r4 * const restrict, const index_type * const restrict, 464 GFC_LOGICAL_4 *, GFC_LOGICAL_4 back); 465 export_proto(sminloc1_16_r4); 466 467 void 468 sminloc1_16_r4 (gfc_array_i16 * const restrict retarray, 469 gfc_array_r4 * const restrict array, 470 const index_type * const restrict pdim, 471 GFC_LOGICAL_4 * mask, GFC_LOGICAL_4 back) 472 { 473 index_type count[GFC_MAX_DIMENSIONS]; 474 index_type extent[GFC_MAX_DIMENSIONS]; 475 index_type dstride[GFC_MAX_DIMENSIONS]; 476 GFC_INTEGER_16 * restrict dest; 477 index_type rank; 478 index_type n; 479 index_type dim; 480 481 482 if (mask == NULL || *mask) 483 { 484 #ifdef HAVE_BACK_ARG 485 minloc1_16_r4 (retarray, array, pdim, back); 486 #else 487 minloc1_16_r4 (retarray, array, pdim); 488 #endif 489 return; 490 } 491 /* Make dim zero based to avoid confusion. */ 492 dim = (*pdim) - 1; 493 rank = GFC_DESCRIPTOR_RANK (array) - 1; 494 495 if (unlikely (dim < 0 || dim > rank)) 496 { 497 runtime_error ("Dim argument incorrect in MINLOC intrinsic: " 498 "is %ld, should be between 1 and %ld", 499 (long int) dim + 1, (long int) rank + 1); 500 } 501 502 for (n = 0; n < dim; n++) 503 { 504 extent[n] = GFC_DESCRIPTOR_EXTENT(array,n); 505 506 if (extent[n] <= 0) 507 extent[n] = 0; 508 } 509 510 for (n = dim; n < rank; n++) 511 { 512 extent[n] = 513 GFC_DESCRIPTOR_EXTENT(array,n + 1); 514 515 if (extent[n] <= 0) 516 extent[n] = 0; 517 } 518 519 if (retarray->base_addr == NULL) 520 { 521 size_t alloc_size, str; 522 523 for (n = 0; n < rank; n++) 524 { 525 if (n == 0) 526 str = 1; 527 else 528 str = GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1]; 529 530 GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str); 531 532 } 533 534 retarray->offset = 0; 535 retarray->dtype.rank = rank; 536 537 alloc_size = GFC_DESCRIPTOR_STRIDE(retarray,rank-1) * extent[rank-1]; 538 539 retarray->base_addr = xmallocarray (alloc_size, sizeof (GFC_INTEGER_16)); 540 if (alloc_size == 0) 541 return; 542 } 543 else 544 { 545 if (rank != GFC_DESCRIPTOR_RANK (retarray)) 546 runtime_error ("rank of return array incorrect in" 547 " MINLOC intrinsic: is %ld, should be %ld", 548 (long int) (GFC_DESCRIPTOR_RANK (retarray)), 549 (long int) rank); 550 551 if (unlikely (compile_options.bounds_check)) 552 { 553 for (n=0; n < rank; n++) 554 { 555 index_type ret_extent; 556 557 ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,n); 558 if (extent[n] != ret_extent) 559 runtime_error ("Incorrect extent in return value of" 560 " MINLOC intrinsic in dimension %ld:" 561 " is %ld, should be %ld", (long int) n + 1, 562 (long int) ret_extent, (long int) extent[n]); 563 } 564 } 565 } 566 567 for (n = 0; n < rank; n++) 568 { 569 count[n] = 0; 570 dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n); 571 } 572 573 dest = retarray->base_addr; 574 575 while(1) 576 { 577 *dest = 0; 578 count[0]++; 579 dest += dstride[0]; 580 n = 0; 581 while (count[n] == extent[n]) 582 { 583 /* When we get to the end of a dimension, reset it and increment 584 the next dimension. */ 585 count[n] = 0; 586 /* We could precalculate these products, but this is a less 587 frequently used path so probably not worth it. */ 588 dest -= dstride[n] * extent[n]; 589 n++; 590 if (n >= rank) 591 return; 592 else 593 { 594 count[n]++; 595 dest += dstride[n]; 596 } 597 } 598 } 599 } 600 601 #endif 602