`syevd!` is reporting more errors than `syevr!`

[VinceNeede]

syevd! chosen as default algorithm for juliav v1.12.0+beta is reporting more errors respect to syevr!, creating
problems when a large number of eigen decompositions must be performed.

The same matrix can give LAPACKError (with info>0) on a device and pass on another and viceversa, furthermore the same matrix can give error or not based on the number of threads used (what I mean is that a specific matrix raises an error, for example, with 1 thread and not for 2, while another matrix can give an error for a different number of threads), meaning that some instability is present (at least in OpenBLAS, MKL seems to be working fine).

test code with matrix saved in HDF5

The matrix is of medium size $658\times 658$ and resulted from a long computation, so it is stored in a `HDF5` file, if you tell me how can I attach the file for the debug I will send it to you. The code is just the following:

using HDF5
using LinearAlgebra

mat = nothing
h5open("mat_error.h5", "r") do file
	global mat = Hermitian(file["mat"][])
end

@info "exectuing with 2 threads"
BLAS.set_num_threads(2)
eigen(mat)

@info "exectuing with 1 thread"
BLAS.set_num_threads(1)
eigen(mat)

Error obtained

The output is

[ Info: exectuing with 2 threads
[ Info: exectuing with 1 thread
ERROR: LoadError: LAPACKException(326121)
Stacktrace:
  [1] chklapackerror_positive(::Int64)
    @ LinearAlgebra.LAPACK ~/.julia/juliaup/julia-1.12.0-beta1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/lapack.jl:44
  [2] chklapackerror
    @ ~/.julia/juliaup/julia-1.12.0-beta1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/lapack.jl:40 [inlined]
  [3] _syevd!(jobz::Char, uplo::Char, A::Matrix{Float64})
    @ LinearAlgebra.LAPACK ~/.julia/juliaup/julia-1.12.0-beta1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/lapack.jl:5463
  [4] syevd!
    @ ~/.julia/juliaup/julia-1.12.0-beta1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/lapack.jl:5324 [inlined]
  [5] #eigen!#251
    @ ~/.julia/juliaup/julia-1.12.0-beta1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/symmetriceigen.jl:14 [inlined]
  [6] eigen!
    @ ~/.julia/juliaup/julia-1.12.0-beta1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/symmetriceigen.jl:12 [inlined]
  [7] #_eigen#253
    @ ~/.julia/juliaup/julia-1.12.0-beta1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/symmetriceigen.jl:55 [inlined]
  [8] _eigen
    @ ~/.julia/juliaup/julia-1.12.0-beta1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/symmetriceigen.jl:53 [inlined]
  [9] #eigen#252
    @ ~/.julia/juliaup/julia-1.12.0-beta1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/symmetriceigen.jl:49 [inlined]
 [10] eigen
    @ ~/.julia/juliaup/julia-1.12.0-beta1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/symmetriceigen.jl:48 [inlined]
 [11] eigen(A::Hermitian{Float64, Matrix{Float64}})
    @ LinearAlgebra ~/.julia/juliaup/julia-1.12.0-beta1+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/symmetriceigen.jl:48
 [12] top-level scope
    @ ~/Tesi/lapack_error/study_from_mat.jl:15
 [13] include(mapexpr::Function, mod::Module, _path::String)
    @ Base ./Base.jl:304
 [14] top-level scope
    @ REPL[2]:1
in expression starting at /home/bisogno/Tesi/lapack_error/study_from_mat.jl:15

[aravindh-krishnamoorthy]

The matrix is of medium size 658× 658 and resulted from a long computation, so it is stored in a HDF5 file, if you tell me how can I attach the file for the debug I will send it to you.

Please just drag and drop the file into the "Add a comment" text box or just click on the "Paste, drop, or click to add files" with a clip icon just at the bottom left of the "Add a comment" text box. Will be nice to have this file to see where exactly the error is originating from...

[VinceNeede]

mat_error.csv I have converted it to a csv file (I couldn't upload hdf5)

modified import

using CSV, DataFrames
using LinearAlgebra

mat = Hermitian(Matrix(CSV.read("mat_error.csv", DataFrame)))

@info "exectuing with 2 threads"
BLAS.set_num_threads(2)
eigen(mat)

@info "exectuing with 1 thread"
BLAS.set_num_threads(1)
eigen(mat)

As previously said, it is unlikely that this will reproduce the error on a different machine since a matrix can give an error on machine A and not on B while another matrix can pass on machine A and give error on B.

This is the hardware on which this matrix gives error:

julia> versioninfo()
Julia Version 1.12.0-beta2
Commit dd74040f22b (2025-04-25 12:03 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 32 × Intel(R) Xeon(R) CPU E5-2670 0 @ 2.60GHz
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, sandybridge)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 32 virtual cores)

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

I also wanted to add that I haven't encountered the error with MKL

[araujoms]

Can't reproduce. Also running 1.12 beta 2 on an Intel chip on Linux x64.

[VinceNeede]

Yes both the machines I used have intel chips. I'm trying to find a way to systematically create a matrix that gives error, but I'm not positive on succeeding.

Still I opened the issue not because I was searching for a solution itself (I'm not sure it even exists given the strong device dependent behavior) but more for reporting my personal experience with DivideAndConquer as default algorithm, since I think this is the meaning in using pre-release versions.

[araujoms]

Well if syevd! gives more errors than syevr! that's a good argument to change back the default, but this needs to be quantifiable and reproducible. Have you done statistics for both?

[VinceNeede]

My code is a montecarlo on mps evolution (quantum jumps) using ITensors, which internally computes eigen decomposition to perform the MPO * MPS operation.

It worked fine on julia v1.11.5, when I changed to the 1.12 beta I started to get this error, which I was able to completely bypass by either setting LinearAlgebra.default_eigen_algorithm() = RobustRepresentation() or moving to MKL for which the error never happened.

To give an idea of the frequency, the mps is composed by $L=100$ sites (which is the number of eigen decompositions performed per step), the state is evolved for around $T=20$ steps and I would get the error on around 100-200 trajectories. So we are talking of a frequency about $1/200k$ which is REALLY rare. But still in computations long enough it can happen, and it is quite unhappy to send a job on a server for a couple of days and find out it stopped for an error after the first day.

I'm sorry I'm not being able to give a reproducible way to get it

[araujoms]

Ok, so you have no errors with syevr!, but an error rate of 5e-6 with syevd!? I think that's a good enough argument to switch back to syevr!. It is in fact supposed to be more reliable.

I'd be happier if I could reproduce your statistics, though. Is your MPS code public? I could leave it running on my machine for a couple of days to see if I get similar error rates (even if for different matrices).

[mtfishman]

Thanks for raising this issue @VinceNeede. I'm one of the ITensor developers, @VinceNeede raised this issue to us originally, I would vote for switching the default back to syevr since it appears to be more reliable (assuming the issue can be reproduced).

[ViralBShah]

Interesting. IIRC, syevr is actually a much more recent implementation. I always thought that syevd was the old reliable battle tested version in LAPACK, but in Julia we had used syevr for a long time until recently. I can't find the issue/PR where we switched and why we did.

cc @antoine-levitt

[ViralBShah]

If we are going to change, we should make sure it makes it into 1.12.

[araujoms]

It was changed in JuliaLang/julia#49355. The reason was that syevd was slightly faster than syevr when using Intel MKL.

[ViralBShah]

So the change requested here is to make RobustRepresentations() the default rather than DivideAndConquer()?