Skip to content

Commit

Permalink
fix test: complex assignment caused accuracy loss
Browse files Browse the repository at this point in the history
  • Loading branch information
@jalvesz
jalvesz committed Jan 4, 2025
1 parent 7cea1fd commit eaffa4a
Showing 1 changed file with 37 additions and 62 deletions.
99 changes: 37 additions & 62 deletions test/intrinsics/test_intrinsics.fypp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,15 +28,15 @@ subroutine test_sum(error)
type(error_type), allocatable, intent(out) :: error

!> Internal parameters and variables
integer, parameter :: n = 1e3, ncalc = 3, niter = 100
integer, parameter :: n = 1e3, ncalc = 3
real(sp) :: u
integer :: iter, i, j
!====================================================================================
#:for k1, t1, s1 in R_KINDS_TYPES
block
${t1}$, allocatable :: x(:)
${t1}$, parameter :: total_sum = 4*atan(1._${k1}$), tolerance = epsilon(1._${k1}$)*100
${t1}$ :: xsum(ncalc), meanval(ncalc), err(ncalc)
${t1}$ :: xsum(ncalc), err(ncalc)
logical, allocatable :: mask(:), nmask(:)

allocate(x(n))
Expand All @@ -54,26 +54,18 @@ subroutine test_sum(error)
call swap( nmask(i), nmask(j) )
end do

err(:) = 0._${k1}$
do iter = 1, niter
xsum(1) = sum(x) ! compiler intrinsic
xsum(2) = fsum_kahan(x) ! chunked Kahan summation
xsum(3) = fsum(x) ! chunked summation
err(1:ncalc) = err(1:ncalc) + abs(1._${k1}$-xsum(1:ncalc)/total_sum)
end do
err(1:ncalc) = err(1:ncalc) / niter
xsum(1) = sum(x) ! compiler intrinsic
xsum(2) = fsum_kahan(x) ! chunked Kahan summation
xsum(3) = fsum(x) ! chunked summation
err(1:ncalc) = abs(1._${k1}$-xsum(1:ncalc)/total_sum)

call check(error, all(err(:)<tolerance) , "real sum is not accurate" )
if (allocated(error)) return

err(:) = 0._${k1}$
do iter = 1, niter
xsum(1) = sum(x,mask)+sum(x,nmask) ! compiler intrinsic
xsum(2) = fsum_kahan(x,mask)+fsum_kahan(x,nmask) ! chunked Kahan summation
xsum(3) = fsum(x,mask)+fsum(x,nmask) ! chunked summation
err(1:ncalc) = err(1:ncalc) + abs(1._${k1}$-xsum(1:ncalc)/total_sum)
end do
err(1:ncalc) = err(1:ncalc) / niter
xsum(1) = sum(x,mask)+sum(x,nmask) ! compiler intrinsic
xsum(2) = fsum_kahan(x,mask)+fsum_kahan(x,nmask) ! chunked Kahan summation
xsum(3) = fsum(x,mask)+fsum(x,nmask) ! chunked summation
err(1:ncalc) = abs(1._${k1}$-xsum(1:ncalc)/total_sum)

call check(error, all(err(:)<tolerance) , "masked real sum is not accurate" )
if (allocated(error)) return
Expand All @@ -85,15 +77,14 @@ subroutine test_sum(error)
${t1}$, allocatable :: x(:)
real(${k1}$), parameter :: total_sum = 4*atan(1._${k1}$), tolerance = epsilon(1._${k1}$)*100
real(${k1}$) :: err(ncalc)
${t1}$ :: xsum(ncalc), meanval(ncalc)
${t1}$ :: xsum(ncalc)
logical, allocatable :: mask(:), nmask(:)

allocate(x(n))
do i = 1, n
x(i) = cmplx(&
8*atan(1._${k1}$)*(real(i,kind=${k1}$)-0.5_${k1}$)/real(2*n,kind=${k1}$)**2,&
8*atan(1._${k1}$)*(real(i+n,kind=${k1}$)-0.5_${k1}$)/real(2*n,kind=${k1}$)**2)
x(i) = (8*atan(1._${k1}$)*(real(i,kind=${k1}$)-0.5_${k1}$)/n**2)*cmplx(1._${k1}$,1._${k1}$)
end do

allocate(mask(n),source=.false.); mask(1:n:2) = .true.
allocate(nmask(n)); nmask = .not.mask
! scramble array
Expand All @@ -105,26 +96,20 @@ subroutine test_sum(error)
call swap( nmask(i), nmask(j) )
end do

err(:) = 0._${k1}$
do iter = 1, niter
xsum(1) = sum(x) ! compiler intrinsic
xsum(2) = fsum_kahan(x) ! chunked Kahan summation
xsum(3) = fsum(x) ! chunked summation
err(1:ncalc) = err(1:ncalc) + abs(1._${k1}$-(xsum(1:ncalc)%re+xsum(1:ncalc)%im)/total_sum)
end do
err(1:ncalc) = err(1:ncalc) / niter
xsum(1) = sum(x) ! compiler intrinsic
xsum(2) = fsum_kahan(x) ! chunked Kahan summation
xsum(3) = fsum(x) ! chunked summation
err(1:ncalc) = abs(1._${k1}$-(xsum(1:ncalc)%re)/total_sum)


call check(error, all(err(:)<tolerance) , "complex sum is not accurate" )
if (allocated(error)) return

