Add syevd LAPACK Hermitian eigensolver

[jaemolihm]

Added syevd! wrapper to LAPACK complex Hermitian eigensolvers zheevd and cheevd.

It is more accurate than the existing wrappers (syev!, syevr), and is also the fastest for large matrices (N > 100) with MKL.

https://netlib.org/lapack/lawnspdf/lawn183.pdf

QR [SYEV] and DC [SYEVD] are the most accurate algorithms, measured both in terms of producing pairwise orthogonal eigenvectors and small residuals norms. MR [SYEVR] is less accurate but still achieves errors of size O(nε), where n is the dimension and ε is machine epsilon

This is useful e.g. in JuliaMolSim/DFTK.jl#842

Benchmark (MKL, N = 100)

Method	Timing (ms)	Allocation (KiB)	Orthogonality error
syev!	1.808	36.36	2.3e-14
syevr!	1.577	394.62	1.5e-13
syevd!	1.641	373.23	1.9e-14

Benchmark (MKL, N = 1000)

Method	Timing (ms)	Allocation (MiB)	Orthogonality error
syev!	793	0.5	2.4e-13
syevr!	547	31.3	5.2e-12
syevd!	504	31.1	9.9e-14

Benchmark code

```julia using LinearAlgebra using BenchmarkTools using MKL

function test_syev(h, tmp, jobz='V')
tmp .= h
LinearAlgebra.LAPACK.syev!(jobz, 'U', tmp)
end
function test_syevr(h, tmp, jobz='V')
tmp .= h
LinearAlgebra.LAPACK.syevr!(jobz, tmp)
end
function test_syevd(h, tmp, jobz='V')
tmp .= h
LinearAlgebra.LAPACK.syevd!(jobz, 'U', tmp)
end

N = 1000
T = Float64
h = rand(Complex{T}, N, N);
h .+= h';
tmp = zeros(Complex{T}, N, N);

@Btime test_syev($h, $tmp);
@Btime test_syevr($h, $tmp);
@Btime test_syevd($h, $tmp);

e, v = LinearAlgebra.LAPACK.syev!('V', 'U', copy(h));
@info norm(v' * v - I)
e, v = LinearAlgebra.LAPACK.syevr!('V', copy(h));
@info norm(v' * v - I)
e, v = LinearAlgebra.LAPACK.syevd!('V', 'U', copy(h));
@info norm(v' * v - I)

</details>

[dkarrasch]

Cool. Are there cases where we want to dispatch to this method by default, or do we just want to provide this?

[ViralBShah]

@jaemolihm The tests introduced for this PR are failing.

[jaemolihm]

Thank you @ViralBShah , I fixed the test by adding the real symmetric eigensolver interface.

Cool. Are there cases where we want to dispatch to this method by default, or do we just want to provide this?

The motivation is just to provide this. Choosing the default is more difficult because the timing also depends on the library used (and possibly on the hardware).

[ViralBShah]

Is it worth considering providing an option to the eigensolver interface to use this? Could be in a separate PR.

[jaemolihm]

Is it worth considering providing an option to the eigensolver interface to use this? Could be in a separate PR.

Yes, it seems useful. Maybe something like below?

using LinearAlgebra.LAPACK: syev!, syevd!, syevr!
using LinearAlgebra: RealHermSymComplexHerm, eigencopy_oftype, BlasReal, Eigen, sorteig!, LAPACK, eigtype
import LinearAlgebra: eigen!, eigen

function eigen(A::RealHermSymComplexHerm, method::Symbol; sortby::Union{Function,Nothing}=nothing)
    S = eigtype(eltype(A))
    eigen!(eigencopy_oftype(A, S), method; sortby)
end

function eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, method::Symbol; sortby::Union{Function,Nothing}=nothing)
    if method == :syevr
      Eigen(sorteig!(LAPACK.syevr!('V', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)..., sortby)...)
    elseif method == :syev
      Eigen(sorteig!(LAPACK.syev!('V', A.uplo, A.data)..., sortby)...)
    elseif method == :syevd
      Eigen(sorteig!(LAPACK.syevd!('V', A.uplo, A.data)..., sortby)...)
    else
      throw(ArgumentError("Wrong method $method. Must be :syevr or :syev or :syevd."))
    end
end

[ViralBShah]

Yes, something along those lines, but perhaps with slightly more descriptive names?

[ViralBShah]

Just noticed this typo in performance after merging.

[jaemolihm]

I fixed this typo in the follow up PR #49355.