BUG: np.dot incorrect for more than 2^30 complex elements

[dchackett]

Describe the issue:

np.dot gives the wrong answer when applied to compute dot products of 1d complex arrays with length greater than 2^30 elements. Notes:

• Based on CPU usage as observed in top, np.dot used multiple threads while (x**2).sum() used only one.
• Error onset always at 2^30+1, for different random seeds.
• Difference between results for 2^30 and 2^30+1 is too large to be round-off error iteratively piling up.
• Complex dtype is important: error does not occur for float dtype for 2^30+1, 2^31+1, or 2^32+1.
• Tested for complex128, have not looked at complex64.
• Besides provided example, similarly fails for two different vectors.

Running on Ubuntu 18.04.3 LTS, reproduced on different machines w/ CPUs Intel(R) Xeon(R) CPU E5-2680 and Intel(R) Xeon(R) Gold 5218 CPU. Probably MKL for multithreaded backend:

>>> np.show_config()
blas_mkl_info:
    libraries = ['mkl_rt', 'pthread']
    library_dirs = # ...
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = # ...
# ... similar for blas_opt_info, lapack_mkl_info, lapack_opt_info ...

Reproduce the code example:

import numpy as np

x = np.random.rand(2**30+1) + 1j * np.random.rand(2**30+1)

# with 2^30 elements
y = x[:-1]
norm_dot = np.dot(y,y)
norm_sqsum = (y**2).sum()
assert np.isclose(norm_dot.real, norm_sqsum.real) and np.isclose(norm_dot.imag, norm_sqsum.imag) # passes

# include one more element
correct = norm_dot + x[-1]**2
norm_sqsum = (x**2).sum()
norm_dot = np.dot(x,x)
assert np.isclose(correct.real, norm_sqsum.real) and np.isclose(correct.imag, norm_sqsum.imag) # passes
assert np.isclose(norm_dot.real, norm_sqsum.real) and np.isclose(norm_dot.imag, norm_sqsum.imag) # fails

Error message:

No response

NumPy/Python version information:

1.21.2 3.9.10 | packaged by conda-forge | (main, Feb 1 2022, 21:24:37)
[GCC 9.4.0]

Context for the issue:

This is a dangerous silent failure.

[MatteoRaso]

This also happens for complex64 when doing 2 ** 16 elements.

[seberg]

Yap, that is a pretty annoying bug in NumPy, thanks for the report! It doesn't +matter here (the bug is clearly in NumPy), but to report your blas information, I would suggest using:

python -m threadpoolctl -i numpy

np.show_config() is currently fairly useless, at least unless you did a local/custom compilation of NumPy.

@MatteoRaso are you sure about the 2**16 observation, it really doesn't make sense to me that this could be related. Maybe it just failed for you due to random fluctuations? I.e. make sure to inspect x[-1]**2 and the difference between norm_dot and correct.

If there is something odd about 2**16, please report your blas version!