Wrapping LAPACK functions using MKL

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.

``````@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

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?

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.

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
``````

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_
I’m out of my depth here. Maybe you should open an issue on `MKL.jl`