Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief sums arrays of real/complex numbers with *much* reduced round-off as compared to
10 : !> a naive implementation (or the one found in most compiler's SUM intrinsic)
11 : !> using an implementation of Kahan's algorithm for summing real numbers
12 : !> that can be used instead of the standard Fortran SUM(array[,mask]).
13 : !>
14 : !> see also http://en.wikipedia.org/wiki/Kahan_summation_algorithm
15 : !> \note
16 : !> if the compiler optimises away the 'tricky' bit, no accuracy is gained,
17 : !> if the compiler uses extended precision inconsistently even worse results might be obtained.
18 : !> This has not been observed.
19 : !> This algorithm is not fast, and thus not recommended for cases where round-off is not a
20 : !> concern but performance is.
21 : !>
22 : !> the standard intrinsic sum can be 'replaced' using the following use statement
23 : !>
24 : !> USE kahan_sum, ONLY: sum => kahan_sum
25 : !> \par History
26 : !> 03.2006 [Joost VandeVondele]
27 : !> \author Joost VandeVondele
28 : ! **************************************************************************************************
29 : MODULE kahan_sum
30 :
31 : IMPLICIT NONE
32 : PRIVATE
33 : PUBLIC :: accurate_dot_product, accurate_sum, accurate_dot_product_2
34 : INTEGER, PARAMETER :: sp = KIND(0.0), dp = KIND(0.0D0)
35 : REAL(KIND=sp), PARAMETER :: szero = 0.0_sp
36 : REAL(KIND=dp), PARAMETER :: dzero = 0.0_dp
37 : COMPLEX(KIND=sp), PARAMETER :: czero = (0.0_sp, 0.0_sp)
38 : COMPLEX(KIND=dp), PARAMETER :: zzero = (0.0_dp, 0.0_dp)
39 : INTEGER, PARAMETER :: dblksize = 8
40 :
41 : INTERFACE accurate_sum
42 : MODULE PROCEDURE &
43 : kahan_sum_s1, kahan_sum_d1, kahan_sum_c1, kahan_sum_z1, &
44 : kahan_sum_s2, kahan_sum_d2, kahan_sum_c2, kahan_sum_z2, &
45 : kahan_sum_s3, kahan_sum_d3, kahan_sum_c3, kahan_sum_z3, &
46 : kahan_sum_s4, kahan_sum_d4, kahan_sum_c4, kahan_sum_z4, &
47 : kahan_sum_s5, kahan_sum_d5, kahan_sum_c5, kahan_sum_z5, &
48 : kahan_sum_s6, kahan_sum_d6, kahan_sum_c6, kahan_sum_z6, &
49 : kahan_sum_s7, kahan_sum_d7, kahan_sum_c7, kahan_sum_z7
50 : END INTERFACE accurate_sum
51 :
52 : INTERFACE accurate_dot_product
53 : MODULE PROCEDURE &
54 : kahan_dot_product_d1, &
55 : kahan_dot_product_s2, kahan_dot_product_d2, kahan_dot_product_z2, &
56 : kahan_dot_product_d3, kahan_dot_product_masked_d3
57 : END INTERFACE accurate_dot_product
58 :
59 : INTERFACE accurate_dot_product_2
60 : MODULE PROCEDURE kahan_blocked_dot_product_d1
61 : END INTERFACE
62 :
63 : CONTAINS
64 : ! **************************************************************************************************
65 : !> \brief ...
66 : !> \param array ...
67 : !> \param mask ...
68 : !> \return ...
69 : ! **************************************************************************************************
70 0 : PURE FUNCTION kahan_sum_s1(array, mask) RESULT(ks)
71 : REAL(KIND=sp), DIMENSION(:), INTENT(IN) :: array
72 : LOGICAL, DIMENSION(:), INTENT(IN), OPTIONAL :: mask
73 : REAL(KIND=sp) :: ks
74 :
75 : INTEGER :: i1
76 : REAL(KIND=sp) :: c, t, y
77 :
78 0 : ks = szero; t = szero; y = szero; c = szero
79 :
80 0 : IF (PRESENT(mask)) THEN
81 0 : DO i1 = 1, SIZE(array, 1)
82 0 : IF (mask(i1)) THEN
83 0 : y = array(i1) - c
84 0 : t = ks + y
85 0 : c = (t - ks) - y
86 0 : ks = t
87 : END IF
88 : END DO
89 : ELSE
90 0 : DO i1 = 1, SIZE(array, 1)
91 0 : y = array(i1) - c
92 0 : t = ks + y
93 0 : c = (t - ks) - y
94 0 : ks = t
95 : END DO
96 : END IF
97 0 : END FUNCTION kahan_sum_s1
98 :
99 : ! **************************************************************************************************
100 : !> \brief ...
101 : !> \param array ...
102 : !> \param mask ...
103 : !> \return ...
104 : ! **************************************************************************************************
105 1497003 : PURE FUNCTION kahan_sum_d1(array, mask) RESULT(ks)
106 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: array
107 : LOGICAL, DIMENSION(:), INTENT(IN), OPTIONAL :: mask
108 : REAL(KIND=dp) :: ks
109 :
110 : INTEGER :: i1
111 : REAL(KIND=dp) :: c, t, y
112 :
113 1497003 : ks = dzero; t = dzero; y = dzero; c = dzero
114 :
115 1497003 : IF (PRESENT(mask)) THEN
116 0 : DO i1 = 1, SIZE(array, 1)
117 0 : IF (mask(i1)) THEN
118 0 : y = array(i1) - c
119 0 : t = ks + y
120 0 : c = (t - ks) - y
121 0 : ks = t
122 : END IF
123 : END DO
124 : ELSE
125 9593683945 : DO i1 = 1, SIZE(array, 1)
126 9592186942 : y = array(i1) - c
127 9592186942 : t = ks + y
128 9592186942 : c = (t - ks) - y
129 9593683945 : ks = t
130 : END DO
131 : END IF
132 1497003 : END FUNCTION kahan_sum_d1
133 :
134 : ! **************************************************************************************************
135 : !> \brief computes the accurate sum of an array that is the element wise product of two input arrays.
136 : !> \param array1 array of real numbers
137 : !> \param array2 another array of real numbers
138 : !> \return dot product
139 : ! **************************************************************************************************
140 28344 : PURE FUNCTION kahan_dot_product_d1(array1, array2) RESULT(ks)
141 : REAL(KIND=dp), DIMENSION(:), INTENT(in) :: array1, array2
142 : REAL(KIND=dp) :: ks
143 :
144 : INTEGER :: i, n
145 : REAL(KIND=dp), DIMENSION(dblksize) :: c, ks_local, t, y
146 :
147 28344 : t = dzero; y = dzero; c = dzero; ks_local = dzero
148 :
149 28344 : n = SIZE(array1)
150 113374 : DO i = 1, MOD(n, dblksize)
151 85030 : y(1) = array1(i)*array2(i) - c(1)
152 85030 : t(1) = ks_local(1) + y(1)
153 85030 : c(1) = (t(1) - ks_local(1)) - y(1)
154 113374 : ks_local(1) = t(1)
155 : END DO
156 28344 : DO i = MOD(n, dblksize) + 1, n, dblksize
157 4895442 : y = array1(i:i + (dblksize - 1))*array2(i:i + (dblksize - 1)) - c
158 4895442 : t = ks_local + y
159 4895442 : c = (t - ks_local) - y
160 546872 : ks_local = t
161 : END DO
162 226752 : DO i = 2, dblksize
163 198408 : y(1) = ks_local(i) - (c(1) + c(i))
164 198408 : t(1) = ks_local(1) + y(1)
165 198408 : c(1) = (t(1) - ks_local(1)) - y(1)
166 226752 : ks_local(1) = t(1)
167 : END DO
168 28344 : ks = ks_local(1)
169 28344 : END FUNCTION kahan_dot_product_d1
170 :
171 : ! **************************************************************************************************
172 : !> \brief ...
173 : !> \param array ...
174 : !> \param mask ...
175 : !> \return ...
176 : ! **************************************************************************************************
177 0 : PURE FUNCTION kahan_sum_c1(array, mask) RESULT(ks)
178 : COMPLEX(KIND=sp), DIMENSION(:), INTENT(IN) :: array
179 : LOGICAL, DIMENSION(:), INTENT(IN), OPTIONAL :: mask
180 : COMPLEX(KIND=sp) :: ks
181 :
182 : COMPLEX(KIND=sp) :: c, t, y
183 : INTEGER :: i1
184 :
185 0 : ks = czero; t = czero; y = czero; c = czero
186 :
187 0 : IF (PRESENT(mask)) THEN
188 0 : DO i1 = 1, SIZE(array, 1)
189 0 : IF (mask(i1)) THEN
190 0 : y = array(i1) - c
191 0 : t = ks + y
192 0 : c = (t - ks) - y
193 0 : ks = t
194 : END IF
195 : END DO
196 : ELSE
197 0 : DO i1 = 1, SIZE(array, 1)
198 0 : y = array(i1) - c
199 0 : t = ks + y
200 0 : c = (t - ks) - y
201 0 : ks = t
202 : END DO
203 : END IF
204 0 : END FUNCTION kahan_sum_c1
205 :
206 : ! **************************************************************************************************
207 : !> \brief ...
208 : !> \param array ...
209 : !> \param mask ...
210 : !> \return ...
211 : ! **************************************************************************************************
212 24064 : PURE FUNCTION kahan_sum_z1(array, mask) RESULT(ks)
213 : COMPLEX(KIND=dp), DIMENSION(:), INTENT(IN) :: array
214 : LOGICAL, DIMENSION(:), INTENT(IN), OPTIONAL :: mask
215 : COMPLEX(KIND=dp) :: ks
216 :
217 : COMPLEX(KIND=dp) :: c, t, y
218 : INTEGER :: i1
219 :
220 24064 : ks = zzero; t = zzero; y = zzero; c = zzero
221 :
222 24064 : IF (PRESENT(mask)) THEN
223 0 : DO i1 = 1, SIZE(array, 1)
224 0 : IF (mask(i1)) THEN
225 0 : y = array(i1) - c
226 0 : t = ks + y
227 0 : c = (t - ks) - y
228 0 : ks = t
229 : END IF
230 : END DO
231 : ELSE
232 621568 : DO i1 = 1, SIZE(array, 1)
233 597504 : y = array(i1) - c
234 597504 : t = ks + y
235 597504 : c = (t - ks) - y
236 621568 : ks = t
237 : END DO
238 : END IF
239 24064 : END FUNCTION kahan_sum_z1
240 :
241 : ! **************************************************************************************************
242 : !> \brief ...
243 : !> \param array ...
244 : !> \param mask ...
245 : !> \return ...
246 : ! **************************************************************************************************
247 0 : PURE FUNCTION kahan_sum_s2(array, mask) RESULT(ks)
248 : REAL(KIND=sp), DIMENSION(:, :), INTENT(IN) :: array
249 : LOGICAL, DIMENSION(:, :), INTENT(IN), OPTIONAL :: mask
250 : REAL(KIND=sp) :: ks
251 :
252 : INTEGER :: i1, i2
253 : REAL(KIND=sp) :: c, t, y
254 :
255 0 : ks = szero; t = szero; y = szero; c = szero
256 :
257 0 : IF (PRESENT(mask)) THEN
258 0 : DO i2 = 1, SIZE(array, 2)
259 0 : DO i1 = 1, SIZE(array, 1)
260 0 : IF (mask(i1, i2)) THEN
261 0 : y = array(i1, i2) - c
262 0 : t = ks + y
263 0 : c = (t - ks) - y
264 0 : ks = t
265 : END IF
266 : END DO
267 : END DO
268 : ELSE
269 0 : DO i2 = 1, SIZE(array, 2)
270 0 : DO i1 = 1, SIZE(array, 1)
271 0 : y = array(i1, i2) - c
272 0 : t = ks + y
273 0 : c = (t - ks) - y
274 0 : ks = t
275 : END DO
276 : END DO
277 : END IF
278 0 : END FUNCTION kahan_sum_s2
279 :
280 : ! **************************************************************************************************
281 : !> \brief computes the accurate sum of an array that is the element wise product of two input arrays.
282 : !> \param array1 2-D array of real numbers
283 : !> \param array2 another 2-D array of real numbers
284 : !> \return dot product
285 : ! **************************************************************************************************
286 0 : PURE FUNCTION kahan_dot_product_s2(array1, array2) RESULT(ks)
287 : REAL(KIND=sp), DIMENSION(:, :), INTENT(in) :: array1, array2
288 : REAL(KIND=dp) :: ks
289 :
290 : INTEGER :: i1, i2, n1, n2
291 : REAL(KIND=dp) :: c, t, y
292 :
293 0 : ks = dzero; t = dzero; y = dzero; c = dzero
294 :
295 0 : n1 = SIZE(array1, 1)
296 0 : n2 = SIZE(array1, 2)
297 0 : DO i2 = 1, n2
298 0 : DO i1 = 1, n1
299 0 : y = REAL(array1(i1, i2), dp)*REAL(array2(i1, i2), dp) - c
300 0 : t = ks + y
301 0 : c = (t - ks) - y
302 0 : ks = t
303 : END DO
304 : END DO
305 0 : END FUNCTION kahan_dot_product_s2
306 :
307 : ! **************************************************************************************************
308 : !> \brief ...
309 : !> \param array ...
310 : !> \param mask ...
311 : !> \return ...
312 : ! **************************************************************************************************
313 20704 : PURE FUNCTION kahan_sum_d2(array, mask) RESULT(ks)
314 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: array
315 : LOGICAL, DIMENSION(:, :), INTENT(IN), OPTIONAL :: mask
316 : REAL(KIND=dp) :: ks
317 :
318 : INTEGER :: i1, i2
319 : REAL(KIND=dp) :: c, t, y
320 :
321 20704 : ks = dzero; t = dzero; y = dzero; c = dzero
322 :
323 20704 : IF (PRESENT(mask)) THEN
324 0 : DO i2 = 1, SIZE(array, 2)
325 0 : DO i1 = 1, SIZE(array, 1)
326 0 : IF (mask(i1, i2)) THEN
327 0 : y = array(i1, i2) - c
328 0 : t = ks + y
329 0 : c = (t - ks) - y
330 0 : ks = t
331 : END IF
332 : END DO
333 : END DO
334 : ELSE
335 1300512 : DO i2 = 1, SIZE(array, 2)
336 44549664 : DO i1 = 1, SIZE(array, 1)
337 43249152 : y = array(i1, i2) - c
338 43249152 : t = ks + y
339 43249152 : c = (t - ks) - y
340 44528960 : ks = t
341 : END DO
342 : END DO
343 : END IF
344 20704 : END FUNCTION kahan_sum_d2
345 :
346 : ! **************************************************************************************************
347 : !> \brief computes the accurate sum of an array that is the element wise product of two input arrays.
348 : !> \param array1 2-D array of real numbers
349 : !> \param array2 another 2-D array of real numbers
350 : !> \return dot product
351 : ! **************************************************************************************************
352 558804 : PURE FUNCTION kahan_dot_product_d2(array1, array2) RESULT(ks)
353 : REAL(KIND=dp), DIMENSION(:, :), INTENT(in) :: array1, array2
354 : REAL(KIND=dp) :: ks
355 :
356 : INTEGER :: i1, i2, n1, n2
357 : REAL(KIND=dp) :: c, t, y
358 :
359 558804 : ks = dzero; t = dzero; y = dzero; c = dzero
360 :
361 558804 : n1 = SIZE(array1, 1)
362 558804 : n2 = SIZE(array1, 2)
363 7480338 : DO i2 = 1, n2
364 188199020 : DO i1 = 1, n1
365 180718682 : y = array1(i1, i2)*array2(i1, i2) - c
366 180718682 : t = ks + y
367 180718682 : c = (t - ks) - y
368 187640216 : ks = t
369 : END DO
370 : END DO
371 558804 : END FUNCTION kahan_dot_product_d2
372 :
373 : ! **************************************************************************************************
374 : !> \brief ...
375 : !> \param array ...
376 : !> \param mask ...
377 : !> \return ...
378 : ! **************************************************************************************************
379 0 : PURE FUNCTION kahan_sum_c2(array, mask) RESULT(ks)
380 : COMPLEX(KIND=sp), DIMENSION(:, :), INTENT(IN) :: array
381 : LOGICAL, DIMENSION(:, :), INTENT(IN), OPTIONAL :: mask
382 : COMPLEX(KIND=sp) :: ks
383 :
384 : COMPLEX(KIND=sp) :: c, t, y
385 : INTEGER :: i1, i2
386 :
387 0 : ks = czero; t = czero; y = czero; c = czero
388 :
389 0 : IF (PRESENT(mask)) THEN
390 0 : DO i2 = 1, SIZE(array, 2)
391 0 : DO i1 = 1, SIZE(array, 1)
392 0 : IF (mask(i1, i2)) THEN
393 0 : y = array(i1, i2) - c
394 0 : t = ks + y
395 0 : c = (t - ks) - y
396 0 : ks = t
397 : END IF
398 : END DO
399 : END DO
400 : ELSE
401 0 : DO i2 = 1, SIZE(array, 2)
402 0 : DO i1 = 1, SIZE(array, 1)
403 0 : y = array(i1, i2) - c
404 0 : t = ks + y
405 0 : c = (t - ks) - y
406 0 : ks = t
407 : END DO
408 : END DO
409 : END IF
410 0 : END FUNCTION kahan_sum_c2
411 :
412 : ! **************************************************************************************************
413 : !> \brief ...
414 : !> \param array ...
415 : !> \param mask ...
416 : !> \return ...
417 : ! **************************************************************************************************
418 0 : PURE FUNCTION kahan_sum_z2(array, mask) RESULT(ks)
419 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: array
420 : LOGICAL, DIMENSION(:, :), INTENT(IN), OPTIONAL :: mask
421 : COMPLEX(KIND=dp) :: ks
422 :
423 : COMPLEX(KIND=dp) :: c, t, y
424 : INTEGER :: i1, i2
425 :
426 0 : ks = zzero; t = zzero; y = zzero; c = zzero
427 :
428 0 : IF (PRESENT(mask)) THEN
429 0 : DO i2 = 1, SIZE(array, 2)
430 0 : DO i1 = 1, SIZE(array, 1)
431 0 : IF (mask(i1, i2)) THEN
432 0 : y = array(i1, i2) - c
433 0 : t = ks + y
434 0 : c = (t - ks) - y
435 0 : ks = t
436 : END IF
437 : END DO
438 : END DO
439 : ELSE
440 0 : DO i2 = 1, SIZE(array, 2)
441 0 : DO i1 = 1, SIZE(array, 1)
442 0 : y = array(i1, i2) - c
443 0 : t = ks + y
444 0 : c = (t - ks) - y
445 0 : ks = t
446 : END DO
447 : END DO
448 : END IF
449 0 : END FUNCTION kahan_sum_z2
450 :
451 : ! **************************************************************************************************
452 : !> \brief computes the accurate sum of an array that is the element wise product of two input arrays.
453 : !> \param array1 2-D array of complex numbers
454 : !> \param array2 another 2-D array of complex numbers
455 : !> \return dot product
456 : ! **************************************************************************************************
457 31315 : PURE FUNCTION kahan_dot_product_z2(array1, array2) RESULT(ks)
458 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(in) :: array1, array2
459 : COMPLEX(KIND=dp) :: ks
460 :
461 : COMPLEX(KIND=dp) :: c, t, y
462 : INTEGER :: i1, i2, n1, n2
463 :
464 31315 : ks = zzero; t = zzero; y = zzero; c = zzero
465 :
466 31315 : n1 = SIZE(array1, 1)
467 31315 : n2 = SIZE(array1, 2)
468 628553 : DO i2 = 1, n2
469 15815599 : DO i1 = 1, n1
470 15187046 : y = array1(i1, i2)*array2(i1, i2) - c
471 15187046 : t = ks + y
472 15187046 : c = (t - ks) - y
473 15784284 : ks = t
474 : END DO
475 : END DO
476 31315 : END FUNCTION kahan_dot_product_z2
477 :
478 : ! **************************************************************************************************
479 : !> \brief ...
480 : !> \param array ...
481 : !> \param mask ...
482 : !> \return ...
483 : ! **************************************************************************************************
484 0 : PURE FUNCTION kahan_sum_s3(array, mask) RESULT(ks)
485 : REAL(KIND=sp), DIMENSION(:, :, :), INTENT(IN) :: array
486 : LOGICAL, DIMENSION(:, :, :), INTENT(IN), OPTIONAL :: mask
487 : REAL(KIND=sp) :: ks
488 :
489 : INTEGER :: i1, i2, i3
490 : REAL(KIND=sp) :: c, t, y
491 :
492 0 : ks = szero; t = szero; y = szero; c = szero
493 :
494 0 : IF (PRESENT(mask)) THEN
495 0 : DO i3 = 1, SIZE(array, 3)
496 0 : DO i2 = 1, SIZE(array, 2)
497 0 : DO i1 = 1, SIZE(array, 1)
498 0 : IF (mask(i1, i2, i3)) THEN
499 0 : y = array(i1, i2, i3) - c
500 0 : t = ks + y
501 0 : c = (t - ks) - y
502 0 : ks = t
503 : END IF
504 : END DO
505 : END DO
506 : END DO
507 : ELSE
508 0 : DO i3 = 1, SIZE(array, 3)
509 0 : DO i2 = 1, SIZE(array, 2)
510 0 : DO i1 = 1, SIZE(array, 1)
511 0 : y = array(i1, i2, i3) - c
512 0 : t = ks + y
513 0 : c = (t - ks) - y
514 0 : ks = t
515 : END DO
516 : END DO
517 : END DO
518 : END IF
519 0 : END FUNCTION kahan_sum_s3
520 :
521 : ! **************************************************************************************************
522 : !> \brief ...
523 : !> \param array ...
524 : !> \param mask ...
525 : !> \return ...
526 : ! **************************************************************************************************
527 357305 : PURE FUNCTION kahan_sum_d3(array, mask) RESULT(ks)
528 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: array
529 : LOGICAL, DIMENSION(:, :, :), INTENT(IN), OPTIONAL :: mask
530 : REAL(KIND=dp) :: ks
531 :
532 : INTEGER :: i1, i2, i3
533 : REAL(KIND=dp) :: c, t, y
534 :
535 357305 : ks = dzero; t = dzero; y = dzero; c = dzero
536 :
537 357305 : IF (PRESENT(mask)) THEN
538 0 : DO i3 = 1, SIZE(array, 3)
539 0 : DO i2 = 1, SIZE(array, 2)
540 0 : DO i1 = 1, SIZE(array, 1)
541 0 : IF (mask(i1, i2, i3)) THEN
542 0 : y = array(i1, i2, i3) - c
543 0 : t = ks + y
544 0 : c = (t - ks) - y
545 0 : ks = t
546 : END IF
547 : END DO
548 : END DO
549 : END DO
550 : ELSE
551 16465980 : DO i3 = 1, SIZE(array, 3)
552 789903635 : DO i2 = 1, SIZE(array, 2)
553 21809179686 : DO i1 = 1, SIZE(array, 1)
554 21019633356 : y = array(i1, i2, i3) - c
555 21019633356 : t = ks + y
556 21019633356 : c = (t - ks) - y
557 21793071011 : ks = t
558 : END DO
559 : END DO
560 : END DO
561 : END IF
562 357305 : END FUNCTION kahan_sum_d3
563 :
564 : ! **************************************************************************************************
565 : !> \brief computes the accurate sum of an array that is the element wise product of two input arrays.
566 : !> \param array1 3-D array of real numbers
567 : !> \param array2 another 3-D array of real numbers
568 : !> \return dot product
569 : ! **************************************************************************************************
570 476392 : PURE FUNCTION kahan_dot_product_d3(array1, array2) RESULT(ks)
571 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(in) :: array1, array2
572 : REAL(KIND=dp) :: ks
573 :
574 : INTEGER :: i1, i2, i3, n1, n2, n3
575 : REAL(KIND=dp) :: c, t, y
576 :
577 476392 : ks = dzero; t = dzero; y = dzero; c = dzero
578 :
579 476392 : n1 = SIZE(array1, 1)
580 476392 : n2 = SIZE(array1, 2)
581 476392 : n3 = SIZE(array1, 3)
582 5696042 : DO i3 = 1, n3
583 138896218 : DO i2 = 1, n2
584 3310383987 : DO i1 = 1, n1
585 3171964161 : y = array1(i1, i2, i3)*array2(i1, i2, i3) - c
586 3171964161 : t = ks + y
587 3171964161 : c = (t - ks) - y
588 3305164337 : ks = t
589 : END DO
590 : END DO
591 : END DO
592 476392 : END FUNCTION kahan_dot_product_d3
593 :
594 : ! **************************************************************************************************
595 : !> \brief computes the accurate sum of an array that is the element wise product of two input arrays.
596 : !> a mask array determines which product array points to include in the sum
597 : !> \param array1 the first input array to compute the product array
598 : !> \param array2 the second input array to compute the product array
599 : !> \param mask the mask array
600 : !> \param th screening threshold: only array points where the value of mask is greater than th are
601 : !> included in the sum
602 : !> \return the result of summation
603 : ! **************************************************************************************************
604 180 : PURE FUNCTION kahan_dot_product_masked_d3(array1, array2, mask, th) RESULT(ks)
605 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN), &
606 : POINTER :: array1, array2, mask
607 : REAL(KIND=dp), INTENT(in) :: th
608 : REAL(KIND=dp) :: ks
609 :
610 : INTEGER :: i1, i2, i3
611 : REAL(KIND=dp) :: c, t, y
612 :
613 180 : ks = dzero; t = dzero; y = dzero; c = dzero
614 8636 : DO i3 = LBOUND(mask, 3), UBOUND(mask, 3)
615 431252 : DO i2 = LBOUND(mask, 2), UBOUND(mask, 2)
616 22808224 : DO i1 = LBOUND(mask, 1), UBOUND(mask, 1)
617 21986560 : IF (mask(i1, i2, i3) .GT. th) THEN
618 7234016 : y = array1(i1, i2, i3)*array2(i1, i2, i3) - c
619 7234016 : t = ks + y
620 7234016 : c = (t - ks) - y
621 7234016 : ks = t
622 : END IF
623 : END DO
624 : END DO
625 : END DO
626 180 : END FUNCTION kahan_dot_product_masked_d3
627 :
628 : ! **************************************************************************************************
629 : !> \brief ...
630 : !> \param array ...
631 : !> \param mask ...
632 : !> \return ...
633 : ! **************************************************************************************************
634 0 : PURE FUNCTION kahan_sum_c3(array, mask) RESULT(ks)
635 : COMPLEX(KIND=sp), DIMENSION(:, :, :), INTENT(IN) :: array
636 : LOGICAL, DIMENSION(:, :, :), INTENT(IN), OPTIONAL :: mask
637 : COMPLEX(KIND=sp) :: ks
638 :
639 : COMPLEX(KIND=sp) :: c, t, y
640 : INTEGER :: i1, i2, i3
641 :
642 0 : ks = czero; t = czero; y = czero; c = czero
643 :
644 0 : IF (PRESENT(mask)) THEN
645 0 : DO i3 = 1, SIZE(array, 3)
646 0 : DO i2 = 1, SIZE(array, 2)
647 0 : DO i1 = 1, SIZE(array, 1)
648 0 : IF (mask(i1, i2, i3)) THEN
649 0 : y = array(i1, i2, i3) - c
650 0 : t = ks + y
651 0 : c = (t - ks) - y
652 0 : ks = t
653 : END IF
654 : END DO
655 : END DO
656 : END DO
657 : ELSE
658 0 : DO i3 = 1, SIZE(array, 3)
659 0 : DO i2 = 1, SIZE(array, 2)
660 0 : DO i1 = 1, SIZE(array, 1)
661 0 : y = array(i1, i2, i3) - c
662 0 : t = ks + y
663 0 : c = (t - ks) - y
664 0 : ks = t
665 : END DO
666 : END DO
667 : END DO
668 : END IF
669 0 : END FUNCTION kahan_sum_c3
670 :
671 : ! **************************************************************************************************
672 : !> \brief ...
673 : !> \param array ...
674 : !> \param mask ...
675 : !> \return ...
676 : ! **************************************************************************************************
677 0 : PURE FUNCTION kahan_sum_z3(array, mask) RESULT(ks)
678 : COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: array
679 : LOGICAL, DIMENSION(:, :, :), INTENT(IN), OPTIONAL :: mask
680 : COMPLEX(KIND=dp) :: ks
681 :
682 : COMPLEX(KIND=dp) :: c, t, y
683 : INTEGER :: i1, i2, i3
684 :
685 0 : ks = zzero; t = zzero; y = zzero; c = zzero
686 :
687 0 : IF (PRESENT(mask)) THEN
688 0 : DO i3 = 1, SIZE(array, 3)
689 0 : DO i2 = 1, SIZE(array, 2)
690 0 : DO i1 = 1, SIZE(array, 1)
691 0 : IF (mask(i1, i2, i3)) THEN
692 0 : y = array(i1, i2, i3) - c
693 0 : t = ks + y
694 0 : c = (t - ks) - y
695 0 : ks = t
696 : END IF
697 : END DO
698 : END DO
699 : END DO
700 : ELSE
701 0 : DO i3 = 1, SIZE(array, 3)
702 0 : DO i2 = 1, SIZE(array, 2)
703 0 : DO i1 = 1, SIZE(array, 1)
704 0 : y = array(i1, i2, i3) - c
705 0 : t = ks + y
706 0 : c = (t - ks) - y
707 0 : ks = t
708 : END DO
709 : END DO
710 : END DO
711 : END IF
712 0 : END FUNCTION kahan_sum_z3
713 :
714 : ! **************************************************************************************************
715 : !> \brief ...
716 : !> \param array ...
717 : !> \param mask ...
718 : !> \return ...
719 : ! **************************************************************************************************
720 0 : PURE FUNCTION kahan_sum_s4(array, mask) RESULT(ks)
721 : REAL(KIND=sp), DIMENSION(:, :, :, :), INTENT(IN) :: array
722 : LOGICAL, DIMENSION(:, :, :, :), INTENT(IN), &
723 : OPTIONAL :: mask
724 : REAL(KIND=sp) :: ks
725 :
726 : INTEGER :: i1, i2, i3, i4
727 : REAL(KIND=sp) :: c, t, y
728 :
729 0 : ks = szero; t = szero; y = szero; c = szero
730 :
731 0 : IF (PRESENT(mask)) THEN
732 0 : DO i4 = 1, SIZE(array, 4)
733 0 : DO i3 = 1, SIZE(array, 3)
734 0 : DO i2 = 1, SIZE(array, 2)
735 0 : DO i1 = 1, SIZE(array, 1)
736 0 : IF (mask(i1, i2, i3, i4)) THEN
737 0 : y = array(i1, i2, i3, i4) - c
738 0 : t = ks + y
739 0 : c = (t - ks) - y
740 0 : ks = t
741 : END IF
742 : END DO
743 : END DO
744 : END DO
745 : END DO
746 : ELSE
747 0 : DO i4 = 1, SIZE(array, 4)
748 0 : DO i3 = 1, SIZE(array, 3)
749 0 : DO i2 = 1, SIZE(array, 2)
750 0 : DO i1 = 1, SIZE(array, 1)
751 0 : y = array(i1, i2, i3, i4) - c
752 0 : t = ks + y
753 0 : c = (t - ks) - y
754 0 : ks = t
755 : END DO
756 : END DO
757 : END DO
758 : END DO
759 : END IF
760 0 : END FUNCTION kahan_sum_s4
761 :
762 : ! **************************************************************************************************
763 : !> \brief ...
764 : !> \param array ...
765 : !> \param mask ...
766 : !> \return ...
767 : ! **************************************************************************************************
768 0 : PURE FUNCTION kahan_sum_d4(array, mask) RESULT(ks)
769 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: array
770 : LOGICAL, DIMENSION(:, :, :, :), INTENT(IN), &
771 : OPTIONAL :: mask
772 : REAL(KIND=dp) :: ks
773 :
774 : INTEGER :: i1, i2, i3, i4
775 : REAL(KIND=dp) :: c, t, y
776 :
777 0 : ks = dzero; t = dzero; y = dzero; c = dzero
778 :
779 0 : IF (PRESENT(mask)) THEN
780 0 : DO i4 = 1, SIZE(array, 4)
781 0 : DO i3 = 1, SIZE(array, 3)
782 0 : DO i2 = 1, SIZE(array, 2)
783 0 : DO i1 = 1, SIZE(array, 1)
784 0 : IF (mask(i1, i2, i3, i4)) THEN
785 0 : y = array(i1, i2, i3, i4) - c
786 0 : t = ks + y
787 0 : c = (t - ks) - y
788 0 : ks = t
789 : END IF
790 : END DO
791 : END DO
792 : END DO
793 : END DO
794 : ELSE
795 0 : DO i4 = 1, SIZE(array, 4)
796 0 : DO i3 = 1, SIZE(array, 3)
797 0 : DO i2 = 1, SIZE(array, 2)
798 0 : DO i1 = 1, SIZE(array, 1)
799 0 : y = array(i1, i2, i3, i4) - c
800 0 : t = ks + y
801 0 : c = (t - ks) - y
802 0 : ks = t
803 : END DO
804 : END DO
805 : END DO
806 : END DO
807 : END IF
808 0 : END FUNCTION kahan_sum_d4
809 :
810 : ! **************************************************************************************************
811 : !> \brief ...
812 : !> \param array ...
813 : !> \param mask ...
814 : !> \return ...
815 : ! **************************************************************************************************
816 0 : PURE FUNCTION kahan_sum_c4(array, mask) RESULT(ks)
817 : COMPLEX(KIND=sp), DIMENSION(:, :, :, :), &
818 : INTENT(IN) :: array
819 : LOGICAL, DIMENSION(:, :, :, :), INTENT(IN), &
820 : OPTIONAL :: mask
821 : COMPLEX(KIND=sp) :: ks
822 :
823 : COMPLEX(KIND=sp) :: c, t, y
824 : INTEGER :: i1, i2, i3, i4
825 :
826 0 : ks = czero; t = czero; y = czero; c = czero
827 :
828 0 : IF (PRESENT(mask)) THEN
829 0 : DO i4 = 1, SIZE(array, 4)
830 0 : DO i3 = 1, SIZE(array, 3)
831 0 : DO i2 = 1, SIZE(array, 2)
832 0 : DO i1 = 1, SIZE(array, 1)
833 0 : IF (mask(i1, i2, i3, i4)) THEN
834 0 : y = array(i1, i2, i3, i4) - c
835 0 : t = ks + y
836 0 : c = (t - ks) - y
837 0 : ks = t
838 : END IF
839 : END DO
840 : END DO
841 : END DO
842 : END DO
843 : ELSE
844 0 : DO i4 = 1, SIZE(array, 4)
845 0 : DO i3 = 1, SIZE(array, 3)
846 0 : DO i2 = 1, SIZE(array, 2)
847 0 : DO i1 = 1, SIZE(array, 1)
848 0 : y = array(i1, i2, i3, i4) - c
849 0 : t = ks + y
850 0 : c = (t - ks) - y
851 0 : ks = t
852 : END DO
853 : END DO
854 : END DO
855 : END DO
856 : END IF
857 0 : END FUNCTION kahan_sum_c4
858 :
859 : ! **************************************************************************************************
860 : !> \brief ...
861 : !> \param array ...
862 : !> \param mask ...
863 : !> \return ...
864 : ! **************************************************************************************************
865 0 : PURE FUNCTION kahan_sum_z4(array, mask) RESULT(ks)
866 : COMPLEX(KIND=dp), DIMENSION(:, :, :, :), &
867 : INTENT(IN) :: array
868 : LOGICAL, DIMENSION(:, :, :, :), INTENT(IN), &
869 : OPTIONAL :: mask
870 : COMPLEX(KIND=dp) :: ks
871 :
872 : COMPLEX(KIND=dp) :: c, t, y
873 : INTEGER :: i1, i2, i3, i4
874 :
875 0 : ks = zzero; t = zzero; y = zzero; c = zzero
876 :
877 0 : IF (PRESENT(mask)) THEN
878 0 : DO i4 = 1, SIZE(array, 4)
879 0 : DO i3 = 1, SIZE(array, 3)
880 0 : DO i2 = 1, SIZE(array, 2)
881 0 : DO i1 = 1, SIZE(array, 1)
882 0 : IF (mask(i1, i2, i3, i4)) THEN
883 0 : y = array(i1, i2, i3, i4) - c
884 0 : t = ks + y
885 0 : c = (t - ks) - y
886 0 : ks = t
887 : END IF
888 : END DO
889 : END DO
890 : END DO
891 : END DO
892 : ELSE
893 0 : DO i4 = 1, SIZE(array, 4)
894 0 : DO i3 = 1, SIZE(array, 3)
895 0 : DO i2 = 1, SIZE(array, 2)
896 0 : DO i1 = 1, SIZE(array, 1)
897 0 : y = array(i1, i2, i3, i4) - c
898 0 : t = ks + y
899 0 : c = (t - ks) - y
900 0 : ks = t
901 : END DO
902 : END DO
903 : END DO
904 : END DO
905 : END IF
906 0 : END FUNCTION kahan_sum_z4
907 :
908 : ! **************************************************************************************************
909 : !> \brief ...
910 : !> \param array ...
911 : !> \param mask ...
912 : !> \return ...
913 : ! **************************************************************************************************
914 0 : PURE FUNCTION kahan_sum_s5(array, mask) RESULT(ks)
915 : REAL(KIND=sp), DIMENSION(:, :, :, :, :), &
916 : INTENT(IN) :: array
917 : LOGICAL, DIMENSION(:, :, :, :, :), INTENT(IN), &
918 : OPTIONAL :: mask
919 : REAL(KIND=sp) :: ks
920 :
921 : INTEGER :: i1, i2, i3, i4, i5
922 : REAL(KIND=sp) :: c, t, y
923 :
924 0 : ks = szero; t = szero; y = szero; c = szero
925 :
926 0 : IF (PRESENT(mask)) THEN
927 0 : DO i5 = 1, SIZE(array, 5)
928 0 : DO i4 = 1, SIZE(array, 4)
929 0 : DO i3 = 1, SIZE(array, 3)
930 0 : DO i2 = 1, SIZE(array, 2)
931 0 : DO i1 = 1, SIZE(array, 1)
932 0 : IF (mask(i1, i2, i3, i4, i5)) THEN
933 0 : y = array(i1, i2, i3, i4, i5) - c
934 0 : t = ks + y
935 0 : c = (t - ks) - y
936 0 : ks = t
937 : END IF
938 : END DO
939 : END DO
940 : END DO
941 : END DO
942 : END DO
943 : ELSE
944 0 : DO i5 = 1, SIZE(array, 5)
945 0 : DO i4 = 1, SIZE(array, 4)
946 0 : DO i3 = 1, SIZE(array, 3)
947 0 : DO i2 = 1, SIZE(array, 2)
948 0 : DO i1 = 1, SIZE(array, 1)
949 0 : y = array(i1, i2, i3, i4, i5) - c
950 0 : t = ks + y
951 0 : c = (t - ks) - y
952 0 : ks = t
953 : END DO
954 : END DO
955 : END DO
956 : END DO
957 : END DO
958 : END IF
959 0 : END FUNCTION kahan_sum_s5
960 :
961 : ! **************************************************************************************************
962 : !> \brief ...
963 : !> \param array ...
964 : !> \param mask ...
965 : !> \return ...
966 : ! **************************************************************************************************
967 0 : PURE FUNCTION kahan_sum_d5(array, mask) RESULT(ks)
968 : REAL(KIND=dp), DIMENSION(:, :, :, :, :), &
969 : INTENT(IN) :: array
970 : LOGICAL, DIMENSION(:, :, :, :, :), INTENT(IN), &
971 : OPTIONAL :: mask
972 : REAL(KIND=dp) :: ks
973 :
974 : INTEGER :: i1, i2, i3, i4, i5
975 : REAL(KIND=dp) :: c, t, y
976 :
977 0 : ks = dzero; t = dzero; y = dzero; c = dzero
978 :
979 0 : IF (PRESENT(mask)) THEN
980 0 : DO i5 = 1, SIZE(array, 5)
981 0 : DO i4 = 1, SIZE(array, 4)
982 0 : DO i3 = 1, SIZE(array, 3)
983 0 : DO i2 = 1, SIZE(array, 2)
984 0 : DO i1 = 1, SIZE(array, 1)
985 0 : IF (mask(i1, i2, i3, i4, i5)) THEN
986 0 : y = array(i1, i2, i3, i4, i5) - c
987 0 : t = ks + y
988 0 : c = (t - ks) - y
989 0 : ks = t
990 : END IF
991 : END DO
992 : END DO
993 : END DO
994 : END DO
995 : END DO
996 : ELSE
997 0 : DO i5 = 1, SIZE(array, 5)
998 0 : DO i4 = 1, SIZE(array, 4)
999 0 : DO i3 = 1, SIZE(array, 3)
1000 0 : DO i2 = 1, SIZE(array, 2)
1001 0 : DO i1 = 1, SIZE(array, 1)
1002 0 : y = array(i1, i2, i3, i4, i5) - c
1003 0 : t = ks + y
1004 0 : c = (t - ks) - y
1005 0 : ks = t
1006 : END DO
1007 : END DO
1008 : END DO
1009 : END DO
1010 : END DO
1011 : END IF
1012 0 : END FUNCTION kahan_sum_d5
1013 :
1014 : ! **************************************************************************************************
1015 : !> \brief ...
1016 : !> \param array ...
1017 : !> \param mask ...
1018 : !> \return ...
1019 : ! **************************************************************************************************
1020 0 : PURE FUNCTION kahan_sum_c5(array, mask) RESULT(ks)
1021 : COMPLEX(KIND=sp), DIMENSION(:, :, :, :, :), &
1022 : INTENT(IN) :: array
1023 : LOGICAL, DIMENSION(:, :, :, :, :), INTENT(IN), &
1024 : OPTIONAL :: mask
1025 : COMPLEX(KIND=sp) :: ks
1026 :
1027 : COMPLEX(KIND=sp) :: c, t, y
1028 : INTEGER :: i1, i2, i3, i4, i5
1029 :
1030 0 : ks = czero; t = czero; y = czero; c = czero
1031 :
1032 0 : IF (PRESENT(mask)) THEN
1033 0 : DO i5 = 1, SIZE(array, 5)
1034 0 : DO i4 = 1, SIZE(array, 4)
1035 0 : DO i3 = 1, SIZE(array, 3)
1036 0 : DO i2 = 1, SIZE(array, 2)
1037 0 : DO i1 = 1, SIZE(array, 1)
1038 0 : IF (mask(i1, i2, i3, i4, i5)) THEN
1039 0 : y = array(i1, i2, i3, i4, i5) - c
1040 0 : t = ks + y
1041 0 : c = (t - ks) - y
1042 0 : ks = t
1043 : END IF
1044 : END DO
1045 : END DO
1046 : END DO
1047 : END DO
1048 : END DO
1049 : ELSE
1050 0 : DO i5 = 1, SIZE(array, 5)
1051 0 : DO i4 = 1, SIZE(array, 4)
1052 0 : DO i3 = 1, SIZE(array, 3)
1053 0 : DO i2 = 1, SIZE(array, 2)
1054 0 : DO i1 = 1, SIZE(array, 1)
1055 0 : y = array(i1, i2, i3, i4, i5) - c
1056 0 : t = ks + y
1057 0 : c = (t - ks) - y
1058 0 : ks = t
1059 : END DO
1060 : END DO
1061 : END DO
1062 : END DO
1063 : END DO
1064 : END IF
1065 0 : END FUNCTION kahan_sum_c5
1066 :
1067 : ! **************************************************************************************************
1068 : !> \brief ...
1069 : !> \param array ...
1070 : !> \param mask ...
1071 : !> \return ...
1072 : ! **************************************************************************************************
1073 0 : PURE FUNCTION kahan_sum_z5(array, mask) RESULT(ks)
1074 : COMPLEX(KIND=dp), DIMENSION(:, :, :, :, :), &
1075 : INTENT(IN) :: array
1076 : LOGICAL, DIMENSION(:, :, :, :, :), INTENT(IN), &
1077 : OPTIONAL :: mask
1078 : COMPLEX(KIND=dp) :: ks
1079 :
1080 : COMPLEX(KIND=dp) :: c, t, y
1081 : INTEGER :: i1, i2, i3, i4, i5
1082 :
1083 0 : ks = zzero; t = zzero; y = zzero; c = zzero
1084 :
1085 0 : IF (PRESENT(mask)) THEN
1086 0 : DO i5 = 1, SIZE(array, 5)
1087 0 : DO i4 = 1, SIZE(array, 4)
1088 0 : DO i3 = 1, SIZE(array, 3)
1089 0 : DO i2 = 1, SIZE(array, 2)
1090 0 : DO i1 = 1, SIZE(array, 1)
1091 0 : IF (mask(i1, i2, i3, i4, i5)) THEN
1092 0 : y = array(i1, i2, i3, i4, i5) - c
1093 0 : t = ks + y
1094 0 : c = (t - ks) - y
1095 0 : ks = t
1096 : END IF
1097 : END DO
1098 : END DO
1099 : END DO
1100 : END DO
1101 : END DO
1102 : ELSE
1103 0 : DO i5 = 1, SIZE(array, 5)
1104 0 : DO i4 = 1, SIZE(array, 4)
1105 0 : DO i3 = 1, SIZE(array, 3)
1106 0 : DO i2 = 1, SIZE(array, 2)
1107 0 : DO i1 = 1, SIZE(array, 1)
1108 0 : y = array(i1, i2, i3, i4, i5) - c
1109 0 : t = ks + y
1110 0 : c = (t - ks) - y
1111 0 : ks = t
1112 : END DO
1113 : END DO
1114 : END DO
1115 : END DO
1116 : END DO
1117 : END IF
1118 0 : END FUNCTION kahan_sum_z5
1119 :
1120 : ! **************************************************************************************************
1121 : !> \brief ...
1122 : !> \param array ...
1123 : !> \param mask ...
1124 : !> \return ...
1125 : ! **************************************************************************************************
1126 0 : PURE FUNCTION kahan_sum_s6(array, mask) RESULT(ks)
1127 : REAL(KIND=sp), DIMENSION(:, :, :, :, :, :), &
1128 : INTENT(IN) :: array
1129 : LOGICAL, DIMENSION(:, :, :, :, :, :), INTENT(IN), &
1130 : OPTIONAL :: mask
1131 : REAL(KIND=sp) :: ks
1132 :
1133 : INTEGER :: i1, i2, i3, i4, i5, i6
1134 : REAL(KIND=sp) :: c, t, y
1135 :
1136 0 : ks = szero; t = szero; y = szero; c = szero
1137 :
1138 0 : IF (PRESENT(mask)) THEN
1139 0 : DO i6 = 1, SIZE(array, 6)
1140 0 : DO i5 = 1, SIZE(array, 5)
1141 0 : DO i4 = 1, SIZE(array, 4)
1142 0 : DO i3 = 1, SIZE(array, 3)
1143 0 : DO i2 = 1, SIZE(array, 2)
1144 0 : DO i1 = 1, SIZE(array, 1)
1145 0 : IF (mask(i1, i2, i3, i4, i5, i6)) THEN
1146 0 : y = array(i1, i2, i3, i4, i5, i6) - c
1147 0 : t = ks + y
1148 0 : c = (t - ks) - y
1149 0 : ks = t
1150 : END IF
1151 : END DO
1152 : END DO
1153 : END DO
1154 : END DO
1155 : END DO
1156 : END DO
1157 : ELSE
1158 0 : DO i6 = 1, SIZE(array, 6)
1159 0 : DO i5 = 1, SIZE(array, 5)
1160 0 : DO i4 = 1, SIZE(array, 4)
1161 0 : DO i3 = 1, SIZE(array, 3)
1162 0 : DO i2 = 1, SIZE(array, 2)
1163 0 : DO i1 = 1, SIZE(array, 1)
1164 0 : y = array(i1, i2, i3, i4, i5, i6) - c
1165 0 : t = ks + y
1166 0 : c = (t - ks) - y
1167 0 : ks = t
1168 : END DO
1169 : END DO
1170 : END DO
1171 : END DO
1172 : END DO
1173 : END DO
1174 : END IF
1175 0 : END FUNCTION kahan_sum_s6
1176 :
1177 : ! **************************************************************************************************
1178 : !> \brief ...
1179 : !> \param array ...
1180 : !> \param mask ...
1181 : !> \return ...
1182 : ! **************************************************************************************************
1183 0 : PURE FUNCTION kahan_sum_d6(array, mask) RESULT(ks)
1184 : REAL(KIND=dp), DIMENSION(:, :, :, :, :, :), &
1185 : INTENT(IN) :: array
1186 : LOGICAL, DIMENSION(:, :, :, :, :, :), INTENT(IN), &
1187 : OPTIONAL :: mask
1188 : REAL(KIND=dp) :: ks
1189 :
1190 : INTEGER :: i1, i2, i3, i4, i5, i6
1191 : REAL(KIND=dp) :: c, t, y
1192 :
1193 0 : ks = dzero; t = dzero; y = dzero; c = dzero
1194 :
1195 0 : IF (PRESENT(mask)) THEN
1196 0 : DO i6 = 1, SIZE(array, 6)
1197 0 : DO i5 = 1, SIZE(array, 5)
1198 0 : DO i4 = 1, SIZE(array, 4)
1199 0 : DO i3 = 1, SIZE(array, 3)
1200 0 : DO i2 = 1, SIZE(array, 2)
1201 0 : DO i1 = 1, SIZE(array, 1)
1202 0 : IF (mask(i1, i2, i3, i4, i5, i6)) THEN
1203 0 : y = array(i1, i2, i3, i4, i5, i6) - c
1204 0 : t = ks + y
1205 0 : c = (t - ks) - y
1206 0 : ks = t
1207 : END IF
1208 : END DO
1209 : END DO
1210 : END DO
1211 : END DO
1212 : END DO
1213 : END DO
1214 : ELSE
1215 0 : DO i6 = 1, SIZE(array, 6)
1216 0 : DO i5 = 1, SIZE(array, 5)
1217 0 : DO i4 = 1, SIZE(array, 4)
1218 0 : DO i3 = 1, SIZE(array, 3)
1219 0 : DO i2 = 1, SIZE(array, 2)
1220 0 : DO i1 = 1, SIZE(array, 1)
1221 0 : y = array(i1, i2, i3, i4, i5, i6) - c
1222 0 : t = ks + y
1223 0 : c = (t - ks) - y
1224 0 : ks = t
1225 : END DO
1226 : END DO
1227 : END DO
1228 : END DO
1229 : END DO
1230 : END DO
1231 : END IF
1232 0 : END FUNCTION kahan_sum_d6
1233 :
1234 : ! **************************************************************************************************
1235 : !> \brief ...
1236 : !> \param array ...
1237 : !> \param mask ...
1238 : !> \return ...
1239 : ! **************************************************************************************************
1240 0 : PURE FUNCTION kahan_sum_c6(array, mask) RESULT(ks)
1241 : COMPLEX(KIND=sp), DIMENSION(:, :, :, :, :, :), &
1242 : INTENT(IN) :: array
1243 : LOGICAL, DIMENSION(:, :, :, :, :, :), INTENT(IN), &
1244 : OPTIONAL :: mask
1245 : COMPLEX(KIND=sp) :: ks
1246 :
1247 : COMPLEX(KIND=sp) :: c, t, y
1248 : INTEGER :: i1, i2, i3, i4, i5, i6
1249 :
1250 0 : ks = czero; t = czero; y = czero; c = czero
1251 :
1252 0 : IF (PRESENT(mask)) THEN
1253 0 : DO i6 = 1, SIZE(array, 6)
1254 0 : DO i5 = 1, SIZE(array, 5)
1255 0 : DO i4 = 1, SIZE(array, 4)
1256 0 : DO i3 = 1, SIZE(array, 3)
1257 0 : DO i2 = 1, SIZE(array, 2)
1258 0 : DO i1 = 1, SIZE(array, 1)
1259 0 : IF (mask(i1, i2, i3, i4, i5, i6)) THEN
1260 0 : y = array(i1, i2, i3, i4, i5, i6) - c
1261 0 : t = ks + y
1262 0 : c = (t - ks) - y
1263 0 : ks = t
1264 : END IF
1265 : END DO
1266 : END DO
1267 : END DO
1268 : END DO
1269 : END DO
1270 : END DO
1271 : ELSE
1272 0 : DO i6 = 1, SIZE(array, 6)
1273 0 : DO i5 = 1, SIZE(array, 5)
1274 0 : DO i4 = 1, SIZE(array, 4)
1275 0 : DO i3 = 1, SIZE(array, 3)
1276 0 : DO i2 = 1, SIZE(array, 2)
1277 0 : DO i1 = 1, SIZE(array, 1)
1278 0 : y = array(i1, i2, i3, i4, i5, i6) - c
1279 0 : t = ks + y
1280 0 : c = (t - ks) - y
1281 0 : ks = t
1282 : END DO
1283 : END DO
1284 : END DO
1285 : END DO
1286 : END DO
1287 : END DO
1288 : END IF
1289 0 : END FUNCTION kahan_sum_c6
1290 :
1291 : ! **************************************************************************************************
1292 : !> \brief ...
1293 : !> \param array ...
1294 : !> \param mask ...
1295 : !> \return ...
1296 : ! **************************************************************************************************
1297 0 : PURE FUNCTION kahan_sum_z6(array, mask) RESULT(ks)
1298 : COMPLEX(KIND=dp), DIMENSION(:, :, :, :, :, :), &
1299 : INTENT(IN) :: array
1300 : LOGICAL, DIMENSION(:, :, :, :, :, :), INTENT(IN), &
1301 : OPTIONAL :: mask
1302 : COMPLEX(KIND=dp) :: ks
1303 :
1304 : COMPLEX(KIND=dp) :: c, t, y
1305 : INTEGER :: i1, i2, i3, i4, i5, i6
1306 :
1307 0 : ks = zzero; t = zzero; y = zzero; c = zzero
1308 :
1309 0 : IF (PRESENT(mask)) THEN
1310 0 : DO i6 = 1, SIZE(array, 6)
1311 0 : DO i5 = 1, SIZE(array, 5)
1312 0 : DO i4 = 1, SIZE(array, 4)
1313 0 : DO i3 = 1, SIZE(array, 3)
1314 0 : DO i2 = 1, SIZE(array, 2)
1315 0 : DO i1 = 1, SIZE(array, 1)
1316 0 : IF (mask(i1, i2, i3, i4, i5, i6)) THEN
1317 0 : y = array(i1, i2, i3, i4, i5, i6) - c
1318 0 : t = ks + y
1319 0 : c = (t - ks) - y
1320 0 : ks = t
1321 : END IF
1322 : END DO
1323 : END DO
1324 : END DO
1325 : END DO
1326 : END DO
1327 : END DO
1328 : ELSE
1329 0 : DO i6 = 1, SIZE(array, 6)
1330 0 : DO i5 = 1, SIZE(array, 5)
1331 0 : DO i4 = 1, SIZE(array, 4)
1332 0 : DO i3 = 1, SIZE(array, 3)
1333 0 : DO i2 = 1, SIZE(array, 2)
1334 0 : DO i1 = 1, SIZE(array, 1)
1335 0 : y = array(i1, i2, i3, i4, i5, i6) - c
1336 0 : t = ks + y
1337 0 : c = (t - ks) - y
1338 0 : ks = t
1339 : END DO
1340 : END DO
1341 : END DO
1342 : END DO
1343 : END DO
1344 : END DO
1345 : END IF
1346 0 : END FUNCTION kahan_sum_z6
1347 :
1348 : ! **************************************************************************************************
1349 : !> \brief ...
1350 : !> \param array ...
1351 : !> \param mask ...
1352 : !> \return ...
1353 : ! **************************************************************************************************
1354 0 : PURE FUNCTION kahan_sum_s7(array, mask) RESULT(ks)
1355 : REAL(KIND=sp), DIMENSION(:, :, :, :, :, :, :), &
1356 : INTENT(IN) :: array
1357 : LOGICAL, DIMENSION(:, :, :, :, :, :, :), &
1358 : INTENT(IN), OPTIONAL :: mask
1359 : REAL(KIND=sp) :: ks
1360 :
1361 : INTEGER :: i1, i2, i3, i4, i5, i6, i7
1362 : REAL(KIND=sp) :: c, t, y
1363 :
1364 0 : ks = szero; t = szero; y = szero; c = szero
1365 :
1366 0 : IF (PRESENT(mask)) THEN
1367 0 : DO i7 = 1, SIZE(array, 7)
1368 0 : DO i6 = 1, SIZE(array, 6)
1369 0 : DO i5 = 1, SIZE(array, 5)
1370 0 : DO i4 = 1, SIZE(array, 4)
1371 0 : DO i3 = 1, SIZE(array, 3)
1372 0 : DO i2 = 1, SIZE(array, 2)
1373 0 : DO i1 = 1, SIZE(array, 1)
1374 0 : IF (mask(i1, i2, i3, i4, i5, i6, i7)) THEN
1375 0 : y = array(i1, i2, i3, i4, i5, i6, i7) - c
1376 0 : t = ks + y
1377 0 : c = (t - ks) - y
1378 0 : ks = t
1379 : END IF
1380 : END DO
1381 : END DO
1382 : END DO
1383 : END DO
1384 : END DO
1385 : END DO
1386 : END DO
1387 : ELSE
1388 0 : DO i7 = 1, SIZE(array, 7)
1389 0 : DO i6 = 1, SIZE(array, 6)
1390 0 : DO i5 = 1, SIZE(array, 5)
1391 0 : DO i4 = 1, SIZE(array, 4)
1392 0 : DO i3 = 1, SIZE(array, 3)
1393 0 : DO i2 = 1, SIZE(array, 2)
1394 0 : DO i1 = 1, SIZE(array, 1)
1395 0 : y = array(i1, i2, i3, i4, i5, i6, i7) - c
1396 0 : t = ks + y
1397 0 : c = (t - ks) - y
1398 0 : ks = t
1399 : END DO
1400 : END DO
1401 : END DO
1402 : END DO
1403 : END DO
1404 : END DO
1405 : END DO
1406 : END IF
1407 0 : END FUNCTION kahan_sum_s7
1408 :
1409 : ! **************************************************************************************************
1410 : !> \brief ...
1411 : !> \param array ...
1412 : !> \param mask ...
1413 : !> \return ...
1414 : ! **************************************************************************************************
1415 0 : PURE FUNCTION kahan_sum_d7(array, mask) RESULT(ks)
1416 : REAL(KIND=dp), DIMENSION(:, :, :, :, :, :, :), &
1417 : INTENT(IN) :: array
1418 : LOGICAL, DIMENSION(:, :, :, :, :, :, :), &
1419 : INTENT(IN), OPTIONAL :: mask
1420 : REAL(KIND=dp) :: ks
1421 :
1422 : INTEGER :: i1, i2, i3, i4, i5, i6, i7
1423 : REAL(KIND=dp) :: c, t, y
1424 :
1425 0 : ks = dzero; t = dzero; y = dzero; c = dzero
1426 :
1427 0 : IF (PRESENT(mask)) THEN
1428 0 : DO i7 = 1, SIZE(array, 7)
1429 0 : DO i6 = 1, SIZE(array, 6)
1430 0 : DO i5 = 1, SIZE(array, 5)
1431 0 : DO i4 = 1, SIZE(array, 4)
1432 0 : DO i3 = 1, SIZE(array, 3)
1433 0 : DO i2 = 1, SIZE(array, 2)
1434 0 : DO i1 = 1, SIZE(array, 1)
1435 0 : IF (mask(i1, i2, i3, i4, i5, i6, i7)) THEN
1436 0 : y = array(i1, i2, i3, i4, i5, i6, i7) - c
1437 0 : t = ks + y
1438 0 : c = (t - ks) - y
1439 0 : ks = t
1440 : END IF
1441 : END DO
1442 : END DO
1443 : END DO
1444 : END DO
1445 : END DO
1446 : END DO
1447 : END DO
1448 : ELSE
1449 0 : DO i7 = 1, SIZE(array, 7)
1450 0 : DO i6 = 1, SIZE(array, 6)
1451 0 : DO i5 = 1, SIZE(array, 5)
1452 0 : DO i4 = 1, SIZE(array, 4)
1453 0 : DO i3 = 1, SIZE(array, 3)
1454 0 : DO i2 = 1, SIZE(array, 2)
1455 0 : DO i1 = 1, SIZE(array, 1)
1456 0 : y = array(i1, i2, i3, i4, i5, i6, i7) - c
1457 0 : t = ks + y
1458 0 : c = (t - ks) - y
1459 0 : ks = t
1460 : END DO
1461 : END DO
1462 : END DO
1463 : END DO
1464 : END DO
1465 : END DO
1466 : END DO
1467 : END IF
1468 0 : END FUNCTION kahan_sum_d7
1469 :
1470 : ! **************************************************************************************************
1471 : !> \brief ...
1472 : !> \param array ...
1473 : !> \param mask ...
1474 : !> \return ...
1475 : ! **************************************************************************************************
1476 0 : PURE FUNCTION kahan_sum_c7(array, mask) RESULT(ks)
1477 : COMPLEX(KIND=sp), DIMENSION(:, :, :, :, :, :, :), &
1478 : INTENT(IN) :: array
1479 : LOGICAL, DIMENSION(:, :, :, :, :, :, :), &
1480 : INTENT(IN), OPTIONAL :: mask
1481 : COMPLEX(KIND=sp) :: ks
1482 :
1483 : COMPLEX(KIND=sp) :: c, t, y
1484 : INTEGER :: i1, i2, i3, i4, i5, i6, i7
1485 :
1486 0 : ks = czero; t = czero; y = czero; c = czero
1487 :
1488 0 : IF (PRESENT(mask)) THEN
1489 0 : DO i7 = 1, SIZE(array, 7)
1490 0 : DO i6 = 1, SIZE(array, 6)
1491 0 : DO i5 = 1, SIZE(array, 5)
1492 0 : DO i4 = 1, SIZE(array, 4)
1493 0 : DO i3 = 1, SIZE(array, 3)
1494 0 : DO i2 = 1, SIZE(array, 2)
1495 0 : DO i1 = 1, SIZE(array, 1)
1496 0 : IF (mask(i1, i2, i3, i4, i5, i6, i7)) THEN
1497 0 : y = array(i1, i2, i3, i4, i5, i6, i7) - c
1498 0 : t = ks + y
1499 0 : c = (t - ks) - y
1500 0 : ks = t
1501 : END IF
1502 : END DO
1503 : END DO
1504 : END DO
1505 : END DO
1506 : END DO
1507 : END DO
1508 : END DO
1509 : ELSE
1510 0 : DO i7 = 1, SIZE(array, 7)
1511 0 : DO i6 = 1, SIZE(array, 6)
1512 0 : DO i5 = 1, SIZE(array, 5)
1513 0 : DO i4 = 1, SIZE(array, 4)
1514 0 : DO i3 = 1, SIZE(array, 3)
1515 0 : DO i2 = 1, SIZE(array, 2)
1516 0 : DO i1 = 1, SIZE(array, 1)
1517 0 : y = array(i1, i2, i3, i4, i5, i6, i7) - c
1518 0 : t = ks + y
1519 0 : c = (t - ks) - y
1520 0 : ks = t
1521 : END DO
1522 : END DO
1523 : END DO
1524 : END DO
1525 : END DO
1526 : END DO
1527 : END DO
1528 : END IF
1529 0 : END FUNCTION kahan_sum_c7
1530 :
1531 : ! **************************************************************************************************
1532 : !> \brief ...
1533 : !> \param array ...
1534 : !> \param mask ...
1535 : !> \return ...
1536 : ! **************************************************************************************************
1537 0 : PURE FUNCTION kahan_sum_z7(array, mask) RESULT(ks)
1538 : COMPLEX(KIND=dp), DIMENSION(:, :, :, :, :, :, :), &
1539 : INTENT(IN) :: array
1540 : LOGICAL, DIMENSION(:, :, :, :, :, :, :), &
1541 : INTENT(IN), OPTIONAL :: mask
1542 : COMPLEX(KIND=dp) :: ks
1543 :
1544 : COMPLEX(KIND=dp) :: c, t, y
1545 : INTEGER :: i1, i2, i3, i4, i5, i6, i7
1546 :
1547 0 : ks = zzero; t = zzero; y = zzero; c = zzero
1548 :
1549 0 : IF (PRESENT(mask)) THEN
1550 0 : DO i7 = 1, SIZE(array, 7)
1551 0 : DO i6 = 1, SIZE(array, 6)
1552 0 : DO i5 = 1, SIZE(array, 5)
1553 0 : DO i4 = 1, SIZE(array, 4)
1554 0 : DO i3 = 1, SIZE(array, 3)
1555 0 : DO i2 = 1, SIZE(array, 2)
1556 0 : DO i1 = 1, SIZE(array, 1)
1557 0 : IF (mask(i1, i2, i3, i4, i5, i6, i7)) THEN
1558 0 : y = array(i1, i2, i3, i4, i5, i6, i7) - c
1559 0 : t = ks + y
1560 0 : c = (t - ks) - y
1561 0 : ks = t
1562 : END IF
1563 : END DO
1564 : END DO
1565 : END DO
1566 : END DO
1567 : END DO
1568 : END DO
1569 : END DO
1570 : ELSE
1571 0 : DO i7 = 1, SIZE(array, 7)
1572 0 : DO i6 = 1, SIZE(array, 6)
1573 0 : DO i5 = 1, SIZE(array, 5)
1574 0 : DO i4 = 1, SIZE(array, 4)
1575 0 : DO i3 = 1, SIZE(array, 3)
1576 0 : DO i2 = 1, SIZE(array, 2)
1577 0 : DO i1 = 1, SIZE(array, 1)
1578 0 : y = array(i1, i2, i3, i4, i5, i6, i7) - c
1579 0 : t = ks + y
1580 0 : c = (t - ks) - y
1581 0 : ks = t
1582 : END DO
1583 : END DO
1584 : END DO
1585 : END DO
1586 : END DO
1587 : END DO
1588 : END DO
1589 : END IF
1590 0 : END FUNCTION kahan_sum_z7
1591 :
1592 : ! **************************************************************************************************
1593 : !> \brief computes the accurate sum of blocks of regular dot products
1594 : !> \param array1 array of real numbers
1595 : !> \param array2 another array of real numbers
1596 : !> \param blksize ...
1597 : !> \return dot product
1598 : ! **************************************************************************************************
1599 6720 : FUNCTION kahan_blocked_dot_product_d1(array1, array2, blksize) RESULT(ks)
1600 : REAL(KIND=dp), DIMENSION(:), INTENT(in) :: array1, array2
1601 : INTEGER, INTENT(IN), OPTIONAL :: blksize
1602 : REAL(KIND=dp) :: ks
1603 :
1604 : INTEGER :: my_blksize
1605 : REAL(KIND=dp) :: DDOT
1606 :
1607 6720 : my_blksize = 32
1608 6720 : IF (PRESENT(blksize)) my_blksize = blksize
1609 :
1610 6720 : IF (my_blksize <= 1) THEN
1611 : ! The original should be faster
1612 0 : ks = accurate_dot_product(array1, array2)
1613 6720 : ELSE IF (my_blksize >= SIZE(array1)) THEN
1614 : ! Just use standard dot product from BLAS for performance
1615 6720 : ks = DDOT(SIZE(array1), array1(1), 1, array2(1), 1)
1616 : ELSE
1617 0 : ks = 0.0_dp
1618 : BLOCK
1619 : INTEGER :: i, n, stripesize
1620 : REAL(KIND=dp) :: c, dotproduct, t, y
1621 0 : t = dzero; y = dzero; c = dzero
1622 :
1623 0 : n = SIZE(array1)
1624 0 : DO i = 1, n, my_blksize
1625 : ! Remove 1 to save an operation in the length
1626 0 : stripesize = MIN(my_blksize, n - i + 1)
1627 : ! Perform dot product on the given stripe
1628 0 : dotproduct = DDOT(stripesize, array1(i), 1, array2(i), 1)
1629 0 : y = dotproduct - c
1630 0 : t = ks + y
1631 0 : c = (t - ks) - y
1632 0 : ks = t
1633 : END DO
1634 : END BLOCK
1635 : END IF
1636 6720 : END FUNCTION kahan_blocked_dot_product_d1
1637 :
1638 : END MODULE kahan_sum
|
