`stdlib_dot_product` and `stdlib_dot_product_kahan` return erroneous results for complex vectors

[loiseaujc]

Description

The implementations of the improved dot_product kernels work like a charm for real vectors but return erroneous results for complex ones. It seems like the complex versions are missing a complex-conjugate operation. See the following code snippet for reproducibility.

program main
   use stdlib_linalg_constants
   use stdlib_intrinsics
   implicit none

   integer(ilp), parameter :: n = 128
   complex(dp), dimension(n) :: x, y
   real(dp), dimension(n, 2) :: z

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

   print *, dot_product(x, y)
   print *, stdlib_dot_product(x, y)
   print *, stdlib_dot_product_kahan(x, y)
end program main

Expected Behaviour

Results between stdlib_dot_product, stdlib_dot_product_kahan and the Fortran intrinsic dot_product should be identical (up to a multiple of the machine precision).

Version of stdlib

latest

Platform and Architecture

Ubuntu 22.04

Additional Information

Ping @perazz, @jalvesz

[jalvesz]

humm indeed ... here is a CE snippet to test it https://godbolt.org/z/Mza5P7bae ... So it would seem that a final conj( ) is missing though it looks strange as X*Y = sum_{i=1,N}( x_i * conj(y_i) ) by definition, right?

[loiseaujc]

It is actually the opposite: <x, y> = x.H @ y

[bjbraams]

If I read it correctly the test program for this Kahan dot product and the Kahan sum only checks that the Kahan implementation gives a result consistent with standard floating point arithmetic, not that it gives a result of the desired high accuracy. (I am looking at test/intrinsics/test_intrinsics.fypp and under the tests for the real sum and dot product I notice the statement tolerance = epsilon(1._${k}$)*100.) In some cases it might be enough to guard against gross errors, as these tests would do, but Kahan summation is vulnerable to subtle errors too. Or maybe my reading is wrong.

[jalvesz]

@bjbraams sorry, I didn't understand your question or proposal. Would you mind reformulating?

Both the chunked and chunked-kahan implementations for sum and dot_product return most of the time results which are at worst 1 fold more accurate than the intrinsics, at best 2 or 3 folds. This comes from the empirical finding that compilers try to respect the stric sequential ordering of the array for the sum. This constrain is purposefully broken here.

Floating-point accuracy is a tricky matter, you could have a blast reading this long thread https://fortran-lang.discourse.group/t/some-intrinsic-sums/5760 and also this short one https://fortran-lang.discourse.group/t/dot-product-in-single-precision-emulating-double-precision/9092/2

If you have a proposal to improve the testing methodology please, feel free to discuss it 🙂

[bjbraams]

Kahan summation (also called compensated summation) is a way of emulating additional precision arithmetic. One expects a result (of a sum or of a dot product) that differs from the exact result by at most macheps times the largest (in magnitude) intermediate result. In particular, for a sequence of terms all of the same sign one expects a result that is off by at most an ulp. My point is that a test program for Kahan summation or Kahan inner product should test that stringent level of accuracy.

[jalvesz]

Sure, would you like to open a PR ? Let's discuss what the test should look like in a separate issue dedicated to your proposal.

[bjbraams]

Sure, would you like to open a PR ? Let's discuss what the test should look like in a separate issue dedicated to your proposal.

I have no time for this now. Anyway, a pragmatic approach for a test of accuracy is to compare the Kahan dot product with BLAS sdsdot or BLAS dqddot. I think that the test may show an unpleasant surprise for at least some implementations of the library because Kahan summation is so sensitive to compiler treatments. For sure the compiler must respect parentheses and while that is demanded by the Fortran standard it may be broken by unfortunate compiler options such as -fast. More subtly, Kahan summation (as implemented in the library codes as they stand now) can be broken if register arithmetic uses higher precision than the storage format, e.g. 80 bits vs 64 bits.