Skip to content

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

@dchackett

Description

@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.

Metadata

Metadata

Assignees

No one assigned

    Labels

    00 - BugPriority: highHigh priority, also add milestones for urgent issues

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions