dgemmjit looks identical to regular MKL, suggesting (a) the JIT isn't working (compiler warnings about ignored preprocessor directives reinforce this), and (b) a way to benchmark both OpenBLAS and MKL.

Author: chriselrod

benchmark/benchmarkflops.jl

@@ -52,7 +52,9 @@ function matmul_bench!(br, C, A, B, i)
     @assert C ≈ Cblas "eigen gemm wrong?"; fill!(C, NaN)
     br[10,i] = n_gflop / @belapsed iegemm!($C, $A, $B)
     @assert C ≈ Cblas "i-eigen gemm wrong?"; fill!(C, NaN)
-    br[11,i] = n_gflop / @belapsed gemmavx!($C, $A, $B)
+    br[11,i] = n_gflop / @belapsed dgemmjit!($C, $A, $B)
+    @assert C ≈ Cblas "MKL JIT gemm wrong?"; fill!(C, NaN)
+    br[12,i] = n_gflop / @belapsed gemmavx!($C, $A, $B)
     @assert C ≈ Cblas "LoopVec gemm wrong?"
 end
 function A_mul_B_bench!(br, s, i)
@@ -97,7 +99,7 @@ const BLASTESTS = [
     "GFortran", "GFort-intrinsic",
     "icc", "ifort", "ifort-intrinsic",
     "Clang++ & Eigen-3", "icpc & Eigen-3",
-    "LoopVectorization"
+    "MKL JIT", "LoopVectorization"
 ]
 
 function benchmark_AmulB(sizes)
@@ -209,7 +211,9 @@ function gemv_bench!(br, x, A, y, i)
     @assert x ≈ xblas "eigen wrong?"; fill!(x, NaN);
     br[10,i] = n_gflop / @belapsed iegemv!($x, $A, $y)
     @assert x ≈ xblas "i-eigen wrong?"; fill!(x, NaN);
-    br[11,i] = n_gflop / @belapsed jgemvavx!($x, $A, $y)
+    br[11,i] = n_gflop / @belapsed dgemmjit!($x, $A, $y)
+    @assert x ≈ xblas "gemmjit wrong?"; fill!(x, NaN);
+    br[12,i] = n_gflop / @belapsed jgemvavx!($x, $A, $y)
     @assert x ≈ xblas "LoopVec wrong?"
 end
 function A_mul_vb_bench!(br, s, i)

benchmark/directcalljit.f90

@@ -0,0 +1,33 @@
+module jitmul
+
+include "/opt/intel/mkl/include/mkl_direct_call.fi"
+
+use ISO_C_BINDING
+implicit none
+
+contains
+
+    ! subroutine dgemmjit(C,A,B,M,K,N,alpha,beta) bind(C, name = "dgemmjit")
+  subroutine dgemmjit(C,A,B,M,K,N,At,Bt) bind(C, name = "dgemmjit")
+    integer(C_int32_t),  intent(in)  :: M, K, N
+    integer(C_int8_t),   intent(in) :: At, Bt
+    real(C_double), parameter :: alpha = 1.0D0, beta = 0.0D0
+    ! real(C_double),                 intent(in)  :: alpha, beta
+    real(C_double), dimension(M,K), intent(in)  :: A
+    real(C_double), dimension(K,N), intent(in)  :: B
+    real(C_double), dimension(M,N), intent(out) :: C
+    character :: Atc, Btc
+    if (At == 1_C_int8_t) then
+       Atc = 'T'
+    else
+       Atc = 'N'
+    end if
+    if (Bt == 1_C_int8_t) then
+       Btc = 'T'
+    else
+       Btc = 'N'
+    end if
+    call dgemm(Atc, Btc, M, N, K, alpha, A, M, B, K, beta, C, M)
+  end subroutine dgemmjit
+
+end module jitmul

benchmark/driver.jl

@@ -51,8 +51,8 @@ save(joinpath(PICTURES, "bench_AplusAt_v$v.$filetype"), plot(AplusAt_bench));
 save(joinpath(PICTURES, "bench_Amulvb_v$v.$filetype"), plot(Amulvb_bench));
 save(joinpath(PICTURES, "bench_Atmulvb_v$v.$filetype"), plot(Atmulvb_bench));
 save(joinpath(PICTURES, "bench_exp_v$v.$filetype"), plot(exp_bench));
-save(joinpath(PICTURES, "bench_logdettriangle_v$v.$filetype"), plot(logdettriangle_bench));
 save(joinpath(PICTURES, "bench_random_access_v$v.$filetype"), plot(randomaccess_bench));
