[ITensors] [BUG] Different results when using BLAS multithreading

[MaryamR99]

Description of bug

When I use BLAS.set_num_threads(1), I get a different result than when I set it to 4 (or a different number other than 1).

Minimal code demonstrating the bug or unexpected behavior

Minimal runnable code

using ITensorMPS, ITensors, LinearAlgebra, Random

function run_test(thread_count::Int, N, p, χ, angle)

    BLAS.set_num_threads(thread_count)
    
    Random.seed!(1234)
    println("Using N = $N, χ = $χ, Threads = $thread_count")
    s = siteinds("S=1/2", N)
    ψ = randomMPS(s; linkdims=χ)

    angles_1 = fill(angle, p)
    angles_2 = fill(angle, p)
    for k in 1:p

        for n1 in 1:N
            
            for n2 in n1+1:N

                ϕ = copy(ψ)
                gate = op("Z", s, n1) 
                ϕ[n1] = noprime( gate* ϕ[n1])

                gate = op("Z", s, n2) 
                ϕ[n2] = noprime( gate* ϕ[n2])
                sine, cose = sincos(angles_1[k])

                ψ_sum = +(cose * ψ, 1im * sine * ϕ; alg="directsum")

                truncate!(ψ_sum; maxdim=χ)
                normalize!(ψ_sum)

                ψ = ψ_sum
            end
        end
            
        for i in 1:N
            gateX = op("X", s, i)
            gatee = cos(angles_2[k]) * op("Id", s, i) - 1im * sin(angles_2[k]) * gateX
            ψ[i] = noprime( gatee * ψ[i])
        end
    end


    H_C = 0.0
    for n1 in 1:N
        for n2 in n1+1:N

            ϕ = copy(ψ)
            gate = op("Z", s, n1) 
            noprime!(ϕ[n1] *= gate)

            gate = op("Z", s, n2) 
            noprime!(ϕ[n2] *= gate)

            # Compute ⟨ψ|Z_i Z_j|ψ⟩
            zz_expectation = ITensorMPS.inner(ψ, ϕ)
            H_C +=  real(zz_expectation)
        end
    end
    return H_C, ψ
end




println("--- Starting Reproducibility Test ---")

N = 16
p = 20
χ = 40
angle = 0.3

println(Threads.nthreads())


println("Running with 1 thread...")
H_ref, ψ_ref = run_test(1, N, p, χ, angle)



println("Running with 4 threads...") 
H_test, ψ_test = run_test(4, N, p, χ, angle) 


# Compare


println("\n------------------------------------------------")
println("Results:")
println("H_ref:  $H_ref")
println("H_test: $H_test")

println("state overlap: ", abs(ITensorMPS.inner(ψ_ref, ψ_test)))

Expected output or behavior

The expectation value of the energy should be the same, between the two cases.
Also, the overlap of the final states, |<ψ_ref|ψ_test>|, should be almost 1.

H_ref:  2.451631441113766
H_test: 2.4615310779333828
state overlap: 0.9984877783164545

Version information

• Output from versioninfo():

julia> versioninfo()
Julia Version 1.12.4
Commit 01a2eadb047 (2026-01-06 16:56 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: macOS (arm64-apple-darwin24.0.0)
  CPU: 8 × Apple M3
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, apple-m3)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 4 virtual cores)

• Output from using Pkg; Pkg.status("ITensors"):

julia> using Pkg; Pkg.status("ITensors")
  [9136182c] ITensors v0.9.18

[mtfishman]

Thanks for the nice issue report!

Multithreading in general (including BLAS and LAPACK multithreading) may not give the exact same floating point results for different numbers of threads and even different runs, because operations are assigned non-deterministically to threads and therefore operations can occur in different orderings, and floating point operations aren't guaranteed commute with each other exactly (i.e. in general there can be numerical roundoff errors).

Having said that, the difference between those results looks quite high for floating point roundoff errors. Do you have a sense for what precision you are expecting to get in your results given the bond dimension/truncations?

