I'm trying to call BLAS in Julia using ccall
like this
ccall((BLAS.@blasfunc(:zgemm_), BLAS.libblas),...other arguments)
For as far as I can tell, this is the same way the LinearAlgebra package calls BLAS (link to source)
I get the following error however:
ccall: could not find function :zgemm_64_ in library libopenblas64_
Anyone have any idea what could be the problem?
EDIT: found out that using :zgemm_64_
directly instead of BLAS.@blasfunc(:zgemm_)
solved the error, but I'd still like to know why.
In case it becomes necessary, here is the full function where I make the BLAS call.
import LinearAlgebra: norm, lmul!, rmul!, BlasInt, BLAS
# Preallocated version of A = A*B
function rmul!(
A::AbstractMatrix{T},
B::AbstractMatrix{T},
workspace::AbstractVector{T}
) where {T<:Number}
m,n,lw = size(A,1), size(B,2), length(workspace)
if(size(A,2) !== size(B,1))
throw(DimensionMismatch("dimensions of A and B don't match"))
end
if(size(B,1) !== n)
throw(DimensionMismatch("A must be square"))
end
if(lw < m*n)
throw(DimensionMismatch("provided workspace is too small"))
end
# Multiplication via direct blas call
ccall((BLAS.@blasfunc(:zgemm_), BLAS.libblas), Cvoid,
(Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
Ref{BlasInt}, Ref{T}, Ptr{T}, Ref{BlasInt},
Ptr{T}, Ref{BlasInt}, Ref{T}, Ptr{T},
Ref{BlasInt}),
'N', 'N', m, n,n, 1.0, A, max(1,stride(A,2)),B, max(1,stride(B,2)), 0.0, workspace, n)
# Copy temp to A
for j=1:n
for i=1:m
A[i,j] = workspace[j+*(i-1)*n]
end
end
end
function test_rmul(m::Integer, n::Integer)
BLAS.set_num_threads(1)
A = rand(ComplexF64, m,n)
Q = rand(ComplexF64, n,n)
workspace = similar(A, m*n)
A_original = copy(A)
Q_original = copy(Q)
rmul!(A,Q,workspace)
@show norm(A_original*Q_original - A)
@show norm(Q_original - Q)
end
test_rmul(100,50)
BLAS.@blasfunc(:zgemm_)
returns Symbol(":zgemm_64_")
, and not :zgemm_64_
, which looks rather strange in the first place... it's hygienic in the technical sense, but admittedly confusing. The reason it works in the original implementation is because there, the symbol with the name is always spliced into @eval
; compare:
julia> @eval begin
BLAS.@blasfunc(:zgemm_)
end
Symbol(":zgemm_64_")
julia> @eval begin
BLAS.@blasfunc($(:zgemm_))
end
:zgemm_64_
So, @blasfunc
expects its argument to be a name (i.e., a symbol in the AST), not a symbol literal (a quoted symbol in the AST). You could equivalently write it like a variable name:
julia> @eval begin
BLAS.@blasfunc zgemm_
end
:zgemm_64_
(without zgemm_
being actually defined in this scope!)