Added benchmarks folder, made a few tweeks to unrolling behavior.

Author: chriselrod

.gitignore

@@ -7,4 +7,8 @@
 *~
 src/#*#
 tests/#*#
+benchmarks/#*#
 *.mem
+*.mod
+*.mod0
+*.so

benchmarks/driver.jl

@@ -0,0 +1,87 @@
+pkgdir(pkg::String) = abspath(joinpath(dirname(Base.find_package(pkg)), ".."))
+const LOOPVECBENCHDIR = joinpath(pkgdir("LoopVectorization"), "benchmarks")
+include(joinpath(LOOPVECBENCHDIR, "looptests.jl"))
+include(joinpath(LOOPVECBENCHDIR, "loadsharedlibs.jl"))
+
+using PrettyTables, BenchmarkTools
+struct SizedResults{V <: AbstractVector} <: AbstractMatrix{String}
+    results::Matrix{Float64}
+    sizes::V
+end
+function Base.size(sr::SizedResults)
+    M, N = size(sr.results)
+    N, M+1
+end
+struct BenchmarkResult{V}
+    tests::Vector{String}
+    sizedresults::SizedResults{V}
+end
+function BenchmarkResult(tests, sizes)
+    ntests = length(tests); nsizes = length(sizes)
+    BenchmarkResult(
+        append!(["Size"], tests),
+        SizedResults(Matrix{Float64}(undef, ntests, nsizes), sizes)
+    )
+end
+function Base.getindex(br::SizedResults, row, col)
+    col == 1 ? string(br.sizes[row]) : string(br.results[col - 1, row])
+end
+Base.setindex!(br::BenchmarkResult, v, i...) = br.sizedresults.results[i...] = v
+function Base.show(io::IO, br::BenchmarkResult)
+    pretty_table(io, br.sizedresults, br.tests)
+end
+function alloc_matrices(s::NTuple{3,Int})
+    M, K, N = s
+    C = Matrix{Float64}(undef, M, N)
+    A = rand(M, K)
+    B = rand(K, N)
+    C, A, B
+end
+alloc_matrices(s::Int) = alloc_matrices((s,s,s))
+gflop(s::Int) = s^3 * 1e-9
+gflop(s::NTuple{3,Int}) = prod(s) * 1e-9
+function benchmark_gemm(sizes)
+    tests = [BLAS.vendor() === :mkl ? "IntelMKL" : "OpenBLAS", "Julia", "Clang-Polly", "GFort-loops", "GFort-intrinsic", "LoopVectorization"]
+    br = BenchmarkResult(tests, sizes)
+    for (i,s) ∈ enumerate(sizes)
+        C, A, B = alloc_matrices(s)
+        n_gflop = gflop(s)
+        br[1,i] = n_gflop / @belapsed mul!($C, $A, $B)
+        Cblas = copy(C)
+        br[2,i] = n_gflop / @belapsed jgemm_nkm!($C, $A, $B)
+        @assert C ≈ Cblas "Julia gemm wrong?"
+        br[3,i] = n_gflop / @belapsed cgemm_nkm!($C, $A, $B)
+        @assert C ≈ Cblas "Polly gemm wrong?"
+        br[4,i] = n_gflop / @belapsed fgemm_nkm!($C, $A, $B)
+        @assert C ≈ Cblas "Fort gemm wrong?"
+        br[5,i] = n_gflop / @belapsed fgemm_builtin!($C, $A, $B)
+        @assert C ≈ Cblas "Fort intrinsic gemm wrong?"
+        br[6,i] = n_gflop / @belapsed gemmavx!($C, $A, $B)
+        @assert C ≈ Cblas "LoopVec gemm wrong?"
+    end
+    br
+end
+
+using VegaLite, IndexedTables
+function plot(br::BenchmarkResult)
+    res = vec(br.sizedresults.results)
+    brsizes = br.sizedresults.sizes
+    sizes = Vector{eltype(brsizes)}(undef, length(res))
+    ntests = length(br.tests) - 1
+    for i ∈ 0:length(brsizes)-1
+        si = brsizes[i+1]
+        for j ∈ 1:ntests
+            sizes[j + i*ntests] = si
+        end
+    end
+    tests = vcat((@view(br.tests[2:end]) for _ ∈ eachindex(brsizes))...)
+    t = table((GFLOPS = res, Size = sizes, Method = tests))
+    t |> @vlplot(
+        :line,
+        x = :Size,
+        y = :GFLOPS,
+        color = :Method
+    )
+end
+
+

benchmarks/loadsharedlibs.jl