+save(joinpath(PICTURES, "bench_logdettriangle_v$v.$filetype"), plot(logdettriangle_bench));
 
 
 

benchmark/loadsharedlibs.jl

@@ -1,4 +1,4 @@
-
+using LinearAlgebra
 using LoopVectorization.VectorizationBase: REGISTER_SIZE
 
 pkgdir(pkg::String) = abspath(joinpath(dirname(Base.find_package(pkg)), ".."))
@@ -11,6 +11,7 @@ const LIBICTEST = joinpath(LOOPVECBENCHDIR, "libictests.so")
 const LIBIFTEST = joinpath(LOOPVECBENCHDIR, "libiftests.so")
 const LIBEIGENTEST = joinpath(LOOPVECBENCHDIR, "libetest.so")
 const LIBIEIGENTEST = joinpath(LOOPVECBENCHDIR, "libietest.so")
+const LIBDIRECTCALLJIT = joinpath(LOOPVECBENCHDIR, "libdcjtest.so")
 
 # requires Clang with polly to build
 cfile = joinpath(LOOPVECBENCHDIR, "looptests.c")
@@ -23,7 +24,7 @@ end
 ffile = joinpath(LOOPVECBENCHDIR, "looptests.f90")
 if !isfile(LIBFTEST) || mtime(ffile) > mtime(LIBFTEST)
     # --param max-unroll-times defaults to ≥8, which is generally excessive
-    run(`gfortran -Ofast -march=native -funroll-loops --param max-unroll-times=4 -floop-nest-optimize -mprefer-vector-width=$(8REGISTER_SIZE) -shared -fPIC $ffile -o $LIBFTEST`)
+    run(`gfortran -Ofast -march=native -funroll-loops -floop-nest-optimize -mprefer-vector-width=$(8REGISTER_SIZE) -shared -fPIC $ffile -o $LIBFTEST`)
 end
 if !isfile(LIBIFTEST) || mtime(ffile) > mtime(LIBIFTEST)
     run(`ifort -fast -qopt-zmm-usage=high -qoverride-limits -shared -fPIC $ffile -o $LIBIFTEST`)
@@ -39,6 +40,26 @@ if !isfile(LIBIEIGENTEST) || mtime(eigenfile) > mtime(LIBIEIGENTEST)
     run(`icpc -fast -qopt-zmm-usage=high -fargument-noalias-global -qoverride-limits -I/usr/include/eigen3 -shared -fPIC $eigenfile -o $LIBIEIGENTEST`)
 end
 
+directcalljitfile = joinpath(LOOPVECBENCHDIR, "directcalljit.f90")
+if !isfile(LIBDIRECTCALLJIT) || mtime(directcalljitfile) > mtime(LIBDIRECTCALLJIT)
+    # run(`ifort -fast -DMKL_DIRECT_CALL_SEQ_JIT -fpp -qopt-zmm-usage=high -shared -fPIC $directcalljitfile -o $LIBDIRECTCALLJIT`)
+    run(`gfortran -Ofast -march=native -DMKL_DIRECT_CALL_SEQ_JIT -cpp -mprefer-vector-width=$(8REGISTER_SIZE) -shared -fPIC $directcalljitfile -o $LIBDIRECTCALLJIT`)
+end
+
+istransposed(x) = false
+istransposed(x::Adjoint) = true
+istransposed(x::Transpose) = true
+function dgemmjit!(C::AbstractVecOrMat{Float64}, A::AbstractVecOrMat{Float64}, B::AbstractVecOrMat{Float64})
+    M, N = size(C); K = size(B, 1)
+    ccall(
+        (:dgemmjit, LIBDIRECTCALLJIT), Cvoid,
+        (Ptr{Float64},Ptr{Float64},Ptr{Float64},Ref{Int},Ref{Int},Ref{Int},Ref{Bool},Ref{Bool}),
+        parent(C), parent(A), parent(B),
+        Ref(M), Ref(K), Ref(N),
+        Ref(istransposed(A)), Ref(istransposed(B))
+    )
+end
+
 for (prefix,Cshared,Fshared,Eshared) ∈ ((Symbol(""),LIBCTEST,LIBFTEST,LIBEIGENTEST), (:i,LIBICTEST,LIBIFTEST,LIBIEIGENTEST))
     for order ∈ (:kmn, :knm, :mkn, :mnk, :nkm, :nmk)
         gemm = Symbol(:gemm_, order)