I'm guessing that you are using OpenBLAS, which is the default BLAS/LAPACK backend that Julia uses (you can check with BLAS.get_config()), and I'm able to reproduce your results on my macOS using OpenBLAS. An interesting thing to note is that when I try Apple Accelerate as a BLAS/LAPACK backend (https://github.com/JuliaLinearAlgebra/AppleAccelerate.jl), the results become reproducible:

julia> include("run.jl")
--- Starting Reproducibility Test ---
1
Running with 1 thread...
Using N = 16, χ = 40, Threads = 1
Running with 4 threads...
Using N = 16, χ = 40, Threads = 4

------------------------------------------------
Results:
H_ref:  2.4516295429352404
H_test: 2.4516295429352404
state overlap: 0.9999999999999967

julia> BLAS.get_config()
LinearAlgebra.BLAS.LBTConfig
Libraries: 
├ [ILP64] libopenblas64_.dylib
├ [ LP64] Accelerate
└ [ILP64] Accelerate

julia> versioninfo()
Julia Version 1.12.4
Commit 01a2eadb047 (2026-01-06 16:56 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: macOS (arm64-apple-darwin24.0.0)
  CPU: 10 × Apple M1 Max
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, apple-m1)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 8 virtual cores)
Environment:
  JULIA_EDITOR = code

Note that you need to control the threads with AppleAccelerate.set_num_threads(thread_count).

You may be hitting a bug in OpenBLAS. Can you reproduce this kind of discrepancy with a more minimal (and faster to run) code?

[MaryamR99]

Thank you for your reply.

About the precision..., for the same set of parameters and different angles, I get the following results:

N = 16
p = 20
χ = 40
angle = 0.2

------------------------------------------------
Results:
H_ref:  16.79735469170988
H_test: 16.79684605444248
state overlap: 0.9999600490701943
------------------------------------------------
------------------------------------------------
N = 16
p = 20
χ = 40
angle = 0.7

------------------------------------------------
Results:
H_ref:  -4.677929029427252
H_test: -4.677927911486912
state overlap: 0.9999999999039759

Also, if I change the bond dimension, I get these results:

N = 16
p = 20
χ = 25
angle = 0.3

------------------------------------------------
Results:
H_ref:  0.7779943967727879
H_test: 0.7779943967727879
state overlap: 0.9999999999999976
------------------------------------------------
------------------------------------------------
N = 16
p = 20
χ = 80
angle = 0.3

------------------------------------------------
Results:
H_ref:  9.394290907711547
H_test: 9.413608009466703
state overlap: 0.9943564741365561

Yes, I am using OpenBLAS,

julia> BLAS.get_config()
LinearAlgebra.BLAS.LBTConfig
Libraries: 
└ [ILP64] libopenblas64_.dylib

I tested the code on Linux, also using OpenBLAS, and there was almost the same amount of difference between the results.

As you suggested, by using AppleAccelerate.set_num_threads(thread_count), the results became the same (but they are different from your results)

N = 16
p = 20
χ = 40
angle = 0.3

------------------------------------------------
Results:
H_ref:  2.461500037657551
H_test: 2.461500037657551
state overlap: 1.0000000000000013


------------------------------------------------
julia> BLAS.get_config()
LinearAlgebra.BLAS.LBTConfig
Libraries: 
├ [ILP64] libopenblas64_.dylib
├ [ LP64] Accelerate
└ [ILP64] Accelerate

Unfortunately, this was my attempt to get this discrepancy in a minimal example from my original code! I will try more to find a better and faster example.

[emstoudenmire]

Based on the different OpenBLAS versus Accelerate behaviors, is it fair to conclude that this is an OpenBLAS issue versus an ITensor issue? I.e. I'm not sure what could be changed inside ITensor to work around a BLAS giving different results for various settings/flags, but the same mathematical inputs.

[mtfishman]

It looks like an OpenBLAS problem, and as you say there is probably not much we can do from our side.

@MaryamR99 as practical next steps for you if this is a blocking issue, can you use the AppleAccelerate.jl backend instead? If you have the time and motivation, you could come up with a more minimal demonstration of the issue and try to file a report of the problem to LinearAlgebra.jl and/or OpenBLAS. That's the only way this could potentially get fixed.

[MaryamR99]

Thanks for the clarification!

I will continue using AppleAccelerate.jl and try to find a more minimal example of the issue to report to OpenBLAS.

Thank you again.