Skip to content

Commit 2552c36

Browse files
jalveszCopilot
andauthored
fix: complex dot_product formulation (fortran-lang#1017)
* fix complex dot_product formulation * smaller array and wider tolerance * Update test/intrinsics/test_intrinsics.fypp Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --------- Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
1 parent a2fadac commit 2552c36

File tree

2 files changed

+25
-4
lines changed

2 files changed

+25
-4
lines changed

src/stdlib_intrinsics_dot_product.fypp

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -34,10 +34,10 @@ pure module function stdlib_dot_product_${s}$(a,b) result(p)
3434
n = size(a,kind=ilp)
3535
r = mod(n,chunk)
3636

37-
abatch(1:r) = a(1:r)*${cnjg(t,'b(1:r)')}$
37+
abatch(1:r) = ${cnjg(t,'a(1:r)')}$*b(1:r)
3838
abatch(r+1:chunk) = zero_${s}$
3939
do i = r+1, n-r, chunk
40-
abatch(1:chunk) = abatch(1:chunk) + a(i:i+chunk-1)*${cnjg(t,'b(i:i+chunk-1)')}$
40+
abatch(1:chunk) = abatch(1:chunk) + ${cnjg(t,'a(i:i+chunk-1)')}$*b(i:i+chunk-1)
4141
end do
4242

4343
p = zero_${s}$
@@ -60,11 +60,11 @@ pure module function stdlib_dot_product_kahan_${s}$(a,b) result(p)
6060
n = size(a,kind=ilp)
6161
r = mod(n,chunk)
6262

63-
abatch(1:r) = a(1:r)*${cnjg(t,'b(1:r)')}$
63+
abatch(1:r) = ${cnjg(t,'a(1:r)')}$*b(1:r)
6464
abatch(r+1:chunk) = zero_${s}$
6565
cbatch = zero_${s}$
6666
do i = r+1, n-r, chunk
67-
call kahan_kernel( a(i:i+chunk-1)*${cnjg(t,'b(i:i+chunk-1)')}$ , abatch(1:chunk) , cbatch(1:chunk) )
67+
call kahan_kernel( ${cnjg(t,'a(i:i+chunk-1)')}$*b(i:i+chunk-1) , abatch(1:chunk) , cbatch(1:chunk) )
6868
end do
6969

7070
p = zero_${s}$

test/intrinsics/test_intrinsics.fypp

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -246,6 +246,27 @@ subroutine test_dot_product(error)
246246
call check(error, all(err(:)<tolerance) , "complex dot_product is not accurate" )
247247
if (allocated(error)) return
248248
end block
249+
250+
block ! test for https://github.com/fortran-lang/stdlib/issues/1016
251+
${t}$ :: x(128), y(128)
252+
real(${k}$) :: z(128,2)
253+
real(${k}$), parameter :: tolerance = epsilon(1._${k}$)*100000
254+
real(${k}$) :: err(2)
255+
${t}$ :: p(3)
256+
257+
call random_number(z)
258+
x%re = z(:, 1); x%im = z(:, 2)
259+
call random_number(z)
260+
y%re = z(:, 1); y%im = z(:, 2)
261+
262+
p(1) = dot_product(x,y) ! compiler intrinsic
263+
p(2) = stdlib_dot_product_kahan(x,y) ! chunked Kahan dot_product
264+
p(3) = stdlib_dot_product(x,y) ! chunked dot_product
265+
err(1:2) = sqrt((p(2:3)%re - p(1)%re)**2 + (p(2:3)%im - p(1)%im)**2)
266+
267+
call check(error, all(err(:)<tolerance) , "complex dot_product does not conform to the standard" )
268+
if (allocated(error)) return
269+
end block
249270
#:endfor
250271

251272
end subroutine

0 commit comments

Comments
 (0)