@@ -59,9 +80,9 @@ for (prefix,Cshared,Fshared,Eshared) ∈ ((Symbol(""),LIBCTEST,LIBFTEST,LIBEIGEN
             )
         end
     end
-    @eval @inline $(Symbol(prefix,:cgemm!))(C, A, B) = $(Symbol(prefix, :cgemm_nkm!))(C, A, B)
-    @eval @inline $(Symbol(prefix,:fgemm!))(C, A, B) = $(Symbol(prefix, :fgemm_nkm!))(C, A, B)
-    @eval @inline function $(Symbol(prefix,:egemm!))(C, A, B)
+    @eval $(Symbol(prefix,:cgemm!))(C, A, B) = $(Symbol(prefix, :cgemm_nkm!))(C, A, B)
+    @eval $(Symbol(prefix,:fgemm!))(C, A, B) = $(Symbol(prefix, :fgemm_nkm!))(C, A, B)
+    @eval function $(Symbol(prefix,:egemm!))(C, A, B)
         M, N = size(C); K = size(B, 1)
         ccall(
             (:AmulB, $Eshared), Cvoid,
@@ -78,7 +99,7 @@ for (prefix,Cshared,Fshared,Eshared) ∈ ((Symbol(""),LIBCTEST,LIBFTEST,LIBEIGEN
         )
     end
 for (p,s) ∈ [(:c,Cshared) (:e,Eshared)]
-    @eval @inline function $(Symbol(prefix,p,:gemm!))(C, A::Adjoint, B)
+    @eval function $(Symbol(prefix,p,:gemm!))(C, A::Adjoint, B)
         M, N = size(C); K = size(B, 1)
         ccall(
             (:AtmulB, $s), Cvoid,
@@ -87,15 +108,15 @@ for (p,s) ∈ [(:c,Cshared) (:e,Eshared)]
         )
     end
 end
-    @eval @inline function $(Symbol(prefix,:fgemm!))(C, A::Adjoint, B)
+    @eval function $(Symbol(prefix,:fgemm!))(C, A::Adjoint, B)
         M, N = size(C); K = size(B, 1)
         ccall(
             (:AtmulB, $Fshared), Cvoid,
             (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}, Ref{Clong}),
             C, parent(A), B, Ref(M), Ref(K), Ref(N)
         )
     end
-    @eval @inline function $(Symbol(prefix,:fgemm_builtin!))(C, A::Adjoint, B)
+    @eval function $(Symbol(prefix,:fgemm_builtin!))(C, A::Adjoint, B)
         M, N = size(C); K = size(B, 1)
         ccall(
             (:AtmulBbuiltin, $Fshared), Cvoid,
@@ -104,7 +125,7 @@ end
         )
     end
 for (p,s) ∈ [(:c,Cshared) (:e,Eshared)]
-    @eval @inline function $(Symbol(prefix,p,:gemm!))(C, A, B::Adjoint)
+    @eval function $(Symbol(prefix,p,:gemm!))(C, A, B::Adjoint)
         M, N = size(C); K = size(B, 1)
         ccall(
             (:AmulBt, $s), Cvoid,
@@ -113,15 +134,15 @@ for (p,s) ∈ [(:c,Cshared) (:e,Eshared)]
         )
     end
 end
-    @eval @inline function $(Symbol(prefix,:fgemm!))(C, A, B::Adjoint)
+    @eval function $(Symbol(prefix,:fgemm!))(C, A, B::Adjoint)
         M, N = size(C); K = size(B, 1)
         ccall(
             (:AmulBt, $Fshared), Cvoid,
             (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}, Ref{Clong}),
             C, A, parent(B), Ref(M), Ref(K), Ref(N)
         )
     end
-    @eval @inline function $(Symbol(prefix,:fgemm_builtin!))(C, A, B::Adjoint)
+    @eval function $(Symbol(prefix,:fgemm_builtin!))(C, A, B::Adjoint)
         M, N = size(C); K = size(B, 1)
         ccall(
             (:AmulBtbuiltin, $Fshared), Cvoid,
@@ -130,7 +151,7 @@ end
         )
     end
 for (p,s) ∈ [(:c,Cshared) (:e,Eshared)]
-    @eval @inline function $(Symbol(prefix,p,:gemm!))(C, A::Adjoint, B::Adjoint)
+    @eval function $(Symbol(prefix,p,:gemm!))(C, A::Adjoint, B::Adjoint)
         M, N = size(C); K = size(B, 1)
         ccall(
             (:AtmulBt, $s), Cvoid,
@@ -139,15 +160,15 @@ for (p,s) ∈ [(:c,Cshared) (:e,Eshared)]
         )
     end
 end
-    @eval @inline function $(Symbol(prefix,:fgemm!))(C, A::Adjoint, B::Adjoint)
+    @eval function $(Symbol(prefix,:fgemm!))(C, A::Adjoint, B::Adjoint)
         M, N = size(C); K = size(B, 1)
         ccall(
             (:AtmulBt, $Fshared), Cvoid,
             (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}, Ref{Clong}),
             C, parent(A), parent(B), Ref(M), Ref(K), Ref(N)
         )
     end