err(:) = 0._${k1}$
do iter = 1, niter
xsum(1) = sum(x,mask)+sum(x,nmask) ! compiler intrinsic
xsum(2) = fsum_kahan(x,mask)+fsum_kahan(x,nmask) ! chunked Kahan summation
xsum(3) = fsum(x,mask)+fsum(x,nmask) ! chunked summation
err(1:ncalc) = err(1:ncalc) + abs(1._${k1}$-(xsum(1:ncalc)%re+xsum(1:ncalc)%im)/total_sum)
end do
err(1:ncalc) = err(1:ncalc) / niter
xsum(1) = sum(x,mask)+sum(x,nmask) ! compiler intrinsic
xsum(2) = fsum_kahan(x,mask)+fsum_kahan(x,nmask) ! chunked Kahan summation
xsum(3) = fsum(x,mask)+fsum(x,nmask) ! chunked summation
err(1:ncalc) = abs(1._${k1}$-(xsum(1:ncalc)%re)/total_sum)


call check(error, all(err(:)<tolerance) , "complex masked sum is not accurate" )
if (allocated(error)) return
Expand All @@ -138,19 +123,19 @@ subroutine test_dot_product(error)
type(error_type), allocatable, intent(out) :: error

!> Internal parameters and variables
integer, parameter :: n = 1e3, ncalc = 3, niter = 100
integer, parameter :: n = 1e3, ncalc = 3
real(sp) :: u
integer :: iter, i, j
!====================================================================================
#:for k1, t1, s1 in R_KINDS_TYPES
block
${t1}$, allocatable :: x(:)
${t1}$, parameter :: total_sum = 4*atan(1._${k1}$), tolerance = epsilon(1._${k1}$)*100
${t1}$ :: xsum(ncalc), meanval(ncalc), err(ncalc)
${t1}$ :: xsum(ncalc), err(ncalc)

allocate(x(n))
do i = 1, n
x(i) = sqrt( 8*atan(1._${k1}$)*(real(i,kind=${k1}$)-0.5_${k1}$)/real(n,kind=${k1}$)**2 )
x(i) = 2*sqrt( 2*atan(1._${k1}$)*(real(i,kind=${k1}$)-0.5_${k1}$) )/n
end do
! scramble array
do i = 1, n
Expand All @@ -159,14 +144,10 @@ subroutine test_dot_product(error)
call swap( x(i), x(j) )
end do

err(:) = 0._${k1}$
do iter = 1, niter
xsum(1) = dot_product(x,x) ! compiler intrinsic
xsum(2) = fprod_kahan(x,x) ! chunked Kahan summation
xsum(3) = fprod(x,x) ! chunked summation
err(1:ncalc) = err(1:ncalc) + abs(1._${k1}$-xsum(1:ncalc)/total_sum)
end do
err(1:ncalc) = err(1:ncalc) / niter
xsum(1) = dot_product(x,x) ! compiler intrinsic
xsum(2) = fprod_kahan(x,x) ! chunked Kahan summation
xsum(3) = fprod(x,x) ! chunked summation
err(1:ncalc) = abs(1._${k1}$-xsum(1:ncalc)/total_sum)

call check(error, all(err(:)<tolerance) , "real dot_product is not accurate" )
if (allocated(error)) return
Expand All @@ -178,13 +159,11 @@ subroutine test_dot_product(error)
${t1}$, allocatable :: x(:)
real(${k1}$), parameter :: total_sum = 4*atan(1._${k1}$), tolerance = epsilon(1._${k1}$)*100
real(${k1}$) :: err(ncalc)
${t1}$ :: xsum(ncalc), meanval(ncalc)
${t1}$ :: xsum(ncalc)

allocate(x(n))
do i = 1, n
x(i) = cmplx(&
sqrt(8*atan(1._${k1}$)*(real(i,kind=${k1}$)-0.5_${k1}$)/real(2*n,kind=${k1}$)**2),&
sqrt(8*atan(1._${k1}$)*(real(i+n,kind=${k1}$)-0.5_${k1}$)/real(2*n,kind=${k1}$)**2))
x(i) = ( 2*sqrt( 2*atan(1._${k1}$)*(real(i,kind=${k1}$)-0.5_${k1}$) ) / n )*cmplx(1._${k1}$,1._${k1}$)
end do
! scramble array
do i = 1, n
Expand All @@ -193,15 +172,11 @@ subroutine test_dot_product(error)
call swap( x(i), x(j) )
end do

err(:) = 0._${k1}$
do iter = 1, niter
xsum(1) = dot_product(x,x) ! compiler intrinsic
xsum(2) = fprod_kahan(x,x) ! chunked Kahan dot_product
xsum(3) = fprod(x,x) ! chunked dot_product
err(1:ncalc) = err(1:ncalc) + abs(1._${k1}$-(xsum(1:ncalc)%re+xsum(1:ncalc)%im)/total_sum)
end do
err(1:ncalc) = err(1:ncalc) / niter

xsum(1) = dot_product(x,x) ! compiler intrinsic
xsum(2) = fprod_kahan(x,x) ! chunked Kahan dot_product
xsum(3) = fprod(x,x) ! chunked dot_product
err(1:ncalc) = abs(1._${k1}$-xsum(1:ncalc)%re/(2*total_sum))

call check(error, all(err(:)<tolerance) , "complex dot_product is not accurate" )
if (allocated(error)) return
end block
Expand Down

0 comments on commit eaffa4a

Please sign in to comment.