Skip to content

Commit eaffa4a

Browse files
committed
fix test: complex assignment caused accuracy loss
1 parent 7cea1fd commit eaffa4a

File tree

1 file changed

+37
-62
lines changed

1 file changed

+37
-62
lines changed

test/intrinsics/test_intrinsics.fypp

Lines changed: 37 additions & 62 deletions
Original file line numberDiff line numberDiff line change
@@ -28,15 +28,15 @@ subroutine test_sum(error)
2828
type(error_type), allocatable, intent(out) :: error
2929

3030
!> Internal parameters and variables
31-
integer, parameter :: n = 1e3, ncalc = 3, niter = 100
31+
integer, parameter :: n = 1e3, ncalc = 3
3232
real(sp) :: u
3333
integer :: iter, i, j
3434
!====================================================================================
3535
#:for k1, t1, s1 in R_KINDS_TYPES
3636
block
3737
${t1}$, allocatable :: x(:)
3838
${t1}$, parameter :: total_sum = 4*atan(1._${k1}$), tolerance = epsilon(1._${k1}$)*100
39-
${t1}$ :: xsum(ncalc), meanval(ncalc), err(ncalc)
39+
${t1}$ :: xsum(ncalc), err(ncalc)
4040
logical, allocatable :: mask(:), nmask(:)
4141

4242
allocate(x(n))
@@ -54,26 +54,18 @@ subroutine test_sum(error)
5454
call swap( nmask(i), nmask(j) )
5555
end do
5656

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

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

69-
err(:) = 0._${k1}$
70-
do iter = 1, niter
71-
xsum(1) = sum(x,mask)+sum(x,nmask) ! compiler intrinsic
72-
xsum(2) = fsum_kahan(x,mask)+fsum_kahan(x,nmask) ! chunked Kahan summation
73-
xsum(3) = fsum(x,mask)+fsum(x,nmask) ! chunked summation
74-
err(1:ncalc) = err(1:ncalc) + abs(1._${k1}$-xsum(1:ncalc)/total_sum)
75-
end do
76-
err(1:ncalc) = err(1:ncalc) / niter
65+
xsum(1) = sum(x,mask)+sum(x,nmask) ! compiler intrinsic
66+
xsum(2) = fsum_kahan(x,mask)+fsum_kahan(x,nmask) ! chunked Kahan summation
67+
xsum(3) = fsum(x,mask)+fsum(x,nmask) ! chunked summation
68+
err(1:ncalc) = abs(1._${k1}$-xsum(1:ncalc)/total_sum)
7769

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

9183
allocate(x(n))
9284
do i = 1, n
93-
x(i) = cmplx(&
94-
8*atan(1._${k1}$)*(real(i,kind=${k1}$)-0.5_${k1}$)/real(2*n,kind=${k1}$)**2,&
95-
8*atan(1._${k1}$)*(real(i+n,kind=${k1}$)-0.5_${k1}$)/real(2*n,kind=${k1}$)**2)
85+
x(i) = (8*atan(1._${k1}$)*(real(i,kind=${k1}$)-0.5_${k1}$)/n**2)*cmplx(1._${k1}$,1._${k1}$)
9686
end do
87+
9788
allocate(mask(n),source=.false.); mask(1:n:2) = .true.
9889
allocate(nmask(n)); nmask = .not.mask
9990
! scramble array
@@ -105,26 +96,20 @@ subroutine test_sum(error)
10596
call swap( nmask(i), nmask(j) )
10697
end do
10798

108-
err(:) = 0._${k1}$
109-
do iter = 1, niter
110-
xsum(1) = sum(x) ! compiler intrinsic
111-
xsum(2) = fsum_kahan(x) ! chunked Kahan summation
112-
xsum(3) = fsum(x) ! chunked summation
113-
err(1:ncalc) = err(1:ncalc) + abs(1._${k1}$-(xsum(1:ncalc)%re+xsum(1:ncalc)%im)/total_sum)
114-
end do
115-
err(1:ncalc) = err(1:ncalc) / niter
99+
xsum(1) = sum(x) ! compiler intrinsic
100+
xsum(2) = fsum_kahan(x) ! chunked Kahan summation
101+
xsum(3) = fsum(x) ! chunked summation
102+
err(1:ncalc) = abs(1._${k1}$-(xsum(1:ncalc)%re)/total_sum)
103+
116104

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

