Skip to content

fix: complex dot_product formulation #1017

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Jul 23, 2025
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 4 additions & 4 deletions src/stdlib_intrinsics_dot_product.fypp
Original file line number Diff line number Diff line change
Expand Up @@ -34,10 +34,10 @@ pure module function stdlib_dot_product_${s}$(a,b) result(p)
n = size(a,kind=ilp)
r = mod(n,chunk)

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

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

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

p = zero_${s}$
Expand Down
21 changes: 21 additions & 0 deletions test/intrinsics/test_intrinsics.fypp
Original file line number Diff line number Diff line change
Expand Up @@ -246,6 +246,27 @@ subroutine test_dot_product(error)
call check(error, all(err(:)<tolerance) , "complex dot_product is not accurate" )
if (allocated(error)) return
end block

block ! test for https://github.com/fortran-lang/stdlib/issues/1016
${t}$ :: x(128), y(128)
real(${k}$) :: z(128,2)
real(${k}$), parameter :: tolerance = epsilon(1._${k}$)*100000
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For a 32-bit precision dot product, this error is ~0.01, which seems a bit large? I believe using sqrt(epsilon(0.0_${k}$)) would be a more established error threshold in numerical analysis. LGTM otherwise.

real(${k}$) :: err(2)
${t}$ :: p(3)

call random_number(z)
x%re = z(:, 1); x%im = z(:, 2)
call random_number(z)
y%re = z(:, 1); y%im = z(:, 2)

p(1) = dot_product(x,y) ! compiler intrinsic
p(2) = stdlib_dot_product_kahan(x,y) ! chunked Kahan dot_product
p(3) = stdlib_dot_product(x,y) ! chunked dot_product
err(1:2) = sqrt((p(2:3)%re - p(1)%re)**2 + (p(2:3)%im - p(1)%im)**2)

call check(error, all(err(:)<tolerance) , "complex dot_product does not conform to the standard" )
if (allocated(error)) return
end block
#:endfor

end subroutine
Expand Down
Loading