@@ -0,0 +1,148 @@
+
+using VectorizationBase: REGISTER_SIZE
+# run(`gfortran `)
+
+pkgdir(pkg::String) = abspath(joinpath(dirname(Base.find_package(pkg)), ".."))
+const LOOPVECBENCHDIR = joinpath(pkgdir("LoopVectorization"), "benchmarks")
+include(joinpath(LOOPVECBENCHDIR, "looptests.jl"))
+
+const LIBCTEST = joinpath(LOOPVECBENCHDIR, "libctests.so")
+const LIBFTEST = joinpath(LOOPVECBENCHDIR, "libftests.so")
+
+# requires Clang with polly to build
+if !isfile(LIBCTEST)
+    cfile = joinpath(LOOPVECBENCHDIR, "looptests.c")
+    run(`clang -Ofast -march=native -mprefer-vector-width=$(8REGISTER_SIZE) -mllvm -polly -mllvm -polly-vectorizer=stripmine -shared -fPIC $cfile -o $LIBCTEST`)
+end
+if !isfile(LIBFTEST)
+    ffile = joinpath(LOOPVECBENCHDIR, "looptests.f90") # --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`)
+end
+
+for order ∈ (:kmn, :knm, :mkn, :mnk, :nkm, :nmk)
+    gemm = Symbol(:gemm_, order)
+    @eval function $(Symbol(:c, gemm, :!))(C, A, B)
+        M, N = size(C); K = size(B, 1)
+        ccall(
+            ($(QuoteNode(gemm)), LIBCTEST), Cvoid,
+            (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong, Clong),
+            C, A, B, M, K, N
+        )
+    end
+    @eval function $(Symbol(:f, gemm, :!))(C, A, B)
+        M, N = size(C); K = size(B, 1)
+        ccall(
+            ($(QuoteNode(gemm)), LIBFTEST), Cvoid,
+            (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}, Ref{Clong}),
+            C, A, B, Ref(M), Ref(K), Ref(N)
+        )
+    end
+end
+
+function fgemm_builtin!(C, A, B)
+    M, N = size(C); K = size(B, 1)
+    ccall(
+        (:gemmbuiltin, LIBFTEST), Cvoid,
+        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}, Ref{Clong}),
+        C, A, B, Ref(M), Ref(K), Ref(N)
+    )
+end
+
+function cdot(a, b)
+    N = length(a)
+    ccall(
+        (:dot, LIBCTEST), Float64,
+        (Ptr{Float64}, Ptr{Float64}, Clong),
+        a, b, N
+    )
+end
+function fdot(a, b)
+    N = length(a)
+    d = Ref{Float64}()
+    ccall(
+        (:dot, LIBFTEST), Cvoid,
+        (Ref{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}),
+        d, a, b, Ref(N)
+    )
+    d[]
+end
+function cselfdot(a)
+    N = length(a)
+    ccall(
+        (:selfdot, LIBCTEST), Float64,
+        (Ptr{Float64}, Clong),
+        a, N
+    )
+end
+function fselfdot(a)
+    N = length(a)
+    d = Ref{Float64}()
+    ccall(
+        (:selfdot, LIBFTEST), Cvoid,
+        (Ref{Float64}, Ptr{Float64}, Ref{Clong}),
+        d, a, Ref(N)
+    )
+    d[]
+end
+
+function cgemv!(y, A, x)
+    M, K = size(A)
+    ccall(
+        (:gemv, LIBCTEST), Cvoid,
+        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong),
+        y, A, x, M, K
+    )
+end
+function fgemv!(y, A, x)
+    M, K = size(A)
+    ccall(
+        (:gemv, LIBFTEST), Cvoid,
+        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}),
+        y, A, x, Ref(M), Ref(K)
+    )
+end
+function fgemv_builtin!(y, A, x)
+    M, K = size(A)
+    ccall(
+        (:gemvbuiltin, LIBFTEST), Cvoid,
+        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}),
+        y, A, x, Ref(M), Ref(K)
+    )
+end
+
+function caplusBc!(D, a, B, c)
+    M, K = size(B)
+    ccall(
+        (:aplusBc, LIBCTEST), Cvoid,
+        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong),
+        y, A, x, M, K
+    )
+end
+function faplusBc!(D, a, B, c)
+    M, K = size(B)
+    ccall(
+        (:aplusBc, LIBFTEST), Cvoid,
+        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}),
+        y, A, x, Ref(M), Ref(K)
+    )
+end
+function cOLSlp(y, X, β)
+    N, P = size(X)
+    ccall(
+        (:OLSlp, LIBCTEST), Float64,
+        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong),
+        y, X, β, N, P
+    )
+end
+function fOLSlp(y, X, β)
+    N, P = size(X)
+    lp = Ref{Float64}()
+    ccall(
+        (:OLSlp, LIBFTEST), Cvoid,
+        (Ref{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}),
+        lp, y, X, β, Ref(N), Ref(P)
+    )
+    lp[]
+end
+
+

