# Call LAPACK's ZHEEV directly

I know that LAPACK wrappers can be found in `Base.LinAlg.LAPACK`. However, for example ZHEEV is missing there. I found zheesv in `lapack.jl`. There, `zheesv!()` is somehow created in an `@eval for` loop that I do not really understand.

How (if possible) do I call ZHEESV from within Julia in an easy manner? I’m really interested in something easy like `Base.LinAlg.LAPACK.zheev` and not manual `ccall`. Also, I am aware that `eigfact` and co are probably already dispatching to the best LAPACK functions for given input. I still want to know how to call the LAPACK routine directly.

(Of course, it would also be really nice if someone could enlighten me what that `@eval for` construction in `lapack.jl` really does)

Ok, AFAIU no `zheesv!()` is created, but instead methods `syev!(jobz::Char, uplo::Char, A::StridedMatrix{elty})` with different `elty`. Depending on `elty` the respective LAPACK routine is called via `ccall` within each of those methods. From this I would conclude that calling LAPACK’s ZHEEV corresponds (via dispatch) to calling `Base.LinAlg.LAPACK.syev!` with complex matrix `A`.

However, then it seems that the documentation isn’t really precise:

``````help?> Base.LinAlg.LAPACK.syev!
syev!(jobz, uplo, A)

Finds the eigenvalues (jobz = N) or eigenvalues and eigenvectors (jobz = V) of a symmetric matrix A. If uplo = U, the upper
triangle of A is used. If uplo = L, the lower triangle of A is used.

julia>
``````

It doesn’t mention that it can also be used to find the eigenvalues of a complex Hermitian matrix.

1 Like

So, this works (calling with usual complex array):

``````julia> H = convert.(Complex128, Hermitian(rand(3,3)))
3×3 Array{Complex{Float64},2}:
0.0816049+0.0im  0.414275+0.0im  0.598706+0.0im
0.414275+0.0im  0.512361+0.0im  0.295375+0.0im
0.598706+0.0im  0.295375+0.0im  0.995971+0.0im

julia> Base.LinAlg.LAPACK.syev!('N', 'L', H)
3-element Array{Float64,1}:
-0.298891
0.386473
1.50236
``````

However, this (unsuprisingly but maybe also unconveniently) doesn’t (calling with a `::Hermitian`):

``````julia> H = Hermitian(convert.(Complex128, Hermitian(rand(3,3))))
3×3 Hermitian{Complex{Float64},Array{Complex{Float64},2}}:
0.122402-0.0im   0.46566+0.0im  0.465734+0.0im
0.46566-0.0im  0.296259-0.0im  0.090323+0.0im
0.465734-0.0im  0.090323-0.0im  0.722005-0.0im

