Update plotting script

Author: YingboMa

perf/lu.jl

@@ -1,5 +1,5 @@
 using BenchmarkTools
-import LinearAlgebra, RecursiveFactorization
+using LinearAlgebra, RecursiveFactorization
 
 BenchmarkTools.DEFAULT_PARAMETERS.seconds = 0.08
 
@@ -17,37 +17,39 @@ function luflop(m, n=m; innerflop=2)
 end
 
 bas_mflops = Float64[]
-rec8_mflops = Float64[]
-rec16_mflops = Float64[]
-rec32_mflops = Float64[]
+rec_mflops = Float64[]
 ref_mflops = Float64[]
-ns = 4:32:500
+ns = 4:8:500
 for n in ns
     @info "$n × $n"
     A = rand(n, n)
     bt = @belapsed LinearAlgebra.lu!($(copy(A)))
     push!(bas_mflops, luflop(n)/bt/1e9)
 
-    rt8 = @belapsed RecursiveFactorization.lu!($(copy(A)); blocksize=8)
-    push!(rec8_mflops, luflop(n)/rt8/1e9)
-
-    rt16 = @belapsed RecursiveFactorization.lu!($(copy(A)); blocksize=16)
-    push!(rec16_mflops, luflop(n)/rt16/1e9)
-
-    rt32 = @belapsed RecursiveFactorization.lu!($(copy(A)); blocksize=32)
-    push!(rec32_mflops, luflop(n)/rt32/1e9)
+    rt = @belapsed RecursiveFactorization.lu!($(copy(A)))
+    push!(rec_mflops, luflop(n)/rt/1e9)
 
     ref = @belapsed LinearAlgebra.generic_lufact!($(copy(A)))
     push!(ref_mflops, luflop(n)/ref/1e9)
 end
 
-using Plots
-plt = plot(ns, bas_mflops, legend=:bottomright, lab="OpenBLAS", title="LU Factorization Benchmark", marker=:auto, dpi=150)
-plot!(plt, ns, rec8_mflops, lab="RF8", marker=:auto)
-plot!(plt, ns, rec16_mflops, lab="RF16", marker=:auto)
-plot!(plt, ns, rec32_mflops, lab="RF32", marker=:auto)
+using DataFrames, VegaLite
+df = DataFrame(Size = ns, RecursiveFactorization = rec_mflops, OpenBLAS = bas_mflops, Reference = ref_mflops)
+df = stack(df, [:RecursiveFactorization, :OpenBLAS, :Reference], variable_name = :Library, value_name = :GFLOPS)
+plt = df |> @vlplot(
+        :line, color = :Library,
+        x = {:Size}, y = {:GFLOPS},
+        width = 2400, height = 600
+    )
+save(joinpath(homedir(), "Pictures", "lu_float64.png"), plt)
+
+#=
+using Plot
+plt = plot(ns, bas_mflops, legend=:bottomright, lab="OpenBLAS", title="LU Factorization Benchmark", marker=:auto, dpi=300)
+plot!(plt, ns, rec_mflops, lab="RecursiveFactorization", marker=:auto)
 plot!(plt, ns, ref_mflops, lab="Reference", marker=:auto)
 xaxis!(plt, "size (N x N)")
 yaxis!(plt, "GFLOPS")
 savefig("lubench.png")
 savefig("lubench.pdf")
+=#

src/lu.jl

@@ -1,5 +1,5 @@
 using LoopVectorization
-using LinearAlgebra: BlasInt, BlasFloat, LU, UnitLowerTriangular, ldiv!, mul!, checknonsingular
+using LinearAlgebra: BlasInt, BlasFloat, LU, UnitLowerTriangular, ldiv!, checknonsingular, BLAS
 
 function lu(A::AbstractMatrix, pivot::Union{Val{false}, Val{true}} = Val(true); kwargs...)
     lu!(copy(A), pivot; kwargs...)
@@ -9,13 +9,16 @@ function lu!(A, pivot::Union{Val{false}, Val{true}} = Val(true); kwargs...)
     lu!(A, Vector{BlasInt}(undef, min(size(A)...)), pivot; kwargs...)
 end
 
+# Use a function here to make sure it gets optimized away
+# OpenBLAS' TRSM isn't very good, we use a higher threshold for recursion
+pick_threshold() = BLAS.vendor() === :mkl ? 48 : 192
+
 function lu!(A::AbstractMatrix{T}, ipiv::AbstractVector{<:Integer},
              pivot::Union{Val{false}, Val{true}} = Val(true);
              check::Bool=true,
              # the performance is not sensitive wrt blocksize, and 16 is a good default
              blocksize::Integer=16,
-             # OpenBLAS' TRSM isn't very good, we use a higher threshold for recursion
-             threshold::Integer=BLAS.vendor() === :mkl ? 48 : 192) where T
+             threshold::Integer=pick_threshold()) where T
     info = Ref(zero(BlasInt))
     m, n = size(A)
     mnmin = min(m, n)