benchmarks/looptests.c

@@ -0,0 +1,142 @@
+
+
+void gemm_mnk(double* restrict C, double* restrict A, double* restrict B, long M, long K, long N){
+  for (long i = 0; i < M*N; i++){
+    C[i] = 0.0;
+  }
+  for (long m = 0; m < M; m++){
+    for (long n = 0; n < N; n++){
+      for (long k = 0; k < K; k++){
+	C[m + n*M] += A[m + k*M] * B[k + n*K];
+      }
+    }
+  }
+  return;
+}
+void gemm_mkn(double* restrict C, double* restrict A, double* restrict B, long M, long K, long N){
+  for (long i = 0; i < M*N; i++){
+    C[i] = 0.0;
+  }
+  for (long m = 0; m < M; m++){
+    for (long k = 0; k < K; k++){
+      for (long n = 0; n < N; n++){
+	C[m + n*M] += A[m + k*M] * B[k + n*K];
+      }
+    }
+  }
+  return;
+}
+void gemm_nmk(double* restrict C, double* restrict A, double* restrict B, long M, long K, long N){
+  for (long i = 0; i < M*N; i++){
+    C[i] = 0.0;
+  }
+  for (long n = 0; n < N; n++){
+    for (long m = 0; m < M; m++){
+      for (long k = 0; k < K; k++){
+	C[m + n*M] += A[m + k*M] * B[k + n*K];
+      }
+    }
+  }
+  return;
+}
+void gemm_nkm(double* restrict C, double* restrict A, double* restrict B, long M, long K, long N){
+  for (long i = 0; i < M*N; i++){
+    C[i] = 0.0;
+  }
+  for (long n = 0; n < N; n++){
+    for (long k = 0; k < K; k++){
+      for (long m = 0; m < M; m++){
+	C[m + n*M] += A[m + k*M] * B[k + n*K];
+      }
+    }
+  }
+  return;
+}
+void gemm_kmn(double* restrict C, double* restrict A, double* restrict B, long M, long K, long N){
+  for (long i = 0; i < M*N; i++){
+    C[i] = 0.0;
+  }
+  for (long k = 0; k < K; k++){
+    for (long m = 0; m < M; m++){
+      for (long n = 0; n < N; n++){
+	C[m + n*M] += A[m + k*M] * B[k + n*K];
+      }
+    }
+  }
+  return;
+}
+void gemm_knm(double* restrict C, double* restrict A, double* restrict B, long M, long K, long N){
+  for (long i = 0; i < M*N; i++){
+    C[i] = 0.0;
+  }
+  for (long k = 0; k < K; k++){
+    for (long n = 0; n < N; n++){
+      for (long m = 0; m < M; m++){
+	C[m + n*M] += A[m + k*M] * B[k + n*K];
+      }
+    }
+  }
+  return;
+}
+double dot(double* restrict a, double* restrict b, long N){
+  double s = 0.0;
+  for (long n = 0; n < N; n++){
+    s += a[n]*b[n];
+  }
+  return s;
+}
+double selfdot(double* restrict a, long N){
+  double s = 0.0;
+  for (long n = 0; n < N; n++){
+    s += a[n]*a[n];
+  }
+  return s;
+}
+
+void gemv(double* restrict y, double* restrict  A, double* restrict x, long M, long K){
+  for (long m = 0; m < M; m++){
+    y[m] = 0.0;
+  }
+  for (long k = 0; k < K; k++){
+    for (long m = 0; m < M; m++){
+      y[m] += A[m + k*M] * x[k]; 
+    }
+  }
+  return;
+}
+
+void unscaledvar(double* restrict s, double* restrict A, double* restrict xb, long M, long N){
+  for (long m = 0; m < M; m++){
+    s[m] = 0.0;
+  }
+  for (long n = 0; n < N; n++){
+    for (long m = 0; m < M; m++){
+      double d = A[m + n*M] - xb[m];
+      s[m] += d*d;
+    }
+  }
+  return;
+}
+
+void aplucBc(double* restrict D, double* restrict a, double* restrict B, double* restrict c, long M, long N){
+  for (long n = 0; n < N; n++){
+    for (long m = 0; m < M; m++){
+      D[m + n*M] = a[m] + B[m + n*M] * c[n];
+    }
+  }
+  return;
+}
+
+double OLSlp(double* restrict y, double* restrict X, double* restrict b, long N, long P){
+  double lp = 0.0;
+  for (long n = 0; n < N; n++){
+    double d = y[n];
+    for (long p = 0; p < P; p++){
+      d -= X[n + p*N] * b[p];
+    }
+    lp += d*d;
+  }
+  return lp;
+}
+
+