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

ITensor Julia Questions
1

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.

2

Does it work with another backend, like OpenBLAS?

3

Also, how far off is it? What is the error?

4

The OpenBLAS works fine, the norm given by OpenBLAS is 1.2327681181461817e-12. The MKL gives 281.4977875674409.

5

Interesting, thanks for the update. I’m not sure if there is much we can do about that since that sounds like a bug in MKL. You could convert that ITensor to a Julia Matrix, save it in a transferable file format (maybe with HDF5.jl or as a delimited file), and report it as a bug to Julia’s MKL wrapper (GitHub - JuliaLinearAlgebra/MKL.jl: Intel MKL linear algebra backend for Julia) or raise it as an issue on Intel’s user forum: Intel® oneAPI Math Kernel Library - Intel Community.

You could also try other versions of MKL as described here: GitHub - JuliaLinearAlgebra/MKL.jl: Intel MKL linear algebra backend for Julia to see if it is an issue particular to a certain version of MKL.

Finally, could you share the output of versioninfo() here so I can see what system and Julia version you are using?

6

Does it change anything if you set ishermitian=false? That would call a different LAPACK routine so it would at least narrow down the issue a little bit.

Also, after you convert the ITensor to a Julia Matrix, you could try calling LinearAlgebra.eigen directly on that Julia Matrix (either wrapping it in Symmetric/Hermitian or not wrapping it to try out different LAPACK eigendecomposition backends), just to rule out the possibility that it is somehow related to the ITensor code wrapping the eigendecomposition.