julia - 直接从 julia 调用 BLAS 的问题(找不到函数 :zgemm_64_ in library libopenblas64_)

Question

我正在尝试使用 ccall 在 Julia 中调用 BLAS像这样

ccall((BLAS.@blasfunc(:zgemm_), BLAS.libblas),...other arguments)

据我所知，这与 LinearAlgebra 包调用 BLAS ( link to source ) 的方式相同

但是，我收到以下错误:

ccall: could not find function :zgemm_64_ in library libopenblas64_

任何人都知道可能是什么问题？

编辑:发现使用 :zgemm_64_直接代替 BLAS.@blasfunc(:zgemm_)解决了错误，但我仍然想知道为什么。

如果有必要，这里是我进行 BLAS 调用的完整函数。

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)

Answer 1

BLAS.@blasfunc(:zgemm_)返回 Symbol(":zgemm_64_") ，而不是 :zgemm_64_ ，首先看起来很奇怪......从技术意义上讲它很卫生，但无可否认地令人困惑。它在原始实现中起作用的原因是因为在那里，带有名称的符号总是拼接成 @eval ;相比:

julia> @eval begin
           BLAS.@blasfunc(:zgemm_)
       end
Symbol(":zgemm_64_")

julia> @eval begin
           BLAS.@blasfunc($(:zgemm_))
       end
:zgemm_64_

所以， @blasfunc期望它的参数是一个名称(即 AST 中的符号)，而不是符号文字(AST 中的带引号的符号)。你可以等效地把它写成一个变量名:

julia> @eval begin
           BLAS.@blasfunc zgemm_
       end
:zgemm_64_

(没有在此范围内实际定义 zgemm_!)