Search code examples
julialinear-algebralapackeigenvalue

Implementing LAPACK routine hseqr() in Julia (creating wrapper error: no method matching Array{Float64,1}(::Int64))


I am trying to implement the hseqr() routine from LAPACK(link:http://www.netlib.org/lapack/explore-html/da/dba/group__double_o_t_h_e_rcomputational_gacb35e85b362ce8ccf9d653cc3f8fb89c.html#gacb35e85b362ce8ccf9d653cc3f8fb89c) to Julia but I get this error when I call the function:

error: no method matching Array{Float64,1}(::Int64)

This is my code:

for (hseqr, elty) in
    ((:dhseqr_,:Float64),
    (:shseqr_,:Float32))
    @eval begin
        """
     JOB    
          JOB is CHARACTER*1
           = 'E':  compute eigenvalues only;
           = 'S':  compute eigenvalues and the Schur form T.
    COMPZ   
          COMPZ is CHARACTER*1
           = 'N':  no Schur vectors are computed;
           = 'I':  Z is initialized to the unit matrix and the matrix Z
                   of Schur vectors of H is returned;
           = 'V':  Z must contain an orthogonal matrix Q on entry, and
                   the product Q*Z is returned.
        """
        function hseqr!(job::Char,compz::Char,ilo::Integer,ihi::Integer,H::StridedMatrix{$elty}, Z::StridedMatrix{$elty})
            N=size(H,1)
            ldh=size(H,1)
            ldz=size(Z,1)
            work = Vector{$elty}(1)
            lwork = BlasInt(-1)
            info = Ref{BlasInt}()
            for i = 1:2  # first call returns lwork as work[1]
                ccall((@blasfunc($hseqr), liblapack),Void,
                    (Ref{UInt8},Ref{UInt8}, Ref{BlasInt},Ref{BlasInt},Ref{BlasInt},
                        Ptr{$elty}, Ref{BlasInt},Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, 
                        Ref{BlasInt}, Ptr{BlasInt}, Ref{BlasInt}, Ref{BlasInt}),
                        job,compz, N,ilo,ihi,H,ldh,wr,wi,Z,ldz,work, 
                        lwork,info) 
                chklapackerror(info[])
                if i == 1
                    lwork = BlasInt(real(work[1]))
                    resize!(work, lwork)
                end
            end
            return wr,wi,Z
        end
    end
end
hseqr!(job::Char,compz::Char,H::StridedMatrix{},Z::StridedMatrix{}) = hseqr!(job,compz,1, size(H, 1), H,Z)

This is my call: hseqr!('E','N',H,Matrix{Float64}(I,5,5)) (H is Hessenberg matrix with dimension 5).

I am not sure if I understand correctly how to make wrappers, so any hint helps.


Solution

  • This is a great example for a question that could benefit from being stripped down quite a bit to isolate the issue.

    If I read your code correctly, the line

    work = Vector{$elty}(1)
    

    is called with Float64 as element type, so the whole thing boils down to:

    julia> Vector{Float64}(1)
    ERROR: MethodError: no method matching Array{Float64,1}(::Int64)
    
    

    You are misusing the constructor for vectors. Now I'm not sure what you're trying to do here, but if you are looking to create a vector of Float64 numbers of length one you might be looking for one of these:

    julia> Array{Float64}(undef, 1)
    1-element Array{Float64,1}:
     6.9274314537094e-310
    
    julia> zeros(1)
    1-element Array{Float64,1}:
     0.0
    

    Note that the first call leaves the array uninitialised, i.e. it is filled with whatever stuff was in memory before so you want to be sure to only use this if you're actually overwriting it later on.