Performance of Hermitian eigensolvers

I usually work with full diagonalisation of Hermitian matrices. After some research I realised that there are two options to retrieve the eigenvalues and eigenvectors of such a matrix.
LinearAlgebra.LAPACK.syev!() and eigen!(). From the source code of eigen!() it seems that it is calling a general eigensolver routine called ggevx().

The LinearAlgebra.LAPACK.syev!() routine required some digging into the source code to figure out that it is also intended for Hermitian matrices ( only talks about symmetric matrices). I realised that it is calling specialised Hermitian routines in LAPACK.

Now after building up a 2000x2000 Hermitian matrix I called both routines the following way and timed their peformance:

H = weyl(10, 10, 10)
Ht = deepcopy(H)

@time eigsys = LAPACK.syev!('V', 'U', H)
@time eigsyst = eigen!(Hermitian(Ht))

To my surprise the eigen! routine is quite a bit faster:

  8.998029 seconds (6 allocations: 1.068 MiB)
  6.128692 seconds (16 allocations: 123.658 MiB, 1.71% gc time)

To me this is quite counterintuitive behaviour. I would have guessed that the specialised routines should work faster.

If you look more closely at LinearAlgebra/src/symmetric.jl you find that eigen!(::Hermitian) calls the syevr routines from LAPACK (note the trailing r). For full diagonalization these use Dhillon and Parlett’s “Relatively Robust Representations” which are more efficient than the traditional syev scheme.