Update documentation and benchmarks.

Author: chriselrod

.travis.yml

@@ -16,7 +16,7 @@ after_success:
 jobs:
   include:
     - stage: Documentation
-      julia: 1.0
+      julia: 1.3
       script: julia --project=docs -e '
           using Pkg;
           Pkg.develop(PackageSpec(path=pwd()));

benchmark/benchmarkflops.jl

@@ -10,59 +10,48 @@ include(joinpath(LOOPVECBENCHDIR, "plotbenchmarks.jl"))
 
 using Distributed
 
-addprocs(13);
+addprocs((Sys.CPU_THREADS >> 1)-1); nprocs()
 
 @everywhere begin
     pkgdir(pkg::String) = abspath(joinpath(dirname(Base.find_package(pkg)), ".."))
     const LOOPVECBENCHDIR = joinpath(pkgdir("LoopVectorization"), "benchmark")
     include(joinpath(LOOPVECBENCHDIR, "benchmarkflops.jl"))
-    BenchmarkTools.DEFAULT_PARAMETERS.seconds = 1
+    # BenchmarkTools.DEFAULT_PARAMETERS.seconds = 1
 end
 
-gemm_future = @spawnat 2 benchmark_gemm(2:256);
-AtmulB_future = @spawnat 3 benchmark_AtmulB(2:256);
-dot_future = @spawnat 4 benchmark_dot(2:256);
-selfdot_future = @spawnat 5 benchmark_selfdot(2:256);
-gemv_future = @spawnat 6 benchmark_gemv(2:256);
-dot3_future = @spawnat 7 benchmark_dot3(2:256);
-sse_future = @spawnat 8 benchmark_sse(2:256);
-exp_future = @spawnat 9 benchmark_exp(2:256);
-aplusBc_future = @spawnat 10 benchmark_aplusBc(2:256);
-AplusAt_future = @spawnat 11 benchmark_AplusAt(2:256);
-randomaccess_future = @spawnat 12 benchmark_random_access(2:256);
-AmulBt_future = @spawnat 13 benchmark_AmulBt(2:256);
-Atmulvb_future = @spawnat 14 benchmark_Atmulvb(2:256);
-
-dot_bench = fetch(dot_future)
-selfdot_bench = fetch(selfdot_future)
-gemv_bench = fetch(gemv_future)
-randomaccess_bench = fetch(randomaccess_future)
-dot3_bench = fetch(dot3_future)
-sse_bench = fetch(sse_future)
-exp_bench = fetch(exp_future)
-AplusAt_bench = fetch(AplusAt_future)
-aplusBc_bench = fetch(aplusBc_future)
-gemm_bench = fetch(gemm_future)
-AtmulB_bench = fetch(AtmulB_future)
-AmulBt_bench = fetch(AmulBt_future)
-Atmulvb_bench = fetch(Atmulvb_future)
-
+AmulB_bench = benchmark_AmulB(2:256)
+AmulBt_bench = benchmark_AmulBt(2:256)
+AtmulB_bench = benchmark_AtmulB(2:256)
+AtmulBt_bench = benchmark_AtmulBt(2:256)
+dot_bench = benchmark_dot(2:256)
+selfdot_bench = benchmark_selfdot(2:256)
+Amulvb_bench = benchmark_Amulvb(2:256)
+Atmulvb_bench = benchmark_Atmulvb(2:256)
+dot3_bench = benchmark_dot3(2:256)
+sse_bench = benchmark_sse(2:256)
+exp_bench = benchmark_exp(2:256)
+aplusBc_bench = benchmark_aplusBc(2:256)
+AplusAt_bench = benchmark_AplusAt(2:256)
+randomaccess_bench = benchmark_random_access(2:256)
+logdettriangle_bench = benchmark_logdettriangle(2:256)
 
 v = 1
 filetype = "svg"
 const PICTURES = joinpath(pkgdir("LoopVectorization"), "docs", "src", "assets")
 save(joinpath(PICTURES, "bench_exp_v$v.$filetype"), plot(exp_bench));
-save(joinpath(PICTURES, "bench_gemm_v$v.$filetype"), plot(gemm_bench));
+save(joinpath(PICTURES, "bench_logdettriangle_v$v.$filetype"), plot(logdettriangle_bench));
+save(joinpath(PICTURES, "bench_AmulB_v$v.$filetype"), plot(AmulB_bench));
+save(joinpath(PICTURES, "bench_AmulBt_v$v.$filetype"), plot(AmulBt_bench));
 save(joinpath(PICTURES, "bench_AtmulB_v$v.$filetype"), plot(AtmulB_bench));
