1/* Implementation of the PRODUCT intrinsic 2 Copyright (C) 2002-2022 Free Software Foundation, Inc. 3 Contributed by Paul Brook <paul@nowt.org> 4 5This file is part of the GNU Fortran runtime library (libgfortran). 6 7Libgfortran is free software; you can redistribute it and/or 8modify it under the terms of the GNU General Public 9License as published by the Free Software Foundation; either 10version 3 of the License, or (at your option) any later version. 11 12Libgfortran is distributed in the hope that it will be useful, 13but WITHOUT ANY WARRANTY; without even the implied warranty of 14MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15GNU General Public License for more details. 16 17Under Section 7 of GPL version 3, you are granted additional 18permissions described in the GCC Runtime Library Exception, version 193.1, as published by the Free Software Foundation. 20 21You should have received a copy of the GNU General Public License and 22a copy of the GCC Runtime Library Exception along with this program; 23see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 24<http://www.gnu.org/licenses/>. */ 25 26#include "libgfortran.h" 27 28 29#if defined (HAVE_GFC_COMPLEX_4) && defined (HAVE_GFC_COMPLEX_4) 30 31 32extern void product_c4 (gfc_array_c4 * const restrict, 33 gfc_array_c4 * const restrict, const index_type * const restrict); 34export_proto(product_c4); 35 36void 37product_c4 (gfc_array_c4 * const restrict retarray, 38 gfc_array_c4 * const restrict array, 39 const index_type * const restrict pdim) 40{ 41 index_type count[GFC_MAX_DIMENSIONS]; 42 index_type extent[GFC_MAX_DIMENSIONS]; 43 index_type sstride[GFC_MAX_DIMENSIONS]; 44 index_type dstride[GFC_MAX_DIMENSIONS]; 45 const GFC_COMPLEX_4 * restrict base; 46 GFC_COMPLEX_4 * restrict dest; 47 index_type rank; 48 index_type n; 49 index_type len; 50 index_type delta; 51 index_type dim; 52 int continue_loop; 53 54 /* Make dim zero based to avoid confusion. */ 55 rank = GFC_DESCRIPTOR_RANK (array) - 1; 56 dim = (*pdim) - 1; 57 58 if (unlikely (dim < 0 || dim > rank)) 59 { 60 runtime_error ("Dim argument incorrect in PRODUCT intrinsic: " 61 "is %ld, should be between 1 and %ld", 62 (long int) dim + 1, (long int) rank + 1); 63 } 64 65 len = GFC_DESCRIPTOR_EXTENT(array,dim); 66 if (len < 0) 67 len = 0; 68 delta = GFC_DESCRIPTOR_STRIDE(array,dim); 69 70 for (n = 0; n < dim; n++) 71 { 72 sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n); 73 extent[n] = GFC_DESCRIPTOR_EXTENT(array,n); 74 75 if (extent[n] < 0) 76 extent[n] = 0; 77 } 78 for (n = dim; n < rank; n++) 79 { 80 sstride[n] = GFC_DESCRIPTOR_STRIDE(array, n + 1); 81 extent[n] = GFC_DESCRIPTOR_EXTENT(array, n + 1); 82 83 if (extent[n] < 0) 84 extent[n] = 0; 85 } 86 87 if (retarray->base_addr == NULL) 88 { 89 size_t alloc_size, str; 90 91 for (n = 0; n < rank; n++) 92 { 93 if (n == 0) 94 str = 1; 95 else 96 str = GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1]; 97 98 GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str); 99 100 } 101 102 retarray->offset = 0; 103 retarray->dtype.rank = rank; 104 105 alloc_size = GFC_DESCRIPTOR_STRIDE(retarray,rank-1) * extent[rank-1]; 106 107 retarray->base_addr = xmallocarray (alloc_size, sizeof (GFC_COMPLEX_4)); 108 if (alloc_size == 0) 109 { 110 /* Make sure we have a zero-sized array. */ 111 GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 1); 112 return; 113 114 } 115 } 116 else 117 { 118 if (rank != GFC_DESCRIPTOR_RANK (retarray)) 119 runtime_error ("rank of return array incorrect in" 120 " PRODUCT intrinsic: is %ld, should be %ld", 121 (long int) (GFC_DESCRIPTOR_RANK (retarray)), 122 (long int) rank); 123 124 if (unlikely (compile_options.bounds_check)) 125 bounds_ifunction_return ((array_t *) retarray, extent, 126 "return value", "PRODUCT"); 127 } 128 129 for (n = 0; n < rank; n++) 130 { 131 count[n] = 0; 132 dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n); 133 if (extent[n] <= 0) 134 return; 135 } 136 137 base = array->base_addr; 138 dest = retarray->base_addr; 139 140 continue_loop = 1; 141 while (continue_loop) 142 { 143 const GFC_COMPLEX_4 * restrict src; 144 GFC_COMPLEX_4 result; 145 src = base; 146 { 147 148 result = 1; 149 if (len <= 0) 150 *dest = 1; 151 else 152 { 153#if ! defined HAVE_BACK_ARG 154 for (n = 0; n < len; n++, src += delta) 155 { 156#endif 157 158 result *= *src; 159 } 160 161 *dest = result; 162 } 163 } 164 /* Advance to the next element. */ 165 count[0]++; 166 base += sstride[0]; 167 dest += dstride[0]; 168 n = 0; 169 while (count[n] == extent[n]) 170 { 171 /* When we get to the end of a dimension, reset it and increment 172 the next dimension. */ 173 count[n] = 0; 174 /* We could precalculate these products, but this is a less 175 frequently used path so probably not worth it. */ 176 base -= sstride[n] * extent[n]; 177 dest -= dstride[n] * extent[n]; 178 n++; 179 if (n >= rank) 180 { 181 /* Break out of the loop. */ 182 continue_loop = 0; 183 break; 184 } 185 else 186 { 187 count[n]++; 188 base += sstride[n]; 189 dest += dstride[n]; 190 } 191 } 192 } 193} 194 195 196extern void mproduct_c4 (gfc_array_c4 * const restrict, 197 gfc_array_c4 * const restrict, const index_type * const restrict, 198 gfc_array_l1 * const restrict); 199export_proto(mproduct_c4); 200 201void 202mproduct_c4 (gfc_array_c4 * const restrict retarray, 203 gfc_array_c4 * const restrict array, 204 const index_type * const restrict pdim, 205 gfc_array_l1 * const restrict mask) 206{ 207 index_type count[GFC_MAX_DIMENSIONS]; 208 index_type extent[GFC_MAX_DIMENSIONS]; 209 index_type sstride[GFC_MAX_DIMENSIONS]; 210 index_type dstride[GFC_MAX_DIMENSIONS]; 211 index_type mstride[GFC_MAX_DIMENSIONS]; 212 GFC_COMPLEX_4 * restrict dest; 213 const GFC_COMPLEX_4 * restrict base; 214 const GFC_LOGICAL_1 * restrict mbase; 215 index_type rank; 216 index_type dim; 217 index_type n; 218 index_type len; 219 index_type delta; 220 index_type mdelta; 221 int mask_kind; 222 223 if (mask == NULL) 224 { 225#ifdef HAVE_BACK_ARG 226 product_c4 (retarray, array, pdim, back); 227#else 228 product_c4 (retarray, array, pdim); 229#endif 230 return; 231 } 232 233 dim = (*pdim) - 1; 234 rank = GFC_DESCRIPTOR_RANK (array) - 1; 235 236 237 if (unlikely (dim < 0 || dim > rank)) 238 { 239 runtime_error ("Dim argument incorrect in PRODUCT intrinsic: " 240 "is %ld, should be between 1 and %ld", 241 (long int) dim + 1, (long int) rank + 1); 242 } 243 244 len = GFC_DESCRIPTOR_EXTENT(array,dim); 245 if (len <= 0) 246 return; 247 248 mbase = mask->base_addr; 249 250 mask_kind = GFC_DESCRIPTOR_SIZE (mask); 251 252 if (mask_kind == 1 || mask_kind == 2 || mask_kind == 4 || mask_kind == 8 253#ifdef HAVE_GFC_LOGICAL_16 254 || mask_kind == 16 255#endif 256 ) 257 mbase = GFOR_POINTER_TO_L1 (mbase, mask_kind); 258 else 259 runtime_error ("Funny sized logical array"); 260 261 delta = GFC_DESCRIPTOR_STRIDE(array,dim); 262 mdelta = GFC_DESCRIPTOR_STRIDE_BYTES(mask,dim); 263 264 for (n = 0; n < dim; n++) 265 { 266 sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n); 267 mstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(mask,n); 268 extent[n] = GFC_DESCRIPTOR_EXTENT(array,n); 269 270 if (extent[n] < 0) 271 extent[n] = 0; 272 273 } 274 for (n = dim; n < rank; n++) 275 { 276 sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n + 1); 277 mstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(mask, n + 1); 278 extent[n] = GFC_DESCRIPTOR_EXTENT(array, n + 1); 279 280 if (extent[n] < 0) 281 extent[n] = 0; 282 } 283 284 if (retarray->base_addr == NULL) 285 { 286 size_t alloc_size, str; 287 288 for (n = 0; n < rank; n++) 289 { 290 if (n == 0) 291 str = 1; 292 else 293 str= GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1]; 294 295 GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str); 296 297 } 298 299 alloc_size = GFC_DESCRIPTOR_STRIDE(retarray,rank-1) * extent[rank-1]; 300 301 retarray->offset = 0; 302 retarray->dtype.rank = rank; 303 304 if (alloc_size == 0) 305 { 306 /* Make sure we have a zero-sized array. */ 307 GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 1); 308 return; 309 } 310 else 311 retarray->base_addr = xmallocarray (alloc_size, sizeof (GFC_COMPLEX_4)); 312 313 } 314 else 315 { 316 if (rank != GFC_DESCRIPTOR_RANK (retarray)) 317 runtime_error ("rank of return array incorrect in PRODUCT intrinsic"); 318 319 if (unlikely (compile_options.bounds_check)) 320 { 321 bounds_ifunction_return ((array_t *) retarray, extent, 322 "return value", "PRODUCT"); 323 bounds_equal_extents ((array_t *) mask, (array_t *) array, 324 "MASK argument", "PRODUCT"); 325 } 326 } 327 328 for (n = 0; n < rank; n++) 329 { 330 count[n] = 0; 331 dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n); 332 if (extent[n] <= 0) 333 return; 334 } 335 336 dest = retarray->base_addr; 337 base = array->base_addr; 338 339 while (base) 340 { 341 const GFC_COMPLEX_4 * restrict src; 342 const GFC_LOGICAL_1 * restrict msrc; 343 GFC_COMPLEX_4 result; 344 src = base; 345 msrc = mbase; 346 { 347 348 result = 1; 349 for (n = 0; n < len; n++, src += delta, msrc += mdelta) 350 { 351 352 if (*msrc) 353 result *= *src; 354 } 355 *dest = result; 356 } 357 /* Advance to the next element. */ 358 count[0]++; 359 base += sstride[0]; 360 mbase += mstride[0]; 361 dest += dstride[0]; 362 n = 0; 363 while (count[n] == extent[n]) 364 { 365 /* When we get to the end of a dimension, reset it and increment 366 the next dimension. */ 367 count[n] = 0; 368 /* We could precalculate these products, but this is a less 369 frequently used path so probably not worth it. */ 370 base -= sstride[n] * extent[n]; 371 mbase -= mstride[n] * extent[n]; 372 dest -= dstride[n] * extent[n]; 373 n++; 374 if (n >= rank) 375 { 376 /* Break out of the loop. */ 377 base = NULL; 378 break; 379 } 380 else 381 { 382 count[n]++; 383 base += sstride[n]; 384 mbase += mstride[n]; 385 dest += dstride[n]; 386 } 387 } 388 } 389} 390 391 392extern void sproduct_c4 (gfc_array_c4 * const restrict, 393 gfc_array_c4 * const restrict, const index_type * const restrict, 394 GFC_LOGICAL_4 *); 395export_proto(sproduct_c4); 396 397void 398sproduct_c4 (gfc_array_c4 * const restrict retarray, 399 gfc_array_c4 * const restrict array, 400 const index_type * const restrict pdim, 401 GFC_LOGICAL_4 * mask) 402{ 403 index_type count[GFC_MAX_DIMENSIONS]; 404 index_type extent[GFC_MAX_DIMENSIONS]; 405 index_type dstride[GFC_MAX_DIMENSIONS]; 406 GFC_COMPLEX_4 * restrict dest; 407 index_type rank; 408 index_type n; 409 index_type dim; 410 411 412 if (mask == NULL || *mask) 413 { 414#ifdef HAVE_BACK_ARG 415 product_c4 (retarray, array, pdim, back); 416#else 417 product_c4 (retarray, array, pdim); 418#endif 419 return; 420 } 421 /* Make dim zero based to avoid confusion. */ 422 dim = (*pdim) - 1; 423 rank = GFC_DESCRIPTOR_RANK (array) - 1; 424 425 if (unlikely (dim < 0 || dim > rank)) 426 { 427 runtime_error ("Dim argument incorrect in PRODUCT intrinsic: " 428 "is %ld, should be between 1 and %ld", 429 (long int) dim + 1, (long int) rank + 1); 430 } 431 432 for (n = 0; n < dim; n++) 433 { 434 extent[n] = GFC_DESCRIPTOR_EXTENT(array,n); 435 436 if (extent[n] <= 0) 437 extent[n] = 0; 438 } 439 440 for (n = dim; n < rank; n++) 441 { 442 extent[n] = 443 GFC_DESCRIPTOR_EXTENT(array,n + 1); 444 445 if (extent[n] <= 0) 446 extent[n] = 0; 447 } 448 449 if (retarray->base_addr == NULL) 450 { 451 size_t alloc_size, str; 452 453 for (n = 0; n < rank; n++) 454 { 455 if (n == 0) 456 str = 1; 457 else 458 str = GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1]; 459 460 GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str); 461 462 } 463 464 retarray->offset = 0; 465 retarray->dtype.rank = rank; 466 467 alloc_size = GFC_DESCRIPTOR_STRIDE(retarray,rank-1) * extent[rank-1]; 468 469 if (alloc_size == 0) 470 { 471 /* Make sure we have a zero-sized array. */ 472 GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 1); 473 return; 474 } 475 else 476 retarray->base_addr = xmallocarray (alloc_size, sizeof (GFC_COMPLEX_4)); 477 } 478 else 479 { 480 if (rank != GFC_DESCRIPTOR_RANK (retarray)) 481 runtime_error ("rank of return array incorrect in" 482 " PRODUCT intrinsic: is %ld, should be %ld", 483 (long int) (GFC_DESCRIPTOR_RANK (retarray)), 484 (long int) rank); 485 486 if (unlikely (compile_options.bounds_check)) 487 { 488 for (n=0; n < rank; n++) 489 { 490 index_type ret_extent; 491 492 ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,n); 493 if (extent[n] != ret_extent) 494 runtime_error ("Incorrect extent in return value of" 495 " PRODUCT intrinsic in dimension %ld:" 496 " is %ld, should be %ld", (long int) n + 1, 497 (long int) ret_extent, (long int) extent[n]); 498 } 499 } 500 } 501 502 for (n = 0; n < rank; n++) 503 { 504 count[n] = 0; 505 dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n); 506 } 507 508 dest = retarray->base_addr; 509 510 while(1) 511 { 512 *dest = 1; 513 count[0]++; 514 dest += dstride[0]; 515 n = 0; 516 while (count[n] == extent[n]) 517 { 518 /* When we get to the end of a dimension, reset it and increment 519 the next dimension. */ 520 count[n] = 0; 521 /* We could precalculate these products, but this is a less 522 frequently used path so probably not worth it. */ 523 dest -= dstride[n] * extent[n]; 524 n++; 525 if (n >= rank) 526 return; 527 else 528 { 529 count[n]++; 530 dest += dstride[n]; 531 } 532 } 533 } 534} 535 536#endif 537