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