Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 3458454 : 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 3458454 : ks = dzero; t = dzero; y = dzero; c = dzero
114 :
115 3458454 : 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 12488346985 : DO i1 = 1, SIZE(array, 1)
126 12484888531 : y = array(i1) - c
127 12484888531 : t = ks + y
128 12484888531 : c = (t - ks) - y
129 12488346985 : ks = t
130 : END DO
131 : END IF
132 3458454 : 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 25426 : 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 25426 : t = dzero; y = dzero; c = dzero; ks_local = dzero
148 :
149 25426 : n = SIZE(array1)
150 101702 : DO i = 1, MOD(n, dblksize)
151 76276 : y(1) = array1(i)*array2(i) - c(1)
152 76276 : t(1) = ks_local(1) + y(1)
153 76276 : c(1) = (t(1) - ks_local(1)) - y(1)
154 101702 : ks_local(1) = t(1)
155 : END DO
156 25426 : DO i = MOD(n, dblksize) + 1, n, dblksize
157 4896234 : y = array1(i:i + (dblksize - 1))*array2(i:i + (dblksize - 1)) - c
158 4896234 : t = ks_local + y
159 4896234 : c = (t - ks_local) - y
160 544026 : ks_local = t
161 : END DO
162 203408 : DO i = 2, dblksize
163 177982 : y(1) = ks_local(i) - (c(1) + c(i))
164 177982 : t(1) = ks_local(1) + y(1)
165 177982 : c(1) = (t(1) - ks_local(1)) - y(1)
166 203408 : ks_local(1) = t(1)
167 : END DO
168 25426 : ks = ks_local(1)
169 25426 : 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 0 : 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 0 : ks = zzero; t = zzero; y = zzero; c = zzero
221 :
222 0 : 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 0 : DO i1 = 1, SIZE(array, 1)
233 0 : y = array(i1) - c
234 0 : t = ks + y
235 0 : c = (t - ks) - y
236 0 : ks = t
237 : END DO
238 : END IF
239 0 : 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 20874 : 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 20874 : ks = dzero; t = dzero; y = dzero; c = dzero
322 :
323 20874 : 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 1311306 : DO i2 = 1, SIZE(array, 2)
336 44934346 : DO i1 = 1, SIZE(array, 1)
337 43623040 : y = array(i1, i2) - c
338 43623040 : t = ks + y
339 43623040 : c = (t - ks) - y
340 44913472 : ks = t
341 : END DO
342 : END DO
343 : END IF
344 20874 : 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 627200 : 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 627200 : ks = dzero; t = dzero; y = dzero; c = dzero
360 :
361 627200 : n1 = SIZE(array1, 1)
362 627200 : n2 = SIZE(array1, 2)
363 8602514 : DO i2 = 1, n2
364 205395574 : DO i1 = 1, n1
365 196793060 : y = array1(i1, i2)*array2(i1, i2) - c
366 196793060 : t = ks + y
367 196793060 : c = (t - ks) - y
368 204768374 : ks = t
369 : END DO
370 : END DO
371 627200 : 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 : !> \param conjugate_a conjugate the first array before multiplication (default: false)
456 : !> \return dot product
457 : ! **************************************************************************************************
458 167358 : PURE FUNCTION kahan_dot_product_z2(array1, array2, conjugate_a) RESULT(ks)
459 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(in) :: array1, array2
460 : LOGICAL, INTENT(IN), OPTIONAL :: conjugate_a
461 : COMPLEX(KIND=dp) :: ks
462 :
463 : COMPLEX(KIND=dp) :: c, t, value, y
464 : INTEGER :: i1, i2, n1, n2
465 : LOGICAL :: conjugate
466 :
467 167358 : ks = zzero; t = zzero; y = zzero; c = zzero
468 167358 : conjugate = .FALSE.
469 167358 : IF (PRESENT(conjugate_a)) conjugate = conjugate_a
470 :
471 167358 : n1 = SIZE(array1, 1)
472 167358 : n2 = SIZE(array1, 2)
473 4010090 : DO i2 = 1, n2
474 160948476 : DO i1 = 1, n1
475 156938386 : value = array1(i1, i2)
476 156938386 : IF (conjugate) value = CONJG(value)
477 156938386 : y = value*array2(i1, i2) - c
478 156938386 : t = ks + y
479 156938386 : c = (t - ks) - y
480 160781118 : ks = t
481 : END DO
482 : END DO
483 167358 : END FUNCTION kahan_dot_product_z2
484 :
485 : ! **************************************************************************************************
486 : !> \brief ...
487 : !> \param array ...
488 : !> \param mask ...
489 : !> \return ...
490 : ! **************************************************************************************************
491 0 : PURE FUNCTION kahan_sum_s3(array, mask) RESULT(ks)
492 : REAL(KIND=sp), DIMENSION(:, :, :), INTENT(IN) :: array
493 : LOGICAL, DIMENSION(:, :, :), INTENT(IN), OPTIONAL :: mask
494 : REAL(KIND=sp) :: ks
495 :
496 : INTEGER :: i1, i2, i3
497 : REAL(KIND=sp) :: c, t, y
498 :
499 0 : ks = szero; t = szero; y = szero; c = szero
500 :
501 0 : IF (PRESENT(mask)) THEN
502 0 : DO i3 = 1, SIZE(array, 3)
503 0 : DO i2 = 1, SIZE(array, 2)
504 0 : DO i1 = 1, SIZE(array, 1)
505 0 : IF (mask(i1, i2, i3)) THEN
506 0 : y = array(i1, i2, i3) - c
507 0 : t = ks + y
508 0 : c = (t - ks) - y
509 0 : ks = t
510 : END IF
511 : END DO
512 : END DO
513 : END DO
514 : ELSE
515 0 : DO i3 = 1, SIZE(array, 3)
516 0 : DO i2 = 1, SIZE(array, 2)
517 0 : DO i1 = 1, SIZE(array, 1)
518 0 : y = array(i1, i2, i3) - c
519 0 : t = ks + y
520 0 : c = (t - ks) - y
521 0 : ks = t
522 : END DO
523 : END DO
524 : END DO
525 : END IF
526 0 : END FUNCTION kahan_sum_s3
527 :
528 : ! **************************************************************************************************
529 : !> \brief ...
530 : !> \param array ...
531 : !> \param mask ...
532 : !> \return ...
533 : ! **************************************************************************************************
534 489673 : PURE FUNCTION kahan_sum_d3(array, mask) RESULT(ks)
535 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: array
536 : LOGICAL, DIMENSION(:, :, :), INTENT(IN), OPTIONAL :: mask
537 : REAL(KIND=dp) :: ks
538 :
539 : INTEGER :: i1, i2, i3
540 : REAL(KIND=dp) :: c, t, y
541 :
542 489673 : ks = dzero; t = dzero; y = dzero; c = dzero
543 :
544 489673 : IF (PRESENT(mask)) THEN
545 0 : DO i3 = 1, SIZE(array, 3)
546 0 : DO i2 = 1, SIZE(array, 2)
547 0 : DO i1 = 1, SIZE(array, 1)
548 0 : IF (mask(i1, i2, i3)) THEN
549 0 : y = array(i1, i2, i3) - c
550 0 : t = ks + y
551 0 : c = (t - ks) - y
552 0 : ks = t
553 : END IF
554 : END DO
555 : END DO
556 : END DO
557 : ELSE
558 22794156 : DO i3 = 1, SIZE(array, 3)
559 1102783473 : DO i2 = 1, SIZE(array, 2)
560 30617283594 : DO i1 = 1, SIZE(array, 1)
561 29514989794 : y = array(i1, i2, i3) - c
562 29514989794 : t = ks + y
563 29514989794 : c = (t - ks) - y
564 30594979111 : ks = t
565 : END DO
566 : END DO
567 : END DO
568 : END IF
569 489673 : END FUNCTION kahan_sum_d3
570 :
571 : ! **************************************************************************************************
572 : !> \brief computes the accurate sum of an array that is the element wise product of two input arrays.
573 : !> \param array1 3-D array of real numbers
574 : !> \param array2 another 3-D array of real numbers
575 : !> \return dot product
576 : ! **************************************************************************************************
577 554306 : PURE FUNCTION kahan_dot_product_d3(array1, array2) RESULT(ks)
578 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(in) :: array1, array2
579 : REAL(KIND=dp) :: ks
580 :
581 : INTEGER :: i1, i2, i3, n1, n2, n3
582 : REAL(KIND=dp) :: c, t, y
583 :
584 554306 : ks = dzero; t = dzero; y = dzero; c = dzero
585 :
586 554306 : n1 = SIZE(array1, 1)
587 554306 : n2 = SIZE(array1, 2)
588 554306 : n3 = SIZE(array1, 3)
589 6829404 : DO i3 = 1, n3
590 182410360 : DO i2 = 1, n2
591 4849012514 : DO i1 = 1, n1
592 4667156460 : y = array1(i1, i2, i3)*array2(i1, i2, i3) - c
593 4667156460 : t = ks + y
594 4667156460 : c = (t - ks) - y
595 4842737416 : ks = t
596 : END DO
597 : END DO
598 : END DO
599 554306 : END FUNCTION kahan_dot_product_d3
600 :
601 : ! **************************************************************************************************
602 : !> \brief computes the accurate sum of an array that is the element wise product of two input arrays.
603 : !> a mask array determines which product array points to include in the sum
604 : !> \param array1 the first input array to compute the product array
605 : !> \param array2 the second input array to compute the product array
606 : !> \param mask the mask array
607 : !> \param th screening threshold: only array points where the value of mask is greater than th are
608 : !> included in the sum
609 : !> \return the result of summation
610 : ! **************************************************************************************************
611 180 : PURE FUNCTION kahan_dot_product_masked_d3(array1, array2, mask, th) RESULT(ks)
612 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN), &
613 : POINTER :: array1, array2, mask
614 : REAL(KIND=dp), INTENT(in) :: th
615 : REAL(KIND=dp) :: ks
616 :
617 : INTEGER :: i1, i2, i3
618 : REAL(KIND=dp) :: c, t, y
619 :
620 180 : ks = dzero; t = dzero; y = dzero; c = dzero
621 8636 : DO i3 = LBOUND(mask, 3), UBOUND(mask, 3)
622 431252 : DO i2 = LBOUND(mask, 2), UBOUND(mask, 2)
623 22808224 : DO i1 = LBOUND(mask, 1), UBOUND(mask, 1)
624 21986560 : IF (mask(i1, i2, i3) > th) THEN
625 7234016 : y = array1(i1, i2, i3)*array2(i1, i2, i3) - c
626 7234016 : t = ks + y
627 7234016 : c = (t - ks) - y
628 7234016 : ks = t
629 : END IF
630 : END DO
631 : END DO
632 : END DO
633 180 : END FUNCTION kahan_dot_product_masked_d3
634 :
635 : ! **************************************************************************************************
636 : !> \brief ...
637 : !> \param array ...
638 : !> \param mask ...
639 : !> \return ...
640 : ! **************************************************************************************************
641 0 : PURE FUNCTION kahan_sum_c3(array, mask) RESULT(ks)
642 : COMPLEX(KIND=sp), DIMENSION(:, :, :), INTENT(IN) :: array
643 : LOGICAL, DIMENSION(:, :, :), INTENT(IN), OPTIONAL :: mask
644 : COMPLEX(KIND=sp) :: ks
645 :
646 : COMPLEX(KIND=sp) :: c, t, y
647 : INTEGER :: i1, i2, i3
648 :
649 0 : ks = czero; t = czero; y = czero; c = czero
650 :
651 0 : IF (PRESENT(mask)) THEN
652 0 : DO i3 = 1, SIZE(array, 3)
653 0 : DO i2 = 1, SIZE(array, 2)
654 0 : DO i1 = 1, SIZE(array, 1)
655 0 : IF (mask(i1, i2, i3)) THEN
656 0 : y = array(i1, i2, i3) - c
657 0 : t = ks + y
658 0 : c = (t - ks) - y
659 0 : ks = t
660 : END IF
661 : END DO
662 : END DO
663 : END DO
664 : ELSE
665 0 : DO i3 = 1, SIZE(array, 3)
666 0 : DO i2 = 1, SIZE(array, 2)
667 0 : DO i1 = 1, SIZE(array, 1)
668 0 : y = array(i1, i2, i3) - c
669 0 : t = ks + y
670 0 : c = (t - ks) - y
671 0 : ks = t
672 : END DO
673 : END DO
674 : END DO
675 : END IF
676 0 : END FUNCTION kahan_sum_c3
677 :
678 : ! **************************************************************************************************
679 : !> \brief ...
680 : !> \param array ...
681 : !> \param mask ...
682 : !> \return ...
683 : ! **************************************************************************************************
684 0 : PURE FUNCTION kahan_sum_z3(array, mask) RESULT(ks)
685 : COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: array
686 : LOGICAL, DIMENSION(:, :, :), INTENT(IN), OPTIONAL :: mask
687 : COMPLEX(KIND=dp) :: ks
688 :
689 : COMPLEX(KIND=dp) :: c, t, y
690 : INTEGER :: i1, i2, i3
691 :
692 0 : ks = zzero; t = zzero; y = zzero; c = zzero
693 :
694 0 : IF (PRESENT(mask)) THEN
695 0 : DO i3 = 1, SIZE(array, 3)
696 0 : DO i2 = 1, SIZE(array, 2)
697 0 : DO i1 = 1, SIZE(array, 1)
698 0 : IF (mask(i1, i2, i3)) THEN
699 0 : y = array(i1, i2, i3) - c
700 0 : t = ks + y
701 0 : c = (t - ks) - y
702 0 : ks = t
703 : END IF
704 : END DO
705 : END DO
706 : END DO
707 : ELSE
708 0 : DO i3 = 1, SIZE(array, 3)
709 0 : DO i2 = 1, SIZE(array, 2)
710 0 : DO i1 = 1, SIZE(array, 1)
711 0 : y = array(i1, i2, i3) - c
712 0 : t = ks + y
713 0 : c = (t - ks) - y
714 0 : ks = t
715 : END DO
716 : END DO
717 : END DO
718 : END IF
719 0 : END FUNCTION kahan_sum_z3
720 :
721 : ! **************************************************************************************************
722 : !> \brief ...
723 : !> \param array ...
724 : !> \param mask ...
725 : !> \return ...
726 : ! **************************************************************************************************
727 0 : PURE FUNCTION kahan_sum_s4(array, mask) RESULT(ks)
728 : REAL(KIND=sp), DIMENSION(:, :, :, :), INTENT(IN) :: array
729 : LOGICAL, DIMENSION(:, :, :, :), INTENT(IN), &
730 : OPTIONAL :: mask
731 : REAL(KIND=sp) :: ks
732 :
733 : INTEGER :: i1, i2, i3, i4
734 : REAL(KIND=sp) :: c, t, y
735 :
736 0 : ks = szero; t = szero; y = szero; c = szero
737 :
738 0 : IF (PRESENT(mask)) THEN
739 0 : DO i4 = 1, SIZE(array, 4)
740 0 : DO i3 = 1, SIZE(array, 3)
741 0 : DO i2 = 1, SIZE(array, 2)
742 0 : DO i1 = 1, SIZE(array, 1)
743 0 : IF (mask(i1, i2, i3, i4)) THEN
744 0 : y = array(i1, i2, i3, i4) - c
745 0 : t = ks + y
746 0 : c = (t - ks) - y
747 0 : ks = t
748 : END IF
749 : END DO
750 : END DO
751 : END DO
752 : END DO
753 : ELSE
754 0 : DO i4 = 1, SIZE(array, 4)
755 0 : DO i3 = 1, SIZE(array, 3)
756 0 : DO i2 = 1, SIZE(array, 2)
757 0 : DO i1 = 1, SIZE(array, 1)
758 0 : y = array(i1, i2, i3, i4) - c
759 0 : t = ks + y
760 0 : c = (t - ks) - y
761 0 : ks = t
762 : END DO
763 : END DO
764 : END DO
765 : END DO
766 : END IF
767 0 : END FUNCTION kahan_sum_s4
768 :
769 : ! **************************************************************************************************
770 : !> \brief ...
771 : !> \param array ...
772 : !> \param mask ...
773 : !> \return ...
774 : ! **************************************************************************************************
775 0 : PURE FUNCTION kahan_sum_d4(array, mask) RESULT(ks)
776 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: array
777 : LOGICAL, DIMENSION(:, :, :, :), INTENT(IN), &
778 : OPTIONAL :: mask
779 : REAL(KIND=dp) :: ks
780 :
781 : INTEGER :: i1, i2, i3, i4
782 : REAL(KIND=dp) :: c, t, y
783 :
784 0 : ks = dzero; t = dzero; y = dzero; c = dzero
785 :
786 0 : IF (PRESENT(mask)) THEN
787 0 : DO i4 = 1, SIZE(array, 4)
788 0 : DO i3 = 1, SIZE(array, 3)
789 0 : DO i2 = 1, SIZE(array, 2)
790 0 : DO i1 = 1, SIZE(array, 1)
791 0 : IF (mask(i1, i2, i3, i4)) THEN
792 0 : y = array(i1, i2, i3, i4) - c
793 0 : t = ks + y
794 0 : c = (t - ks) - y
795 0 : ks = t
796 : END IF
797 : END DO
798 : END DO
799 : END DO
800 : END DO
801 : ELSE
802 0 : DO i4 = 1, SIZE(array, 4)
803 0 : DO i3 = 1, SIZE(array, 3)
804 0 : DO i2 = 1, SIZE(array, 2)
805 0 : DO i1 = 1, SIZE(array, 1)
806 0 : y = array(i1, i2, i3, i4) - c
807 0 : t = ks + y
808 0 : c = (t - ks) - y
809 0 : ks = t
810 : END DO
811 : END DO
812 : END DO
813 : END DO
814 : END IF
815 0 : END FUNCTION kahan_sum_d4
816 :
817 : ! **************************************************************************************************
818 : !> \brief ...
819 : !> \param array ...
820 : !> \param mask ...
821 : !> \return ...
822 : ! **************************************************************************************************
823 0 : PURE FUNCTION kahan_sum_c4(array, mask) RESULT(ks)
824 : COMPLEX(KIND=sp), DIMENSION(:, :, :, :), &
825 : INTENT(IN) :: array
826 : LOGICAL, DIMENSION(:, :, :, :), INTENT(IN), &
827 : OPTIONAL :: mask
828 : COMPLEX(KIND=sp) :: ks
829 :
830 : COMPLEX(KIND=sp) :: c, t, y
831 : INTEGER :: i1, i2, i3, i4
832 :
833 0 : ks = czero; t = czero; y = czero; c = czero
834 :
835 0 : IF (PRESENT(mask)) THEN
836 0 : DO i4 = 1, SIZE(array, 4)
837 0 : DO i3 = 1, SIZE(array, 3)
838 0 : DO i2 = 1, SIZE(array, 2)
839 0 : DO i1 = 1, SIZE(array, 1)
840 0 : IF (mask(i1, i2, i3, i4)) THEN
841 0 : y = array(i1, i2, i3, i4) - c
842 0 : t = ks + y
843 0 : c = (t - ks) - y
844 0 : ks = t
845 : END IF
846 : END DO
847 : END DO
848 : END DO
849 : END DO
850 : ELSE
851 0 : DO i4 = 1, SIZE(array, 4)
852 0 : DO i3 = 1, SIZE(array, 3)
853 0 : DO i2 = 1, SIZE(array, 2)
854 0 : DO i1 = 1, SIZE(array, 1)
855 0 : y = array(i1, i2, i3, i4) - c
856 0 : t = ks + y
857 0 : c = (t - ks) - y
858 0 : ks = t
859 : END DO
860 : END DO
861 : END DO
862 : END DO
863 : END IF
864 0 : END FUNCTION kahan_sum_c4
865 :
866 : ! **************************************************************************************************
867 : !> \brief ...
868 : !> \param array ...
869 : !> \param mask ...
870 : !> \return ...
871 : ! **************************************************************************************************
872 0 : PURE FUNCTION kahan_sum_z4(array, mask) RESULT(ks)
873 : COMPLEX(KIND=dp), DIMENSION(:, :, :, :), &
874 : INTENT(IN) :: array
875 : LOGICAL, DIMENSION(:, :, :, :), INTENT(IN), &
876 : OPTIONAL :: mask
877 : COMPLEX(KIND=dp) :: ks
878 :
879 : COMPLEX(KIND=dp) :: c, t, y
880 : INTEGER :: i1, i2, i3, i4
881 :
882 0 : ks = zzero; t = zzero; y = zzero; c = zzero
883 :
884 0 : IF (PRESENT(mask)) THEN
885 0 : DO i4 = 1, SIZE(array, 4)
886 0 : DO i3 = 1, SIZE(array, 3)
887 0 : DO i2 = 1, SIZE(array, 2)
888 0 : DO i1 = 1, SIZE(array, 1)
889 0 : IF (mask(i1, i2, i3, i4)) THEN
890 0 : y = array(i1, i2, i3, i4) - c
891 0 : t = ks + y
892 0 : c = (t - ks) - y
893 0 : ks = t
894 : END IF
895 : END DO
896 : END DO
897 : END DO
898 : END DO
899 : ELSE
900 0 : DO i4 = 1, SIZE(array, 4)
901 0 : DO i3 = 1, SIZE(array, 3)
902 0 : DO i2 = 1, SIZE(array, 2)
903 0 : DO i1 = 1, SIZE(array, 1)
904 0 : y = array(i1, i2, i3, i4) - c
905 0 : t = ks + y
906 0 : c = (t - ks) - y
907 0 : ks = t
908 : END DO
909 : END DO
910 : END DO
911 : END DO
912 : END IF
913 0 : END FUNCTION kahan_sum_z4
914 :
915 : ! **************************************************************************************************
916 : !> \brief ...
917 : !> \param array ...
918 : !> \param mask ...
919 : !> \return ...
920 : ! **************************************************************************************************
921 0 : PURE FUNCTION kahan_sum_s5(array, mask) RESULT(ks)
922 : REAL(KIND=sp), DIMENSION(:, :, :, :, :), &
923 : INTENT(IN) :: array
924 : LOGICAL, DIMENSION(:, :, :, :, :), INTENT(IN), &
925 : OPTIONAL :: mask
926 : REAL(KIND=sp) :: ks
927 :
928 : INTEGER :: i1, i2, i3, i4, i5
929 : REAL(KIND=sp) :: c, t, y
930 :
931 0 : ks = szero; t = szero; y = szero; c = szero
932 :
933 0 : IF (PRESENT(mask)) THEN
934 0 : DO i5 = 1, SIZE(array, 5)
935 0 : DO i4 = 1, SIZE(array, 4)
936 0 : DO i3 = 1, SIZE(array, 3)
937 0 : DO i2 = 1, SIZE(array, 2)
938 0 : DO i1 = 1, SIZE(array, 1)
939 0 : IF (mask(i1, i2, i3, i4, i5)) THEN
940 0 : y = array(i1, i2, i3, i4, i5) - c
941 0 : t = ks + y
942 0 : c = (t - ks) - y
943 0 : ks = t
944 : END IF
945 : END DO
946 : END DO
947 : END DO
948 : END DO
949 : END DO
950 : ELSE
951 0 : DO i5 = 1, SIZE(array, 5)
952 0 : DO i4 = 1, SIZE(array, 4)
953 0 : DO i3 = 1, SIZE(array, 3)
954 0 : DO i2 = 1, SIZE(array, 2)
955 0 : DO i1 = 1, SIZE(array, 1)
956 0 : y = array(i1, i2, i3, i4, i5) - c
957 0 : t = ks + y
958 0 : c = (t - ks) - y
959 0 : ks = t
960 : END DO
961 : END DO
962 : END DO
963 : END DO
964 : END DO
965 : END IF
966 0 : END FUNCTION kahan_sum_s5
967 :
968 : ! **************************************************************************************************
969 : !> \brief ...
970 : !> \param array ...
971 : !> \param mask ...
972 : !> \return ...
973 : ! **************************************************************************************************
974 0 : PURE FUNCTION kahan_sum_d5(array, mask) RESULT(ks)
975 : REAL(KIND=dp), DIMENSION(:, :, :, :, :), &
976 : INTENT(IN) :: array
977 : LOGICAL, DIMENSION(:, :, :, :, :), INTENT(IN), &
978 : OPTIONAL :: mask
979 : REAL(KIND=dp) :: ks
980 :
981 : INTEGER :: i1, i2, i3, i4, i5
982 : REAL(KIND=dp) :: c, t, y
983 :
984 0 : ks = dzero; t = dzero; y = dzero; c = dzero
985 :
986 0 : IF (PRESENT(mask)) THEN
987 0 : DO i5 = 1, SIZE(array, 5)
988 0 : DO i4 = 1, SIZE(array, 4)
989 0 : DO i3 = 1, SIZE(array, 3)
990 0 : DO i2 = 1, SIZE(array, 2)
991 0 : DO i1 = 1, SIZE(array, 1)
992 0 : IF (mask(i1, i2, i3, i4, i5)) THEN
993 0 : y = array(i1, i2, i3, i4, i5) - c
994 0 : t = ks + y
995 0 : c = (t - ks) - y
996 0 : ks = t
997 : END IF
998 : END DO
999 : END DO
1000 : END DO
1001 : END DO
1002 : END DO
1003 : ELSE
1004 0 : DO i5 = 1, SIZE(array, 5)
1005 0 : DO i4 = 1, SIZE(array, 4)
1006 0 : DO i3 = 1, SIZE(array, 3)
1007 0 : DO i2 = 1, SIZE(array, 2)
1008 0 : DO i1 = 1, SIZE(array, 1)
1009 0 : y = array(i1, i2, i3, i4, i5) - c
1010 0 : t = ks + y
1011 0 : c = (t - ks) - y
1012 0 : ks = t
1013 : END DO
1014 : END DO
1015 : END DO
1016 : END DO
1017 : END DO
1018 : END IF
1019 0 : END FUNCTION kahan_sum_d5
1020 :
1021 : ! **************************************************************************************************
1022 : !> \brief ...
1023 : !> \param array ...
1024 : !> \param mask ...
1025 : !> \return ...
1026 : ! **************************************************************************************************
1027 0 : PURE FUNCTION kahan_sum_c5(array, mask) RESULT(ks)
1028 : COMPLEX(KIND=sp), DIMENSION(:, :, :, :, :), &
1029 : INTENT(IN) :: array
1030 : LOGICAL, DIMENSION(:, :, :, :, :), INTENT(IN), &
1031 : OPTIONAL :: mask
1032 : COMPLEX(KIND=sp) :: ks
1033 :
1034 : COMPLEX(KIND=sp) :: c, t, y
1035 : INTEGER :: i1, i2, i3, i4, i5
1036 :
1037 0 : ks = czero; t = czero; y = czero; c = czero
1038 :
1039 0 : IF (PRESENT(mask)) THEN
1040 0 : DO i5 = 1, SIZE(array, 5)
1041 0 : DO i4 = 1, SIZE(array, 4)
1042 0 : DO i3 = 1, SIZE(array, 3)
1043 0 : DO i2 = 1, SIZE(array, 2)
1044 0 : DO i1 = 1, SIZE(array, 1)
1045 0 : IF (mask(i1, i2, i3, i4, i5)) THEN
1046 0 : y = array(i1, i2, i3, i4, i5) - c
1047 0 : t = ks + y
1048 0 : c = (t - ks) - y
1049 0 : ks = t
1050 : END IF
1051 : END DO
1052 : END DO
1053 : END DO
1054 : END DO
1055 : END DO
1056 : ELSE
1057 0 : DO i5 = 1, SIZE(array, 5)
1058 0 : DO i4 = 1, SIZE(array, 4)
1059 0 : DO i3 = 1, SIZE(array, 3)
1060 0 : DO i2 = 1, SIZE(array, 2)
1061 0 : DO i1 = 1, SIZE(array, 1)
1062 0 : y = array(i1, i2, i3, i4, i5) - c
1063 0 : t = ks + y
1064 0 : c = (t - ks) - y
1065 0 : ks = t
1066 : END DO
1067 : END DO
1068 : END DO
1069 : END DO
1070 : END DO
1071 : END IF
1072 0 : END FUNCTION kahan_sum_c5
1073 :
1074 : ! **************************************************************************************************
1075 : !> \brief ...
1076 : !> \param array ...
1077 : !> \param mask ...
1078 : !> \return ...
1079 : ! **************************************************************************************************
1080 0 : PURE FUNCTION kahan_sum_z5(array, mask) RESULT(ks)
1081 : COMPLEX(KIND=dp), DIMENSION(:, :, :, :, :), &
1082 : INTENT(IN) :: array
1083 : LOGICAL, DIMENSION(:, :, :, :, :), INTENT(IN), &
1084 : OPTIONAL :: mask
1085 : COMPLEX(KIND=dp) :: ks
1086 :
1087 : COMPLEX(KIND=dp) :: c, t, y
1088 : INTEGER :: i1, i2, i3, i4, i5
1089 :
1090 0 : ks = zzero; t = zzero; y = zzero; c = zzero
1091 :
1092 0 : IF (PRESENT(mask)) THEN
1093 0 : DO i5 = 1, SIZE(array, 5)
1094 0 : DO i4 = 1, SIZE(array, 4)
1095 0 : DO i3 = 1, SIZE(array, 3)
1096 0 : DO i2 = 1, SIZE(array, 2)
1097 0 : DO i1 = 1, SIZE(array, 1)
1098 0 : IF (mask(i1, i2, i3, i4, i5)) THEN
1099 0 : y = array(i1, i2, i3, i4, i5) - c
1100 0 : t = ks + y
1101 0 : c = (t - ks) - y
1102 0 : ks = t
1103 : END IF
1104 : END DO
1105 : END DO
1106 : END DO
1107 : END DO
1108 : END DO
1109 : ELSE
1110 0 : DO i5 = 1, SIZE(array, 5)
1111 0 : DO i4 = 1, SIZE(array, 4)
1112 0 : DO i3 = 1, SIZE(array, 3)
1113 0 : DO i2 = 1, SIZE(array, 2)
1114 0 : DO i1 = 1, SIZE(array, 1)
1115 0 : y = array(i1, i2, i3, i4, i5) - c
1116 0 : t = ks + y
1117 0 : c = (t - ks) - y
1118 0 : ks = t
1119 : END DO
1120 : END DO
1121 : END DO
1122 : END DO
1123 : END DO
1124 : END IF
1125 0 : END FUNCTION kahan_sum_z5
1126 :
1127 : ! **************************************************************************************************
1128 : !> \brief ...
1129 : !> \param array ...
1130 : !> \param mask ...
1131 : !> \return ...
1132 : ! **************************************************************************************************
1133 0 : PURE FUNCTION kahan_sum_s6(array, mask) RESULT(ks)
1134 : REAL(KIND=sp), DIMENSION(:, :, :, :, :, :), &
1135 : INTENT(IN) :: array
1136 : LOGICAL, DIMENSION(:, :, :, :, :, :), INTENT(IN), &
1137 : OPTIONAL :: mask
1138 : REAL(KIND=sp) :: ks
1139 :
1140 : INTEGER :: i1, i2, i3, i4, i5, i6
1141 : REAL(KIND=sp) :: c, t, y
1142 :
1143 0 : ks = szero; t = szero; y = szero; c = szero
1144 :
1145 0 : IF (PRESENT(mask)) THEN
1146 0 : DO i6 = 1, SIZE(array, 6)
1147 0 : DO i5 = 1, SIZE(array, 5)
1148 0 : DO i4 = 1, SIZE(array, 4)
1149 0 : DO i3 = 1, SIZE(array, 3)
1150 0 : DO i2 = 1, SIZE(array, 2)
1151 0 : DO i1 = 1, SIZE(array, 1)
1152 0 : IF (mask(i1, i2, i3, i4, i5, i6)) THEN
1153 0 : y = array(i1, i2, i3, i4, i5, i6) - c
1154 0 : t = ks + y
1155 0 : c = (t - ks) - y
1156 0 : ks = t
1157 : END IF
1158 : END DO
1159 : END DO
1160 : END DO
1161 : END DO
1162 : END DO
1163 : END DO
1164 : ELSE
1165 0 : DO i6 = 1, SIZE(array, 6)
1166 0 : DO i5 = 1, SIZE(array, 5)
1167 0 : DO i4 = 1, SIZE(array, 4)
1168 0 : DO i3 = 1, SIZE(array, 3)
1169 0 : DO i2 = 1, SIZE(array, 2)
1170 0 : DO i1 = 1, SIZE(array, 1)
1171 0 : y = array(i1, i2, i3, i4, i5, i6) - c
1172 0 : t = ks + y
1173 0 : c = (t - ks) - y
1174 0 : ks = t
1175 : END DO
1176 : END DO
1177 : END DO
1178 : END DO
1179 : END DO
1180 : END DO
1181 : END IF
1182 0 : END FUNCTION kahan_sum_s6
1183 :
1184 : ! **************************************************************************************************
1185 : !> \brief ...
1186 : !> \param array ...
1187 : !> \param mask ...
1188 : !> \return ...
1189 : ! **************************************************************************************************
1190 0 : PURE FUNCTION kahan_sum_d6(array, mask) RESULT(ks)
1191 : REAL(KIND=dp), DIMENSION(:, :, :, :, :, :), &
1192 : INTENT(IN) :: array
1193 : LOGICAL, DIMENSION(:, :, :, :, :, :), INTENT(IN), &
1194 : OPTIONAL :: mask
1195 : REAL(KIND=dp) :: ks
1196 :
1197 : INTEGER :: i1, i2, i3, i4, i5, i6
1198 : REAL(KIND=dp) :: c, t, y
1199 :
1200 0 : ks = dzero; t = dzero; y = dzero; c = dzero
1201 :
1202 0 : IF (PRESENT(mask)) THEN
1203 0 : DO i6 = 1, SIZE(array, 6)
1204 0 : DO i5 = 1, SIZE(array, 5)
1205 0 : DO i4 = 1, SIZE(array, 4)
1206 0 : DO i3 = 1, SIZE(array, 3)
1207 0 : DO i2 = 1, SIZE(array, 2)
1208 0 : DO i1 = 1, SIZE(array, 1)
1209 0 : IF (mask(i1, i2, i3, i4, i5, i6)) THEN
1210 0 : y = array(i1, i2, i3, i4, i5, i6) - c
1211 0 : t = ks + y
1212 0 : c = (t - ks) - y
1213 0 : ks = t
1214 : END IF
1215 : END DO
1216 : END DO
1217 : END DO
1218 : END DO
1219 : END DO
1220 : END DO
1221 : ELSE
1222 0 : DO i6 = 1, SIZE(array, 6)
1223 0 : DO i5 = 1, SIZE(array, 5)
1224 0 : DO i4 = 1, SIZE(array, 4)
1225 0 : DO i3 = 1, SIZE(array, 3)
1226 0 : DO i2 = 1, SIZE(array, 2)
1227 0 : DO i1 = 1, SIZE(array, 1)
1228 0 : y = array(i1, i2, i3, i4, i5, i6) - c
1229 0 : t = ks + y
1230 0 : c = (t - ks) - y
1231 0 : ks = t
1232 : END DO
1233 : END DO
1234 : END DO
1235 : END DO
1236 : END DO
1237 : END DO
1238 : END IF
1239 0 : END FUNCTION kahan_sum_d6
1240 :
1241 : ! **************************************************************************************************
1242 : !> \brief ...
1243 : !> \param array ...
1244 : !> \param mask ...
1245 : !> \return ...
1246 : ! **************************************************************************************************
1247 0 : PURE FUNCTION kahan_sum_c6(array, mask) RESULT(ks)
1248 : COMPLEX(KIND=sp), DIMENSION(:, :, :, :, :, :), &
1249 : INTENT(IN) :: array
1250 : LOGICAL, DIMENSION(:, :, :, :, :, :), INTENT(IN), &
1251 : OPTIONAL :: mask
1252 : COMPLEX(KIND=sp) :: ks
1253 :
1254 : COMPLEX(KIND=sp) :: c, t, y
1255 : INTEGER :: i1, i2, i3, i4, i5, i6
1256 :
1257 0 : ks = czero; t = czero; y = czero; c = czero
1258 :
1259 0 : IF (PRESENT(mask)) THEN
1260 0 : DO i6 = 1, SIZE(array, 6)
1261 0 : DO i5 = 1, SIZE(array, 5)
1262 0 : DO i4 = 1, SIZE(array, 4)
1263 0 : DO i3 = 1, SIZE(array, 3)
1264 0 : DO i2 = 1, SIZE(array, 2)
1265 0 : DO i1 = 1, SIZE(array, 1)
1266 0 : IF (mask(i1, i2, i3, i4, i5, i6)) THEN
1267 0 : y = array(i1, i2, i3, i4, i5, i6) - c
1268 0 : t = ks + y
1269 0 : c = (t - ks) - y
1270 0 : ks = t
1271 : END IF
1272 : END DO
1273 : END DO
1274 : END DO
1275 : END DO
1276 : END DO
1277 : END DO
1278 : ELSE
1279 0 : DO i6 = 1, SIZE(array, 6)
1280 0 : DO i5 = 1, SIZE(array, 5)
1281 0 : DO i4 = 1, SIZE(array, 4)
1282 0 : DO i3 = 1, SIZE(array, 3)
1283 0 : DO i2 = 1, SIZE(array, 2)
1284 0 : DO i1 = 1, SIZE(array, 1)
1285 0 : y = array(i1, i2, i3, i4, i5, i6) - c
1286 0 : t = ks + y
1287 0 : c = (t - ks) - y
1288 0 : ks = t
1289 : END DO
1290 : END DO
1291 : END DO
1292 : END DO
1293 : END DO
1294 : END DO
1295 : END IF
1296 0 : END FUNCTION kahan_sum_c6
1297 :
1298 : ! **************************************************************************************************
1299 : !> \brief ...
1300 : !> \param array ...
1301 : !> \param mask ...
1302 : !> \return ...
1303 : ! **************************************************************************************************
1304 0 : PURE FUNCTION kahan_sum_z6(array, mask) RESULT(ks)
1305 : COMPLEX(KIND=dp), DIMENSION(:, :, :, :, :, :), &
1306 : INTENT(IN) :: array
1307 : LOGICAL, DIMENSION(:, :, :, :, :, :), INTENT(IN), &
1308 : OPTIONAL :: mask
1309 : COMPLEX(KIND=dp) :: ks
1310 :
1311 : COMPLEX(KIND=dp) :: c, t, y
1312 : INTEGER :: i1, i2, i3, i4, i5, i6
1313 :
1314 0 : ks = zzero; t = zzero; y = zzero; c = zzero
1315 :
1316 0 : IF (PRESENT(mask)) THEN
1317 0 : DO i6 = 1, SIZE(array, 6)
1318 0 : DO i5 = 1, SIZE(array, 5)
1319 0 : DO i4 = 1, SIZE(array, 4)
1320 0 : DO i3 = 1, SIZE(array, 3)
1321 0 : DO i2 = 1, SIZE(array, 2)
1322 0 : DO i1 = 1, SIZE(array, 1)
1323 0 : IF (mask(i1, i2, i3, i4, i5, i6)) THEN
1324 0 : y = array(i1, i2, i3, i4, i5, i6) - c
1325 0 : t = ks + y
1326 0 : c = (t - ks) - y
1327 0 : ks = t
1328 : END IF
1329 : END DO
1330 : END DO
1331 : END DO
1332 : END DO
1333 : END DO
1334 : END DO
1335 : ELSE
1336 0 : DO i6 = 1, SIZE(array, 6)
1337 0 : DO i5 = 1, SIZE(array, 5)
1338 0 : DO i4 = 1, SIZE(array, 4)
1339 0 : DO i3 = 1, SIZE(array, 3)
1340 0 : DO i2 = 1, SIZE(array, 2)
1341 0 : DO i1 = 1, SIZE(array, 1)
1342 0 : y = array(i1, i2, i3, i4, i5, i6) - c
1343 0 : t = ks + y
1344 0 : c = (t - ks) - y
1345 0 : ks = t
1346 : END DO
1347 : END DO
1348 : END DO
1349 : END DO
1350 : END DO
1351 : END DO
1352 : END IF
1353 0 : END FUNCTION kahan_sum_z6
1354 :
1355 : ! **************************************************************************************************
1356 : !> \brief ...
1357 : !> \param array ...
1358 : !> \param mask ...
1359 : !> \return ...
1360 : ! **************************************************************************************************
1361 0 : PURE FUNCTION kahan_sum_s7(array, mask) RESULT(ks)
1362 : REAL(KIND=sp), DIMENSION(:, :, :, :, :, :, :), &
1363 : INTENT(IN) :: array
1364 : LOGICAL, DIMENSION(:, :, :, :, :, :, :), &
1365 : INTENT(IN), OPTIONAL :: mask
1366 : REAL(KIND=sp) :: ks
1367 :
1368 : INTEGER :: i1, i2, i3, i4, i5, i6, i7
1369 : REAL(KIND=sp) :: c, t, y
1370 :
1371 0 : ks = szero; t = szero; y = szero; c = szero
1372 :
1373 0 : IF (PRESENT(mask)) THEN
1374 0 : DO i7 = 1, SIZE(array, 7)
1375 0 : DO i6 = 1, SIZE(array, 6)
1376 0 : DO i5 = 1, SIZE(array, 5)
1377 0 : DO i4 = 1, SIZE(array, 4)
1378 0 : DO i3 = 1, SIZE(array, 3)
1379 0 : DO i2 = 1, SIZE(array, 2)
1380 0 : DO i1 = 1, SIZE(array, 1)
1381 0 : IF (mask(i1, i2, i3, i4, i5, i6, i7)) THEN
1382 0 : y = array(i1, i2, i3, i4, i5, i6, i7) - c
1383 0 : t = ks + y
1384 0 : c = (t - ks) - y
1385 0 : ks = t
1386 : END IF
1387 : END DO
1388 : END DO
1389 : END DO
1390 : END DO
1391 : END DO
1392 : END DO
1393 : END DO
1394 : ELSE
1395 0 : DO i7 = 1, SIZE(array, 7)
1396 0 : DO i6 = 1, SIZE(array, 6)
1397 0 : DO i5 = 1, SIZE(array, 5)
1398 0 : DO i4 = 1, SIZE(array, 4)
1399 0 : DO i3 = 1, SIZE(array, 3)
1400 0 : DO i2 = 1, SIZE(array, 2)
1401 0 : DO i1 = 1, SIZE(array, 1)
1402 0 : y = array(i1, i2, i3, i4, i5, i6, i7) - c
1403 0 : t = ks + y
1404 0 : c = (t - ks) - y
1405 0 : ks = t
1406 : END DO
1407 : END DO
1408 : END DO
1409 : END DO
1410 : END DO
1411 : END DO
1412 : END DO
1413 : END IF
1414 0 : END FUNCTION kahan_sum_s7
1415 :
1416 : ! **************************************************************************************************
1417 : !> \brief ...
1418 : !> \param array ...
1419 : !> \param mask ...
1420 : !> \return ...
1421 : ! **************************************************************************************************
1422 0 : PURE FUNCTION kahan_sum_d7(array, mask) RESULT(ks)
1423 : REAL(KIND=dp), DIMENSION(:, :, :, :, :, :, :), &
1424 : INTENT(IN) :: array
1425 : LOGICAL, DIMENSION(:, :, :, :, :, :, :), &
1426 : INTENT(IN), OPTIONAL :: mask
1427 : REAL(KIND=dp) :: ks
1428 :
1429 : INTEGER :: i1, i2, i3, i4, i5, i6, i7
1430 : REAL(KIND=dp) :: c, t, y
1431 :
1432 0 : ks = dzero; t = dzero; y = dzero; c = dzero
1433 :
1434 0 : IF (PRESENT(mask)) THEN
1435 0 : DO i7 = 1, SIZE(array, 7)
1436 0 : DO i6 = 1, SIZE(array, 6)
1437 0 : DO i5 = 1, SIZE(array, 5)
1438 0 : DO i4 = 1, SIZE(array, 4)
1439 0 : DO i3 = 1, SIZE(array, 3)
1440 0 : DO i2 = 1, SIZE(array, 2)
1441 0 : DO i1 = 1, SIZE(array, 1)
1442 0 : IF (mask(i1, i2, i3, i4, i5, i6, i7)) THEN
1443 0 : y = array(i1, i2, i3, i4, i5, i6, i7) - c
1444 0 : t = ks + y
1445 0 : c = (t - ks) - y
1446 0 : ks = t
1447 : END IF
1448 : END DO
1449 : END DO
1450 : END DO
1451 : END DO
1452 : END DO
1453 : END DO
1454 : END DO
1455 : ELSE
1456 0 : DO i7 = 1, SIZE(array, 7)
1457 0 : DO i6 = 1, SIZE(array, 6)
1458 0 : DO i5 = 1, SIZE(array, 5)
1459 0 : DO i4 = 1, SIZE(array, 4)
1460 0 : DO i3 = 1, SIZE(array, 3)
1461 0 : DO i2 = 1, SIZE(array, 2)
1462 0 : DO i1 = 1, SIZE(array, 1)
1463 0 : y = array(i1, i2, i3, i4, i5, i6, i7) - c
1464 0 : t = ks + y
1465 0 : c = (t - ks) - y
1466 0 : ks = t
1467 : END DO
1468 : END DO
1469 : END DO
1470 : END DO
1471 : END DO
1472 : END DO
1473 : END DO
1474 : END IF
1475 0 : END FUNCTION kahan_sum_d7
1476 :
1477 : ! **************************************************************************************************
1478 : !> \brief ...
1479 : !> \param array ...
1480 : !> \param mask ...
1481 : !> \return ...
1482 : ! **************************************************************************************************
1483 0 : PURE FUNCTION kahan_sum_c7(array, mask) RESULT(ks)
1484 : COMPLEX(KIND=sp), DIMENSION(:, :, :, :, :, :, :), &
1485 : INTENT(IN) :: array
1486 : LOGICAL, DIMENSION(:, :, :, :, :, :, :), &
1487 : INTENT(IN), OPTIONAL :: mask
1488 : COMPLEX(KIND=sp) :: ks
1489 :
1490 : COMPLEX(KIND=sp) :: c, t, y
1491 : INTEGER :: i1, i2, i3, i4, i5, i6, i7
1492 :
1493 0 : ks = czero; t = czero; y = czero; c = czero
1494 :
1495 0 : IF (PRESENT(mask)) THEN
1496 0 : DO i7 = 1, SIZE(array, 7)
1497 0 : DO i6 = 1, SIZE(array, 6)
1498 0 : DO i5 = 1, SIZE(array, 5)
1499 0 : DO i4 = 1, SIZE(array, 4)
1500 0 : DO i3 = 1, SIZE(array, 3)
1501 0 : DO i2 = 1, SIZE(array, 2)
1502 0 : DO i1 = 1, SIZE(array, 1)
1503 0 : IF (mask(i1, i2, i3, i4, i5, i6, i7)) THEN
1504 0 : y = array(i1, i2, i3, i4, i5, i6, i7) - c
1505 0 : t = ks + y
1506 0 : c = (t - ks) - y
1507 0 : ks = t
1508 : END IF
1509 : END DO
1510 : END DO
1511 : END DO
1512 : END DO
1513 : END DO
1514 : END DO
1515 : END DO
1516 : ELSE
1517 0 : DO i7 = 1, SIZE(array, 7)
1518 0 : DO i6 = 1, SIZE(array, 6)
1519 0 : DO i5 = 1, SIZE(array, 5)
1520 0 : DO i4 = 1, SIZE(array, 4)
1521 0 : DO i3 = 1, SIZE(array, 3)
1522 0 : DO i2 = 1, SIZE(array, 2)
1523 0 : DO i1 = 1, SIZE(array, 1)
1524 0 : y = array(i1, i2, i3, i4, i5, i6, i7) - c
1525 0 : t = ks + y
1526 0 : c = (t - ks) - y
1527 0 : ks = t
1528 : END DO
1529 : END DO
1530 : END DO
1531 : END DO
1532 : END DO
1533 : END DO
1534 : END DO
1535 : END IF
1536 0 : END FUNCTION kahan_sum_c7
1537 :
1538 : ! **************************************************************************************************
1539 : !> \brief ...
1540 : !> \param array ...
1541 : !> \param mask ...
1542 : !> \return ...
1543 : ! **************************************************************************************************
1544 0 : PURE FUNCTION kahan_sum_z7(array, mask) RESULT(ks)
1545 : COMPLEX(KIND=dp), DIMENSION(:, :, :, :, :, :, :), &
1546 : INTENT(IN) :: array
1547 : LOGICAL, DIMENSION(:, :, :, :, :, :, :), &
1548 : INTENT(IN), OPTIONAL :: mask
1549 : COMPLEX(KIND=dp) :: ks
1550 :
1551 : COMPLEX(KIND=dp) :: c, t, y
1552 : INTEGER :: i1, i2, i3, i4, i5, i6, i7
1553 :
1554 0 : ks = zzero; t = zzero; y = zzero; c = zzero
1555 :
1556 0 : IF (PRESENT(mask)) THEN
1557 0 : DO i7 = 1, SIZE(array, 7)
1558 0 : DO i6 = 1, SIZE(array, 6)
1559 0 : DO i5 = 1, SIZE(array, 5)
1560 0 : DO i4 = 1, SIZE(array, 4)
1561 0 : DO i3 = 1, SIZE(array, 3)
1562 0 : DO i2 = 1, SIZE(array, 2)
1563 0 : DO i1 = 1, SIZE(array, 1)
1564 0 : IF (mask(i1, i2, i3, i4, i5, i6, i7)) THEN
1565 0 : y = array(i1, i2, i3, i4, i5, i6, i7) - c
1566 0 : t = ks + y
1567 0 : c = (t - ks) - y
1568 0 : ks = t
1569 : END IF
1570 : END DO
1571 : END DO
1572 : END DO
1573 : END DO
1574 : END DO
1575 : END DO
1576 : END DO
1577 : ELSE
1578 0 : DO i7 = 1, SIZE(array, 7)
1579 0 : DO i6 = 1, SIZE(array, 6)
1580 0 : DO i5 = 1, SIZE(array, 5)
1581 0 : DO i4 = 1, SIZE(array, 4)
1582 0 : DO i3 = 1, SIZE(array, 3)
1583 0 : DO i2 = 1, SIZE(array, 2)
1584 0 : DO i1 = 1, SIZE(array, 1)
1585 0 : y = array(i1, i2, i3, i4, i5, i6, i7) - c
1586 0 : t = ks + y
1587 0 : c = (t - ks) - y
1588 0 : ks = t
1589 : END DO
1590 : END DO
1591 : END DO
1592 : END DO
1593 : END DO
1594 : END DO
1595 : END DO
1596 : END IF
1597 0 : END FUNCTION kahan_sum_z7
1598 :
1599 : ! **************************************************************************************************
1600 : !> \brief computes the accurate sum of blocks of regular dot products
1601 : !> \param array1 array of real numbers
1602 : !> \param array2 another array of real numbers
1603 : !> \param blksize ...
1604 : !> \return dot product
1605 : ! **************************************************************************************************
1606 6720 : FUNCTION kahan_blocked_dot_product_d1(array1, array2, blksize) RESULT(ks)
1607 : REAL(KIND=dp), DIMENSION(:), INTENT(in) :: array1, array2
1608 : INTEGER, INTENT(IN), OPTIONAL :: blksize
1609 : REAL(KIND=dp) :: ks
1610 :
1611 : INTEGER :: my_blksize
1612 : REAL(KIND=dp) :: DDOT
1613 :
1614 6720 : my_blksize = 32
1615 6720 : IF (PRESENT(blksize)) my_blksize = blksize
1616 :
1617 6720 : IF (my_blksize <= 1) THEN
1618 : ! The original should be faster
1619 0 : ks = accurate_dot_product(array1, array2)
1620 6720 : ELSE IF (my_blksize >= SIZE(array1)) THEN
1621 : ! Just use standard dot product from BLAS for performance
1622 6720 : ks = DDOT(SIZE(array1), array1(1), 1, array2(1), 1)
1623 : ELSE
1624 0 : ks = 0.0_dp
1625 : BLOCK
1626 : INTEGER :: i, n, stripesize
1627 : REAL(KIND=dp) :: c, dotproduct, t, y
1628 0 : t = dzero; y = dzero; c = dzero
1629 :
1630 0 : n = SIZE(array1)
1631 0 : DO i = 1, n, my_blksize
1632 : ! Remove 1 to save an operation in the length
1633 0 : stripesize = MIN(my_blksize, n - i + 1)
1634 : ! Perform dot product on the given stripe
1635 0 : dotproduct = DDOT(stripesize, array1(i), 1, array2(i), 1)
1636 0 : y = dotproduct - c
1637 0 : t = ks + y
1638 0 : c = (t - ks) - y
1639 0 : ks = t
1640 : END DO
1641 : END BLOCK
1642 : END IF
1643 6720 : END FUNCTION kahan_blocked_dot_product_d1
1644 :
1645 : END MODULE kahan_sum
|
