Call LAPACK's ZHEEV directly

General Usage
1

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)

2

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
3

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?
4

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.

5

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).

6

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

1 Like
7

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
8

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.

Thanks for your answers!

1 Like
9

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.