120-
err(:) = 0._${k1}$
121-
do iter = 1, niter
122-
xsum(1) = sum(x,mask)+sum(x,nmask) ! compiler intrinsic
123-
xsum(2) = fsum_kahan(x,mask)+fsum_kahan(x,nmask) ! chunked Kahan summation
124-
xsum(3) = fsum(x,mask)+fsum(x,nmask) ! chunked summation
125-
err(1:ncalc) = err(1:ncalc) + abs(1._${k1}$-(xsum(1:ncalc)%re+xsum(1:ncalc)%im)/total_sum)
126-
end do
127-
err(1:ncalc) = err(1:ncalc) / niter
108+
xsum(1) = sum(x,mask)+sum(x,nmask) ! compiler intrinsic
109+
xsum(2) = fsum_kahan(x,mask)+fsum_kahan(x,nmask) ! chunked Kahan summation
110+
xsum(3) = fsum(x,mask)+fsum(x,nmask) ! chunked summation
111+
err(1:ncalc) = abs(1._${k1}$-(xsum(1:ncalc)%re)/total_sum)
112+
128113

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

140125
!> Internal parameters and variables
141-
integer, parameter :: n = 1e3, ncalc = 3, niter = 100
126+
integer, parameter :: n = 1e3, ncalc = 3
142127
real(sp) :: u
143128
integer :: iter, i, j
144129
!====================================================================================
145130
#:for k1, t1, s1 in R_KINDS_TYPES
146131
block
147132
${t1}$, allocatable :: x(:)
148133
${t1}$, parameter :: total_sum = 4*atan(1._${k1}$), tolerance = epsilon(1._${k1}$)*100
149-
${t1}$ :: xsum(ncalc), meanval(ncalc), err(ncalc)
134+
${t1}$ :: xsum(ncalc), err(ncalc)
150135

151136
allocate(x(n))
152137
do i = 1, n
153-
x(i) = sqrt( 8*atan(1._${k1}$)*(real(i,kind=${k1}$)-0.5_${k1}$)/real(n,kind=${k1}$)**2 )
138+
x(i) = 2*sqrt( 2*atan(1._${k1}$)*(real(i,kind=${k1}$)-0.5_${k1}$) )/n
154139
end do
155140
! scramble array
156141
do i = 1, n
@@ -159,14 +144,10 @@ subroutine test_dot_product(error)
159144
call swap( x(i), x(j) )
160145
end do
161146

162-
err(:) = 0._${k1}$
163-
do iter = 1, niter
164-
xsum(1) = dot_product(x,x) ! compiler intrinsic
165-
xsum(2) = fprod_kahan(x,x) ! chunked Kahan summation
166-
xsum(3) = fprod(x,x) ! chunked summation
167-
err(1:ncalc) = err(1:ncalc) + abs(1._${k1}$-xsum(1:ncalc)/total_sum)
168-
end do
169-
err(1:ncalc) = err(1:ncalc) / niter
147+
xsum(1) = dot_product(x,x) ! compiler intrinsic
148+
xsum(2) = fprod_kahan(x,x) ! chunked Kahan summation
149+
xsum(3) = fprod(x,x) ! chunked summation
150+
err(1:ncalc) = abs(1._${k1}$-xsum(1:ncalc)/total_sum)
170151

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

183164
allocate(x(n))
184165
do i = 1, n
185-
x(i) = cmplx(&
186-
sqrt(8*atan(1._${k1}$)*(real(i,kind=${k1}$)-0.5_${k1}$)/real(2*n,kind=${k1}$)**2),&
187-
sqrt(8*atan(1._${k1}$)*(real(i+n,kind=${k1}$)-0.5_${k1}$)/real(2*n,kind=${k1}$)**2))
166+
x(i) = ( 2*sqrt( 2*atan(1._${k1}$)*(real(i,kind=${k1}$)-0.5_${k1}$) ) / n )*cmplx(1._${k1}$,1._${k1}$)
188167
end do
189168
! scramble array
190169
do i = 1, n
@@ -193,15 +172,11 @@ subroutine test_dot_product(error)
193172
call swap( x(i), x(j) )
194173
end do
195174

196-
err(:) = 0._${k1}$
197-
do iter = 1, niter
198-
xsum(1) = dot_product(x,x) ! compiler intrinsic
199-
xsum(2) = fprod_kahan(x,x) ! chunked Kahan dot_product
200-
xsum(3) = fprod(x,x) ! chunked dot_product
201-
err(1:ncalc) = err(1:ncalc) + abs(1._${k1}$-(xsum(1:ncalc)%re+xsum(1:ncalc)%im)/total_sum)
202-
end do
203-
err(1:ncalc) = err(1:ncalc) / niter
204-
175+
xsum(1) = dot_product(x,x) ! compiler intrinsic
176+
xsum(2) = fprod_kahan(x,x) ! chunked Kahan dot_product
177+
xsum(3) = fprod(x,x) ! chunked dot_product
178+
err(1:ncalc) = abs(1._${k1}$-xsum(1:ncalc)%re/(2*total_sum))
179+
205180
call check(error, all(err(:)<tolerance) , "complex dot_product is not accurate" )
206181
if (allocated(error)) return
207182
end block

0 commit comments

Comments
 (0)