julia> Base.LinAlg.LAPACK.syev!('N', 'L', H)
ERROR: MethodError: no method matching syev!(::Char, ::Char, ::Hermitian{Complex{Float64},Array{Complex{Float64},2}})
Closest candidates are:
syev!(::Char, ::Char, ::Union{Base.ReshapedArray{Float64,2,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{DenseArray, SubArray{T,N,P,I,true} where I<:Tuple{Union{Base.Slice, UnitRange},Vararg{Any,N} where N} where P where N where T}, DenseArray{Float64,2}, SubArray{Float64,2,A,I,L} where L} where I<:Tuple{Vararg{Union{Base.AbstractCartesianIndex, Int64, Range{Int64}},N} where N} where A<:Union{Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{DenseArray, SubArray{T,N,P,I,true} where I<:Tuple{Union{Base.Slice, UnitRange},Vararg{Any,N} where N} where P where N where T} where N where T, DenseArray}) at linalg\lapack.jl:4785
syev!(::Char, ::Char, ::Union{Base.ReshapedArray{Float32,2,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{DenseArray, SubArray{T,N,P,I,true} where I<:Tuple{Union{Base.Slice, UnitRange},Vararg{Any,N} where N} where P where N where T}, DenseArray{Float32,2}, SubArray{Float32,2,A,I,L} where L} where I<:Tuple{Vararg{Union{Base.AbstractCartesianIndex, Int64, Range{Int64}},N} where N} where A<:Union{Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{DenseArray, SubArray{T,N,P,I,true} where I<:Tuple{Union{Base.Slice, UnitRange},Vararg{Any,N} where N} where P where N where T} where N where T, DenseArray}) at linalg\lapack.jl:4785
syev!(::Char, ::Char, ::Union{Base.ReshapedArray{Complex{Float64},2,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{DenseArray, SubArray{T,N,P,I,true} where I<:Tuple{Union{Base.Slice, UnitRange},Vararg{Any,N} where N} where P where N where T}, DenseArray{Complex{Float64},2}, SubArray{Complex{Float64},2,A,I,L} where L} where I<:Tuple{Vararg{Union{Base.AbstractCartesianIndex, Int64, Range{Int64}},N} where N} where A<:Union{Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{DenseArray, SubArray{T,N,P,I,true} where I<:Tuple{Union{Base.Slice, UnitRange},Vararg{Any,N} where N} where P where N where T} where N where T, DenseArray}) at linalg\lapack.jl:4929
...

``````
• Is there a way to test if ZHEEV is actually called?
• Should the documentation be adjusted? (A Hermitian matrix is, of course, also symmetric, so the doc isn’t wrong but maybe misleading)
• Should the second example also work (i.e. add more methods)?
• Why not make ZHEEV and others explicitly available as `Base.LinAlg.LAPACK.zheev`?

In Julia, it is more natural and convenient to use dispatch on the argument types rather than having multiple function names ala Fortran 77.

Moreover, it’s not clear to me what purpose is served by making easier to call low-level LAPACK internals if they don’t have any functionality not available in `eigfact`.

In Julia, it is more natural and convenient to use dispatch on the argument types rather than having multiple function names ala Fortran 77.

Sure, but `Base.LinAlg.LAPACK` provides wrappers to precisely those old Fortran routines and it might by highly unintuitive to not have those methods in there explicitly but hidden in dispatch logic (although more natural to Julia). If I search for ZHEEV I will probably not assume that it might be a differently dispatched version of SYEV which is only for real matrices in Fortran.

However, I agree that in general there shouldn’t really be a reason to call those functions. Nonetheless, I have to call `Base.LinAlg.LAPACK.gesvd!` explicitly in one of my codes, because `svdfact` is based on `Base.LinAlg.LAPACK.gesdd!` which is unstable in this particular usecase. Moreover, people switching to Julia might like having direct access to routines they are used to (but this is of course a thin argument).

The right thing to do would be to expose that option in `svdfact` if it is needed in some cases.

1 Like

Certainly, the LAPACK documentation is rather lacking. In particular, we should explicitly mention in the docs which routines are called (we could also link to the LAPACK docs).

This is a bit harder to say. My personal interpretation of the LAPACK module is that it is basically LAPACK + error handling and dispatch on matrix element types.

All the other functions drop the element type prefix (since we use dispatch), so it would be `LAPACK.heev!`. I agree that there is an inconsistency here.

2 Likes

So I guess that is the way to go for now. If Julia does not provide more or less one-to-one mappings to the LAPACK functions, which is probably fine/preferable as per @stevengj’s comment, then we should at least make the documentation of the current state clearer. As I just noticed, the docs at least indicate dispatching with one sentence.

If I find the time, I will prepare a PR that tries to clarify the docstrings of the wrapper functions.

LAPACK has two set of eigensolver methods, `[sd]syev` for real-Hermitian (i.e. real-symmetric) problems and `[cz]heev` for complex-Hermitian problems. From Julia’s perspective, these are all the “same” thing, just dispatching on the method type, so they should all have the same name, either `syev!` or `heev!`. We had to pick one, so the former was chosen.