Search code examples
juliametaprogrammingopenblas

Problem calling BLAS from julia directly (could not find function :zgemm_64_ in library libopenblas64_)


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)

Solution

  • 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!)