1 /* Implementation of the PRODUCT intrinsic
2 Copyright 2002, 2007, 2009 Free Software Foundation, Inc.
3 Contributed by Paul Brook <paul@nowt.org>
5 This file is part of the GNU Fortran 95 runtime library (libgfortran).
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.
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.
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.
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/>. */
26 #include "libgfortran.h"
31 #if defined (HAVE_GFC_REAL_16) && defined (HAVE_GFC_REAL_16)
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
);
39 product_r16 (gfc_array_r16
* const restrict retarray
,
40 gfc_array_r16
* const restrict array
,
41 const index_type
* const restrict pdim
)
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
;
56 /* Make dim zero based to avoid confusion. */
58 rank
= GFC_DESCRIPTOR_RANK (array
) - 1;
60 len
= array
->dim
[dim
].ubound
+ 1 - array
->dim
[dim
].lbound
;
63 delta
= array
->dim
[dim
].stride
;
65 for (n
= 0; n
< dim
; n
++)
67 sstride
[n
] = array
->dim
[n
].stride
;
68 extent
[n
] = array
->dim
[n
].ubound
+ 1 - array
->dim
[n
].lbound
;
73 for (n
= dim
; n
< rank
; n
++)
75 sstride
[n
] = array
->dim
[n
+ 1].stride
;
77 array
->dim
[n
+ 1].ubound
+ 1 - array
->dim
[n
+ 1].lbound
;
83 if (retarray
->data
== NULL
)
87 for (n
= 0; n
< rank
; n
++)
89 retarray
->dim
[n
].lbound
= 0;
90 retarray
->dim
[n
].ubound
= extent
[n
]-1;
92 retarray
->dim
[n
].stride
= 1;
94 retarray
->dim
[n
].stride
= retarray
->dim
[n
-1].stride
* extent
[n
-1];
98 retarray
->dtype
= (array
->dtype
& ~GFC_DTYPE_RANK_MASK
) | rank
;
100 alloc_size
= sizeof (GFC_REAL_16
) * retarray
->dim
[rank
-1].stride
105 /* Make sure we have a zero-sized array. */
106 retarray
->dim
[0].lbound
= 0;
107 retarray
->dim
[0].ubound
= -1;
111 retarray
->data
= internal_malloc_size (alloc_size
);
115 if (rank
!= GFC_DESCRIPTOR_RANK (retarray
))
116 runtime_error ("rank of return array incorrect in"
117 " PRODUCT intrinsic: is %ld, should be %ld",
118 (long int) (GFC_DESCRIPTOR_RANK (retarray
)),
121 if (unlikely (compile_options
.bounds_check
))
123 for (n
=0; n
< rank
; n
++)
125 index_type ret_extent
;
127 ret_extent
= retarray
->dim
[n
].ubound
+ 1
128 - retarray
->dim
[n
].lbound
;
129 if (extent
[n
] != ret_extent
)
130 runtime_error ("Incorrect extent in return value of"
131 " PRODUCT intrinsic in dimension %ld:"
132 " is %ld, should be %ld", (long int) n
+ 1,
133 (long int) ret_extent
, (long int) extent
[n
]);
138 for (n
= 0; n
< rank
; n
++)
141 dstride
[n
] = retarray
->dim
[n
].stride
;
147 dest
= retarray
->data
;
150 while (continue_loop
)
152 const GFC_REAL_16
* restrict src
;
162 for (n
= 0; n
< len
; n
++, src
+= delta
)
170 /* Advance to the next element. */
175 while (count
[n
] == extent
[n
])
177 /* When we get to the end of a dimension, reset it and increment
178 the next dimension. */
180 /* We could precalculate these products, but this is a less
181 frequently used path so probably not worth it. */
182 base
-= sstride
[n
] * extent
[n
];
183 dest
-= dstride
[n
] * extent
[n
];
187 /* Break out of the look. */
202 extern void mproduct_r16 (gfc_array_r16
* const restrict
,
203 gfc_array_r16
* const restrict
, const index_type
* const restrict
,
204 gfc_array_l1
* const restrict
);
205 export_proto(mproduct_r16
);
208 mproduct_r16 (gfc_array_r16
* const restrict retarray
,
209 gfc_array_r16
* const restrict array
,
210 const index_type
* const restrict pdim
,
211 gfc_array_l1
* const restrict mask
)
213 index_type count
[GFC_MAX_DIMENSIONS
];
214 index_type extent
[GFC_MAX_DIMENSIONS
];
215 index_type sstride
[GFC_MAX_DIMENSIONS
];
216 index_type dstride
[GFC_MAX_DIMENSIONS
];
217 index_type mstride
[GFC_MAX_DIMENSIONS
];
218 GFC_REAL_16
* restrict dest
;
219 const GFC_REAL_16
* restrict base
;
220 const GFC_LOGICAL_1
* restrict mbase
;
230 rank
= GFC_DESCRIPTOR_RANK (array
) - 1;
232 len
= array
->dim
[dim
].ubound
+ 1 - array
->dim
[dim
].lbound
;
238 mask_kind
= GFC_DESCRIPTOR_SIZE (mask
);
240 if (mask_kind
== 1 || mask_kind
== 2 || mask_kind
== 4 || mask_kind
== 8
241 #ifdef HAVE_GFC_LOGICAL_16
245 mbase
= GFOR_POINTER_TO_L1 (mbase
, mask_kind
);
247 runtime_error ("Funny sized logical array");
249 delta
= array
->dim
[dim
].stride
;
250 mdelta
= mask
->dim
[dim
].stride
* mask_kind
;
252 for (n
= 0; n
< dim
; n
++)
254 sstride
[n
] = array
->dim
[n
].stride
;
255 mstride
[n
] = mask
->dim
[n
].stride
* mask_kind
;
256 extent
[n
] = array
->dim
[n
].ubound
+ 1 - array
->dim
[n
].lbound
;
262 for (n
= dim
; n
< rank
; n
++)
264 sstride
[n
] = array
->dim
[n
+ 1].stride
;
265 mstride
[n
] = mask
->dim
[n
+ 1].stride
* mask_kind
;
267 array
->dim
[n
+ 1].ubound
+ 1 - array
->dim
[n
+ 1].lbound
;
273 if (retarray
->data
== NULL
)
277 for (n
= 0; n
< rank
; n
++)
279 retarray
->dim
[n
].lbound
= 0;
280 retarray
->dim
[n
].ubound
= extent
[n
]-1;
282 retarray
->dim
[n
].stride
= 1;
284 retarray
->dim
[n
].stride
= retarray
->dim
[n
-1].stride
* extent
[n
-1];
287 alloc_size
= sizeof (GFC_REAL_16
) * retarray
->dim
[rank
-1].stride
290 retarray
->offset
= 0;
291 retarray
->dtype
= (array
->dtype
& ~GFC_DTYPE_RANK_MASK
) | rank
;
295 /* Make sure we have a zero-sized array. */
296 retarray
->dim
[0].lbound
= 0;
297 retarray
->dim
[0].ubound
= -1;
301 retarray
->data
= internal_malloc_size (alloc_size
);
306 if (rank
!= GFC_DESCRIPTOR_RANK (retarray
))
307 runtime_error ("rank of return array incorrect in PRODUCT intrinsic");
309 if (unlikely (compile_options
.bounds_check
))
311 for (n
=0; n
< rank
; n
++)
313 index_type ret_extent
;
315 ret_extent
= retarray
->dim
[n
].ubound
+ 1
316 - retarray
->dim
[n
].lbound
;
317 if (extent
[n
] != ret_extent
)
318 runtime_error ("Incorrect extent in return value of"
319 " PRODUCT intrinsic in dimension %ld:"
320 " is %ld, should be %ld", (long int) n
+ 1,
321 (long int) ret_extent
, (long int) extent
[n
]);
323 for (n
=0; n
<= rank
; n
++)
325 index_type mask_extent
, array_extent
;
327 array_extent
= array
->dim
[n
].ubound
+ 1 - array
->dim
[n
].lbound
;
328 mask_extent
= mask
->dim
[n
].ubound
+ 1 - mask
->dim
[n
].lbound
;
329 if (array_extent
!= mask_extent
)
330 runtime_error ("Incorrect extent in MASK argument of"
331 " PRODUCT intrinsic in dimension %ld:"
332 " is %ld, should be %ld", (long int) n
+ 1,
333 (long int) mask_extent
, (long int) array_extent
);
338 for (n
= 0; n
< rank
; n
++)
341 dstride
[n
] = retarray
->dim
[n
].stride
;
346 dest
= retarray
->data
;
351 const GFC_REAL_16
* restrict src
;
352 const GFC_LOGICAL_1
* restrict msrc
;
363 for (n
= 0; n
< len
; n
++, src
+= delta
, msrc
+= mdelta
)
372 /* Advance to the next element. */
378 while (count
[n
] == extent
[n
])
380 /* When we get to the end of a dimension, reset it and increment
381 the next dimension. */
383 /* We could precalculate these products, but this is a less
384 frequently used path so probably not worth it. */
385 base
-= sstride
[n
] * extent
[n
];
386 mbase
-= mstride
[n
] * extent
[n
];
387 dest
-= dstride
[n
] * extent
[n
];
391 /* Break out of the look. */
407 extern void sproduct_r16 (gfc_array_r16
* const restrict
,
408 gfc_array_r16
* const restrict
, const index_type
* const restrict
,
410 export_proto(sproduct_r16
);
413 sproduct_r16 (gfc_array_r16
* const restrict retarray
,
414 gfc_array_r16
* const restrict array
,
415 const index_type
* const restrict pdim
,
416 GFC_LOGICAL_4
* mask
)
418 index_type count
[GFC_MAX_DIMENSIONS
];
419 index_type extent
[GFC_MAX_DIMENSIONS
];
420 index_type sstride
[GFC_MAX_DIMENSIONS
];
421 index_type dstride
[GFC_MAX_DIMENSIONS
];
422 GFC_REAL_16
* restrict dest
;
430 product_r16 (retarray
, array
, pdim
);
433 /* Make dim zero based to avoid confusion. */
435 rank
= GFC_DESCRIPTOR_RANK (array
) - 1;
437 for (n
= 0; n
< dim
; n
++)
439 sstride
[n
] = array
->dim
[n
].stride
;
440 extent
[n
] = array
->dim
[n
].ubound
+ 1 - array
->dim
[n
].lbound
;
446 for (n
= dim
; n
< rank
; n
++)
448 sstride
[n
] = array
->dim
[n
+ 1].stride
;
450 array
->dim
[n
+ 1].ubound
+ 1 - array
->dim
[n
+ 1].lbound
;
456 if (retarray
->data
== NULL
)
460 for (n
= 0; n
< rank
; n
++)
462 retarray
->dim
[n
].lbound
= 0;
463 retarray
->dim
[n
].ubound
= extent
[n
]-1;
465 retarray
->dim
[n
].stride
= 1;
467 retarray
->dim
[n
].stride
= retarray
->dim
[n
-1].stride
* extent
[n
-1];
470 retarray
->offset
= 0;
471 retarray
->dtype
= (array
->dtype
& ~GFC_DTYPE_RANK_MASK
) | rank
;
473 alloc_size
= sizeof (GFC_REAL_16
) * retarray
->dim
[rank
-1].stride
478 /* Make sure we have a zero-sized array. */
479 retarray
->dim
[0].lbound
= 0;
480 retarray
->dim
[0].ubound
= -1;
484 retarray
->data
= internal_malloc_size (alloc_size
);
488 if (rank
!= GFC_DESCRIPTOR_RANK (retarray
))
489 runtime_error ("rank of return array incorrect in"
490 " PRODUCT intrinsic: is %ld, should be %ld",
491 (long int) (GFC_DESCRIPTOR_RANK (retarray
)),
494 if (unlikely (compile_options
.bounds_check
))
496 for (n
=0; n
< rank
; n
++)
498 index_type ret_extent
;
500 ret_extent
= retarray
->dim
[n
].ubound
+ 1
501 - retarray
->dim
[n
].lbound
;
502 if (extent
[n
] != ret_extent
)
503 runtime_error ("Incorrect extent in return value of"
504 " PRODUCT intrinsic in dimension %ld:"
505 " is %ld, should be %ld", (long int) n
+ 1,
506 (long int) ret_extent
, (long int) extent
[n
]);
511 for (n
= 0; n
< rank
; n
++)
514 dstride
[n
] = retarray
->dim
[n
].stride
;
517 dest
= retarray
->data
;
525 while (count
[n
] == extent
[n
])
527 /* When we get to the end of a dimension, reset it and increment
528 the next dimension. */
530 /* We could precalculate these products, but this is a less
531 frequently used path so probably not worth it. */
532 dest
-= dstride
[n
] * extent
[n
];