Started working on documentation, stop tiler from choosing 0 size, update precompile.

Author: chriselrod

.gitignore

@@ -5,11 +5,9 @@
 *.jl.*.cov
 *.jl.mem
 *~
-src/#*#
-tests/#*#
-benchmark/#*#
 *.mem
 *.mod
 *.mod0
 *.so
 *.s
+*#

benchmark/driver.jl

@@ -35,33 +35,35 @@ Atmulvb_future = @spawnat 14 benchmark_Atmulvb(2:256);
 
 dot_bench = fetch(dot_future)
 selfdot_bench = fetch(selfdot_future)
-AplusAt_bench = fetch(AplusAt_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)
 
-v = 6
-const PICTURES = "/home/chriselrod/Pictures"
-save(joinpath(PICTURES, "bench_gemm_v$v.png"), plot(gemm_bench));
-save(joinpath(PICTURES, "bench_AtmulB_v$v.png"), plot(AtmulB_bench));
-save(joinpath(PICTURES, "bench_dot_v$v.png"), plot(dot_bench));
-save(joinpath(PICTURES, "bench_selfdot_v$v.png"), plot(selfdot_bench));
-save(joinpath(PICTURES, "bench_gemv_v$v.png"), plot(gemv_bench));
-save(joinpath(PICTURES, "bench_dot3_v$v.png"), plot(dot3_bench));
-save(joinpath(PICTURES, "bench_sse_v$v.png"), plot(sse_bench));
-save(joinpath(PICTURES, "bench_exp_v$v.png"), plot(exp_bench));
-save(joinpath(PICTURES, "bench_aplusBc_v$v.png"), plot(aplusBc_bench));
-save(joinpath(PICTURES, "bench_AplusAt_v$v.png"), plot(AplusAt_bench));
-save(joinpath(PICTURES, "bench_random_access_v$v.png"), plot(randomaccess_bench));
-save(joinpath(PICTURES, "bench_AmulBt_v$v.png"), plot(AmulBt_bench));
-save(joinpath(PICTURES, "bench_Atmulvb_v$v.png"), plot(Atmulvb_bench));
+
+v = 1
+filetype = "svg"
+const PICTURES = joinpath(pkgdir("LoopVectorization"), "docs", "src", "assets")
+save(joinpath(PICTURES, "bench_gemm_v$v.$filetype"), plot(gemm_bench));
+save(joinpath(PICTURES, "bench_AtmulB_v$v.$filetype"), plot(AtmulB_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_Atmulvb_v$v.$filetype"), plot(Atmulvb_bench));
+save(joinpath(PICTURES, "bench_exp_v$v.$filetype"), plot(exp_bench));
 
 
 

benchmark/looptests.jl

@@ -1,67 +1,67 @@
 using LoopVectorization, LinearAlgebra
 BLAS.set_num_threads(1)
 
-function jgemm!(C, A, B)
-    C .= 0
-    M, N = size(C); K = size(B,1)
+function jgemm!(𝐂, 𝐀, 𝐁)
+    𝐂 .= 0
+    M, N = size(𝐂); K = size(𝐁,1)
     @inbounds for n ∈ 1:N, k ∈ 1:K
         @simd ivdep for m ∈ 1:M
-            C[m,n] += A[m,k] * B[k,n]
+            𝐂[m,n] += 𝐀[m,k] * 𝐁[k,n]
         end
     end
 end
-@inline function jgemm!(C, Aᵀ::Adjoint, B)
-    A = parent(Aᵀ)
-    @inbounds for n ∈ 1:size(C,2), m ∈ 1:size(C,1)
-        Cₘₙ = zero(eltype(C))
-        @simd ivdep for k ∈ 1:size(A,1)
-            Cₘₙ += A[k,m] * B[k,n]
+@inline function jgemm!(𝐂, 𝐀ᵀ::Adjoint, 𝐁)
+    𝐀 = parent(𝐀ᵀ)
+    @inbounds for n ∈ 1:size(𝐂,2), m ∈ 1:size(𝐂,1)
+        𝐂ₘₙ = zero(eltype(𝐂))
+        @simd ivdep for k ∈ 1:size(𝐀,1)
+            𝐂ₘₙ += 𝐀[k,m] * 𝐁[k,n]
         end
-        C[m,n] = Cₘₙ
+        𝐂[m,n] = 𝐂ₘₙ
     end
 end
-@inline function jgemm!(C, A, Bᵀ::Adjoint)
-    C .= 0
-    B = parent(Bᵀ)
-    M, N = size(C); K = size(B,1)
+@inline function jgemm!(𝐂, 𝐀, 𝐁ᵀ::Adjoint)
+    𝐂 .= 0
+    𝐁 = parent(𝐁ᵀ)
+    M, N = size(𝐂); K = size(𝐁,1)
     @inbounds for k ∈ 1:K, n ∈ 1:N
         @simd ivdep for m ∈ 1:M
-            C[m,n] += A[m,k] * B[n,k]
+            𝐂[m,n] += 𝐀[m,k] * 𝐁[n,k]
         end
     end
 end
-@inline function gemmavx!(C, A, B)
-    @avx for i ∈ 1:size(A,1), j ∈ 1:size(B,2)
-        Cᵢⱼ = zero(eltype(C))
-        for k ∈ 1:size(A,2)
-            Cᵢⱼ += A[i,k] * B[k,j]
+@inline function gemmavx!(𝐂, 𝐀, 𝐁)
+    @avx for m ∈ 1:size(𝐀,1), n ∈ 1:size(𝐁,2)
+        𝐂ₘₙ = zero(eltype(𝐂))
+        for k ∈ 1:size(𝐀,2)
+            𝐂ₘₙ += 𝐀[m,k] * 𝐁[k,n]
         end
-        C[i,j] = Cᵢⱼ
+        𝐂[m,n] = 𝐂ₘₙ
     end
 end
 function jdot(a, b)
-    s = 0.0
+    s = zero(eltype(a))
     @inbounds @simd ivdep for i ∈ eachindex(a, b)
         s += a[i] * b[i]
     end
     s
 end
 function jdotavx(a, b)
-    s = 0.0
+    s = zero(eltype(a))
     @avx for i ∈ eachindex(a, b)
         s += a[i] * b[i]
     end
     s
 end
 function jselfdot(a)
-    s = 0.0
+    s = zero(eltype(a))
     @inbounds @simd ivdep for i ∈ eachindex(a)
         s += a[i] * a[i]
     end
     s
 end
 function jselfdotavx(a)
-    s = 0.0
+    s = zero(eltype(a))
     @avx for i ∈ eachindex(a)
         s += a[i] * a[i]
     end
@@ -96,71 +96,71 @@ function jvexpavx!(b, a)
     end
 end
 function jsvexp(a)
-    s = 0.0
+    s = zero(eltype(a))
     @inbounds for i ∈ eachindex(a)
         s += exp(a[i])
     end
     s
 end
 function jsvexpavx(a)
-    s = 0.0
+    s = zero(eltype(a))
     @avx for i ∈ eachindex(a)
         s += exp(a[i])
     end
     s
 end
-function jgemv!(y, A, x)
-    y .= 0.0
+function jgemv!(y, 𝐀, x)
+    y .= zero(eltype(y))
     @inbounds for j ∈ eachindex(x)
         @simd ivdep for i ∈ eachindex(y)
-            y[i] += A[i,j] * x[j]
+            y[i] += 𝐀[i,j] * x[j]
         end
     end
 end
-@inline function jgemv!(y, Aᵀ::Adjoint, x)
-    A = parent(Aᵀ)
-    @inbounds for i ∈ eachindex(y)
-        yᵢ = 0.0
-        @simd ivdep for j ∈ eachindex(x)
-            yᵢ += A[j,i] * x[j]
+@inline function jgemv!(𝐲, 𝐀ᵀ::Adjoint, 𝐱)
+    𝐀 = parent(𝐀ᵀ)
+    @inbounds for i ∈ eachindex(𝐲)
+        𝐲ᵢ = zero(eltype(𝐲))
+        @simd ivdep for j ∈ eachindex(𝐱)
+            𝐲ᵢ += 𝐀[j,i] * 𝐱[j]
         end
-        y[i] = yᵢ
+        𝐲[i] = 𝐲ᵢ
     end
 end
-@inline function jgemvavx!(y, A, x)
-    @avx for i ∈ eachindex(y)
-        yᵢ = 0.0
-        for j ∈ eachindex(x)
-            yᵢ += A[i,j] * x[j]
+@inline function jgemvavx!(𝐲, 𝐀, 𝐱)
+    @avx for i ∈ eachindex(𝐲)
+        𝐲ᵢ = zero(eltype(𝐲))
+        for j ∈ eachindex(𝐱)
+            𝐲ᵢ += 𝐀[i,j] * 𝐱[j]
         end
-        y[i] = yᵢ
+        𝐲[i] = 𝐲ᵢ
     end
 end
-function jvar!(s², A, x̄)
-    @. s² = 0
-    @inbounds for i ∈ 1:size(A,2)
-        @simd for j ∈ eachindex(s²)
-            δ = A[j,i] - x̄[j]
-            s²[j] += δ*δ
+function jvar!(𝐬², 𝐀, x̄)
+    @. s² = zero(eltype(𝐬²))
+    @inbounds for i ∈ 1:size(𝐀,2)
+        @simd for j ∈ eachindex(𝐬²)
+            δ = 𝐀[j,i] - x̄[j]
+            𝐬²[j] += δ*δ
         end
     end
 end
-function jvaravx!(s², A, x̄)
-    @avx for j ∈ eachindex(s²)
-        s²ⱼ = 0.0
+function jvaravx!(𝐬², 𝐀, x̄)
+    @avx for j ∈ eachindex(𝐬²)
+        𝐬²ⱼ = zero(eltype(𝐬²))
         x̄ⱼ = x̄[j]
-        for i ∈ 1:size(A,2)
-            δ = A[j,i] - x̄ⱼ
-            s²ⱼ += δ*δ
+        for i ∈ 1:size(𝐀,2)
+            δ = 𝐀[j,i] - x̄ⱼ
+            𝐬²ⱼ += δ*δ
         end
-        s²[j] = s²ⱼ
+        𝐬²[j] = 𝐬²ⱼ
     end
 end
 japlucBc!(d, a, B, c) =      @. d = a + B * c';
 japlucBcavx!(d, a, B, c) = @avx @. d = a + B * c';
 
 function jOLSlp(y, X, β)
-    lp = 0.0
+    lp = zero(eltype(y))
     @inbounds @fastmath for i ∈ eachindex(y)
         δ = y[i]
         @simd for j ∈ eachindex(β)
@@ -171,7 +171,7 @@ function jOLSlp(y, X, β)
     lp
 end
 function jOLSlp_avx(y, X, β)
-    lp = 0.0
+    lp = zero(eltype(y))
     @avx for i ∈ eachindex(y)
         δ = y[i]
         for j ∈ eachindex(β)

docs/make.jl

@@ -5,6 +5,15 @@ makedocs(;
     format=Documenter.HTML(),
     pages=[
         "Home" => "index.md",
+        "Getting Started" => "getting_started.md",
+        "Examples" => Any[
+            "examples/matrix_multiplication.md",
+            "examples/matrix_vector_ops.md",
+            "examples/dot_product.md",
+            "examples/sum_of_squared_error.md"
+        ],
+        "Vectorized Convenience Functions" => "vectorized_convenience_functions.md",
+        "Future Work" => "future_work.md"
     ],
     repo="https://github.com/chriselrod/LoopVectorization.jl/blob/{commit}{path}#L{line}",
     sitename="LoopVectorization.jl",

docs/src/assets/bench_AmulBt_v1.svg

@@ -0,0 +1,27 @@
+# Dot Products
+
+Dot products are so simple, it is almost surprising that compilers leave any performance on the table.
+
+```julia
+function jdotavx(a, b)
+    s = zero(eltype(a))
+    @avx for i ∈ eachindex(a, b)
+        s += a[i] * b[i]
+    end
+    s
+end
+```
+![dot](../assets/bench_dot_v1.svg)
+
+
+```julia
+function jselfdotavx(a)
+    s = zero(eltype(a))
+    @avx for i ∈ eachindex(a)
+        s += a[i] * a[i]
+    end
+    s
+end
+```
+![selfdot](../assets/bench_selfdot_v1.svg)
+

docs/src/assets/bench_AplusAt_v1.svg

@@ -0,0 +1,25 @@
+# Matrix Multiplication
+
+```julia
+@inline function A_mul_B!(𝐂, 𝐀, 𝐁)
+    @avx for m ∈ 1:size(𝐀,1), n ∈ 1:size(𝐁,2)
+        𝐂ₘₙ = zero(eltype(𝐂))
+        for k ∈ 1:size(𝐀,2)
+            𝐂ₘₙ += 𝐀[m,k] * 𝐁[k,n]
+        end
+        𝐂[m,n] = 𝐂ₘₙ
+    end
+end
+```
+
+Letting all three matrices be square and `Size` x `Size`, we attain the following benchmark results:
+
+![AmulB](../assets/bench_gemm_v1.svg)
+
+
+![AtmulB](../assets/bench_AtmulB_v1.svg)
+
+
+![AmulBt](../assets/bench_AmulBt_v1.svg)
+
+

docs/src/assets/bench_AtmulB_v1.svg

@@ -0,0 +1,27 @@
+# Matrix-Vector Operations
+
+Here I'll discuss a variety of Matrix-vector operations, naturally starting with matrix-vector multiplication.
+
+```julia
+@inline function jgemvavx!(𝐲, 𝐀, 𝐱)
+    @avx for i ∈ eachindex(𝐲)
+        𝐲ᵢ = zero(eltype(𝐲))
+        for j ∈ eachindex(𝐱)
+            𝐲ᵢ += 𝐀[i,j] * 𝐱[j]
+        end
+        𝐲[i] = 𝐲ᵢ
+    end
+end
+```
+
+Using a square `Size` x `Size` matrix `A`, we find the following results.
+![Amulvb](../assets/bench_gemv_v1.svg)
+
+
+![Atmulvb](../assets/bench_Atmulvb_v1.svg)
+
+
+![dot3](../assets/bench_dot3_v1.svg)
+
+
+

docs/src/assets/bench_Atmulvb_v1.svg

@@ -0,0 +1,26 @@
+# Sum of squared error
+
+To calculate `(y - X * β)'(y - X * β)`, we can use the following loop.
+```julia
+function sse_avx(y, X, β)
+    lp = zero(eltype(y))
+    @avx for i ∈ eachindex(y)
+        δ = y[i]
+        for j ∈ eachindex(β)
+            δ -= X[i,j] * β[j]
+        end
+        lp += δ * δ
+    end
+    lp
+end
+```
+
+LoopVectorization does not model memory access yet.
+That seems important for this example, where performance starts to decline for sizes larger than 60.
+Letting `N` be the size, `X` was a `3N/2`x `N/2` matrix. Therefore, performance started to suffer
+when `X` had more than about 30 columns (performance is much less sensitive to the number of rows).
+
+![sse](../assets/bench_sse_v1.svg)
+
+
+