+save(joinpath(PICTURES, "bench_AtmulBt_v$v.$filetype"), plot(AtmulBt_bench));
 save(joinpath(PICTURES, "bench_dot_v$v.$filetype"), plot(dot_bench));
 save(joinpath(PICTURES, "bench_selfdot_v$v.$filetype"), plot(selfdot_bench));
-save(joinpath(PICTURES, "bench_gemv_v$v.$filetype"), plot(gemv_bench));
 save(joinpath(PICTURES, "bench_dot3_v$v.$filetype"), plot(dot3_bench));
 save(joinpath(PICTURES, "bench_sse_v$v.$filetype"), plot(sse_bench));
 save(joinpath(PICTURES, "bench_aplusBc_v$v.$filetype"), plot(aplusBc_bench));
 save(joinpath(PICTURES, "bench_AplusAt_v$v.$filetype"), plot(AplusAt_bench));
 save(joinpath(PICTURES, "bench_random_access_v$v.$filetype"), plot(randomaccess_bench));
-save(joinpath(PICTURES, "bench_AmulBt_v$v.$filetype"), plot(AmulBt_bench));
+save(joinpath(PICTURES, "bench_Amulvb_v$v.$filetype"), plot(Amulvb_bench));
 save(joinpath(PICTURES, "bench_Atmulvb_v$v.$filetype"), plot(Atmulvb_bench));
 
 

benchmark/driver.jl

@@ -16,15 +16,15 @@ if !isfile(LIBCTEST) || mtime(cfile) > mtime(LIBCTEST)
     run(`clang -Ofast -march=native -mprefer-vector-width=$(8REGISTER_SIZE) -lm -mllvm -polly -mllvm -polly-vectorizer=stripmine -shared -fPIC $cfile -o $LIBCTEST`)
 end
 if !isfile(LIBICTEST) || mtime(cfile) > mtime(LIBICTEST)
