Проблема вызова BLAS напрямую от Юлии (не удалось найти функцию: zgemm_64_ в библиотеке libopenblas64_) - PullRequest
Новая
       36

Проблема вызова BLAS напрямую от Юлии (не удалось найти функцию: zgemm_64_ в библиотеке libopenblas64_)

4 голосов
/ 10 января 2020

Я пытаюсь вызвать BLAS в Julia, используя ccall, например, 

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

Насколько я могу судить, это так же, как пакет LinearAlgebra вызывает BLAS ( ссылка к источнику )

Я получаю следующую ошибку, однако: 

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)

1 Ответ

4 голосов
/ 10 января 2020

BLAS.@blasfunc(:zgemm_) возвращает Symbol(":zgemm_64_"), а не :zgemm_64_, что выглядит довольно странно, во-первых ... это гигиенично c в техническом смысле, но, по общему признанию, сбивает с толку. Причина, по которой он работает в оригинальной реализации, заключается в том, что там символ с именем всегда вставляется в @eval; Сравните: 

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

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

Таким образом, @blasfunc ожидает, что его аргумент будет имя (т. е. символ в AST), а не символьный литерал ( в кавычках *) 1012 * символ в AST). Вы могли бы эквивалентно написать его как имя переменной: 

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

(без фактического определения zgemm_ в этой области!)

...