30 use mpp_mod,
only :
mpp_error, fatal, warning, note
35 implicit none ;
private
38 public :: mpp_efp_plus, mpp_efp_minus, mpp_efp_to_real, mpp_real_to_efp, mpp_efp_real_diff
39 public ::
operator(+),
operator(-),
assignment(=)
40 public :: mpp_query_efp_overflow_error, mpp_reset_efp_overflow_error
42 integer,
parameter :: NUMBIT = 46
46 integer(i8_kind),
parameter ::
prec=2_8**numbit
47 real(r8_kind),
parameter ::
r_prec=2.0_8**numbit
48 real(r8_kind),
parameter ::
i_prec=1.0_8/(2.0_8**numbit)
53 real(r8_kind),
parameter,
dimension(NUMINT) :: &
55 real(r8_kind),
parameter,
dimension(NUMINT) :: &
58 logical :: overflow_error = .false., nan_error = .false.
81 integer(i8_kind),
dimension(NUMINT) :: v
87 interface operator (+); module
procedure mpp_efp_plus ; end interface
90 interface operator (-); module
procedure mpp_efp_minus ; end interface
93 interface assignment(=); module
procedure mpp_efp_assign ; end interface
102 overflow_check, err)
result(sum)
103 real(r8_kind),
dimension(:,:),
intent(in) :: array
104 integer,
optional,
intent(in) :: isr, ier, jsr, jer
106 logical,
optional,
intent(in) :: reproducing
107 logical,
optional,
intent(in) :: overflow_check
108 integer,
optional,
intent(out) :: err
111 integer(i8_kind),
dimension(NUMINT) :: ints_sum
112 integer(i8_kind) :: ival, prec_error
113 real(r8_kind) :: rsum(1), rs
114 real(r8_kind) :: max_mag_term
115 logical :: repro, over_check
116 character(len=256) :: mesg
117 integer :: i, j, n, is, ie, js, je, sgn
120 "mpp_reproducing_sum: Too many processors are being used for the value of "//&
121 "prec. Reduce prec to (2^63-1)/mpp_npes.")
123 prec_error = (2_8**62 + (2_8**62 - 1)) /
mpp_npes()
125 is = 1 ; ie =
size(array,1) ; js = 1 ; je =
size(array,2 )
126 if (
present(isr))
then
128 "Value of isr too small in mpp_reproducing_sum_2d.")
131 if (
present(ier))
then
133 "Value of ier too large in mpp_reproducing_sum_2d.")
136 if (
present(jsr))
then
138 "Value of jsr too small in mpp_reproducing_sum_2d.")
141 if (
present(jer))
then
143 "Value of jer too large in mpp_reproducing_sum_2d.")
147 repro = .true. ;
if (
present(reproducing)) repro = reproducing
148 over_check = .true. ;
if (
present(overflow_check)) over_check = overflow_check
151 overflow_error = .false. ; nan_error = .false. ; max_mag_term = 0.0
155 do j=js,je ;
do i=is,ie
167 do j=js,je ;
do i=is,ie
173 do j=js,je ;
do i=is,ie
174 sgn = 1 ;
if (array(i,j)<0.0) sgn = -1
177 ival = int(rs*i_pr(n), 8)
179 ints_sum(n) = ints_sum(n) + sgn*ival
185 if (
present(err))
then
187 if (overflow_error) &
191 if (err > 0)
then ;
do n=1,
numint ; ints_sum(n) = 0 ;
enddo ;
endif
194 call mpp_error(fatal,
"NaN in input field of mpp_reproducing_sum(_2d), this indicates numerical instability")
196 if (abs(max_mag_term) >= prec_error*pr(1))
then
197 write(mesg,
'(ES13.5)') max_mag_term
198 call mpp_error(fatal,
"Overflow in mpp_reproducing_sum(_2d) conversion of "//trim(mesg))
200 if (overflow_error)
then
201 call mpp_error(fatal,
"Overflow in mpp_reproducing_sum(_2d).")
211 do j=js,je ;
do i=is,ie
212 rsum(1) = rsum(1) + array(i,j);
217 if (
present(err))
then ; err = 0 ;
endif
219 if (
debug .or.
present(efp_sum))
then
220 overflow_error = .false.
221 ints_sum =
real_to_ints(sum, prec_error, overflow_error)
222 if (overflow_error)
then
223 if (
present(err))
then
226 write(mesg,
'(ES13.5)') sum
227 call mpp_error(fatal,
"Repro_sum_2d: Overflow in real_to_ints conversion of "//trim(mesg))
233 if (
present(efp_sum)) efp_sum%v(:) = ints_sum(:)
236 write(mesg,
'("2d RS: ", ES24.16, 6 Z17.16)') sum, ints_sum(1:
numint)
243 overflow_check, err)
result(sum)
244 real(r4_kind),
dimension(:,:),
intent(in) :: array
245 integer,
optional,
intent(in) :: isr, ier, jsr, jer
247 logical,
optional,
intent(in) :: reproducing
248 logical,
optional,
intent(in) :: overflow_check
249 integer,
optional,
intent(out) :: err
252 real(r8_kind) :: array_r8(size(array,1), size(array,2))
257 overflow_check, err), r4_kind)
270 real(r8_kind),
dimension(:,:,:),
intent(in) :: array
271 integer,
optional,
intent(in) :: isr, ier, jsr, jer
272 real(r8_kind),
dimension(:),
optional,
intent(out) :: sums
274 integer,
optional,
intent(out) :: err
277 real(r8_kind) :: max_mag_term
278 integer(i8_kind),
dimension(NUMINT) :: ints_sum
279 integer(i8_kind),
dimension(NUMINT,size(array,3)) :: ints_sums
280 integer(i8_kind) :: prec_error
281 character(len=256) :: mesg
282 integer :: i, j, k, is, ie, js, je, ke, isz, jsz, n
285 "mpp_reproducing_sum: Too many processors are being used for the value of "//&
286 "prec. Reduce prec to (2^63-1)/mpp_npes.")
288 prec_error = (2_8**62 + (2_8**62 - 1)) /
mpp_npes()
291 is = 1 ; ie =
size(array,1) ; js = 1 ; je =
size(array,2) ; ke =
size(array,3)
292 if (
present(isr))
then
294 "Value of isr too small in mpp_reproducing_sum(_3d).")
297 if (
present(ier))
then
299 "Value of ier too large in mpp_reproducing_sum(_3d).")
302 if (
present(jsr))
then
304 "Value of jsr too small in mpp_reproducing_sum(_3d).")
307 if (
present(jer))
then
309 "Value of jer too large in mpp_reproducing_sum(_3d).")
312 jsz = je+1-js; isz = ie+1-is
314 if (
present(sums))
then
315 if (
size(sums) > ke)
call mpp_error(fatal,
"Sums is smaller than "//&
316 "the vertical extent of array in mpp_reproducing_sum(_3d).")
318 overflow_error = .false. ; nan_error = .false. ; max_mag_term = 0.0
321 do j=js,je ;
do i=is,ie
327 do k=1,ke ;
do j=js,je
334 do k=1,ke ;
do j=js,je ;
do i=is,ie
337 enddo ;
enddo ;
enddo
339 if (
present(err))
then
341 if (abs(max_mag_term) >= prec_error*pr(1)) err = err+1
342 if (overflow_error) err = err+2
343 if (nan_error) err = err+2
344 if (err > 0)
then ;
do k=1,ke ;
do n=1,
numint ; ints_sums(n,k) = 0 ;
enddo ;
enddo ;
endif
347 "NaN in input field of mpp_reproducing_sum(_3d), this indicates numerical instability")
348 if (abs(max_mag_term) >= prec_error*pr(1))
then
349 write(mesg,
'(ES13.5)') max_mag_term
350 call mpp_error(fatal,
"Overflow in mpp_reproducing_sum(_3d) conversion of "//trim(mesg))
352 if (overflow_error)
call mpp_error(fatal,
"Overflow in mpp_reproducing_sum(_3d).")
364 if (
present(efp_sum))
then
366 do k=1,ke ;
call increment_ints(efp_sum%v(:), ints_sums(:,k)) ;
enddo
370 do n=1,
numint ; ints_sum(n) = 0 ;
enddo
371 do k=1,ke ;
do n=1,
numint ; ints_sum(n) = ints_sum(n) + ints_sums(n,k) ;
enddo ;
enddo
372 write(mesg,
'("3D RS: ", ES24.16, 6 Z17.16)') sum, ints_sum(1:
numint)
377 overflow_error = .false. ; nan_error = .false. ; max_mag_term = 0.0
380 do j=js,je ;
do i=is,ie
386 do k=1,ke ;
do j=js,je
393 do k=1,ke ;
do j=js,je ;
do i=is,ie
396 enddo ;
enddo ;
enddo
398 if (
present(err))
then
400 if (abs(max_mag_term) >= prec_error*pr(1)) err = err+1
401 if (overflow_error) err = err+2
402 if (nan_error) err = err+2
403 if (err > 0)
then ;
do n=1,
numint ; ints_sum(n) = 0 ;
enddo ;
endif
406 "NaN in input field of mpp_reproducing_sum(_3d), this indicates numerical instability")
407 if (abs(max_mag_term) >= prec_error*pr(1))
then
408 write(mesg,
'(ES13.5)') max_mag_term
409 call mpp_error(fatal,
"Overflow in mpp_reproducing_sum(_3d) conversion of "//trim(mesg))
411 if (overflow_error)
call mpp_error(fatal,
"Overflow in mpp_reproducing_sum(_3d).")
419 if (
present(efp_sum)) efp_sum%v(:) = ints_sum(:)
422 write(mesg,
'("3d RS: ", ES24.16, 6 Z17.16)') sum, ints_sum(1:
numint)
432 real(r8_kind),
intent(in) :: r
433 integer(i8_kind),
optional,
intent(in) :: prec_error
434 logical,
optional,
intent(inout) :: overflow
435 integer(i8_kind),
dimension(NUMINT) :: ints
438 character(len=80) :: mesg
439 integer(i8_kind) :: ival, prec_err
442 prec_err =
prec ;
if (
present(prec_error)) prec_err = prec_error
444 if ((r >= 1e30) .eqv. (r < 1e30))
then ; nan_error = .true. ;
return ;
endif
446 sgn = 1 ;
if (r<0.0) sgn = -1
449 if (
present(overflow))
then
450 if (.not.(rs < prec_err*pr(1))) overflow = .true.
451 if ((r >= 1e30) .eqv. (r < 1e30)) overflow = .true.
452 elseif (.not.(rs < prec_err*pr(1)))
then
453 write(mesg,
'(ES13.5)') r
454 call mpp_error(fatal,
"Overflow in real_to_ints conversion of "//trim(mesg))
458 ival = int(rs*i_pr(i), 8)
467 integer(i8_kind),
dimension(NUMINT),
intent(in) :: ints
473 do i=1,
numint ; r = r + pr(i)*ints(i) ;
enddo
479 integer(i8_kind),
dimension(NUMINT),
intent(inout) :: int_sum
480 integer(i8_kind),
dimension(NUMINT),
intent(in) :: int2
481 integer(i8_kind),
optional,
intent(in) :: prec_error
486 int_sum(i) = int_sum(i) + int2(i)
488 if (int_sum(i) >
prec)
then
489 int_sum(i) = int_sum(i) -
prec
490 int_sum(i-1) = int_sum(i-1) + 1
491 elseif (int_sum(i) < -
prec)
then
492 int_sum(i) = int_sum(i) +
prec
493 int_sum(i-1) = int_sum(i-1) - 1
496 int_sum(1) = int_sum(1) + int2(1)
497 if (
present(prec_error))
then
498 if (abs(int_sum(1)) > prec_error) overflow_error = .true.
500 if (abs(int_sum(1)) >
prec) overflow_error = .true.
509 integer(i8_kind),
dimension(NUMINT),
intent(inout) :: int_sum
510 real(r8_kind),
intent(in) :: r
511 real(r8_kind),
intent(inout) :: max_mag_term
514 integer(i8_kind) :: ival
517 if ((r >= 1e30) .eqv. (r < 1e30))
then ; nan_error = .true. ;
return ;
endif
518 sgn = 1 ;
if (r<0.0) sgn = -1
520 if (rs > abs(max_mag_term)) max_mag_term = r
523 ival = int(rs*i_pr(i), 8)
525 int_sum(i) = int_sum(i) + sgn*ival
532 integer(i8_kind),
dimension(NUMINT),
intent(inout) :: int_sum
533 integer(i8_kind),
intent(in) :: prec_error
535 integer :: i, num_carry
537 do i=
numint,2,-1 ;
if (abs(int_sum(i)) >
prec)
then
538 num_carry = int(int_sum(i) *
i_prec)
539 int_sum(i) = int_sum(i) - num_carry*
prec
540 int_sum(i-1) = int_sum(i-1) + num_carry
542 if (abs(int_sum(1)) > prec_error)
then
543 overflow_error = .true.
551 integer(i8_kind),
dimension(NUMINT),
intent(inout) :: int_sum
554 integer :: i, num_carry
556 do i=
numint,2,-1 ;
if (abs(int_sum(i)) >
prec)
then
557 num_carry = int(int_sum(i) *
i_prec)
558 int_sum(i) = int_sum(i) - num_carry*
prec
559 int_sum(i-1) = int_sum(i-1) + num_carry
565 if (abs(int_sum(i)) > 0)
then
566 if (int_sum(i) < 0) positive = .false.
572 do i=
numint,2,-1 ;
if (int_sum(i) < 0)
then
573 int_sum(i) = int_sum(i) +
prec
574 int_sum(i-1) = int_sum(i-1) - 1
577 do i=
numint,2,-1 ;
if (int_sum(i) > 0)
then
578 int_sum(i) = int_sum(i) -
prec
579 int_sum(i-1) = int_sum(i-1) + 1
585 function mpp_query_efp_overflow_error()
586 logical :: mpp_query_efp_overflow_error
587 mpp_query_efp_overflow_error = overflow_error
588 end function mpp_query_efp_overflow_error
590 subroutine mpp_reset_efp_overflow_error()
591 overflow_error = .false.
592 end subroutine mpp_reset_efp_overflow_error
594 function mpp_efp_plus(EFP1, EFP2)
601 end function mpp_efp_plus
603 function mpp_efp_minus(EFP1, EFP2)
608 do i=1,
numint ; mpp_efp_minus%v(i) = -1*efp2%v(i) ;
enddo
611 end function mpp_efp_minus
616 subroutine mpp_efp_assign(EFP1, EFP2)
621 do i=1,
numint ; efp1%v(i) = efp2%v(i) ;
enddo
622 end subroutine mpp_efp_assign
624 function mpp_efp_to_real(EFP1)
626 real(r8_kind) :: mpp_efp_to_real
630 end function mpp_efp_to_real
632 function mpp_efp_real_diff(EFP1, EFP2)
634 real(r8_kind) :: mpp_efp_real_diff
638 efp_diff = efp1 - efp2
639 mpp_efp_real_diff = mpp_efp_to_real(efp_diff)
641 end function mpp_efp_real_diff
643 function mpp_real_to_efp(val, overflow)
644 real(r8_kind),
intent(in) :: val
645 logical,
optional,
intent(inout) :: overflow
649 character(len=80) :: mesg
651 if (
present(overflow))
then
652 mpp_real_to_efp%v(:) =
real_to_ints(val, overflow=overflow)
657 write(mesg,
'(ES13.5)') val
658 call mpp_error(fatal,
"Overflow in mpp_real_to_efp conversion of "//trim(mesg))
662 end function mpp_real_to_efp
668 integer,
intent(in) :: nval
669 logical,
dimension(:),
optional,
intent(out) :: errors
671 integer(i8_kind),
dimension(NUMINT,nval) :: ints
672 integer(i8_kind) :: prec_error
673 logical :: error_found
674 character(len=256) :: mesg
678 "mpp_efp_list_sum_across_PEs: Too many processors are being used for the value of "//&
679 "prec. Reduce prec to (2^63-1)/mpp_npes.")
681 prec_error = (2_8**62 + (2_8**62 - 1)) /
mpp_npes()
683 overflow_error = .false. ; error_found = .false.
685 do i=1,nval ;
do n=1,
numint ; ints(n,i) = efps(i)%v(n) ;
enddo ;
enddo
689 if (
present(errors)) errors(:) = .false.
691 overflow_error = .false.
693 do n=1,
numint ; efps(i)%v(n) = ints(n,i) ;
enddo
694 if (
present(errors)) errors(i) = overflow_error
695 if (overflow_error)
then
696 write (mesg,
'("mpp_efp_list_sum_across_PEs error at ",i6," val was ",ES13.6, ", prec_error = ",ES13.6)') &
697 i, mpp_efp_to_real(efps(i)), real(prec_error)
700 error_found = error_found .or. overflow_error
702 if (error_found .and. .not.(
present(errors)))
then
703 call mpp_error(fatal,
"Overflow in mpp_efp_list_sum_across_PEs.")
708 end module mpp_efp_mod
integer, parameter max_count_prec
The number of values that can be added together with the current value of prec before there will be r...
subroutine increment_ints_faster(int_sum, r, max_mag_term)
This subroutine increments a number with another, both using the integer representation in real_to_in...
logical debug
Making this true enables debugging output.
subroutine carry_overflow(int_sum, prec_error)
This subroutine handles carrying of the overflow.
real(r8_kind), parameter i_prec
The inverse of prec.
integer(i8_kind) function, dimension(numint) real_to_ints(r, prec_error, overflow)
This function converts a real number to an equivalent representation using several long integers.
real(r8_kind) function mpp_reproducing_sum_r8_2d(array, isr, ier, jsr, jer, EFP_sum, reproducing, overflow_check, err)
Calculates a reproducing sum for a 2D, 8-byte real array.
real(r8_kind), parameter r_prec
A real version of prec.
integer, parameter numint
The number of long integers to use to represent a real number.
integer(i8_kind), parameter prec
The precision of each integer.
real(r8_kind) function mpp_reproducing_sum_r8_3d(array, isr, ier, jsr, jer, sums, EFP_sum, err)
Reproducing sum for 3d arrays of 8-bit reals.
subroutine increment_ints(int_sum, int2, prec_error)
This subroutine increments a number with another, both using the integer representation in real_to_in...
real(r4_kind) function mpp_reproducing_sum_r4_2d(array, isr, ier, jsr, jer, EFP_sum, reproducing, overflow_check, err)
real(r8_kind) function ints_to_real(ints)
This function reverses the conversion in real_to_ints.
subroutine regularize_ints(int_sum)
This subroutine carries the overflow, and then makes sure that all integers are of the same sign as t...
subroutine, public mpp_efp_list_sum_across_pes(EFPs, nval, errors)
This subroutine does a sum across PEs of a list of EFP variables, returning the sums in place,...
This interface uses a conversion to an integer representation of real numbers to give order-invariant...
The Extended Fixed Point (mpp_efp) type provides a public interface for doing sums and taking differe...
integer function mpp_npes()
Returns processor count for current pelist.
integer function mpp_pe()
Returns processor ID.