1 /* Implementation of the PRODUCT intrinsic
2 Copyright (C) 2002-2014 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 <stdlib.h>
28 #include <assert.h>
29
30
31 #if defined (HAVE_GFC_REAL_16) && defined (HAVE_GFC_REAL_16)
32
33
34 extern void product_r16 (gfc_array_r16 * const restrict,
35 gfc_array_r16 * const restrict, const index_type * const restrict);
36 export_proto(product_r16);
37
38 void
product_r16(gfc_array_r16 * const restrict retarray,gfc_array_r16 * const restrict array,const index_type * const restrict pdim)39 product_r16 (gfc_array_r16 * const restrict retarray,
40 gfc_array_r16 * const restrict array,
41 const index_type * const restrict pdim)
42 {
43 index_type count[GFC_MAX_DIMENSIONS];
44 index_type extent[GFC_MAX_DIMENSIONS];
45 index_type sstride[GFC_MAX_DIMENSIONS];
46 index_type dstride[GFC_MAX_DIMENSIONS];
47 const GFC_REAL_16 * restrict base;
48 GFC_REAL_16 * restrict dest;
49 index_type rank;
50 index_type n;
51 index_type len;
52 index_type delta;
53 index_type dim;
54 int continue_loop;
55
56 /* Make dim zero based to avoid confusion. */
57 dim = (*pdim) - 1;
58 rank = GFC_DESCRIPTOR_RANK (array) - 1;
59
60 len = GFC_DESCRIPTOR_EXTENT(array,dim);
61 if (len < 0)
62 len = 0;
63 delta = GFC_DESCRIPTOR_STRIDE(array,dim);
64
65 for (n = 0; n < dim; n++)
66 {
67 sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n);
68 extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);
69
70 if (extent[n] < 0)
71 extent[n] = 0;
72 }
73 for (n = dim; n < rank; n++)
74 {
75 sstride[n] = GFC_DESCRIPTOR_STRIDE(array, n + 1);
76 extent[n] = GFC_DESCRIPTOR_EXTENT(array, n + 1);
77
78 if (extent[n] < 0)
79 extent[n] = 0;
80 }
81
82 if (retarray->base_addr == NULL)
83 {
84 size_t alloc_size, str;
85
86 for (n = 0; n < rank; n++)
87 {
88 if (n == 0)
89 str = 1;
90 else
91 str = GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1];
92
93 GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str);
94
95 }
96
97 retarray->offset = 0;
98 retarray->dtype = (array->dtype & ~GFC_DTYPE_RANK_MASK) | rank;
99
100 alloc_size = GFC_DESCRIPTOR_STRIDE(retarray,rank-1) * extent[rank-1];
101
102 retarray->base_addr = xmallocarray (alloc_size, sizeof (GFC_REAL_16));
103 if (alloc_size == 0)
104 {
105 /* Make sure we have a zero-sized array. */
106 GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 1);
107 return;
108
109 }
110 }
111 else
112 {
113 if (rank != GFC_DESCRIPTOR_RANK (retarray))
114 runtime_error ("rank of return array incorrect in"
115 " PRODUCT intrinsic: is %ld, should be %ld",
116 (long int) (GFC_DESCRIPTOR_RANK (retarray)),
117 (long int) rank);
118
119 if (unlikely (compile_options.bounds_check))
120 bounds_ifunction_return ((array_t *) retarray, extent,
121 "return value", "PRODUCT");
122 }
123
124 for (n = 0; n < rank; n++)
125 {
126 count[n] = 0;
127 dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);
128 if (extent[n] <= 0)
129 return;
130 }
131
132 base = array->base_addr;
133 dest = retarray->base_addr;
134
135 continue_loop = 1;
136 while (continue_loop)
137 {
138 const GFC_REAL_16 * restrict src;
139 GFC_REAL_16 result;
140 src = base;
141 {
142
143 result = 1;
144 if (len <= 0)
145 *dest = 1;
146 else
147 {
148 for (n = 0; n < len; n++, src += delta)
149 {
150
151 result *= *src;
152 }
153
154 *dest = result;
155 }
156 }
157 /* Advance to the next element. */
158 count[0]++;
159 base += sstride[0];
160 dest += dstride[0];
161 n = 0;
162 while (count[n] == extent[n])
163 {
164 /* When we get to the end of a dimension, reset it and increment
165 the next dimension. */
166 count[n] = 0;
167 /* We could precalculate these products, but this is a less
168 frequently used path so probably not worth it. */
169 base -= sstride[n] * extent[n];
170 dest -= dstride[n] * extent[n];
171 n++;
172 if (n == rank)
173 {
174 /* Break out of the look. */
175 continue_loop = 0;
176 break;
177 }
178 else
179 {
180 count[n]++;
181 base += sstride[n];
182 dest += dstride[n];
183 }
184 }
185 }
186 }
187
188
189 extern void mproduct_r16 (gfc_array_r16 * const restrict,
190 gfc_array_r16 * const restrict, const index_type * const restrict,
191 gfc_array_l1 * const restrict);
192 export_proto(mproduct_r16);
193
194 void
mproduct_r16(gfc_array_r16 * const restrict retarray,gfc_array_r16 * const restrict array,const index_type * const restrict pdim,gfc_array_l1 * const restrict mask)195 mproduct_r16 (gfc_array_r16 * const restrict retarray,
196 gfc_array_r16 * const restrict array,
197 const index_type * const restrict pdim,
198 gfc_array_l1 * const restrict mask)
199 {
200 index_type count[GFC_MAX_DIMENSIONS];
201 index_type extent[GFC_MAX_DIMENSIONS];
202 index_type sstride[GFC_MAX_DIMENSIONS];
203 index_type dstride[GFC_MAX_DIMENSIONS];
204 index_type mstride[GFC_MAX_DIMENSIONS];
205 GFC_REAL_16 * restrict dest;
206 const GFC_REAL_16 * restrict base;
207 const GFC_LOGICAL_1 * restrict mbase;
208 int rank;
209 int dim;
210 index_type n;
211 index_type len;
212 index_type delta;
213 index_type mdelta;
214 int mask_kind;
215
216 dim = (*pdim) - 1;
217 rank = GFC_DESCRIPTOR_RANK (array) - 1;
218
219 len = GFC_DESCRIPTOR_EXTENT(array,dim);
220 if (len <= 0)
221 return;
222
223 mbase = mask->base_addr;
224
225 mask_kind = GFC_DESCRIPTOR_SIZE (mask);
226
227 if (mask_kind == 1 || mask_kind == 2 || mask_kind == 4 || mask_kind == 8
228 #ifdef HAVE_GFC_LOGICAL_16
229 || mask_kind == 16
230 #endif
231 )
232 mbase = GFOR_POINTER_TO_L1 (mbase, mask_kind);
233 else
234 runtime_error ("Funny sized logical array");
235
236 delta = GFC_DESCRIPTOR_STRIDE(array,dim);
237 mdelta = GFC_DESCRIPTOR_STRIDE_BYTES(mask,dim);
238
239 for (n = 0; n < dim; n++)
240 {
241 sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n);
242 mstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(mask,n);
243 extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);
244
245 if (extent[n] < 0)
246 extent[n] = 0;
247
248 }
249 for (n = dim; n < rank; n++)
250 {
251 sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n + 1);
252 mstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(mask, n + 1);
253 extent[n] = GFC_DESCRIPTOR_EXTENT(array, n + 1);
254
255 if (extent[n] < 0)
256 extent[n] = 0;
257 }
258
259 if (retarray->base_addr == NULL)
260 {
261 size_t alloc_size, str;
262
263 for (n = 0; n < rank; n++)
264 {
265 if (n == 0)
266 str = 1;
267 else
268 str= GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1];
269
270 GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str);
271
272 }
273
274 alloc_size = GFC_DESCRIPTOR_STRIDE(retarray,rank-1) * extent[rank-1];
275
276 retarray->offset = 0;
277 retarray->dtype = (array->dtype & ~GFC_DTYPE_RANK_MASK) | rank;
278
279 if (alloc_size == 0)
280 {
281 /* Make sure we have a zero-sized array. */
282 GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 1);
283 return;
284 }
285 else
286 retarray->base_addr = xmallocarray (alloc_size, sizeof (GFC_REAL_16));
287
288 }
289 else
290 {
291 if (rank != GFC_DESCRIPTOR_RANK (retarray))
292 runtime_error ("rank of return array incorrect in PRODUCT intrinsic");
293
294 if (unlikely (compile_options.bounds_check))
295 {
296 bounds_ifunction_return ((array_t *) retarray, extent,
297 "return value", "PRODUCT");
298 bounds_equal_extents ((array_t *) mask, (array_t *) array,
299 "MASK argument", "PRODUCT");
300 }
301 }
302
303 for (n = 0; n < rank; n++)
304 {
305 count[n] = 0;
306 dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);
307 if (extent[n] <= 0)
308 return;
309 }
310
311 dest = retarray->base_addr;
312 base = array->base_addr;
313
314 while (base)
315 {
316 const GFC_REAL_16 * restrict src;
317 const GFC_LOGICAL_1 * restrict msrc;
318 GFC_REAL_16 result;
319 src = base;
320 msrc = mbase;
321 {
322
323 result = 1;
324 for (n = 0; n < len; n++, src += delta, msrc += mdelta)
325 {
326
327 if (*msrc)
328 result *= *src;
329 }
330 *dest = result;
331 }
332 /* Advance to the next element. */
333 count[0]++;
334 base += sstride[0];
335 mbase += mstride[0];
336 dest += dstride[0];
337 n = 0;
338 while (count[n] == extent[n])
339 {
340 /* When we get to the end of a dimension, reset it and increment
341 the next dimension. */
342 count[n] = 0;
343 /* We could precalculate these products, but this is a less
344 frequently used path so probably not worth it. */
345 base -= sstride[n] * extent[n];
346 mbase -= mstride[n] * extent[n];
347 dest -= dstride[n] * extent[n];
348 n++;
349 if (n == rank)
350 {
351 /* Break out of the look. */
352 base = NULL;
353 break;
354 }
355 else
356 {
357 count[n]++;
358 base += sstride[n];
359 mbase += mstride[n];
360 dest += dstride[n];
361 }
362 }
363 }
364 }
365
366
367 extern void sproduct_r16 (gfc_array_r16 * const restrict,
368 gfc_array_r16 * const restrict, const index_type * const restrict,
369 GFC_LOGICAL_4 *);
370 export_proto(sproduct_r16);
371
372 void
sproduct_r16(gfc_array_r16 * const restrict retarray,gfc_array_r16 * const restrict array,const index_type * const restrict pdim,GFC_LOGICAL_4 * mask)373 sproduct_r16 (gfc_array_r16 * const restrict retarray,
374 gfc_array_r16 * const restrict array,
375 const index_type * const restrict pdim,
376 GFC_LOGICAL_4 * mask)
377 {
378 index_type count[GFC_MAX_DIMENSIONS];
379 index_type extent[GFC_MAX_DIMENSIONS];
380 index_type dstride[GFC_MAX_DIMENSIONS];
381 GFC_REAL_16 * restrict dest;
382 index_type rank;
383 index_type n;
384 index_type dim;
385
386
387 if (*mask)
388 {
389 product_r16 (retarray, array, pdim);
390 return;
391 }
392 /* Make dim zero based to avoid confusion. */
393 dim = (*pdim) - 1;
394 rank = GFC_DESCRIPTOR_RANK (array) - 1;
395
396 for (n = 0; n < dim; n++)
397 {
398 extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);
399
400 if (extent[n] <= 0)
401 extent[n] = 0;
402 }
403
404 for (n = dim; n < rank; n++)
405 {
406 extent[n] =
407 GFC_DESCRIPTOR_EXTENT(array,n + 1);
408
409 if (extent[n] <= 0)
410 extent[n] = 0;
411 }
412
413 if (retarray->base_addr == NULL)
414 {
415 size_t alloc_size, str;
416
417 for (n = 0; n < rank; n++)
418 {
419 if (n == 0)
420 str = 1;
421 else
422 str = GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1];
423
424 GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str);
425
426 }
427
428 retarray->offset = 0;
429 retarray->dtype = (array->dtype & ~GFC_DTYPE_RANK_MASK) | rank;
430
431 alloc_size = GFC_DESCRIPTOR_STRIDE(retarray,rank-1) * extent[rank-1];
432
433 if (alloc_size == 0)
434 {
435 /* Make sure we have a zero-sized array. */
436 GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 1);
437 return;
438 }
439 else
440 retarray->base_addr = xmallocarray (alloc_size, sizeof (GFC_REAL_16));
441 }
442 else
443 {
444 if (rank != GFC_DESCRIPTOR_RANK (retarray))
445 runtime_error ("rank of return array incorrect in"
446 " PRODUCT intrinsic: is %ld, should be %ld",
447 (long int) (GFC_DESCRIPTOR_RANK (retarray)),
448 (long int) rank);
449
450 if (unlikely (compile_options.bounds_check))
451 {
452 for (n=0; n < rank; n++)
453 {
454 index_type ret_extent;
455
456 ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,n);
457 if (extent[n] != ret_extent)
458 runtime_error ("Incorrect extent in return value of"
459 " PRODUCT intrinsic in dimension %ld:"
460 " is %ld, should be %ld", (long int) n + 1,
461 (long int) ret_extent, (long int) extent[n]);
462 }
463 }
464 }
465
466 for (n = 0; n < rank; n++)
467 {
468 count[n] = 0;
469 dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);
470 }
471
472 dest = retarray->base_addr;
473
474 while(1)
475 {
476 *dest = 1;
477 count[0]++;
478 dest += dstride[0];
479 n = 0;
480 while (count[n] == extent[n])
481 {
482 /* When we get to the end of a dimension, reset it and increment
483 the next dimension. */
484 count[n] = 0;
485 /* We could precalculate these products, but this is a less
486 frequently used path so probably not worth it. */
487 dest -= dstride[n] * extent[n];
488 n++;
489 if (n == rank)
490 return;
491 else
492 {
493 count[n]++;
494 dest += dstride[n];
495 }
496 }
497 }
498 }
499
500 #endif
501