Wrapping LAPACK functions using MKL

General Usage
1

The following implements a wrapper to the LAPACK functions dlanv2/slanv2:

import LinearAlgebra.BLAS.@blasfunc

using LinearAlgebra
import LinearAlgebra: libblastrampoline, BlasFloat, BlasReal, BlasComplex, BlasInt


const liblapack = Base.liblapack_name


for (fn, elty) in ((:dlanv2_, :Float64),
                   (:slanv2_, :Float32))
                   
    @eval begin
        function lanv2(A::$elty, B::$elty, C::$elty, D::$elty)
           """
           SUBROUTINE DLANV2( A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN )

           DOUBLE PRECISION A, B, C, CS, D, RT1I, RT1R, RT2I, RT2R, SN
           """
           RT1R = Ref{$elty}(1.0)
           RT1I = Ref{$elty}(1.0)
           RT2R = Ref{$elty}(1.0)
           RT2I = Ref{$elty}(1.0)
           CS = Ref{$elty}(1.0)
           SN = Ref{$elty}(1.0)
           ccall((@blasfunc($fn), liblapack), Cvoid,
                 (Ref{$elty},Ref{$elty},Ref{$elty},Ref{$elty},
                 Ref{$elty},Ref{$elty},Ref{$elty},Ref{$elty},Ref{$elty},Ref{$elty}),
                 A, B, C, D,
                 RT1R, RT1I, RT2R, RT2I, CS, SN)
           return RT1R[], RT1I[], RT2R[], RT2I[], CS[], SN[]
        end
    end
end
"""
    lanv2(A, B, C, D) -> (RT1R, RT1I, RT2R, RT2I, CS, SN)

Compute the Schur factorization of a real 2-by-2 nonsymmetric matrix `[A,B;C,D]` in
standard form. `A`, `B`, `C`, `D` are overwritten on output by the corresponding elements of the
standardised Schur form. `RT1R+im*RT1I` and `RT2R+im*RT2I` are the resulting eigenvalues.
`CS` and `SN` are the parameters of the rotation matrix.
Interface to the LAPACK subroutines DLANV2/SLANV2.
"""
lanv2(A::BlasReal, B::BlasReal, C::BlasReal, D::BlasReal)

I can call lanv2 and get the correct result

julia> lanv2(1.,2.,3.,4.)
(-0.3722813232690143, 0.0, 5.372281323269014, 0.0, -0.8245648401323938, 0.5657674649689923)

After loading MKL, I get an error message and wrong results:

julia> using MKL

julia> lanv2(1.,2.,3.,4.)
Error: no BLAS/LAPACK library loaded!
(1.0, 1.0, 1.0, 1.0, 1.0, 1.0)

I wonder what is wrong with my interface. It would be desirable that this routine works with both OpenBLAS and MKL.

2

In GitHub - DynareJulia/FastLapackInterface.jl we use

@static if VERSION < v"1.7"
    using LinearAlgebra.LAPACK: liblapack
elseif VERSION < v"1.9"
    const liblapack = "libblastrampoline"
else
    const liblapack = LinearAlgebra.libblastrampoline
end

and it works with both OpenBlas and MKL

3

Thanks for hinting to this. However, I already tried the above code, but I doesn’t work for me (for some reason). In both cases, I obtained

liblapack = "libblastrampoline-5.dll"

Is this the expected output?

4

Yes, I get something similar under Linux. It is the libblastrampoline that loads OpenBlas or MKL at run time. See GitHub - JuliaLinearAlgebra/libblastrampoline: Using PLT trampolines to provide a BLAS and LAPACK demuxing library.

5

I still don’t understand the exact situation. Could it be that no interface to the LAPACK auxiliary function dlanv2 (to compute the Schur form and eigenvalues of a 2x2 matrix) is provided in MKL or the interface is different from that provided by OpenBLAS?. I assume it is somehow possible to check this (perhaps can somebody execute the above tests in Linux!). I included the code with the proposed fix (not working for me), just in case somebody has time for a short check. Many thanks in advance.


import LinearAlgebra.BLAS.@blasfunc

using LinearAlgebra
import LinearAlgebra: libblastrampoline, BlasFloat, BlasReal, BlasComplex, BlasInt

using Base: iszero, has_offset_axes
using LinearAlgebra.LAPACK

#const liblapack = Base.liblapack_name

@static if VERSION < v"1.7"
    using LinearAlgebra.LAPACK: liblapack
elseif VERSION < v"1.9"
    const liblapack = "libblastrampoline"
else
    const liblapack = LinearAlgebra.libblastrampoline
end

