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_COMPLEX_10) && defined (HAVE_GFC_COMPLEX_10)
39 extern void product_c10 (gfc_array_c10
* const restrict
,
40 gfc_array_c10
* const restrict
, const index_type
* const restrict
);
41 export_proto(product_c10
);
44 product_c10 (gfc_array_c10
* const restrict retarray
,
45 gfc_array_c10
* 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_COMPLEX_10
* restrict base
;
53 GFC_COMPLEX_10
* 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_COMPLEX_10
) * 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 %ld, should be %ld",
120 (long int) (GFC_DESCRIPTOR_RANK (retarray
)),
123 if (compile_options
.bounds_check
)
125 for (n
=0; n
< rank
; n
++)
127 index_type ret_extent
;
129 ret_extent
= retarray
->dim
[n
].ubound
+ 1
130 - retarray
->dim
[n
].lbound
;
131 if (extent
[n
] != ret_extent
)
132 runtime_error ("Incorrect extent in return value of"
133 " PRODUCT intrinsic in dimension %ld:"
134 " is %ld, should be %ld", (long int) n
+ 1,
135 (long int) ret_extent
, (long int) extent
[n
]);
140 for (n
= 0; n
< rank
; n
++)
143 dstride
[n
] = retarray
->dim
[n
].stride
;
149 dest
= retarray
->data
;
153 const GFC_COMPLEX_10
* restrict src
;
154 GFC_COMPLEX_10 result
;
163 for (n
= 0; n
< len
; n
++, src
+= delta
)
171 /* Advance to the next element. */
176 while (count
[n
] == extent
[n
])
178 /* When we get to the end of a dimension, reset it and increment
179 the next dimension. */
181 /* We could precalculate these products, but this is a less
182 frequently used path so probably not worth it. */
183 base
-= sstride
[n
] * extent
[n
];
184 dest
-= dstride
[n
] * extent
[n
];
188 /* Break out of the look. */
203 extern void mproduct_c10 (gfc_array_c10
* const restrict
,
204 gfc_array_c10
* const restrict
, const index_type
* const restrict
,
205 gfc_array_l1
* const restrict
);
206 export_proto(mproduct_c10
);
209 mproduct_c10 (gfc_array_c10
* const restrict retarray
,
210 gfc_array_c10
* const restrict array
,
211 const index_type
* const restrict pdim
,
212 gfc_array_l1
* const restrict mask
)
214 index_type count
[GFC_MAX_DIMENSIONS
];
215 index_type extent
[GFC_MAX_DIMENSIONS
];
216 index_type sstride
[GFC_MAX_DIMENSIONS
];
217 index_type dstride
[GFC_MAX_DIMENSIONS
];
218 index_type mstride
[GFC_MAX_DIMENSIONS
];
219 GFC_COMPLEX_10
* restrict dest
;
220 const GFC_COMPLEX_10
* restrict base
;
221 const GFC_LOGICAL_1
* restrict mbase
;
231 rank
= GFC_DESCRIPTOR_RANK (array
) - 1;
233 len
= array
->dim
[dim
].ubound
+ 1 - array
->dim
[dim
].lbound
;
239 mask_kind
= GFC_DESCRIPTOR_SIZE (mask
);
241 if (mask_kind
== 1 || mask_kind
== 2 || mask_kind
== 4 || mask_kind
== 8
242 #ifdef HAVE_GFC_LOGICAL_16
246 mbase
= GFOR_POINTER_TO_L1 (mbase
, mask_kind
);
248 runtime_error ("Funny sized logical array");
250 delta
= array
->dim
[dim
].stride
;
251 mdelta
= mask
->dim
[dim
].stride
* mask_kind
;
253 for (n
= 0; n
< dim
; n
++)
255 sstride
[n
] = array
->dim
[n
].stride
;
256 mstride
[n
] = mask
->dim
[n
].stride
* mask_kind
;
257 extent
[n
] = array
->dim
[n
].ubound
+ 1 - array
->dim
[n
].lbound
;
263 for (n
= dim
; n
< rank
; n
++)
265 sstride
[n
] = array
->dim
[n
+ 1].stride
;
266 mstride
[n
] = mask
->dim
[n
+ 1].stride
* mask_kind
;
268 array
->dim
[n
+ 1].ubound
+ 1 - array
->dim
[n
+ 1].lbound
;
274 if (retarray
->data
== NULL
)
278 for (n
= 0; n
< rank
; n
++)
280 retarray
->dim
[n
].lbound
= 0;
281 retarray
->dim
[n
].ubound
= extent
[n
]-1;
283 retarray
->dim
[n
].stride
= 1;
285 retarray
->dim
[n
].stride
= retarray
->dim
[n
-1].stride
* extent
[n
-1];
288 alloc_size
= sizeof (GFC_COMPLEX_10
) * retarray
->dim
[rank
-1].stride
291 retarray
->offset
= 0;
292 retarray
->dtype
= (array
->dtype
& ~GFC_DTYPE_RANK_MASK
) | rank
;
296 /* Make sure we have a zero-sized array. */
297 retarray
->dim
[0].lbound
= 0;
298 retarray
->dim
[0].ubound
= -1;
302 retarray
->data
= internal_malloc_size (alloc_size
);
307 if (rank
!= GFC_DESCRIPTOR_RANK (retarray
))
308 runtime_error ("rank of return array incorrect in PRODUCT intrinsic");
310 if (compile_options
.bounds_check
)
312 for (n
=0; n
< rank
; n
++)
314 index_type ret_extent
;
316 ret_extent
= retarray
->dim
[n
].ubound
+ 1
317 - retarray
->dim
[n
].lbound
;
318 if (extent
[n
] != ret_extent
)
319 runtime_error ("Incorrect extent in return value of"
320 " PRODUCT intrinsic in dimension %ld:"
321 " is %ld, should be %ld", (long int) n
+ 1,
322 (long int) ret_extent
, (long int) extent
[n
]);
324 for (n
=0; n
<= rank
; n
++)
326 index_type mask_extent
, array_extent
;
328 array_extent
= array
->dim
[n
].ubound
+ 1 - array
->dim
[n
].lbound
;
329 mask_extent
= mask
->dim
[n
].ubound
+ 1 - mask
->dim
[n
].lbound
;
330 if (array_extent
!= mask_extent
)
331 runtime_error ("Incorrect extent in MASK argument of"
332 " PRODUCT intrinsic in dimension %ld:"
333 " is %ld, should be %ld", (long int) n
+ 1,
334 (long int) mask_extent
, (long int) array_extent
);
339 for (n
= 0; n
< rank
; n
++)
342 dstride
[n
] = retarray
->dim
[n
].stride
;
347 dest
= retarray
->data
;
352 const GFC_COMPLEX_10
* restrict src
;
353 const GFC_LOGICAL_1
* restrict msrc
;
354 GFC_COMPLEX_10 result
;
364 for (n
= 0; n
< len
; n
++, src
+= delta
, msrc
+= mdelta
)
373 /* Advance to the next element. */
379 while (count
[n
] == extent
[n
])
381 /* When we get to the end of a dimension, reset it and increment
382 the next dimension. */
384 /* We could precalculate these products, but this is a less
385 frequently used path so probably not worth it. */
386 base
-= sstride
[n
] * extent
[n
];
387 mbase
-= mstride
[n
] * extent
[n
];
388 dest
-= dstride
[n
] * extent
[n
];
392 /* Break out of the look. */
408 extern void sproduct_c10 (gfc_array_c10
* const restrict
,
409 gfc_array_c10
* const restrict
, const index_type
* const restrict
,
411 export_proto(sproduct_c10
);
414 sproduct_c10 (gfc_array_c10
* const restrict retarray
,
415 gfc_array_c10
* const restrict array
,
416 const index_type
* const restrict pdim
,
417 GFC_LOGICAL_4
* mask
)
422 GFC_COMPLEX_10
*dest
;
426 product_c10 (retarray
, array
, pdim
);
429 rank
= GFC_DESCRIPTOR_RANK (array
);
431 runtime_error ("Rank of array needs to be > 0");
433 if (retarray
->data
== NULL
)
435 retarray
->dim
[0].lbound
= 0;
436 retarray
->dim
[0].ubound
= rank
-1;
437 retarray
->dim
[0].stride
= 1;
438 retarray
->dtype
= (retarray
->dtype
& ~GFC_DTYPE_RANK_MASK
) | 1;
439 retarray
->offset
= 0;
440 retarray
->data
= internal_malloc_size (sizeof (GFC_COMPLEX_10
) * rank
);
444 if (compile_options
.bounds_check
)
447 index_type ret_extent
;
449 ret_rank
= GFC_DESCRIPTOR_RANK (retarray
);
451 runtime_error ("rank of return array in PRODUCT intrinsic"
452 " should be 1, is %ld", (long int) ret_rank
);
454 ret_extent
= retarray
->dim
[0].ubound
+ 1 - retarray
->dim
[0].lbound
;
455 if (ret_extent
!= rank
)
456 runtime_error ("dimension of return array incorrect");
459 dstride
= retarray
->dim
[0].stride
;
460 dest
= retarray
->data
;
462 for (n
= 0; n
< rank
; n
++)
463 dest
[n
* dstride
] = 1 ;