4

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)
phipsgabler
  • 20,535
  • 4
  • 40
  • 60
Thijs Steel
  • 1,190
  • 7
  • 16

1 Answers1

4

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

phipsgabler
  • 20,535
  • 4
  • 40
  • 60
  • I mean -- if it's only used inside `@eval`, it's actually less confusing than the alternative of having to splice `$(Meta.quot(:zgemm_))`. – phipsgabler Jan 10 '20 at 16:31