Add benchmark suite

.JuliaFormatter.toml

@@ -1,2 +1,3 @@
 # See https://domluna.github.io/JuliaFormatter.jl/stable/ for a list of options
 style = "blue"
+margin = 140

.gitignore

@@ -1,6 +1,10 @@
 *.jl.*.cov
 *.jl.cov
 *.jl.mem
-/Manifest.toml
+Manifest.toml
 /docs/Manifest.toml
 /docs/build/
+
+.vscode
+benchmark
+benchmarks/**/exports

Makefile

@@ -0,0 +1,22 @@
+SHELL = /bin/bash
+.DEFAULT_GOAL = help
+
+
+.PHONY: test
+
+test: ## Run tests
+	julia -e 'import Pkg; Pkg.activate("."); Pkg.test();'
+
+.PHONY: benchmark
+
+benchmark: ## Run benchmarks
+	julia -e '\
+    import Pkg; \
+    Pkg.activate("benchmarks/speed"); \
+    Pkg.develop(Pkg.PackageSpec(path=pwd())); \
+    Pkg.instantiate(); \
+    using PkgBenchmark, FastCholesky, Dates; \
+    Base.Filesystem.mkpath("benchmarks/speed/exports/"); \
+    results = benchmarkpkg("FastCholesky"; script="benchmarks/speed/benchmarks.jl"); \
+    export_markdown("benchmarks/speed/exports/benchmark-"*Dates.format(Dates.now(), "yyyy-mm-ddTHH-MM")*".md", results); \
+    Base.Filesystem.rm("benchmark"; force=true, recursive=true)'

benchmarks/speed/Project.toml

