Small fixes and benchmark for README

Andy Ferris · Andy Ferris · commit 9562d0958ffb · 2016-11-10T15:53:12.000+10:00
Probably should generate a full benchmark when I have time
diff --git a/README.md b/README.md
@@ -29,30 +29,63 @@ than `Base.Array`. See this sample benchmark (or see the full results [here](htt
 
 ```
 =====================================
-   Benchmarks for 3×3 matrices
+    Benchmarks for 3×3 matrices
 =====================================
 
-Matrix multiplication
----------------------
-Array  ->  3.973188 seconds (74.07 M allocations: 6.623 GB, 12.92% gc time)
-SArray ->  0.326989 seconds (5 allocations: 240 bytes)
-MArray ->  2.248258 seconds (37.04 M allocations: 2.759 GB, 14.06% gc time)
-
-Matrix multiplication (mutating)
---------------------------------
-Array  ->  2.237091 seconds (6 allocations: 480 bytes)
-MArray ->  0.795372 seconds (6 allocations: 320 bytes)
-
-Matrix addition
----------------
-Array  ->  2.610709 seconds (44.44 M allocations: 3.974 GB, 11.81% gc time)
-SArray ->  0.073024 seconds (5 allocations: 240 bytes)
-MArray ->  0.896849 seconds (22.22 M allocations: 1.656 GB, 21.33% gc time)
-
-Matrix addition (mutating)
---------------------------
-Array  ->  0.872791 seconds (6 allocations: 480 bytes)
-MArray ->  0.145895 seconds (5 allocations: 240 bytes)
+Matrix multiplication (8.2x speedup)
+------------------------------------
+Array               ->  2.024568 seconds (74.07 M allocations: 6.623 GB, 9.37% gc time)
+SArray              ->  0.247364 seconds (5 allocations: 240 bytes)
+MArray              ->  1.603798 seconds (37.04 M allocations: 2.759 GB, 13.61% gc time)
+SizedArray          ->  2.223853 seconds (74.07 M allocations: 6.071 GB, 10.05% gc time)
+
+Matrix multiplication (mutating) (3.1x speedup)
+-----------------------------------------------
+Array               ->  1.360940 seconds (6 allocations: 480 bytes)
+MArray              ->  0.443528 seconds (7 allocations: 400 bytes)
+SizedArray          ->  0.681896 seconds (7 allocations: 416 bytes)
+
+Matrix addition (45x speedup)
+-----------------------------
+Array               ->  1.458899 seconds (44.44 M allocations: 3.974 GB, 7.96% gc time)
+SArray              ->  0.032043 seconds (5 allocations: 240 bytes)
+MArray              ->  0.682318 seconds (22.22 M allocations: 1.656 GB, 18.70% gc time)
+SizedArray          ->  1.111785 seconds (44.44 M allocations: 3.643 GB, 12.02% gc time)
+
+Matrix addition (mutating) (5.1x speedup)
+-----------------------------------------
+Array      ->  0.493796 seconds (5 allocations: 320 bytes)
+MArray     ->  0.096303 seconds (5 allocations: 240 bytes)
+SizedArray ->  0.135803 seconds (6 allocations: 336 bytes)
+
+Matrix determinant (170x speedup)
+---------------------------------
+Array      -> 15.291557 seconds (222.22 M allocations: 12.694 GB, 16.97% gc time)
+SArray     ->  0.094409 seconds (4 allocations: 160 bytes)
+MArray     ->  0.089569 seconds (4 allocations: 160 bytes)
+SizedArray ->  0.114134 seconds (4 allocations: 160 bytes)
+
+Matrix inverse (125x speedup)
+-----------------------------
+Array      -> 47.704314 seconds (407.41 M allocations: 82.232 GB, 23.67% gc time)
+SArray     ->  0.379657 seconds (4 allocations: 160 bytes)
+MArray     ->  1.294672 seconds (37.04 M allocations: 2.759 GB, 19.77% gc time)
+SizedArray ->  2.136363 seconds (74.07 M allocations: 6.071 GB, 8.61% gc time)
+
+Matrix symmetric eigenvalue (105x speedup)
+------------------------------------------
+Array      -> 418.304283 seconds (740.74 M allocations: 89.407 GB, 1.22% gc time)
+SArray     ->   3.963118 seconds (5 allocations: 256 bytes)
+MArray     ->   3.964029 seconds (6 allocations: 272 bytes)
+SizedArray ->   4.028497 seconds (6 allocations: 208 bytes)
+
+Matrix Cholesky (23.6x speedup)
+-------------------------------
+Array      ->  8.139431 seconds (222.22 M allocations: 9.934 GB, 6.28% gc time)
+SArray     ->  0.344283 seconds (5 allocations: 256 bytes)
+MArray     ->  0.812532 seconds (37.04 M allocations: 2.759 GB, 6.19% gc time)
+SizedArray ->  2.225999 seconds (74.07 M allocations: 6.071 GB, 11.41% gc time)
+
 ```
 
 (Run with `julia -O3` for even faster SIMD code with immutable static arrays!)
diff --git a/perf/benchmark2.jl b/perf/benchmark2.jl
@@ -16,8 +16,10 @@ M_g = div(2*10^8, N^2)
 # Size
 
 A = rand(Float64,N,N)
+A = A*A'
 As = SMatrix{N,N}(A)
 Am = MMatrix{N,N}(A)
+Az = Size(N,N)(copy(A))
 @static if fsa
     Af = Mat(ntuple(j -> ntuple(i->A[i,j], N), N)) # there is a bug in FixedSizeArrays Mat constructor (13 July 2016)
 end
@@ -36,9 +38,19 @@ if !isdefined(:f_mut_marray) || !isdefined(:benchmark_suite) || benchmark_suite
     @generated f_blas_marray(n::Integer, A) = :(@inbounds (C = similar(A); C[:] = A[:]; tmp = similar(A); for i = 1:n; StaticArrays.A_mul_B_blas!(tmp, C, A); C.data = tmp.data; end; return C))
 
     @generated g(n::Integer, A) = :(@inbounds (C = A; for i = 1:n; C = C + A; end; return C))
-    @generated g_mut(n::Integer, A) = :(@inbounds (C = similar(A); C[:] = A[:]; for i = 1:n; @inbounds map!(+, C, C, A); end; return C))
+    @generated g_mut(n::Integer, A) = :(@inbounds (C = copy(A); for i = 1:n; @inbounds map!(+, C, C, A); end; return C))
     @generated g_via_sarray{M}(n::Integer, A::MMatrix{M,M}) = :(@inbounds (C = similar(A); C[:] = A[:]; for i = 1:n; C = MMatrix{M,M}(SMatrix{M,M}(C) + SMatrix{M,M}(A)); end; return C))
 
+    @noinline _det(x) = det(x)
+    @noinline _inv(x) = inv(x)
+    @noinline _eig(x) = eig(x)
+    @noinline _chol(x) = chol(x)
+
+    f_det(n::Int, A) = (for i = 1:n; _det(A); end)
+    f_inv(n::Int, A) = (for i = 1:n; _inv(A); end)
+    f_eig(n::Int, A) = (for i = 1:n; _eig(A); end)
+    f_chol(n::Int, A) = (for i = 1:n; _chol(A); end)
+
     # Notes: - A[:] = B[:] is allocating in Base, unlike `map!`
     #        - Also, the same goes for Base's implementation of broadcast!(f, A, B, C) (but map! is OK).
     #        - I really need to implement copy() in StaticArrays... (and maybe a special instance of setindex!(C, :))
@@ -122,11 +134,14 @@ end
 # Warmup and some checks
 Cs = f(2, As)
 Cm = f(2, Am)
+Cz = f(2, Az)
 
 Cs::SMatrix
 if N <= 4; @assert Cs ≈ C; end
 Cm::MMatrix
 if N <= 4; @assert Cm ≈ C; end
+Cz::SizedMatrix
+if N <= 4; @assert Cz ≈ C; end
 
 @static if fsa
     @static if all_methods
@@ -149,6 +164,10 @@ Cm_mut = f_mut_marray(2, Am)
 Cm_mut::MMatrix
 if N <= 4; @assert Cm_mut ≈ C; end
 
+Cz_mut = f_mut_array(2, Az)
+Cz_mut::SizedMatrix
+if N <= 4; @assert Cz_mut ≈ C; end
+
 @static if all_methods
     println()
     print("A_mul_B!(MMatrix, MMatrix) compilation time (unrolled):")
@@ -196,13 +215,19 @@ if N <= 4; @assert C_mut ≈ C; end
 Cs = g(2, As)
 Cm = g(2, Am)
 Cm_mut = g_mut(2, Am)
+Cz = g(2, Az)
+Cz_mut = g_mut(2, Az)
 
 Cs::SMatrix
 if N <= 4; @assert Cs == C; end
 Cm::MMatrix
 if N <= 4; @assert Cm == C; end
 Cm_mut::MMatrix
 if N <= 4; @assert Cm_mut == C; end
+Cz::SizedMatrix
+if N <= 4; @assert Cz == C; end
+Cz_mut::SizedMatrix
+if N <= 4; @assert Cz_mut == C; end
 
 @static if all_methods
     Cm_via_sarray = g_via_sarray(2, Am)
@@ -216,6 +241,29 @@ end
     if N <= 4; @assert Cf == C; end
 end
 
+if N <= 3
+    # det, eig etc
+    C = f_det(2, A)
+    Cs = f_det(2, As)
+    Cm = f_det(2, Am)
+    Cz = f_det(2, Az)
+
+    C = f_inv(2, A)
+    Cs = f_inv(2, As)
+    Cm = f_inv(2, Am)
+    Cz = f_inv(2, Az)
+
+    C = f_eig(2, Symmetric(A))
+    Cs = f_eig(2, Symmetric(As))
+    Cm = f_eig(2, Symmetric(Am))
+    Cz = f_eig(2, Symmetric(Az))
+
+    C = f_chol(2, Symmetric(A))
+    Cs = f_chol(2, Symmetric(As))
+    Cm = f_chol(2, Symmetric(Am))
+    Cz = f_chol(2, Symmetric(Az))
+end
+
 println()
 
 # Do the performance tests
@@ -229,6 +277,7 @@ begin
    end
    print("SArray              ->"); @time f(M_f, As)
    print("MArray              ->"); @time f(M_f, Am)
+   print("SizedArray          ->"); @time f(M_f, Az)
    @static if all_methods
        print("SArray (unrolled)   ->"); @time f_unrolled(M_f, As)
        print("SArray (chunks)     ->"); @time f_unrolled_chunks(M_f, As)
@@ -246,6 +295,7 @@ println("--------------------------------")
 begin
    print("Array               ->"); @time f_mut_array(M_f, A)
    print("MArray              ->"); @time f_mut_marray(M_f, Am)
+   print("SizedArray          ->"); @time f_mut_array(M_f, Az)
    @static if all_methods
        print("MArray (unrolled)   ->"); @time f_mut_unrolled(M_f, Am)
        print("MArray (chunks)     ->"); @time f_mut_chunks(M_f, Am)
@@ -263,6 +313,7 @@ begin
    end
    print("SArray              ->"); @time g(M_g, As)
    print("MArray              ->"); @time g(M_g, Am)
+   print("SizedArray          ->"); @time g(M_g, Az)
    @static if all_methods
        print("MArray (via SArray) ->"); @time g_via_sarray(M_g, Am)
    end
@@ -272,7 +323,50 @@ println()
 println("Matrix addition (mutating)")
 println("--------------------------")
 begin
-   print("Array  ->"); @time g_mut(M_g, A) # broadcast! seems to be broken!
-   print("MArray ->"); @time g_mut(M_g, Am)
+   print("Array      ->"); @time g_mut(M_g, A) # broadcast! seems to be broken!
+   print("MArray     ->"); @time g_mut(M_g, Am)
+   print("SizedArray ->"); @time g_mut(M_g, Az)
 end
 println()
+
+if N <= 3
+    println("Matrix determinant")
+    println("------------------")
+    begin
+       print("Array      ->"); @time f_det(M_f, A)
+       print("SArray     ->"); @time f_det(M_f, As)
+       print("MArray     ->"); @time f_det(M_f, Am)
+       print("SizedArray ->"); @time f_det(M_f, Az)
+    end
+    println()
+
+    println("Matrix inverse")
+    println("--------------")
+    begin
+       print("Array      ->"); @time f_inv(M_f, A)
+       print("SArray     ->"); @time f_inv(M_f, As)
+       print("MArray     ->"); @time f_inv(M_f, Am)
+       print("SizedArray ->"); @time f_inv(M_f, Az)
+    end
+    println()
+
+    println("Matrix symmetric eigenvalue")
+    println("---------------------------")
+    begin
+       print("Array      ->"); @time f_eig(M_f, Symmetric(A))
+       print("SArray     ->"); @time f_eig(M_f, Symmetric(As))
+       print("MArray     ->"); @time f_eig(M_f, Symmetric(Am))
+       print("SizedArray ->"); @time f_eig(M_f, Symmetric(Az))
+    end
+    println()
+
+    println("Matrix Cholesky")
+    println("---------------")
+    begin
+       print("Array      ->"); @time f_chol(M_f, Symmetric(A))
+       print("SArray     ->"); @time f_chol(M_f, Symmetric(As))
+       print("MArray     ->"); @time f_chol(M_f, Symmetric(Am))
+       print("SizedArray ->"); @time f_chol(M_f, Symmetric(Az))
+    end
+    println()
+end #if N <= 3
diff --git a/src/cholesky.jl b/src/cholesky.jl
@@ -5,14 +5,12 @@
 end
 
 @inline function Base.chol{T<:Real, SM <: StaticMatrix}(A::Base.LinAlg.RealHermSymComplexHerm{T,SM})
-    ishermitian(A) || Base.LinAlg.non_hermitian_error("chol")
-    _chol(Size(A), A)
-end
-
-@inline function Base.chol{SM<:StaticMatrix}(A::Symmetric{SM})
-    eltype(A) <: Real && (ishermitian(A) || Base.LinAlg.non_hermitian_error("chol"))
-    _chol(Size(A), A)
+    _chol(Size(A), A.data)
 end
+#=
+@inline function Base.chol{T<:Real,SM<:StaticMatrix}(A::Symmetric{T,SM})
+    _chol(Size(A), A.data)
+end=#
 
 @generated function _chol(::Size{(1,1)}, A::StaticMatrix)
     @assert size(A) == (1,1)
diff --git a/src/linalg.jl b/src/linalg.jl