-    @eval @inline function $(Symbol(prefix,:fgemm_builtin!))(C, A::Adjoint, B::Adjoint)
+    @eval function $(Symbol(prefix,:fgemm_builtin!))(C, A::Adjoint, B::Adjoint)
         M, N = size(C); K = size(B, 1)
         ccall(
             (:AtmulBtbuiltin, $Fshared), Cvoid,
@@ -242,7 +263,7 @@ end
         )
     end
 for (p,s) ∈ [(:c,Cshared) (:e,Eshared)]
-    @eval @inline function $(Symbol(prefix,p,:gemv!))(y, A::Adjoint, x)
+    @eval function $(Symbol(prefix,p,:gemv!))(y, A::Adjoint, x)
         M, K = size(A)
         ccall(
             (:Atmulvb, $s), Cvoid,
@@ -251,7 +272,7 @@ for (p,s) ∈ [(:c,Cshared) (:e,Eshared)]
         )
     end
 end
-@eval @inline function $(Symbol(prefix,:fgemv!))(y, A::Adjoint, x)
+@eval function $(Symbol(prefix,:fgemv!))(y, A::Adjoint, x)
         M, K = size(A)
         ccall(
             (:Atmulvb, $Fshared), Cvoid,

benchmark/looptests.jl

@@ -10,7 +10,7 @@ function jgemm!(𝐂, 𝐀, 𝐁)
         end
     end
 end
-@inline function jgemm!(𝐂, 𝐀ᵀ::Adjoint, 𝐁)
+function jgemm!(𝐂, 𝐀ᵀ::Adjoint, 𝐁)
     𝐀 = parent(𝐀ᵀ)
     @inbounds for n ∈ 1:size(𝐂,2), m ∈ 1:size(𝐂,1)
         𝐂ₘₙ = zero(eltype(𝐂))
@@ -20,7 +20,7 @@ end
         𝐂[m,n] = 𝐂ₘₙ
     end
 end
-@inline function jgemm!(𝐂, 𝐀, 𝐁ᵀ::Adjoint)
+function jgemm!(𝐂, 𝐀, 𝐁ᵀ::Adjoint)
     𝐂 .= 0
     𝐁 = parent(𝐁ᵀ)
     M, N = size(𝐂); K = size(𝐁ᵀ,1)
@@ -30,7 +30,7 @@ end
         end
     end
 end
-@inline function jgemm!(𝐂, 𝐀ᵀ::Adjoint, 𝐁ᵀ::Adjoint)
+function jgemm!(𝐂, 𝐀ᵀ::Adjoint, 𝐁ᵀ::Adjoint)
     𝐂 .= 0
     𝐀 = parent(𝐀ᵀ)
     𝐁 = parent(𝐁ᵀ)
@@ -41,7 +41,7 @@ end
         end
     end
 end
-@inline function gemmavx!(𝐂, 𝐀, 𝐁)
+function gemmavx!(𝐂, 𝐀, 𝐁)
     @avx for m ∈ 1:size(𝐀,1), n ∈ 1:size(𝐁,2)
         𝐂ₘₙ = zero(eltype(𝐂))
         for k ∈ 1:size(𝐀,2)
@@ -128,7 +128,7 @@ function jgemv!(y, 𝐀, x)
         end
     end
 end
-@inline function jgemv!(𝐲, 𝐀ᵀ::Adjoint, 𝐱)
+function jgemv!(𝐲, 𝐀ᵀ::Adjoint, 𝐱)
     𝐀 = parent(𝐀ᵀ)
     @inbounds for i ∈ eachindex(𝐲)
         𝐲ᵢ = zero(eltype(𝐲))
@@ -138,7 +138,7 @@ end
         𝐲[i] = 𝐲ᵢ
     end
 end
-@inline function jgemvavx!(𝐲, 𝐀, 𝐱)
+function jgemvavx!(𝐲, 𝐀, 𝐱)
     @avx for i ∈ eachindex(𝐲)
         𝐲ᵢ = zero(eltype(𝐲))
         for j ∈ eachindex(𝐱)