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

[Huangzhiqiang]

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]

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]

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

[ViralBShah]

Similar to #129

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