@@ -0,0 +1,5 @@
+[deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+FastCholesky = "2d5283b6-8564-42b6-bb00-83ed8e915756"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+PkgBenchmark = "32113eaa-f34f-5b0d-bd6c-c81e245fc73d"

benchmarks/speed/benchmarks.jl

@@ -0,0 +1,74 @@
+using BenchmarkTools, LinearAlgebra
+using FastCholesky
+
+const SUITE = BenchmarkGroup()
+const BenchmarkFloatTypes = (Float32, Float64, BigFloat)
+
+for T in BenchmarkFloatTypes
+    SUITE[string(T)] = BenchmarkGroup()
+    SUITE[string(T)]["Number"] = BenchmarkGroup()
+    SUITE[string(T)]["Matrix"] = BenchmarkGroup()
+    SUITE[string(T)]["Diagonal"] = BenchmarkGroup()
+end
+
+# general benchmarks
+for T in BenchmarkFloatTypes
+
+    # number
+    a = rand(T)
+
+    SUITE[string(T)]["Number"]["fastcholesky"] = @benchmarkable fastcholesky($a)
+    SUITE[string(T)]["Number"]["fastcholesky!"] = @benchmarkable fastcholesky!($a)
+    SUITE[string(T)]["Number"]["cholinv"] = @benchmarkable cholinv($a)
+    SUITE[string(T)]["Number"]["cholsqrt"] = @benchmarkable cholsqrt($a)
+    SUITE[string(T)]["Number"]["chollogdet"] = @benchmarkable chollogdet($a)
+    SUITE[string(T)]["Number"]["cholinv_logdet"] = @benchmarkable cholinv_logdet($a)
+
+    for dim in (2, 5, 10, 20, 50, 100, 200, 500)
+
+        SUITE[string(T)]["Matrix"]["dim:" * string(dim)] = BenchmarkGroup()
+
+        # Skip `BigFloat` benchmarks for large dimensions because they are extremely slow
+        if (T == BigFloat) && (dim >= 100)
+            continue
+        end
+
+        # generate positive-definite matrix of specified dimensions
+        A = randn(T, dim, dim)
+        B = A * A' + dim * I(dim)
+        C = similar(B)
+
+        # The inplace version works only for `Matrix{<:BlasFloat}`
+        if T <: LinearAlgebra.BlasFloat
+            SUITE[string(T)]["Matrix"]["dim:" * string(dim)]["fastcholesky!"] = @benchmarkable fastcholesky!($C)
+        end
+
+        # `sqrt` does not work for `BigFloat` inputs
+        if T != BigFloat
+            SUITE[string(T)]["Matrix"]["dim:" * string(dim)]["cholsqrt"] = @benchmarkable cholsqrt($B)
+        end
+
+        # define benchmarks
+        SUITE[string(T)]["Matrix"]["dim:" * string(dim)]["fastcholesky"] = @benchmarkable fastcholesky($B)
+        SUITE[string(T)]["Matrix"]["dim:" * string(dim)]["cholinv"] = @benchmarkable cholinv($B)
+        SUITE[string(T)]["Matrix"]["dim:" * string(dim)]["chollogdet"] = @benchmarkable chollogdet($B)
+        SUITE[string(T)]["Matrix"]["dim:" * string(dim)]["cholinv_logdet"] = @benchmarkable cholinv_logdet($B)
+
+        SUITE[string(T)]["Diagonal"]["dim:" * string(dim)] = BenchmarkGroup()
+
+        # Diagonal{T} benchmarks
+        SUITE[string(T)]["Diagonal"]["dim:" * string(dim)]["fastcholesky"]= @benchmarkable fastcholesky($(Diagonal(ones(T, dim))))
+        SUITE[string(T)]["Diagonal"]["dim:" * string(dim)]["cholinv"] = @benchmarkable cholinv($(Diagonal(ones(T, dim))))
+        SUITE[string(T)]["Diagonal"]["dim:" * string(dim)]["cholsqrt"] = @benchmarkable cholsqrt($(Diagonal(ones(T, dim))))
+        SUITE[string(T)]["Diagonal"]["dim:" * string(dim)]["chollogdet"] = @benchmarkable chollogdet($(Diagonal(ones(T, dim))))
+        SUITE[string(T)]["Diagonal"]["dim:" * string(dim)]["cholinv_logdet"] = @benchmarkable cholinv_logdet($(Diagonal(ones(T, dim))))
+    end
+end
+
+SUITE["UniformScaling"] = BenchmarkGroup()
+
+# Uniformscaling benchmarks
+SUITE["UniformScaling"]["fastcholesky"] = @benchmarkable fastcholesky($(3.0I))
+SUITE["UniformScaling"]["fastcholesky!"] = @benchmarkable fastcholesky!($(3.0I))
+SUITE["UniformScaling"]["cholinv"] = @benchmarkable cholinv($(3.0I))
+SUITE["UniformScaling"]["cholsqrt"] = @benchmarkable cholsqrt($(3.0I))

src/FastCholesky.jl

@@ -32,9 +32,10 @@ function fastcholesky(input::AbstractMatrix)
     # The `PositiveFactorizations.default_δ` should small enough in majority of the cases
     return cholesky(PositiveFactorizations.Positive, Hermitian(input), tol = PositiveFactorizations.default_δ(input))
 end
+
+fastcholesky(input::Number) = cholesky(input)
 fastcholesky(input::Diagonal) = cholesky(input)
 fastcholesky(input::Hermitian) = cholesky(PositiveFactorizations.Positive, input)
-fastcholesky(input::Number) = input
 
 function fastcholesky(x::UniformScaling)
     return error(

test/fastcholesky_tests.jl

@@ -28,6 +28,10 @@
                     @test all(cholinv_logdet(input) .≈ (inv(input), logdet(input)))
                     @test cholsqrt(input) * cholsqrt(input)' ≈ input
 
+                    if Type <: LinearAlgebra.BlasFloat && input isa Matrix 
+                        @test collect(fastcholesky(input).L) ≈ collect(fastcholesky!(deepcopy(input)).L)
+                    end
+
                     # `sqrt` does not work on `BigFloat` matrices
                     if Type !== BigFloat
                         @test cholsqrt(input) * cholsqrt(input)' ≈ sqrt(input) * sqrt(input)'
@@ -69,7 +73,8 @@ end
 
     for Type in SupportedTypes
         let number = rand(Type)
-            @test fastcholesky(number) === number
+            @test size(fastcholesky(number).L) == (1, 1)
+            @test all(fastcholesky(number).L .≈ sqrt(number))
             @test cholinv(number) ≈ inv(number)
             @test cholsqrt(number) ≈ sqrt(number)
             @test chollogdet(number) ≈ logdet(number)
@@ -91,4 +96,4 @@ end
     @test inv(fastcholesky(F)) * F ≈ Diagonal(ones(4)) rtol=1e-4
     @test cholinv(F) * F ≈ Diagonal(ones(4)) rtol=1e-4
     @test fastcholesky(F).L ≈ cholesky(F).L
-end
+end