Writing custom LAPACK bindings

New to Julia
1

Hi,
I’m relatively new to Julia. I’ve been using it to prototype some linear algebra methods. I now have need for some LAPACK methods which are not included in LinearAlgebra.LAPACK, (specifically dtpqrt and dtpmqrt). I think I must be missing something basic, because even just copying an existing binding isn’t working.

Here is a minimal example, which I’ve tried in both 1.6 and 1.7:

using LinearAlgebra
using LinearAlgebra: BlasInt, require_one_based_indexing
using LinearAlgebra.LAPACK: liblapack, chkstride1, chklapackerror
using LinearAlgebra.BLAS: @blasfunc

function geqrt!(A::AbstractMatrix{Float64}, T::AbstractMatrix{Float64})
    require_one_based_indexing(A, T)
    chkstride1(A)
    m, n = size(A)
    minmn = min(m, n)
    nb = size(T, 1)
    if nb > minmn
        throw(ArgumentError("block size $nb > $minmn too large"))
    end
    lda = max(1, stride(A,2))
    work = Vector{Float64}(undef, nb*n)
    if n > 0
        info = Ref{BlasInt}()
        ccall((@blasfunc(:dgeqrt_), liblapack), Cvoid,
            (Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ptr{Float64},
                Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}, Ptr{Float64},
                Ptr{BlasInt}),
                m, n, nb, A,
                lda, T, max(1,stride(T,2)), work,
                info)
        chklapackerror(info[])
    end
    A, T
end

Calling this gives

ERROR: could not load symbol ":dgeqrt_64_":
dlsym(0x7fb7a1d04d50, :dgeqrt_64_): symbol not found

Thanks in advance for any assistance.

2

I’m not sure if I fully understand what you want, but what is the problem with just LinearAlgebra.LAPACK.geqrt!(A, T)?

3

You’ve stumbled on a subtlety of the way symbols are treated in macro expressions. Instead of @blasfunc(:dgeqrt_) you want @blasfunc(dgeqrt_).

2 Likes
4

Ah, that was it. Thank you!

For a bit more context, I’m implementing a method where I need to update an existing QR factorization by either removing a few columns, or permuting a block of columns on the left to the right, so update Q R = [A1, A2] to Q' R' = [A2, A1] (' isn’t transpose here, just denotes the modified factorization). tpqrt can be used to do this more cheaply than computing the new factorization from scratch.