for (fn, elty) in ((:dlanv2_, :Float64),
                   (:slanv2_, :Float32))
                   
    @eval begin
        function lanv2(A::$elty, B::$elty, C::$elty, D::$elty)
           """
           SUBROUTINE DLANV2( A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN )

           DOUBLE PRECISION A, B, C, CS, D, RT1I, RT1R, RT2I, RT2R, SN
           """
           RT1R = Ref{$elty}(1.0)
           RT1I = Ref{$elty}(1.0)
           RT2R = Ref{$elty}(1.0)
           RT2I = Ref{$elty}(1.0)
           CS = Ref{$elty}(1.0)
           SN = Ref{$elty}(1.0)
           ccall((@blasfunc($fn), liblapack), Cvoid,
                 (Ref{$elty},Ref{$elty},Ref{$elty},Ref{$elty},
                 Ref{$elty},Ref{$elty},Ref{$elty},Ref{$elty},Ref{$elty},Ref{$elty}),
                 A, B, C, D,
                 RT1R, RT1I, RT2R, RT2I, CS, SN)
           return RT1R[], RT1I[], RT2R[], RT2I[], CS[], SN[]
        end
    end
end
6

I can reproduce the problem in Linux. It seems specific to lanv2 as I can compute gqref with a similar code. However, the way lanv2 is packaged in MKL_jll seems different from gqref
If I look at the symbol in libmkl_rt.so.2, I get:

michel@21HFCTO1WW:~/.julia/artifacts/347e4bf25d69805922225ce6bf819ef0b8715426/lib$ nm libmkl_rt.so.2 |grep dgeqrf
0000000000454840 T dgeqrf
0000000000454840 T dgeqrf_
0000000000454910 T dgeqrf_64
0000000000454910 T dgeqrf_64_
0000000000454b00 T dgeqrfp
0000000000454b00 T dgeqrfp_
0000000000454bd0 T dgeqrfp_64
0000000000454bd0 T dgeqrfp_64_
0000000000508310 T lapacke_dgeqrf
0000000000508310 T lapacke_dgeqrf_
0000000000508310 T LAPACKE_dgeqrf
0000000000508270 T lapacke_dgeqrf_64
0000000000508270 T lapacke_dgeqrf_64_
0000000000508270 T LAPACKE_dgeqrf_64
0000000000508610 T lapacke_dgeqrfp
0000000000508610 T lapacke_dgeqrfp_
0000000000508610 T LAPACKE_dgeqrfp
0000000000508570 T lapacke_dgeqrfp_64
0000000000508570 T lapacke_dgeqrfp_64_
0000000000508570 T LAPACKE_dgeqrfp_64
0000000000508790 T lapacke_dgeqrfp_work
0000000000508790 T lapacke_dgeqrfp_work_
0000000000508790 T LAPACKE_dgeqrfp_work
00000000005086b0 T lapacke_dgeqrfp_work_64
00000000005086b0 T lapacke_dgeqrfp_work_64_
00000000005086b0 T LAPACKE_dgeqrfp_work_64
0000000000508490 T lapacke_dgeqrf_work
0000000000508490 T lapacke_dgeqrf_work_
0000000000508490 T LAPACKE_dgeqrf_work
00000000005083b0 T lapacke_dgeqrf_work_64
00000000005083b0 T lapacke_dgeqrf_work_64_
00000000005083b0 T LAPACKE_dgeqrf_work_64
00000000004549e0 T mkl_dgeqrf_compact
00000000004549e0 T mkl_dgeqrf_compact_
0000000000454840 T mkl_lapack__dgeqrf_
0000000000454910 T mkl_lapack__dgeqrf_64_
00000000006b89c0 T mkl_lapack_dgeqrf_omp_offload_ilp64
00000000006b8ac0 T mkl_lapack_dgeqrf_omp_offload_ilp64_
00000000006b8dc0 T mkl_lapack_dgeqrf_omp_offload_lp64
00000000006b8ec0 T mkl_lapack_dgeqrf_omp_offload_lp64_
0000000000454b00 T mkl_lapack__dgeqrfp_
0000000000454bd0 T mkl_lapack__dgeqrfp_64_
00000000004549e0 T mkl_lapack__mkl_dgeqrf_compact_
michel@21HFCTO1WW:~/.julia/artifacts/347e4bf25d69805922225ce6bf819ef0b8715426/lib$ nm libmkl_rt.so.2 |grep dlanv2
0000000000466fe0 T dlanv2
0000000000466fe0 T dlanv2_
0000000000466fe0 T mkl_lapack__dlanv2_

I’m out of my depth here. Maybe you should open an issue on MKL.jl

7

Thank you for your time. The problem of non-uniform interface in MKL has been already addressed in this issue and a PR is underway to fix it JuliaLinearAlgebra/MKL.jl#140. (not yet merged).

8

:+1: