@@ -8,8 +8,6 @@ BenchmarkTools.DEFAULT_PARAMETERS.seconds = 0.5
88# Why do I need to set this ?
99BenchmarkTools. DEFAULT_PARAMETERS. samples = 10
1010
11-
12-
1311# Sparse matrix generation on a n-dimensional rectangular grid. After
1412# https://discourse.julialang.org/t/seven-lines-of-julia-examples-sought/50416/135
1513# by A. Braunstein.
@@ -29,15 +27,19 @@ lattice(L...; Tv = Float64) = lattice(L[1]; Tv) ⊕ lattice(L[2:end]...; Tv)
2927# Create a matrix similar to that of a finite difference discretization in a `dim`-dimensional
3028# unit cube of ``-Δu + δu`` with approximately N unknowns. It is strictly diagonally dominant.
3129#
32- function fdmatrix (N; dim= 2 , Tv = Float64, δ = 1.0e-2 )
30+ function fdmatrix (N; dim = 2 , Tv = Float64, δ = 1.0e-2 )
3331 n = N^ (1 / dim) |> ceil |> Int
3432 lattice ([n for i in 1 : dim]. .. ; Tv) + Tv (δ) * I
3533end
3634
37-
38- algs = [UMFPACKFactorization (),KLUFactorization (),MKLPardisoFactorize (),SparspakFactorization ()]
39- cols= [:red , :blue , :green , :magenta ,:turqoise ] # one color per alg
40- lst= [:dash ,:solid , :dashdot ] # one line style per dim
35+ algs = [
36+ UMFPACKFactorization (),
37+ KLUFactorization (),
38+ MKLPardisoFactorize (),
39+ SparspakFactorization (),
40+ ]
41+ cols = [:red , :blue , :green , :magenta , :turqoise ] # one color per alg
42+ lst = [:dash , :solid , :dashdot ] # one line style per dim
4143
4244__parameterless_type (T) = Base. typename (T). wrapper
4345parameterless_type (x) = __parameterless_type (typeof (x))
@@ -48,46 +50,48 @@ parameterless_type(::Type{T}) where {T} = __parameterless_type(T)
4850# kmax=15 gives ≈ 328_000 unknows, you can go make a coffee.
4951# Main culprit is KLU factorization in 3D.
5052#
51- function run_and_plot (;dims= [1 ,2 ,3 ], kmax= 12 )
52-
53+ function run_and_plot (; dims = [1 , 2 , 3 ], kmax = 12 )
5354 ns = [10 * 2 ^ k for k in 0 : kmax]
54-
55- res = [ [Float64[] for i in 1 : length (algs)] for dim in dims ]
56-
55+
56+ res = [[Float64[] for i in 1 : length (algs)] for dim in dims]
57+
5758 for dim in dims
5859 for i in 1 : length (ns)
5960 rng = MersenneTwister (123 )
60- A = fdmatrix (ns[i];dim)
61- n = size (A,1 )
61+ A = fdmatrix (ns[i]; dim)
62+ n = size (A, 1 )
6263 @info " dim=$(dim) : $n × $n "
6364 b = rand (rng, n)
64- u0= rand (rng, n)
65-
65+ u0 = rand (rng, n)
66+
6667 for j in 1 : length (algs)
67- bt = @belapsed solve (prob, $ (algs[j])). u setup= (prob = LinearProblem (copy ($ A), copy ($ b); u0 = copy ($ u0), alias_A= true , alias_b= true ))
68+ bt = @belapsed solve (prob, $ (algs[j])). u setup= (prob = LinearProblem (copy ($ A),
69+ copy ($ b);
70+ u0 = copy ($ u0),
71+ alias_A = true ,
72+ alias_b = true ))
6873 push! (res[dim][j], bt)
6974 end
7075 end
7176 end
72-
77+
7378 p = plot (;
74- ylabel = " Time/s" ,
75- xlabel = " N" ,
76- yscale= :log10 ,
77- xscale= :log10 ,
78- title = " Time for NxN sparse LU Factorization" ,
79- label = string (Symbol (parameterless_type (algs[1 ]))),
80- legend= :outertopright )
79+ ylabel = " Time/s" ,
80+ xlabel = " N" ,
81+ yscale = :log10 ,
82+ xscale = :log10 ,
83+ title = " Time for NxN sparse LU Factorization" ,
84+ label = string (Symbol (parameterless_type (algs[1 ]))),
85+ legend = :outertopright )
8186
8287 for dim in dims
8388 for i in 1 : length (algs)
8489 plot! (p, ns, res[dim][i];
85- linecolor= cols[i],
86- linestyle= lst[dim],
87- label = " $(string (Symbol (parameterless_type (algs[i])))) $(dim) D" )
90+ linecolor = cols[i],
91+ linestyle = lst[dim],
92+ label = " $(string (Symbol (parameterless_type (algs[i])))) $(dim) D" )
8893 end
8994 end
9095 savefig (" sparselubench.png" )
9196 savefig (" sparselubench.pdf" )
9297end
93-
0 commit comments