-    run(`icc -fast -qopt-zmm-usage=high -shared -fPIC $cfile -o $LIBICTEST`)
+    run(`icc -fast -qopt-zmm-usage=high -fargument-noalias-global -qoverride-limits -shared -fPIC $cfile -o $LIBICTEST`)
 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`)
 end
 if !isfile(LIBIFTEST) || mtime(ffile) > mtime(LIBIFTEST)
-    run(`ifort -fast -qopt-zmm-usage=high -shared -fPIC $ffile -o $LIBIFTEST`)
+    run(`ifort -fast -qopt-zmm-usage=high -qoverride-limits -shared -fPIC $ffile -o $LIBIFTEST`)
 end
 
 for (prefix,Cshared,Fshared) ∈ ((Symbol(""),LIBCTEST,LIBFTEST), (:i,LIBICTEST,LIBIFTEST))
@@ -105,6 +105,30 @@ for (prefix,Cshared,Fshared) ∈ ((Symbol(""),LIBCTEST,LIBFTEST), (:i,LIBICTEST,
             C, A, parent(B), Ref(M), Ref(K), Ref(N)
         )
     end
+    @eval @inline function $(Symbol(prefix,:cgemm!))(C, A::Adjoint, B::Adjoint)
+        M, N = size(C); K = size(B, 1)
+        ccall(
+            (:AtmulBt, $Cshared), Cvoid,
+            (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong, Clong),
+            C, parent(A), parent(B), M, K, N
+        )
+    end
+    @eval @inline 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)
+        M, N = size(C); K = size(B, 1)
+        ccall(
+            (:AtmulBtbuiltin, $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 function $(Symbol(prefix,:cdot))(a, b)
         N = length(a)
         ccall(
@@ -312,5 +336,24 @@ for (prefix,Cshared,Fshared) ∈ ((Symbol(""),LIBCTEST,LIBFTEST), (:i,LIBICTEST,
         )
         p[]
     end
-    
+    @eval function $(Symbol(prefix,:clogdettriangle))(T::Union{LowerTriangular,UpperTriangular})
+        N = size(T,1)
+        Tp = parent(T)
+        ccall(
+            (:logdettriangle, $Cshared), Float64,
+            (Ptr{Float64}, Clong),
+            Tp, N
+        )
+    end
+    @eval function $(Symbol(prefix,:flogdettriangle))(T::Union{LowerTriangular,UpperTriangular})
+        N = size(T,1)
+        Tp = parent(T)
+        ld = Ref{Float64}()
+        ccall(
+            (:logdettriangle, $Fshared), Cvoid,
+            (Ref{Float64}, Ptr{Float64}, Ref{Clong}),
+            ld, Tp, Ref(N)
+        )
+        ld[]
+    end
 end

benchmark/loadsharedlibs.jl

@@ -104,6 +104,19 @@ void AmulBt(double* restrict C, double* restrict A, double* restrict Bt, long M,
   }
   return;
 }
+void AtmulBt(double* restrict C, double* restrict A, double* restrict Bt, 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[k + K*m] * Bt[n + N*k];
+      }
+    }
+  }
+  return;
+}
 double dot(double* restrict a, double* restrict b, long N){
   double s = 0.0;
   for (long n = 0; n < N; n++){
@@ -212,5 +225,11 @@ double randomaccess(double* restrict P, long* restrict basis, double* restrict c
   }
   return p;
 }
-
+double logdettriangle(double* T, long N){
+  double ld = 0;
+  for (long n = 0; n < N; n++){
+    ld += log(T[n + n*N]);
+  }
+  return ld;
+}
 

benchmark/looptests.c

@@ -147,6 +147,28 @@ subroutine AmulBtbuiltin(C, A, B, M, K, N) BIND(C, name="AmulBtbuiltin")
       real(C_double), dimension(N, K), intent(in) :: B
       C = matmul(A, transpose(B))
     end subroutine AmulBtbuiltin
+    subroutine AtmulBt(C, A, B, M, K, N) BIND(C, name="AtmulBt")
+      integer(C_long), intent(in) :: M, K, N
+      real(C_double), dimension(M, N), intent(out) :: C
+      real(C_double), dimension(K, M), intent(in) :: A
+      real(C_double), dimension(N, K), intent(in) :: B
+      integer(C_long) :: mm, kk, nn
+      C = 0.0
+      do concurrent(nn = 1:N)
+         do concurrent(kk = 1:K)
+            do concurrent(mm = 1:M)
+               C(mm,nn) = C(mm,nn) + A(kk,mm) * B(nn,kk)
+            end do
+         end do
+      end do
+    end subroutine AtmulBt
+    subroutine AtmulBtbuiltin(C, A, B, M, K, N) BIND(C, name="AtmulBtbuiltin")
+      integer(C_long), intent(in) :: M, K, N
+      real(C_double), dimension(M, N), intent(out) :: C
+      real(C_double), dimension(K, M), intent(in) :: A
+      real(C_double), dimension(N, K), intent(in) :: B
+      C = transpose(matmul(B, A))
+    end subroutine AtmulBtbuiltin
     subroutine dot(s, a, b, N) BIND(C, name="dot")
       integer(C_long), intent(in) :: N
       real(C_double), dimension(N), intent(in) :: a, b
@@ -309,4 +331,14 @@ subroutine randomaccess(pp, P, basis, coefs, A, C) BIND(C, name="randomaccess")
          pp = pp + pc
       end do
     end subroutine randomaccess
+    subroutine logdettriangle(ld, T, N) BIND(C, name="logdettriangle")
+      integer(C_long), intent(in) :: N
+      real(C_double), intent(in) :: T(N,N)
+      real(C_double), intent(out) :: ld
+      integer(C_long) :: nn
+      ld = 0
+      do concurrent(nn = 1:N)
+         ld = ld + log(T(nn,nn))
+      end do
+    end subroutine logdettriangle
   end module looptests

benchmark/looptests.f90

@@ -6,7 +6,7 @@ function jgemm!(𝐂, 𝐀, 𝐁)
     M, N = size(𝐂); K = size(𝐁,1)
     @inbounds for n ∈ 1:N, k ∈ 1:K
         @simd ivdep for m ∈ 1:M
-            𝐂[m,n] += 𝐀[m,k] * 𝐁[k,n]
+            @fastmath 𝐂[m,n] += 𝐀[m,k] * 𝐁[k,n]
         end
     end
 end
@@ -15,18 +15,29 @@ end
     @inbounds for n ∈ 1:size(𝐂,2), m ∈ 1:size(𝐂,1)
         𝐂ₘₙ = zero(eltype(𝐂))
         @simd ivdep for k ∈ 1:size(𝐀,1)
-            𝐂ₘₙ += 𝐀[k,m] * 𝐁[k,n]
+            @fastmath 𝐂ₘₙ += 𝐀[k,m] * 𝐁[k,n]
         end
         𝐂[m,n] = 𝐂ₘₙ
     end
 end
 @inline function jgemm!(𝐂, 𝐀, 𝐁ᵀ::Adjoint)
     𝐂 .= 0
     𝐁 = parent(𝐁ᵀ)
-    M, N = size(𝐂); K = size(𝐁,1)
+    M, N = size(𝐂); K = size(𝐁ᵀ,1)
     @inbounds for k ∈ 1:K, n ∈ 1:N
         @simd ivdep for m ∈ 1:M
-            𝐂[m,n] += 𝐀[m,k] * 𝐁[n,k]
+            @fastmath 𝐂[m,n] += 𝐀[m,k] * 𝐁[n,k]
+        end
+    end
+end
+@inline function jgemm!(𝐂, 𝐀ᵀ::Adjoint, 𝐁ᵀ::Adjoint)
+    𝐂 .= 0
+    𝐀 = parent(𝐀ᵀ)
+    𝐁 = parent(𝐁ᵀ)
+    M, N = size(𝐂); K = size(𝐁ᵀ,1)
+    @inbounds for n ∈ 1:N, k ∈ 1:K
+        @simd ivdep for m ∈ 1:M
+            @fastmath 𝐂[m,n] += 𝐀[k,m] * 𝐁[n,k]
         end
     end
 end
@@ -72,7 +83,7 @@ function jdot3(x, A, y)
     s = zero(promote_type(eltype(x), eltype(A), eltype(y)))
     @inbounds for n ∈ 1:N
         @simd ivdep for m ∈ 1:M
-            s += x[m] * A[m,n] * y[n]
+            @fastmath s += x[m] * A[m,n] * y[n]
         end
     end
     s
@@ -113,7 +124,7 @@ function jgemv!(y, 𝐀, x)
     y .= zero(eltype(y))
     @inbounds for j ∈ eachindex(x)
         @simd ivdep for i ∈ eachindex(y)
-            y[i] += 𝐀[i,j] * x[j]
+            @fastmath y[i] += 𝐀[i,j] * x[j]
         end
     end
 end
@@ -122,7 +133,7 @@ end
     @inbounds for i ∈ eachindex(𝐲)
         𝐲ᵢ = zero(eltype(𝐲))
         @simd ivdep for j ∈ eachindex(𝐱)
-            𝐲ᵢ += 𝐀[j,i] * 𝐱[j]
+            @fastmath 𝐲ᵢ += 𝐀[j,i] * 𝐱[j]
         end
         𝐲[i] = 𝐲ᵢ
     end
@@ -138,7 +149,7 @@ end
 end
 function jvar!(𝐬², 𝐀, x̄)
     @. s² = zero(eltype(𝐬²))
-    @inbounds for i ∈ 1:size(𝐀,2)
+    @inbounds @fastmath for i ∈ 1:size(𝐀,2)
         @simd for j ∈ eachindex(𝐬²)
             δ = 𝐀[j,i] - x̄[j]
             𝐬²[j] += δ*δ
@@ -207,6 +218,19 @@ function randomaccessavx(P, basis, coeffs::Vector{T}) where {T}
     end
     return p
 end
-
+function jlogdettriangle(T::Union{LowerTriangular,UpperTriangular})
+    ld = 0.0
+    @inbounds for n ∈ 1:size(T,1)
+        ld += log(T[n,n])
+    end
+    ld
+end
+function jlogdettriangleavx(T::Union{LowerTriangular,UpperTriangular})
+    ld = 0.0
+    @avx for n ∈ 1:size(T,1)
+        ld += log(T[n,n])
+    end
+    ld
+end
 
 

benchmark/looptests.jl

@@ -15,10 +15,10 @@ makedocs(;
         "Vectorized Convenience Functions" => "vectorized_convenience_functions.md",
         "Future Work" => "future_work.md"
     ],
-    repo="https://github.com/chriselrod/LoopVectorization.jl/blob/{commit}{path}#L{line}",
+    # repo="https://github.com/chriselrod/LoopVectorization.jl/blob/{commit}{path}#L{line}",
     sitename="LoopVectorization.jl",
-    authors="Chris Elrod",
-    assets=[],
+    authors="Chris Elrod"
+    # assets=[],
 )
 
 deploydocs(;