From 0859d674582a938820909e27b3cfdb8b1410b93d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Mon, 14 Aug 2023 07:58:19 +0200 Subject: [PATCH 1/2] Add sparselu benchmark So far with direct solvers. --- benchmarks/sparselu.jl | 93 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 93 insertions(+) create mode 100644 benchmarks/sparselu.jl diff --git a/benchmarks/sparselu.jl b/benchmarks/sparselu.jl new file mode 100644 index 000000000..4b9ab049c --- /dev/null +++ b/benchmarks/sparselu.jl @@ -0,0 +1,93 @@ +using BenchmarkTools, Random, VectorizationBase +using LinearAlgebra, SparseArrays, LinearSolve, Sparspak +import Pardiso +using Plots + +BenchmarkTools.DEFAULT_PARAMETERS.seconds = 0.5 + +# Why do I need to set this ? +BenchmarkTools.DEFAULT_PARAMETERS.samples = 10 + + + +# Sparse matrix generation on a n-dimensional rectangular grid. After +# https://discourse.julialang.org/t/seven-lines-of-julia-examples-sought/50416/135 +# by A. Braunstein. + +A ⊕ B = kron(I(size(B, 1)), A) + kron(B, I(size(A, 1))) + +function lattice(n; Tv = Float64) + d = fill(2 * one(Tv), n) + d[1] = one(Tv) + d[end] = one(Tv) + spdiagm(1 => -ones(Tv, n - 1), 0 => d, -1 => -ones(Tv, n - 1)) +end + +lattice(L...; Tv = Float64) = lattice(L[1]; Tv) ⊕ lattice(L[2:end]...; Tv) + +# +# Create a matrix similar to that of a finite difference discretization in a `dim`-dimensional +# unit cube of ``-Δu + δu`` with approximately N unknowns. It is strictly diagonally dominant. +# +function fdmatrix(N; dim=2, Tv = Float64, δ = 1.0e-2) + n = N^(1 / dim) |> ceil |> Int + lattice([n for i in 1:dim]...; Tv) + Tv(δ) * I +end + + +algs = [UMFPACKFactorization(),KLUFactorization(),MKLPardisoFactorize(),SparspakFactorization()] +cols= [:red, :blue, :green, :magenta,:turqoise] # one color per alg +lst=[:dash,:solid, :dashdot] # one line style per dim + +__parameterless_type(T) = Base.typename(T).wrapper +parameterless_type(x) = __parameterless_type(typeof(x)) +parameterless_type(::Type{T}) where {T} = __parameterless_type(T) + +# +# kmax=12 gives ≈ 40_000 unknowns max, can be watched in real time +# kmax=15 gives ≈ 328_000 unknows, you can go make a coffee. +# Main culprit is KLU factorization in 3D. +# +function run_and_plot(;dims=[1,2,3], kmax=12) + + ns = [10 * 2^k for k in 0:kmax] + + res = [ [Float64[] for i in 1:length(algs)] for dim in dims ] + + for dim in dims + for i in 1:length(ns) + rng = MersenneTwister(123) + A = fdmatrix(ns[i];dim) + n = size(A,1) + @info "dim=$(dim): $n × $n" + b = rand(rng, n) + u0= rand(rng, n) + + for j in 1:length(algs) + bt = @belapsed solve(prob, $(algs[j])).u setup=(prob = LinearProblem(copy($A), copy($b); u0 = copy($u0), alias_A=true, alias_b=true)) + push!(res[dim][j], bt) + end + end + end + + p = plot(; + ylabel = "Time/s", + xlabel = "N", + yscale=:log10, + xscale=:log10, + title = "Time for NxN sparse LU Factorization", + label = string(Symbol(parameterless_type(algs[1]))), + legend=:outertopright) + + for dim in dims + for i in 1:length(algs) + plot!(p, ns, res[dim][i]; + linecolor=cols[i], + linestyle= lst[dim], + label = "$(string(Symbol(parameterless_type(algs[i])))) $(dim)D" ) + end + end + savefig("sparselubench.png") + savefig("sparselubench.pdf") +end + From 8c1c1755940a1b1a29bc0ae0de40bfd4e26ae31c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Mon, 14 Aug 2023 08:37:17 +0200 Subject: [PATCH 2/2] format --- benchmarks/sparselu.jl | 62 ++++++++++++++++++++++-------------------- 1 file changed, 33 insertions(+), 29 deletions(-) diff --git a/benchmarks/sparselu.jl b/benchmarks/sparselu.jl index 4b9ab049c..db917b1cc 100644 --- a/benchmarks/sparselu.jl +++ b/benchmarks/sparselu.jl @@ -8,8 +8,6 @@ BenchmarkTools.DEFAULT_PARAMETERS.seconds = 0.5 # Why do I need to set this ? BenchmarkTools.DEFAULT_PARAMETERS.samples = 10 - - # Sparse matrix generation on a n-dimensional rectangular grid. After # https://discourse.julialang.org/t/seven-lines-of-julia-examples-sought/50416/135 # by A. Braunstein. @@ -29,15 +27,19 @@ lattice(L...; Tv = Float64) = lattice(L[1]; Tv) ⊕ lattice(L[2:end]...; Tv) # Create a matrix similar to that of a finite difference discretization in a `dim`-dimensional # unit cube of ``-Δu + δu`` with approximately N unknowns. It is strictly diagonally dominant. # -function fdmatrix(N; dim=2, Tv = Float64, δ = 1.0e-2) +function fdmatrix(N; dim = 2, Tv = Float64, δ = 1.0e-2) n = N^(1 / dim) |> ceil |> Int lattice([n for i in 1:dim]...; Tv) + Tv(δ) * I end - -algs = [UMFPACKFactorization(),KLUFactorization(),MKLPardisoFactorize(),SparspakFactorization()] -cols= [:red, :blue, :green, :magenta,:turqoise] # one color per alg -lst=[:dash,:solid, :dashdot] # one line style per dim +algs = [ + UMFPACKFactorization(), + KLUFactorization(), + MKLPardisoFactorize(), + SparspakFactorization(), +] +cols = [:red, :blue, :green, :magenta, :turqoise] # one color per alg +lst = [:dash, :solid, :dashdot] # one line style per dim __parameterless_type(T) = Base.typename(T).wrapper parameterless_type(x) = __parameterless_type(typeof(x)) @@ -48,46 +50,48 @@ parameterless_type(::Type{T}) where {T} = __parameterless_type(T) # kmax=15 gives ≈ 328_000 unknows, you can go make a coffee. # Main culprit is KLU factorization in 3D. # -function run_and_plot(;dims=[1,2,3], kmax=12) - +function run_and_plot(; dims = [1, 2, 3], kmax = 12) ns = [10 * 2^k for k in 0:kmax] - - res = [ [Float64[] for i in 1:length(algs)] for dim in dims ] - + + res = [[Float64[] for i in 1:length(algs)] for dim in dims] + for dim in dims for i in 1:length(ns) rng = MersenneTwister(123) - A = fdmatrix(ns[i];dim) - n = size(A,1) + A = fdmatrix(ns[i]; dim) + n = size(A, 1) @info "dim=$(dim): $n × $n" b = rand(rng, n) - u0= rand(rng, n) - + u0 = rand(rng, n) + for j in 1:length(algs) - bt = @belapsed solve(prob, $(algs[j])).u setup=(prob = LinearProblem(copy($A), copy($b); u0 = copy($u0), alias_A=true, alias_b=true)) + bt = @belapsed solve(prob, $(algs[j])).u setup=(prob = LinearProblem(copy($A), + copy($b); + u0 = copy($u0), + alias_A = true, + alias_b = true)) push!(res[dim][j], bt) end end end - + p = plot(; - ylabel = "Time/s", - xlabel = "N", - yscale=:log10, - xscale=:log10, - title = "Time for NxN sparse LU Factorization", - label = string(Symbol(parameterless_type(algs[1]))), - legend=:outertopright) + ylabel = "Time/s", + xlabel = "N", + yscale = :log10, + xscale = :log10, + title = "Time for NxN sparse LU Factorization", + label = string(Symbol(parameterless_type(algs[1]))), + legend = :outertopright) for dim in dims for i in 1:length(algs) plot!(p, ns, res[dim][i]; - linecolor=cols[i], - linestyle= lst[dim], - label = "$(string(Symbol(parameterless_type(algs[i])))) $(dim)D" ) + linecolor = cols[i], + linestyle = lst[dim], + label = "$(string(Symbol(parameterless_type(algs[i])))) $(dim)D") end end savefig("sparselubench.png") savefig("sparselubench.pdf") end -