1 /* Implementation of the MINLOC intrinsic
2 Copyright (C) 2002-2023 Free Software Foundation, Inc.
3 Contributed by Paul Brook <paul@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_10) && defined (HAVE_GFC_INTEGER_4)
31
32 #define HAVE_BACK_ARG 1
33
34
35 extern void minloc1_4_r10 (gfc_array_i4 * const restrict,
36 gfc_array_r10 * const restrict, const index_type * const restrict, GFC_LOGICAL_4 back);
37 export_proto(minloc1_4_r10);
38
39 void
40 minloc1_4_r10 (gfc_array_i4 * const restrict retarray,
41 gfc_array_r10 * 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_10 * restrict base;
49 GFC_INTEGER_4 * 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_4));
111 if (alloc_size == 0)
112 {
113 /* Make sure we have a zero-sized array. */
114 GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 1);
115 return;
116
117 }
118 }
119 else
120 {
121 if (rank != GFC_DESCRIPTOR_RANK (retarray))
122 runtime_error ("rank of return array incorrect in"
123 " MINLOC intrinsic: is %ld, should be %ld",
124 (long int) (GFC_DESCRIPTOR_RANK (retarray)),
125 (long int) rank);
126
127 if (unlikely (compile_options.bounds_check))
128 bounds_ifunction_return ((array_t *) retarray, extent,
129 "return value", "MINLOC");
130 }
131
132 for (n = 0; n < rank; n++)
133 {
134 count[n] = 0;
135 dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);
136 if (extent[n] <= 0)
137 return;
138 }
139
140 base = array->base_addr;
141 dest = retarray->base_addr;
142
143 continue_loop = 1;
144 while (continue_loop)
145 {
146 const GFC_REAL_10 * restrict src;
147 GFC_INTEGER_4 result;
148 src = base;
149 {
150
151 GFC_REAL_10 minval;
152 #if defined (GFC_REAL_10_INFINITY)
153 minval = GFC_REAL_10_INFINITY;
154 #else
155 minval = GFC_REAL_10_HUGE;
156 #endif
157 result = 1;
158 if (len <= 0)
159 *dest = 0;
160 else
161 {
162 #if ! defined HAVE_BACK_ARG
163 for (n = 0; n < len; n++, src += delta)
164 {
165 #endif
166
167 #if defined (GFC_REAL_10_QUIET_NAN)
168 for (n = 0; n < len; n++, src += delta)
169 {
170 if (*src <= minval)
171 {
172 minval = *src;
173 result = (GFC_INTEGER_4)n + 1;
174 break;
175 }
176 }
177 #else
178 n = 0;
179 #endif
180 if (back)
181 for (; n < len; n++, src += delta)
182 {
183 if (unlikely (*src <= minval))
184 {
185 minval = *src;
186 result = (GFC_INTEGER_4)n + 1;
187 }
188 }
189 else
190 for (; n < len; n++, src += delta)
191 {
192 if (unlikely (*src < minval))
193 {
194 minval = *src;
195 result = (GFC_INTEGER_4) n + 1;
196 }
197 }
198
199 *dest = result;
200 }
201 }
202 /* Advance to the next element. */
203 count[0]++;
204 base += sstride[0];
205 dest += dstride[0];
206 n = 0;
207 while (count[n] == extent[n])
208 {
209 /* When we get to the end of a dimension, reset it and increment
210 the next dimension. */
211 count[n] = 0;
212 /* We could precalculate these products, but this is a less
213 frequently used path so probably not worth it. */
214 base -= sstride[n] * extent[n];
215 dest -= dstride[n] * extent[n];
216 n++;
217 if (n >= rank)
218 {
219 /* Break out of the loop. */
220 continue_loop = 0;
221 break;
222 }
223 else
224 {
225 count[n]++;
226 base += sstride[n];
227 dest += dstride[n];
228 }
229 }
230 }
231 }
232
233
234 extern void mminloc1_4_r10 (gfc_array_i4 * const restrict,
235 gfc_array_r10 * const restrict, const index_type * const restrict,
236 gfc_array_l1 * const restrict, GFC_LOGICAL_4 back);
237 export_proto(mminloc1_4_r10);
238
239 void
240 mminloc1_4_r10 (gfc_array_i4 * const restrict retarray,
241 gfc_array_r10 * const restrict array,
242 const index_type * const restrict pdim,
243 gfc_array_l1 * const restrict mask, GFC_LOGICAL_4 back)
244 {
245 index_type count[GFC_MAX_DIMENSIONS];
246 index_type extent[GFC_MAX_DIMENSIONS];
247 index_type sstride[GFC_MAX_DIMENSIONS];
248 index_type dstride[GFC_MAX_DIMENSIONS];
249 index_type mstride[GFC_MAX_DIMENSIONS];
250 GFC_INTEGER_4 * restrict dest;
251 const GFC_REAL_10 * restrict base;
252 const GFC_LOGICAL_1 * restrict mbase;
253 index_type rank;
254 index_type dim;
255 index_type n;
256 index_type len;
257 index_type delta;
258 index_type mdelta;
259 int mask_kind;
260
261 if (mask == NULL)
262 {
263 #ifdef HAVE_BACK_ARG
264 minloc1_4_r10 (retarray, array, pdim, back);
265 #else
266 minloc1_4_r10 (retarray, array, pdim);
267 #endif
268 return;
269 }
270
271 dim = (*pdim) - 1;
272 rank = GFC_DESCRIPTOR_RANK (array) - 1;
273
274
275 if (unlikely (dim < 0 || dim > rank))
276 {
277 runtime_error ("Dim argument incorrect in MINLOC intrinsic: "
278 "is %ld, should be between 1 and %ld",
279 (long int) dim + 1, (long int) rank + 1);
280 }
281
282 len = GFC_DESCRIPTOR_EXTENT(array,dim);
283 if (len <= 0)
284 return;
285
286 mbase = mask->base_addr;
287
288 mask_kind = GFC_DESCRIPTOR_SIZE (mask);
289
290 if (mask_kind == 1 || mask_kind == 2 || mask_kind == 4 || mask_kind == 8
291 #ifdef HAVE_GFC_LOGICAL_16
292 || mask_kind == 16
293 #endif
294 )
295 mbase = GFOR_POINTER_TO_L1 (mbase, mask_kind);
296 else
297 runtime_error ("Funny sized logical array");
298
299 delta = GFC_DESCRIPTOR_STRIDE(array,dim);
300 mdelta = GFC_DESCRIPTOR_STRIDE_BYTES(mask,dim);
301
302 for (n = 0; n < dim; n++)
303 {
304 sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n);
305 mstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(mask,n);
306 extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);
307
308 if (extent[n] < 0)
309 extent[n] = 0;
310
311 }
312 for (n = dim; n < rank; n++)
313 {
314 sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n + 1);
315 mstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(mask, n + 1);
316 extent[n] = GFC_DESCRIPTOR_EXTENT(array, n + 1);
317
318 if (extent[n] < 0)
319 extent[n] = 0;
320 }
321
322 if (retarray->base_addr == NULL)
323 {
324 size_t alloc_size, str;
325
326 for (n = 0; n < rank; n++)
327 {
328 if (n == 0)
329 str = 1;
330 else
331 str= GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1];
332
333 GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str);
334
335 }
336
337 alloc_size = GFC_DESCRIPTOR_STRIDE(retarray,rank-1) * extent[rank-1];
338
339 retarray->offset = 0;
340 retarray->dtype.rank = rank;
341
342 if (alloc_size == 0)
343 {
344 /* Make sure we have a zero-sized array. */
345 GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 1);
346 return;
347 }
348 else
349 retarray->base_addr = xmallocarray (alloc_size, sizeof (GFC_INTEGER_4));
350
351 }
352 else
353 {
354 if (rank != GFC_DESCRIPTOR_RANK (retarray))
355 runtime_error ("rank of return array incorrect in MINLOC intrinsic");
356
357 if (unlikely (compile_options.bounds_check))
358 {
359 bounds_ifunction_return ((array_t *) retarray, extent,
360 "return value", "MINLOC");
361 bounds_equal_extents ((array_t *) mask, (array_t *) array,
362 "MASK argument", "MINLOC");
363 }
364 }
365
366 for (n = 0; n < rank; n++)
367 {
368 count[n] = 0;
369 dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);
370 if (extent[n] <= 0)
371 return;
372 }
373
374 dest = retarray->base_addr;
375 base = array->base_addr;
376
377 while (base)
378 {
379 const GFC_REAL_10 * restrict src;
380 const GFC_LOGICAL_1 * restrict msrc;
381 GFC_INTEGER_4 result;
382 src = base;
383 msrc = mbase;
384 {
385
386 GFC_REAL_10 minval;
387 #if defined (GFC_REAL_10_INFINITY)
388 minval = GFC_REAL_10_INFINITY;
389 #else
390 minval = GFC_REAL_10_HUGE;
391 #endif
392 #if defined (GFC_REAL_10_QUIET_NAN)
393 GFC_INTEGER_4 result2 = 0;
394 #endif
395 result = 0;
396 for (n = 0; n < len; n++, src += delta, msrc += mdelta)
397 {
398
399 if (*msrc)
400 {
401 #if defined (GFC_REAL_10_QUIET_NAN)
402 if (!result2)
403 result2 = (GFC_INTEGER_4)n + 1;
404 if (*src <= minval)
405 #endif
406 {
407 minval = *src;
408 result = (GFC_INTEGER_4)n + 1;
409 break;
410 }
411 }
412 }
413 #if defined (GFC_REAL_10_QUIET_NAN)
414 if (unlikely (n >= len))
415 result = result2;
416 else
417 #endif
418 if (back)
419 for (; n < len; n++, src += delta, msrc += mdelta)
420 {
421 if (*msrc && unlikely (*src <= minval))
422 {
423 minval = *src;
424 result = (GFC_INTEGER_4)n + 1;
425 }
426 }
427 else
428 for (; n < len; n++, src += delta, msrc += mdelta)
429 {
430 if (*msrc && unlikely (*src < minval))
431 {
432 minval = *src;
433 result = (GFC_INTEGER_4) n + 1;
434 }
435 }
436 *dest = result;
437 }
438 /* Advance to the next element. */
439 count[0]++;
440 base += sstride[0];
441 mbase += mstride[0];
442 dest += dstride[0];
443 n = 0;
444 while (count[n] == extent[n])
445 {
446 /* When we get to the end of a dimension, reset it and increment
447 the next dimension. */
448 count[n] = 0;
449 /* We could precalculate these products, but this is a less
450 frequently used path so probably not worth it. */
451 base -= sstride[n] * extent[n];
452 mbase -= mstride[n] * extent[n];
453 dest -= dstride[n] * extent[n];
454 n++;
455 if (n >= rank)
456 {
457 /* Break out of the loop. */
458 base = NULL;
459 break;
460 }
461 else
462 {
463 count[n]++;
464 base += sstride[n];
465 mbase += mstride[n];
466 dest += dstride[n];
467 }
468 }
469 }
470 }
471
472
473 extern void sminloc1_4_r10 (gfc_array_i4 * const restrict,
474 gfc_array_r10 * const restrict, const index_type * const restrict,
475 GFC_LOGICAL_4 *, GFC_LOGICAL_4 back);
476 export_proto(sminloc1_4_r10);
477
478 void
479 sminloc1_4_r10 (gfc_array_i4 * const restrict retarray,
480 gfc_array_r10 * const restrict array,
481 const index_type * const restrict pdim,
482 GFC_LOGICAL_4 * mask, GFC_LOGICAL_4 back)
483 {
484 index_type count[GFC_MAX_DIMENSIONS];
485 index_type extent[GFC_MAX_DIMENSIONS];
486 index_type dstride[GFC_MAX_DIMENSIONS];
487 GFC_INTEGER_4 * restrict dest;
488 index_type rank;
489 index_type n;
490 index_type dim;
491
492
493 if (mask == NULL || *mask)
494 {
495 #ifdef HAVE_BACK_ARG
496 minloc1_4_r10 (retarray, array, pdim, back);
497 #else
498 minloc1_4_r10 (retarray, array, pdim);
499 #endif
500 return;
501 }
502 /* Make dim zero based to avoid confusion. */
503 dim = (*pdim) - 1;
504 rank = GFC_DESCRIPTOR_RANK (array) - 1;
505
506 if (unlikely (dim < 0 || dim > rank))
507 {
508 runtime_error ("Dim argument incorrect in MINLOC intrinsic: "
509 "is %ld, should be between 1 and %ld",
510 (long int) dim + 1, (long int) rank + 1);
511 }
512
513 for (n = 0; n < dim; n++)
514 {
515 extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);
516
517 if (extent[n] <= 0)
518 extent[n] = 0;
519 }
520
521 for (n = dim; n < rank; n++)
522 {
523 extent[n] =
524 GFC_DESCRIPTOR_EXTENT(array,n + 1);
525
526 if (extent[n] <= 0)
527 extent[n] = 0;
528 }
529
530 if (retarray->base_addr == NULL)
531 {
532 size_t alloc_size, str;
533
534 for (n = 0; n < rank; n++)
535 {
536 if (n == 0)
537 str = 1;
538 else
539 str = GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1];
540
541 GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str);
542
543 }
544
545 retarray->offset = 0;
546 retarray->dtype.rank = rank;
547
548 alloc_size = GFC_DESCRIPTOR_STRIDE(retarray,rank-1) * extent[rank-1];
549
550 if (alloc_size == 0)
551 {
552 /* Make sure we have a zero-sized array. */
553 GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 1);
554 return;
555 }
556 else
557 retarray->base_addr = xmallocarray (alloc_size, sizeof (GFC_INTEGER_4));
558 }
559 else
560 {
561 if (rank != GFC_DESCRIPTOR_RANK (retarray))
562 runtime_error ("rank of return array incorrect in"
563 " MINLOC intrinsic: is %ld, should be %ld",
564 (long int) (GFC_DESCRIPTOR_RANK (retarray)),
565 (long int) rank);
566
567 if (unlikely (compile_options.bounds_check))
568 {
569 for (n=0; n < rank; n++)
570 {
571 index_type ret_extent;
572
573 ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,n);
574 if (extent[n] != ret_extent)
575 runtime_error ("Incorrect extent in return value of"
576 " MINLOC intrinsic in dimension %ld:"
577 " is %ld, should be %ld", (long int) n + 1,
578 (long int) ret_extent, (long int) extent[n]);
579 }
580 }
581 }
582
583 for (n = 0; n < rank; n++)
584 {
585 count[n] = 0;
586 dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);
587 }
588
589 dest = retarray->base_addr;
590
591 while(1)
592 {
593 *dest = 0;
594 count[0]++;
595 dest += dstride[0];
596 n = 0;
597 while (count[n] == extent[n])
598 {
599 /* When we get to the end of a dimension, reset it and increment
600 the next dimension. */
601 count[n] = 0;
602 /* We could precalculate these products, but this is a less
603 frequently used path so probably not worth it. */
604 dest -= dstride[n] * extent[n];
605 n++;
606 if (n >= rank)
607 return;
608 else
609 {
610 count[n]++;
611 dest += dstride[n];
612 }
613 }
614 }
615 }
616
617 #endif