JuliaLinearAlgebra / MKL.jl

Intel MKL linear algebra backend for Julia
Other
205 stars 32 forks source link

Problems when using the MKL library to perform eigenvalue decomposition on the transverse field ising model #170

Closed Huangzhiqiang closed 2 months ago

Huangzhiqiang commented 2 months ago

As the title says, when running the following code in Julia:

using ITensors, ITensorMPS, MKL N = 12 g = 1.05 sites = siteinds("S=1/2", N) os = OpSum() for j=1:N os += -4,"Sz",j,"Sz",j%N+1 os += 2*g,"Sx",j; end H = MPO(os, sites) ITensors.disable_warn_order() H_full = contract(H); evt, Vt = eigen(H_full; ishermitian=true);

The program will give the wrong decomposition result, and norm(H_full-Vt’evtVt) !≈0. I repeated the above error in computers using AMD 7840hs and 5900x. After I convert the ITensor to a Julia Matrix, and try calling LinearAlgebra.eigen directly on that Julia Matrix. The result is still wrong.

ViralBShah commented 2 months ago

Try disabling MKL threading and see if that is better. Not sure what we can do about MKL doing things incorrectly. If there is a minimum reproducer, we can send it to Intel.

Huangzhiqiang commented 2 months ago

Thank you for your reply, the problem indeed disappeared when closing the MKL thread with BLAS.set_num_threads(1).

ViralBShah commented 2 months ago

Similar to https://github.com/JuliaLinearAlgebra/MKL.jl/issues/129

cc @mkrainiuk. Closing since there seems to be some issue with threading in MKL upstream.