1 /* Implementation of the PRODUCT intrinsic
2 Copyright 2002, 2007 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 2 of the License, or (at your option) any later version.
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
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.
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. */
31 #include "libgfortran.h"
36 #if defined (HAVE_GFC_REAL_16) && defined (HAVE_GFC_REAL_16)
39 extern void product_r16 (gfc_array_r16
* const restrict
,
40 gfc_array_r16
* const restrict
, const index_type
* const restrict
);
41 export_proto(product_r16
);
44 product_r16 (gfc_array_r16
* const restrict retarray
,
45 gfc_array_r16
* const restrict array
,
46 const index_type
* const restrict pdim
)
48 index_type count
[GFC_MAX_DIMENSIONS
];
49 index_type extent
[GFC_MAX_DIMENSIONS
];
50 index_type sstride
[GFC_MAX_DIMENSIONS
];
51 index_type dstride
[GFC_MAX_DIMENSIONS
];
52 const GFC_REAL_16
* restrict base
;
53 GFC_REAL_16
* restrict dest
;
60 /* Make dim zero based to avoid confusion. */
62 rank
= GFC_DESCRIPTOR_RANK (array
) - 1;
64 len
= array
->dim
[dim
].ubound
+ 1 - array
->dim
[dim
].lbound
;
65 delta
= array
->dim
[dim
].stride
;
67 for (n
= 0; n
< dim
; n
++)
69 sstride
[n
] = array
->dim
[n
].stride
;
70 extent
[n
] = array
->dim
[n
].ubound
+ 1 - array
->dim
[n
].lbound
;
75 for (n
= dim
; n
< rank
; n
++)
77 sstride
[n
] = array
->dim
[n
+ 1].stride
;
79 array
->dim
[n
+ 1].ubound
+ 1 - array
->dim
[n
+ 1].lbound
;
85 if (retarray
->data
== NULL
)
89 for (n
= 0; n
< rank
; n
++)
91 retarray
->dim
[n
].lbound
= 0;
92 retarray
->dim
[n
].ubound
= extent
[n
]-1;
94 retarray
->dim
[n
].stride
= 1;
96 retarray
->dim
[n
].stride
= retarray
->dim
[n
-1].stride
* extent
[n
-1];
100 retarray
->dtype
= (array
->dtype
& ~GFC_DTYPE_RANK_MASK
) | rank
;
102 alloc_size
= sizeof (GFC_REAL_16
) * retarray
->dim
[rank
-1].stride
107 /* Make sure we have a zero-sized array. */
108 retarray
->dim
[0].lbound
= 0;
109 retarray
->dim
[0].ubound
= -1;
113 retarray
->data
= internal_malloc_size (alloc_size
);
117 if (rank
!= GFC_DESCRIPTOR_RANK (retarray
))
118 runtime_error ("rank of return array incorrect in"
119 " PRODUCT intrinsic: is %d, should be %d",
120 GFC_DESCRIPTOR_RANK (retarray
), rank
);
122 if (compile_options
.bounds_check
)
124 for (n
=0; n
< rank
; n
++)
126 index_type ret_extent
;
128 ret_extent
= retarray
->dim
[n
].ubound
+ 1
129 - retarray
->dim
[n
].lbound
;
130 if (extent
[n
] != ret_extent
)
131 runtime_error ("Incorrect extent in return value of"
132 " PRODUCT intrinsic in dimension %d:"
133 " is %ld, should be %ld", n
+ 1,
134 (long int) ret_extent
, (long int) extent
[n
]);
139 for (n
= 0; n
< rank
; n
++)
142 dstride
[n
] = retarray
->dim
[n
].stride
;
148 dest
= retarray
->data
;
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 (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 %d:"
320 " is %ld, should be %ld", 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 %d:"
332 " is %ld, should be %ld", 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
)
425 product_r16 (retarray
, array
, pdim
);
428 rank
= GFC_DESCRIPTOR_RANK (array
);
430 runtime_error ("Rank of array needs to be > 0");
432 if (retarray
->data
== NULL
)
434 retarray
->dim
[0].lbound
= 0;
435 retarray
->dim
[0].ubound
= rank
-1;
436 retarray
->dim
[0].stride
= 1;
437 retarray
->dtype
= (retarray
->dtype
& ~GFC_DTYPE_RANK_MASK
) | 1;
438 retarray
->offset
= 0;
439 retarray
->data
= internal_malloc_size (sizeof (GFC_REAL_16
) * rank
);
443 if (compile_options
.bounds_check
)
446 index_type ret_extent
;
448 ret_rank
= GFC_DESCRIPTOR_RANK (retarray
);
450 runtime_error ("rank of return array in PRODUCT intrinsic"
451 " should be 1, is %d", ret_rank
);
453 ret_extent
= retarray
->dim
[0].ubound
+ 1 - retarray
->dim
[0].lbound
;
454 if (ret_extent
!= rank
)
455 runtime_error ("dimension of return array incorrect");
458 dstride
= retarray
->dim
[0].stride
;
459 dest
= retarray
->data
;
461 for (n
= 0; n
< rank
; n
++)
462 dest
[n
* dstride
] = 1 ;