`cov(A, A)` test failure on neoverse-v2 CPU

[giordano]

julia> Base.runtests(["Statistics"])
Running parallel tests with:
  getpid() = 4055299
  nworkers() = 1
  nthreads() = 1
  Sys.CPU_THREADS = 36
  Sys.total_memory() = 477.556 GiB
  Sys.free_memory() = 317.340 GiB

Test   (Worker) | Time (s) | GC (s) | GC % | Alloc (MB) | RSS (MB)
Statistics  (1) |        started at 2025-09-26T14:18:15.898
Statistics  (1) |         failed at 2025-09-26T14:18:58.077
Test Failed at /data/cceamgi/julia-depot/juliaup/julia-1.11.7+0.aarch64.linux.gnu/share/julia/stdlib/v1.11/Statistics/test/runtests.jl:477
  Expression: cov(A) ≈ cov(A, A) ≈ (cov(reshape(A, :, 1)))[1] ≈ (cov(reshape(A, :, 1)))[1] ≈ var(A)
   Evaluated: 400.0748f0 ≈ 399.82156f0 ≈ 400.07483f0 ≈ 400.07483f0 ≈ 400.0748f0



Test Summary:  | Pass  Fail  Total   Time
  Overall      |  838     1    839  43.6s
    Statistics |  838     1    839  43.0s
    FAILURE

The global RNG seed was 0x57895fed49e65e1f3d0e86c30f372f74.

Error in testset Statistics:
Test Failed at /data/cceamgi/julia-depot/juliaup/julia-1.11.7+0.aarch64.linux.gnu/share/julia/stdlib/v1.11/Statistics/test/runtests.jl:477
  Expression: cov(A) ≈ cov(A, A) ≈ (cov(reshape(A, :, 1)))[1] ≈ (cov(reshape(A, :, 1)))[1] ≈ var(A)
   Evaluated: 400.0748f0 ≈ 399.82156f0 ≈ 400.07483f0 ≈ 400.07483f0 ≈ 400.0748f0

ERROR: LoadError: Test run finished with errors
in expression starting at /data/cceamgi/julia-depot/juliaup/julia-1.11.7+0.aarch64.linux.gnu/share/julia/test/runtests.jl:100
ERROR: A test has failed. Please submit a bug report (https://github.com/JuliaLang/julia/issues)
including error messages above and the output of versioninfo():
Julia Version 1.11.7
Commit f2b3dbda30a (2025-09-08 12:10 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (aarch64-linux-gnu)
  CPU: 72 × unknown
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, neoverse-v2)
Threads: 1 default, 0 interactive, 1 GC (on 72 virtual cores)

First reported in JuliaLang/julia#59664.

[vchuravy]

Taking a random guess here, cov looks relatively harmless until you get to unscaled_covzm which does adjoint(y) * x which then ends up calling dot/BLAS.dot.

So there is a remote chance that this is somewhere in cblas_sdot.

[giordano]

This is fixed on master, I believe by #166 which makes cov(A::Vector, A::Vector) do cov(A), which is what's being tested here.

[vchuravy]

Yeah, but it shouldn't be wrong though?

[giordano]

Yes, but I'm having a hard time reproducing the error outside of the tests. I replaced

Statistics.jl/test/runtests.jl

Lines 476 to 479 in 68869af

	A = 20*randn(Float32, 10_000_000) .+ 100
	@test cov(A) ≈ cov(A, A) ≈
	cov(reshape(A, :, 1))[1] ≈ cov(reshape(A, :, 1))[1] ≈
	var(A)

with

        # Large Vector{Float32}
        A = 20*randn(Float32, 10_000_000) .+ 100
        write("/home/cceamgi/tmp/repro-statistics-194/A.dat", A)
        @test cov(A) ≈ cov(A, A) ≈
            cov(reshape(A, :, 1))[1] ≈ cov(reshape(A, :, 1))[1] ≈
            var(A)

to exfiltrate the offending data, but then

julia> A = collect(reinterpret(Float32, read("A.dat")));

julia> using Statistics

julia> cov(A), cov(A, A)
(400.0748f0, 400.07352f0)

julia> cov(A) ≈ cov(A, A)
true

Somehow in the REPL cov(A, A) is more accurate.

[giordano]

I finally managed to reproduce the issue: setting a single BLAS thread (what julia's runtests does) is important. When using multiple BLAS threads (the default in the REPL), the inaccuracy goes away. So there's indeed something wrong in OpenBLAS's single-threaded cblas_sdot, and this is specific to the neoversev1 kernel:

$ OPENBLAS_VERBOSE=2 OPENBLAS_NUM_THREADS=1 julia -q
Core: neoversev1
julia> using Statistics

julia> A = 20*randn(Float32, 10_000_000) .+ 100;

julia> cov(A), cov(A, A)
(400.16788f0, 399.9151f0)

julia> cov(A) ≈ cov(A, A)
false

julia> 
$ OPENBLAS_VERBOSE=2 OPENBLAS_CORETYPE=armv8 OPENBLAS_NUM_THREADS=1 julia -q
Core: armv8
julia> using Statistics

julia> A = 20*randn(Float32, 10_000_000) .+ 100;

julia> cov(A), cov(A, A)
(400.22498f0, 400.1347f0)

julia> cov(A) ≈ cov(A, A)
true

julia>
$ OPENBLAS_VERBOSE=2 OPENBLAS_CORETYPE=neoversen1 OPENBLAS_NUM_THREADS=1 julia -q
Core: neoversen1
julia> using Statistics

julia> A = 20*randn(Float32, 10_000_000) .+ 100;

julia> cov(A), cov(A, A)
(399.82593f0, 399.7934f0)

julia> cov(A) ≈ cov(A, A)
true