(git:374b731)
Loading...
Searching...
No Matches
kahan_sum.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 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! **************************************************************************************************
29MODULE 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
63CONTAINS
64! **************************************************************************************************
65!> \brief ...
66!> \param array ...
67!> \param mask ...
68!> \return ...
69! **************************************************************************************************
70 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 ks = szero; t = szero; y = szero; c = szero
79
80 IF (PRESENT(mask)) THEN
81 DO i1 = 1, SIZE(array, 1)
82 IF (mask(i1)) THEN
83 y = array(i1) - c
84 t = ks + y
85 c = (t - ks) - y
86 ks = t
87 END IF
88 END DO
89 ELSE
90 DO i1 = 1, SIZE(array, 1)
91 y = array(i1) - c
92 t = ks + y
93 c = (t - ks) - y
94 ks = t
95 END DO
96 END IF
97 END FUNCTION kahan_sum_s1
98
99! **************************************************************************************************
100!> \brief ...
101!> \param array ...
102!> \param mask ...
103!> \return ...
104! **************************************************************************************************
105 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 ks = dzero; t = dzero; y = dzero; c = dzero
114
115 IF (PRESENT(mask)) THEN
116 DO i1 = 1, SIZE(array, 1)
117 IF (mask(i1)) THEN
118 y = array(i1) - c
119 t = ks + y
120 c = (t - ks) - y
121 ks = t
122 END IF
123 END DO
124 ELSE
125 DO i1 = 1, SIZE(array, 1)
126 y = array(i1) - c
127 t = ks + y
128 c = (t - ks) - y
129 ks = t
130 END DO
131 END IF
132 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 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 t = dzero; y = dzero; c = dzero; ks_local = dzero
148
149 n = SIZE(array1)
150 DO i = 1, mod(n, dblksize)
151 y(1) = array1(i)*array2(i) - c(1)
152 t(1) = ks_local(1) + y(1)
153 c(1) = (t(1) - ks_local(1)) - y(1)
154 ks_local(1) = t(1)
155 END DO
156 DO i = mod(n, dblksize) + 1, n, dblksize
157 y = array1(i:i + (dblksize - 1))*array2(i:i + (dblksize - 1)) - c
158 t = ks_local + y
159 c = (t - ks_local) - y
160 ks_local = t
161 END DO
162 DO i = 2, dblksize
163 y(1) = ks_local(i) - (c(1) + c(i))
164 t(1) = ks_local(1) + y(1)
165 c(1) = (t(1) - ks_local(1)) - y(1)
166 ks_local(1) = t(1)
167 END DO
168 ks = ks_local(1)
169 END FUNCTION kahan_dot_product_d1
170
171! **************************************************************************************************
172!> \brief ...
173!> \param array ...
174!> \param mask ...
175!> \return ...
176! **************************************************************************************************
177 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 ks = czero; t = czero; y = czero; c = czero
186
187 IF (PRESENT(mask)) THEN
188 DO i1 = 1, SIZE(array, 1)
189 IF (mask(i1)) THEN
190 y = array(i1) - c
191 t = ks + y
192 c = (t - ks) - y
193 ks = t
194 END IF
195 END DO
196 ELSE
197 DO i1 = 1, SIZE(array, 1)
198 y = array(i1) - c
199 t = ks + y
200 c = (t - ks) - y
201 ks = t
202 END DO
203 END IF
204 END FUNCTION kahan_sum_c1
205
206! **************************************************************************************************
207!> \brief ...
208!> \param array ...
209!> \param mask ...
210!> \return ...
211! **************************************************************************************************
212 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 ks = zzero; t = zzero; y = zzero; c = zzero
221
222 IF (PRESENT(mask)) THEN
223 DO i1 = 1, SIZE(array, 1)
224 IF (mask(i1)) THEN
225 y = array(i1) - c
226 t = ks + y
227 c = (t - ks) - y
228 ks = t
229 END IF
230 END DO
231 ELSE
232 DO i1 = 1, SIZE(array, 1)
233 y = array(i1) - c
234 t = ks + y
235 c = (t - ks) - y
236 ks = t
237 END DO
238 END IF
239 END FUNCTION kahan_sum_z1
240
241! **************************************************************************************************
242!> \brief ...
243!> \param array ...
244!> \param mask ...
245!> \return ...
246! **************************************************************************************************
247 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 ks = szero; t = szero; y = szero; c = szero
256
257 IF (PRESENT(mask)) THEN
258 DO i2 = 1, SIZE(array, 2)
259 DO i1 = 1, SIZE(array, 1)
260 IF (mask(i1, i2)) THEN
261 y = array(i1, i2) - c
262 t = ks + y
263 c = (t - ks) - y
264 ks = t
265 END IF
266 END DO
267 END DO
268 ELSE
269 DO i2 = 1, SIZE(array, 2)
270 DO i1 = 1, SIZE(array, 1)
271 y = array(i1, i2) - c
272 t = ks + y
273 c = (t - ks) - y
274 ks = t
275 END DO
276 END DO
277 END IF
278 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 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 ks = dzero; t = dzero; y = dzero; c = dzero
294
295 n1 = SIZE(array1, 1)
296 n2 = SIZE(array1, 2)
297 DO i2 = 1, n2
298 DO i1 = 1, n1
299 y = real(array1(i1, i2), dp)*real(array2(i1, i2), dp) - c
300 t = ks + y
301 c = (t - ks) - y
302 ks = t
303 END DO
304 END DO
305 END FUNCTION kahan_dot_product_s2
306
307! **************************************************************************************************
308!> \brief ...
309!> \param array ...
310!> \param mask ...
311!> \return ...
312! **************************************************************************************************
313 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 ks = dzero; t = dzero; y = dzero; c = dzero
322
323 IF (PRESENT(mask)) THEN
324 DO i2 = 1, SIZE(array, 2)
325 DO i1 = 1, SIZE(array, 1)
326 IF (mask(i1, i2)) THEN
327 y = array(i1, i2) - c
328 t = ks + y
329 c = (t - ks) - y
330 ks = t
331 END IF
332 END DO
333 END DO
334 ELSE
335 DO i2 = 1, SIZE(array, 2)
336 DO i1 = 1, SIZE(array, 1)
337 y = array(i1, i2) - c
338 t = ks + y
339 c = (t - ks) - y
340 ks = t
341 END DO
342 END DO
343 END IF
344 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 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 ks = dzero; t = dzero; y = dzero; c = dzero
360
361 n1 = SIZE(array1, 1)
362 n2 = SIZE(array1, 2)
363 DO i2 = 1, n2
364 DO i1 = 1, n1
365 y = array1(i1, i2)*array2(i1, i2) - c
366 t = ks + y
367 c = (t - ks) - y
368 ks = t
369 END DO
370 END DO
371 END FUNCTION kahan_dot_product_d2
372
373! **************************************************************************************************
374!> \brief ...
375!> \param array ...
376!> \param mask ...
377!> \return ...
378! **************************************************************************************************
379 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 ks = czero; t = czero; y = czero; c = czero
388
389 IF (PRESENT(mask)) THEN
390 DO i2 = 1, SIZE(array, 2)
391 DO i1 = 1, SIZE(array, 1)
392 IF (mask(i1, i2)) THEN
393 y = array(i1, i2) - c
394 t = ks + y
395 c = (t - ks) - y
396 ks = t
397 END IF
398 END DO
399 END DO
400 ELSE
401 DO i2 = 1, SIZE(array, 2)
402 DO i1 = 1, SIZE(array, 1)
403 y = array(i1, i2) - c
404 t = ks + y
405 c = (t - ks) - y
406 ks = t
407 END DO
408 END DO
409 END IF
410 END FUNCTION kahan_sum_c2
411
412! **************************************************************************************************
413!> \brief ...
414!> \param array ...
415!> \param mask ...
416!> \return ...
417! **************************************************************************************************
418 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 ks = zzero; t = zzero; y = zzero; c = zzero
427
428 IF (PRESENT(mask)) THEN
429 DO i2 = 1, SIZE(array, 2)
430 DO i1 = 1, SIZE(array, 1)
431 IF (mask(i1, i2)) THEN
432 y = array(i1, i2) - c
433 t = ks + y
434 c = (t - ks) - y
435 ks = t
436 END IF
437 END DO
438 END DO
439 ELSE
440 DO i2 = 1, SIZE(array, 2)
441 DO i1 = 1, SIZE(array, 1)
442 y = array(i1, i2) - c
443 t = ks + y
444 c = (t - ks) - y
445 ks = t
446 END DO
447 END DO
448 END IF
449 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 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 ks = zzero; t = zzero; y = zzero; c = zzero
465
466 n1 = SIZE(array1, 1)
467 n2 = SIZE(array1, 2)
468 DO i2 = 1, n2
469 DO i1 = 1, n1
470 y = array1(i1, i2)*array2(i1, i2) - c
471 t = ks + y
472 c = (t - ks) - y
473 ks = t
474 END DO
475 END DO
476 END FUNCTION kahan_dot_product_z2
477
478! **************************************************************************************************
479!> \brief ...
480!> \param array ...
481!> \param mask ...
482!> \return ...
483! **************************************************************************************************
484 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 ks = szero; t = szero; y = szero; c = szero
493
494 IF (PRESENT(mask)) THEN
495 DO i3 = 1, SIZE(array, 3)
496 DO i2 = 1, SIZE(array, 2)
497 DO i1 = 1, SIZE(array, 1)
498 IF (mask(i1, i2, i3)) THEN
499 y = array(i1, i2, i3) - c
500 t = ks + y
501 c = (t - ks) - y
502 ks = t
503 END IF
504 END DO
505 END DO
506 END DO
507 ELSE
508 DO i3 = 1, SIZE(array, 3)
509 DO i2 = 1, SIZE(array, 2)
510 DO i1 = 1, SIZE(array, 1)
511 y = array(i1, i2, i3) - c
512 t = ks + y
513 c = (t - ks) - y
514 ks = t
515 END DO
516 END DO
517 END DO
518 END IF
519 END FUNCTION kahan_sum_s3
520
521! **************************************************************************************************
522!> \brief ...
523!> \param array ...
524!> \param mask ...
525!> \return ...
526! **************************************************************************************************
527 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 ks = dzero; t = dzero; y = dzero; c = dzero
536
537 IF (PRESENT(mask)) THEN
538 DO i3 = 1, SIZE(array, 3)
539 DO i2 = 1, SIZE(array, 2)
540 DO i1 = 1, SIZE(array, 1)
541 IF (mask(i1, i2, i3)) THEN
542 y = array(i1, i2, i3) - c
543 t = ks + y
544 c = (t - ks) - y
545 ks = t
546 END IF
547 END DO
548 END DO
549 END DO
550 ELSE
551 DO i3 = 1, SIZE(array, 3)
552 DO i2 = 1, SIZE(array, 2)
553 DO i1 = 1, SIZE(array, 1)
554 y = array(i1, i2, i3) - c
555 t = ks + y
556 c = (t - ks) - y
557 ks = t
558 END DO
559 END DO
560 END DO
561 END IF
562 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 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 ks = dzero; t = dzero; y = dzero; c = dzero
578
579 n1 = SIZE(array1, 1)
580 n2 = SIZE(array1, 2)
581 n3 = SIZE(array1, 3)
582 DO i3 = 1, n3
583 DO i2 = 1, n2
584 DO i1 = 1, n1
585 y = array1(i1, i2, i3)*array2(i1, i2, i3) - c
586 t = ks + y
587 c = (t - ks) - y
588 ks = t
589 END DO
590 END DO
591 END DO
592 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 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 ks = dzero; t = dzero; y = dzero; c = dzero
614 DO i3 = lbound(mask, 3), ubound(mask, 3)
615 DO i2 = lbound(mask, 2), ubound(mask, 2)
616 DO i1 = lbound(mask, 1), ubound(mask, 1)
617 IF (mask(i1, i2, i3) .GT. th) THEN
618 y = array1(i1, i2, i3)*array2(i1, i2, i3) - c
619 t = ks + y
620 c = (t - ks) - y
621 ks = t
622 END IF
623 END DO
624 END DO
625 END DO
626 END FUNCTION kahan_dot_product_masked_d3
627
628! **************************************************************************************************
629!> \brief ...
630!> \param array ...
631!> \param mask ...
632!> \return ...
633! **************************************************************************************************
634 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 ks = czero; t = czero; y = czero; c = czero
643
644 IF (PRESENT(mask)) THEN
645 DO i3 = 1, SIZE(array, 3)
646 DO i2 = 1, SIZE(array, 2)
647 DO i1 = 1, SIZE(array, 1)
648 IF (mask(i1, i2, i3)) THEN
649 y = array(i1, i2, i3) - c
650 t = ks + y
651 c = (t - ks) - y
652 ks = t
653 END IF
654 END DO
655 END DO
656 END DO
657 ELSE
658 DO i3 = 1, SIZE(array, 3)
659 DO i2 = 1, SIZE(array, 2)
660 DO i1 = 1, SIZE(array, 1)
661 y = array(i1, i2, i3) - c
662 t = ks + y
663 c = (t - ks) - y
664 ks = t
665 END DO
666 END DO
667 END DO
668 END IF
669 END FUNCTION kahan_sum_c3
670
671! **************************************************************************************************
672!> \brief ...
673!> \param array ...
674!> \param mask ...
675!> \return ...
676! **************************************************************************************************
677 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 ks = zzero; t = zzero; y = zzero; c = zzero
686
687 IF (PRESENT(mask)) THEN
688 DO i3 = 1, SIZE(array, 3)
689 DO i2 = 1, SIZE(array, 2)
690 DO i1 = 1, SIZE(array, 1)
691 IF (mask(i1, i2, i3)) THEN
692 y = array(i1, i2, i3) - c
693 t = ks + y
694 c = (t - ks) - y
695 ks = t
696 END IF
697 END DO
698 END DO
699 END DO
700 ELSE
701 DO i3 = 1, SIZE(array, 3)
702 DO i2 = 1, SIZE(array, 2)
703 DO i1 = 1, SIZE(array, 1)
704 y = array(i1, i2, i3) - c
705 t = ks + y
706 c = (t - ks) - y
707 ks = t
708 END DO
709 END DO
710 END DO
711 END IF
712 END FUNCTION kahan_sum_z3
713
714! **************************************************************************************************
715!> \brief ...
716!> \param array ...
717!> \param mask ...
718!> \return ...
719! **************************************************************************************************
720 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 ks = szero; t = szero; y = szero; c = szero
730
731 IF (PRESENT(mask)) THEN
732 DO i4 = 1, SIZE(array, 4)
733 DO i3 = 1, SIZE(array, 3)
734 DO i2 = 1, SIZE(array, 2)
735 DO i1 = 1, SIZE(array, 1)
736 IF (mask(i1, i2, i3, i4)) THEN
737 y = array(i1, i2, i3, i4) - c
738 t = ks + y
739 c = (t - ks) - y
740 ks = t
741 END IF
742 END DO
743 END DO
744 END DO
745 END DO
746 ELSE
747 DO i4 = 1, SIZE(array, 4)
748 DO i3 = 1, SIZE(array, 3)
749 DO i2 = 1, SIZE(array, 2)
750 DO i1 = 1, SIZE(array, 1)
751 y = array(i1, i2, i3, i4) - c
752 t = ks + y
753 c = (t - ks) - y
754 ks = t
755 END DO
756 END DO
757 END DO
758 END DO
759 END IF
760 END FUNCTION kahan_sum_s4
761
762! **************************************************************************************************
763!> \brief ...
764!> \param array ...
765!> \param mask ...
766!> \return ...
767! **************************************************************************************************
768 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 ks = dzero; t = dzero; y = dzero; c = dzero
778
779 IF (PRESENT(mask)) THEN
780 DO i4 = 1, SIZE(array, 4)
781 DO i3 = 1, SIZE(array, 3)
782 DO i2 = 1, SIZE(array, 2)
783 DO i1 = 1, SIZE(array, 1)
784 IF (mask(i1, i2, i3, i4)) THEN
785 y = array(i1, i2, i3, i4) - c
786 t = ks + y
787 c = (t - ks) - y
788 ks = t
789 END IF
790 END DO
791 END DO
792 END DO
793 END DO
794 ELSE
795 DO i4 = 1, SIZE(array, 4)
796 DO i3 = 1, SIZE(array, 3)
797 DO i2 = 1, SIZE(array, 2)
798 DO i1 = 1, SIZE(array, 1)
799 y = array(i1, i2, i3, i4) - c
800 t = ks + y
801 c = (t - ks) - y
802 ks = t
803 END DO
804 END DO
805 END DO
806 END DO
807 END IF
808 END FUNCTION kahan_sum_d4
809
810! **************************************************************************************************
811!> \brief ...
812!> \param array ...
813!> \param mask ...
814!> \return ...
815! **************************************************************************************************
816 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 ks = czero; t = czero; y = czero; c = czero
827
828 IF (PRESENT(mask)) THEN
829 DO i4 = 1, SIZE(array, 4)
830 DO i3 = 1, SIZE(array, 3)
831 DO i2 = 1, SIZE(array, 2)
832 DO i1 = 1, SIZE(array, 1)
833 IF (mask(i1, i2, i3, i4)) THEN
834 y = array(i1, i2, i3, i4) - c
835 t = ks + y
836 c = (t - ks) - y
837 ks = t
838 END IF
839 END DO
840 END DO
841 END DO
842 END DO
843 ELSE
844 DO i4 = 1, SIZE(array, 4)
845 DO i3 = 1, SIZE(array, 3)
846 DO i2 = 1, SIZE(array, 2)
847 DO i1 = 1, SIZE(array, 1)
848 y = array(i1, i2, i3, i4) - c
849 t = ks + y
850 c = (t - ks) - y
851 ks = t
852 END DO
853 END DO
854 END DO
855 END DO
856 END IF
857 END FUNCTION kahan_sum_c4
858
859! **************************************************************************************************
860!> \brief ...
861!> \param array ...
862!> \param mask ...
863!> \return ...
864! **************************************************************************************************
865 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 ks = zzero; t = zzero; y = zzero; c = zzero
876
877 IF (PRESENT(mask)) THEN
878 DO i4 = 1, SIZE(array, 4)
879 DO i3 = 1, SIZE(array, 3)
880 DO i2 = 1, SIZE(array, 2)
881 DO i1 = 1, SIZE(array, 1)
882 IF (mask(i1, i2, i3, i4)) THEN
883 y = array(i1, i2, i3, i4) - c
884 t = ks + y
885 c = (t - ks) - y
886 ks = t
887 END IF
888 END DO
889 END DO
890 END DO
891 END DO
892 ELSE
893 DO i4 = 1, SIZE(array, 4)
894 DO i3 = 1, SIZE(array, 3)
895 DO i2 = 1, SIZE(array, 2)
896 DO i1 = 1, SIZE(array, 1)
897 y = array(i1, i2, i3, i4) - c
898 t = ks + y
899 c = (t - ks) - y
900 ks = t
901 END DO
902 END DO
903 END DO
904 END DO
905 END IF
906 END FUNCTION kahan_sum_z4
907
908! **************************************************************************************************
909!> \brief ...
910!> \param array ...
911!> \param mask ...
912!> \return ...
913! **************************************************************************************************
914 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 ks = szero; t = szero; y = szero; c = szero
925
926 IF (PRESENT(mask)) THEN
927 DO i5 = 1, SIZE(array, 5)
928 DO i4 = 1, SIZE(array, 4)
929 DO i3 = 1, SIZE(array, 3)
930 DO i2 = 1, SIZE(array, 2)
931 DO i1 = 1, SIZE(array, 1)
932 IF (mask(i1, i2, i3, i4, i5)) THEN
933 y = array(i1, i2, i3, i4, i5) - c
934 t = ks + y
935 c = (t - ks) - y
936 ks = t
937 END IF
938 END DO
939 END DO
940 END DO
941 END DO
942 END DO
943 ELSE
944 DO i5 = 1, SIZE(array, 5)
945 DO i4 = 1, SIZE(array, 4)
946 DO i3 = 1, SIZE(array, 3)
947 DO i2 = 1, SIZE(array, 2)
948 DO i1 = 1, SIZE(array, 1)
949 y = array(i1, i2, i3, i4, i5) - c
950 t = ks + y
951 c = (t - ks) - y
952 ks = t
953 END DO
954 END DO
955 END DO
956 END DO
957 END DO
958 END IF
959 END FUNCTION kahan_sum_s5
960
961! **************************************************************************************************
962!> \brief ...
963!> \param array ...
964!> \param mask ...
965!> \return ...
966! **************************************************************************************************
967 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 ks = dzero; t = dzero; y = dzero; c = dzero
978
979 IF (PRESENT(mask)) THEN
980 DO i5 = 1, SIZE(array, 5)
981 DO i4 = 1, SIZE(array, 4)
982 DO i3 = 1, SIZE(array, 3)
983 DO i2 = 1, SIZE(array, 2)
984 DO i1 = 1, SIZE(array, 1)
985 IF (mask(i1, i2, i3, i4, i5)) THEN
986 y = array(i1, i2, i3, i4, i5) - c
987 t = ks + y
988 c = (t - ks) - y
989 ks = t
990 END IF
991 END DO
992 END DO
993 END DO
994 END DO
995 END DO
996 ELSE
997 DO i5 = 1, SIZE(array, 5)
998 DO i4 = 1, SIZE(array, 4)
999 DO i3 = 1, SIZE(array, 3)
1000 DO i2 = 1, SIZE(array, 2)
1001 DO i1 = 1, SIZE(array, 1)
1002 y = array(i1, i2, i3, i4, i5) - c
1003 t = ks + y
1004 c = (t - ks) - y
1005 ks = t
1006 END DO
1007 END DO
1008 END DO
1009 END DO
1010 END DO
1011 END IF
1012 END FUNCTION kahan_sum_d5
1013
1014! **************************************************************************************************
1015!> \brief ...
1016!> \param array ...
1017!> \param mask ...
1018!> \return ...
1019! **************************************************************************************************
1020 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 ks = czero; t = czero; y = czero; c = czero
1031
1032 IF (PRESENT(mask)) THEN
1033 DO i5 = 1, SIZE(array, 5)
1034 DO i4 = 1, SIZE(array, 4)
1035 DO i3 = 1, SIZE(array, 3)
1036 DO i2 = 1, SIZE(array, 2)
1037 DO i1 = 1, SIZE(array, 1)
1038 IF (mask(i1, i2, i3, i4, i5)) THEN
1039 y = array(i1, i2, i3, i4, i5) - c
1040 t = ks + y
1041 c = (t - ks) - y
1042 ks = t
1043 END IF
1044 END DO
1045 END DO
1046 END DO
1047 END DO
1048 END DO
1049 ELSE
1050 DO i5 = 1, SIZE(array, 5)
1051 DO i4 = 1, SIZE(array, 4)
1052 DO i3 = 1, SIZE(array, 3)
1053 DO i2 = 1, SIZE(array, 2)
1054 DO i1 = 1, SIZE(array, 1)
1055 y = array(i1, i2, i3, i4, i5) - c
1056 t = ks + y
1057 c = (t - ks) - y
1058 ks = t
1059 END DO
1060 END DO
1061 END DO
1062 END DO
1063 END DO
1064 END IF
1065 END FUNCTION kahan_sum_c5
1066
1067! **************************************************************************************************
1068!> \brief ...
1069!> \param array ...
1070!> \param mask ...
1071!> \return ...
1072! **************************************************************************************************
1073 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 ks = zzero; t = zzero; y = zzero; c = zzero
1084
1085 IF (PRESENT(mask)) THEN
1086 DO i5 = 1, SIZE(array, 5)
1087 DO i4 = 1, SIZE(array, 4)
1088 DO i3 = 1, SIZE(array, 3)
1089 DO i2 = 1, SIZE(array, 2)
1090 DO i1 = 1, SIZE(array, 1)
1091 IF (mask(i1, i2, i3, i4, i5)) THEN
1092 y = array(i1, i2, i3, i4, i5) - c
1093 t = ks + y
1094 c = (t - ks) - y
1095 ks = t
1096 END IF
1097 END DO
1098 END DO
1099 END DO
1100 END DO
1101 END DO
1102 ELSE
1103 DO i5 = 1, SIZE(array, 5)
1104 DO i4 = 1, SIZE(array, 4)
1105 DO i3 = 1, SIZE(array, 3)
1106 DO i2 = 1, SIZE(array, 2)
1107 DO i1 = 1, SIZE(array, 1)
1108 y = array(i1, i2, i3, i4, i5) - c
1109 t = ks + y
1110 c = (t - ks) - y
1111 ks = t
1112 END DO
1113 END DO
1114 END DO
1115 END DO
1116 END DO
1117 END IF
1118 END FUNCTION kahan_sum_z5
1119
1120! **************************************************************************************************
1121!> \brief ...
1122!> \param array ...
1123!> \param mask ...
1124!> \return ...
1125! **************************************************************************************************
1126 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 ks = szero; t = szero; y = szero; c = szero
1137
1138 IF (PRESENT(mask)) THEN
1139 DO i6 = 1, SIZE(array, 6)
1140 DO i5 = 1, SIZE(array, 5)
1141 DO i4 = 1, SIZE(array, 4)
1142 DO i3 = 1, SIZE(array, 3)
1143 DO i2 = 1, SIZE(array, 2)
1144 DO i1 = 1, SIZE(array, 1)
1145 IF (mask(i1, i2, i3, i4, i5, i6)) THEN
1146 y = array(i1, i2, i3, i4, i5, i6) - c
1147 t = ks + y
1148 c = (t - ks) - y
1149 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 DO i6 = 1, SIZE(array, 6)
1159 DO i5 = 1, SIZE(array, 5)
1160 DO i4 = 1, SIZE(array, 4)
1161 DO i3 = 1, SIZE(array, 3)
1162 DO i2 = 1, SIZE(array, 2)
1163 DO i1 = 1, SIZE(array, 1)
1164 y = array(i1, i2, i3, i4, i5, i6) - c
1165 t = ks + y
1166 c = (t - ks) - y
1167 ks = t
1168 END DO
1169 END DO
1170 END DO
1171 END DO
1172 END DO
1173 END DO
1174 END IF
1175 END FUNCTION kahan_sum_s6
1176
1177! **************************************************************************************************
1178!> \brief ...
1179!> \param array ...
1180!> \param mask ...
1181!> \return ...
1182! **************************************************************************************************
1183 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 ks = dzero; t = dzero; y = dzero; c = dzero
1194
1195 IF (PRESENT(mask)) THEN
1196 DO i6 = 1, SIZE(array, 6)
1197 DO i5 = 1, SIZE(array, 5)
1198 DO i4 = 1, SIZE(array, 4)
1199 DO i3 = 1, SIZE(array, 3)
1200 DO i2 = 1, SIZE(array, 2)
1201 DO i1 = 1, SIZE(array, 1)
1202 IF (mask(i1, i2, i3, i4, i5, i6)) THEN
1203 y = array(i1, i2, i3, i4, i5, i6) - c
1204 t = ks + y
1205 c = (t - ks) - y
1206 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 DO i6 = 1, SIZE(array, 6)
1216 DO i5 = 1, SIZE(array, 5)
1217 DO i4 = 1, SIZE(array, 4)
1218 DO i3 = 1, SIZE(array, 3)
1219 DO i2 = 1, SIZE(array, 2)
1220 DO i1 = 1, SIZE(array, 1)
1221 y = array(i1, i2, i3, i4, i5, i6) - c
1222 t = ks + y
1223 c = (t - ks) - y
1224 ks = t
1225 END DO
1226 END DO
1227 END DO
1228 END DO
1229 END DO
1230 END DO
1231 END IF
1232 END FUNCTION kahan_sum_d6
1233
1234! **************************************************************************************************
1235!> \brief ...
1236!> \param array ...
1237!> \param mask ...
1238!> \return ...
1239! **************************************************************************************************
1240 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 ks = czero; t = czero; y = czero; c = czero
1251
1252 IF (PRESENT(mask)) THEN
1253 DO i6 = 1, SIZE(array, 6)
1254 DO i5 = 1, SIZE(array, 5)
1255 DO i4 = 1, SIZE(array, 4)
1256 DO i3 = 1, SIZE(array, 3)
1257 DO i2 = 1, SIZE(array, 2)
1258 DO i1 = 1, SIZE(array, 1)
1259 IF (mask(i1, i2, i3, i4, i5, i6)) THEN
1260 y = array(i1, i2, i3, i4, i5, i6) - c
1261 t = ks + y
1262 c = (t - ks) - y
1263 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 DO i6 = 1, SIZE(array, 6)
1273 DO i5 = 1, SIZE(array, 5)
1274 DO i4 = 1, SIZE(array, 4)
1275 DO i3 = 1, SIZE(array, 3)
1276 DO i2 = 1, SIZE(array, 2)
1277 DO i1 = 1, SIZE(array, 1)
1278 y = array(i1, i2, i3, i4, i5, i6) - c
1279 t = ks + y
1280 c = (t - ks) - y
1281 ks = t
1282 END DO
1283 END DO
1284 END DO
1285 END DO
1286 END DO
1287 END DO
1288 END IF
1289 END FUNCTION kahan_sum_c6
1290
1291! **************************************************************************************************
1292!> \brief ...
1293!> \param array ...
1294!> \param mask ...
1295!> \return ...
1296! **************************************************************************************************
1297 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 ks = zzero; t = zzero; y = zzero; c = zzero
1308
1309 IF (PRESENT(mask)) THEN
1310 DO i6 = 1, SIZE(array, 6)
1311 DO i5 = 1, SIZE(array, 5)
1312 DO i4 = 1, SIZE(array, 4)
1313 DO i3 = 1, SIZE(array, 3)
1314 DO i2 = 1, SIZE(array, 2)
1315 DO i1 = 1, SIZE(array, 1)
1316 IF (mask(i1, i2, i3, i4, i5, i6)) THEN
1317 y = array(i1, i2, i3, i4, i5, i6) - c
1318 t = ks + y
1319 c = (t - ks) - y
1320 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 DO i6 = 1, SIZE(array, 6)
1330 DO i5 = 1, SIZE(array, 5)
1331 DO i4 = 1, SIZE(array, 4)
1332 DO i3 = 1, SIZE(array, 3)
1333 DO i2 = 1, SIZE(array, 2)
1334 DO i1 = 1, SIZE(array, 1)
1335 y = array(i1, i2, i3, i4, i5, i6) - c
1336 t = ks + y
1337 c = (t - ks) - y
1338 ks = t
1339 END DO
1340 END DO
1341 END DO
1342 END DO
1343 END DO
1344 END DO
1345 END IF
1346 END FUNCTION kahan_sum_z6
1347
1348! **************************************************************************************************
1349!> \brief ...
1350!> \param array ...
1351!> \param mask ...
1352!> \return ...
1353! **************************************************************************************************
1354 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 ks = szero; t = szero; y = szero; c = szero
1365
1366 IF (PRESENT(mask)) THEN
1367 DO i7 = 1, SIZE(array, 7)
1368 DO i6 = 1, SIZE(array, 6)
1369 DO i5 = 1, SIZE(array, 5)
1370 DO i4 = 1, SIZE(array, 4)
1371 DO i3 = 1, SIZE(array, 3)
1372 DO i2 = 1, SIZE(array, 2)
1373 DO i1 = 1, SIZE(array, 1)
1374 IF (mask(i1, i2, i3, i4, i5, i6, i7)) THEN
1375 y = array(i1, i2, i3, i4, i5, i6, i7) - c
1376 t = ks + y
1377 c = (t - ks) - y
1378 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 DO i7 = 1, SIZE(array, 7)
1389 DO i6 = 1, SIZE(array, 6)
1390 DO i5 = 1, SIZE(array, 5)
1391 DO i4 = 1, SIZE(array, 4)
1392 DO i3 = 1, SIZE(array, 3)
1393 DO i2 = 1, SIZE(array, 2)
1394 DO i1 = 1, SIZE(array, 1)
1395 y = array(i1, i2, i3, i4, i5, i6, i7) - c
1396 t = ks + y
1397 c = (t - ks) - y
1398 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 END FUNCTION kahan_sum_s7
1408
1409! **************************************************************************************************
1410!> \brief ...
1411!> \param array ...
1412!> \param mask ...
1413!> \return ...
1414! **************************************************************************************************
1415 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 ks = dzero; t = dzero; y = dzero; c = dzero
1426
1427 IF (PRESENT(mask)) THEN
1428 DO i7 = 1, SIZE(array, 7)
1429 DO i6 = 1, SIZE(array, 6)
1430 DO i5 = 1, SIZE(array, 5)
1431 DO i4 = 1, SIZE(array, 4)
1432 DO i3 = 1, SIZE(array, 3)
1433 DO i2 = 1, SIZE(array, 2)
1434 DO i1 = 1, SIZE(array, 1)
1435 IF (mask(i1, i2, i3, i4, i5, i6, i7)) THEN
1436 y = array(i1, i2, i3, i4, i5, i6, i7) - c
1437 t = ks + y
1438 c = (t - ks) - y
1439 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 DO i7 = 1, SIZE(array, 7)
1450 DO i6 = 1, SIZE(array, 6)
1451 DO i5 = 1, SIZE(array, 5)
1452 DO i4 = 1, SIZE(array, 4)
1453 DO i3 = 1, SIZE(array, 3)
1454 DO i2 = 1, SIZE(array, 2)
1455 DO i1 = 1, SIZE(array, 1)
1456 y = array(i1, i2, i3, i4, i5, i6, i7) - c
1457 t = ks + y
1458 c = (t - ks) - y
1459 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 END FUNCTION kahan_sum_d7
1469
1470! **************************************************************************************************
1471!> \brief ...
1472!> \param array ...
1473!> \param mask ...
1474!> \return ...
1475! **************************************************************************************************
1476 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 ks = czero; t = czero; y = czero; c = czero
1487
1488 IF (PRESENT(mask)) THEN
1489 DO i7 = 1, SIZE(array, 7)
1490 DO i6 = 1, SIZE(array, 6)
1491 DO i5 = 1, SIZE(array, 5)
1492 DO i4 = 1, SIZE(array, 4)
1493 DO i3 = 1, SIZE(array, 3)
1494 DO i2 = 1, SIZE(array, 2)
1495 DO i1 = 1, SIZE(array, 1)
1496 IF (mask(i1, i2, i3, i4, i5, i6, i7)) THEN
1497 y = array(i1, i2, i3, i4, i5, i6, i7) - c
1498 t = ks + y
1499 c = (t - ks) - y
1500 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 DO i7 = 1, SIZE(array, 7)
1511 DO i6 = 1, SIZE(array, 6)
1512 DO i5 = 1, SIZE(array, 5)
1513 DO i4 = 1, SIZE(array, 4)
1514 DO i3 = 1, SIZE(array, 3)
1515 DO i2 = 1, SIZE(array, 2)
1516 DO i1 = 1, SIZE(array, 1)
1517 y = array(i1, i2, i3, i4, i5, i6, i7) - c
1518 t = ks + y
1519 c = (t - ks) - y
1520 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 END FUNCTION kahan_sum_c7
1530
1531! **************************************************************************************************
1532!> \brief ...
1533!> \param array ...
1534!> \param mask ...
1535!> \return ...
1536! **************************************************************************************************
1537 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 ks = zzero; t = zzero; y = zzero; c = zzero
1548
1549 IF (PRESENT(mask)) THEN
1550 DO i7 = 1, SIZE(array, 7)
1551 DO i6 = 1, SIZE(array, 6)
1552 DO i5 = 1, SIZE(array, 5)
1553 DO i4 = 1, SIZE(array, 4)
1554 DO i3 = 1, SIZE(array, 3)
1555 DO i2 = 1, SIZE(array, 2)
1556 DO i1 = 1, SIZE(array, 1)
1557 IF (mask(i1, i2, i3, i4, i5, i6, i7)) THEN
1558 y = array(i1, i2, i3, i4, i5, i6, i7) - c
1559 t = ks + y
1560 c = (t - ks) - y
1561 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 DO i7 = 1, SIZE(array, 7)
1572 DO i6 = 1, SIZE(array, 6)
1573 DO i5 = 1, SIZE(array, 5)
1574 DO i4 = 1, SIZE(array, 4)
1575 DO i3 = 1, SIZE(array, 3)
1576 DO i2 = 1, SIZE(array, 2)
1577 DO i1 = 1, SIZE(array, 1)
1578 y = array(i1, i2, i3, i4, i5, i6, i7) - c
1579 t = ks + y
1580 c = (t - ks) - y
1581 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 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 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 my_blksize = 32
1608 IF (PRESENT(blksize)) my_blksize = blksize
1609
1610 IF (my_blksize <= 1) THEN
1611 ! The original should be faster
1612 ks = accurate_dot_product(array1, array2)
1613 ELSE IF (my_blksize >= SIZE(array1)) THEN
1614 ! Just use standard dot product from BLAS for performance
1615 ks = ddot(SIZE(array1), array1(1), 1, array2(1), 1)
1616 ELSE
1617 ks = 0.0_dp
1618 block
1619 INTEGER :: i, n, stripesize
1620 REAL(kind=dp) :: c, dotproduct, t, y
1621 t = dzero; y = dzero; c = dzero
1622
1623 n = SIZE(array1)
1624 DO i = 1, n, my_blksize
1625 ! Remove 1 to save an operation in the length
1626 stripesize = min(my_blksize, n - i + 1)
1627 ! Perform dot product on the given stripe
1628 dotproduct = ddot(stripesize, array1(i), 1, array2(i), 1)
1629 y = dotproduct - c
1630 t = ks + y
1631 c = (t - ks) - y
1632 ks = t
1633 END DO
1634 END block
1635 END IF
1636 END FUNCTION kahan_blocked_dot_product_d1
1637
1638END MODULE kahan_sum
kahan_sum::accurate_dot_product_2
Definition kahan_sum.F:59
kahan_sum::accurate_dot_product
Definition kahan_sum.F:52
kahan_sum::accurate_sum
Definition kahan_sum.F:41
kahan_sum
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Definition kahan_sum.F:29