@@ -8,8 +8,6 @@ BenchmarkTools.DEFAULT_PARAMETERS.seconds = 0.5
8
8
# Why do I need to set this ?
9
9
BenchmarkTools. DEFAULT_PARAMETERS. samples = 10
10
10
11
-
12
-
13
11
# Sparse matrix generation on a n-dimensional rectangular grid. After
14
12
# https://discourse.julialang.org/t/seven-lines-of-julia-examples-sought/50416/135
15
13
# by A. Braunstein.
@@ -29,15 +27,19 @@ lattice(L...; Tv = Float64) = lattice(L[1]; Tv) ⊕ lattice(L[2:end]...; Tv)
29
27
# Create a matrix similar to that of a finite difference discretization in a `dim`-dimensional
30
28
# unit cube of ``-Δu + δu`` with approximately N unknowns. It is strictly diagonally dominant.
31
29
#
32
- function fdmatrix (N; dim= 2 , Tv = Float64, δ = 1.0e-2 )
30
+ function fdmatrix (N; dim = 2 , Tv = Float64, δ = 1.0e-2 )
33
31
n = N^ (1 / dim) |> ceil |> Int
34
32
lattice ([n for i in 1 : dim]. .. ; Tv) + Tv (δ) * I
35
33
end
36
34
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
41
43
42
44
__parameterless_type (T) = Base. typename (T). wrapper
43
45
parameterless_type (x) = __parameterless_type (typeof (x))
@@ -48,46 +50,48 @@ parameterless_type(::Type{T}) where {T} = __parameterless_type(T)
48
50
# kmax=15 gives ≈ 328_000 unknows, you can go make a coffee.
49
51
# Main culprit is KLU factorization in 3D.
50
52
#
51
- function run_and_plot (;dims= [1 ,2 ,3 ], kmax= 12 )
52
-
53
+ function run_and_plot (; dims = [1 , 2 , 3 ], kmax = 12 )
53
54
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
+
57
58
for dim in dims
58
59
for i in 1 : length (ns)
59
60
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 )
62
63
@info " dim=$(dim) : $n × $n "
63
64
b = rand (rng, n)
64
- u0= rand (rng, n)
65
-
65
+ u0 = rand (rng, n)
66
+
66
67
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 ))
68
73
push! (res[dim][j], bt)
69
74
end
70
75
end
71
76
end
72
-
77
+
73
78
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 )
81
86
82
87
for dim in dims
83
88
for i in 1 : length (algs)
84
89
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" )
88
93
end
89
94
end
90
95
savefig (" sparselubench.png" )
91
96
savefig (" sparselubench.pdf" )
92
97
end
93
-
0 commit comments