Timing for banded matrices¶
We’ll use a large banded matrix to observe the speedup possible in LU factorization. We’ll also need to load in a (standard) package for sparse matrices.
using FundamentalsNumericalComputation
If we use an ordinary “dense” matrix, then there’s no way to exploit a banded structure such as tridiagonality.
n = 10000
A = diagm(0=>1:n,1=>n-1:-1:1,-1=>ones(n-1))
10000×10000 Array{Float64,2}:
1.0 9999.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0
1.0 2.0 9998.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 1.0 3.0 9997.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 1.0 4.0 9996.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 1.0 5.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 1.0 … 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
⋮ ⋱
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 … 4.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 9997.0 3.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 1.0 9998.0 2.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 1.0 9999.0 1.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 10000.0
@time lu(A);
6.691195 seconds (108 allocations: 763.022 MiB, 2.07% gc time)
If instead we construct a proper “sparse” matrix, though, the speedup can be dramatic.
A = spdiagm(0=>1:n,1=>n-1:-1:1,-1=>ones(n-1))
10000×10000 SparseMatrixCSC{Float64,Int64} with 29998 stored entries:
[1 , 1] = 1.0
[2 , 1] = 1.0
[1 , 2] = 9999.0
[2 , 2] = 2.0
[3 , 2] = 1.0
[2 , 3] = 9998.0
[3 , 3] = 3.0
[4 , 3] = 1.0
[3 , 4] = 9997.0
[4 , 4] = 4.0
[5 , 4] = 1.0
[4 , 5] = 9996.0
⋮
[9996 , 9996] = 9996.0
[9997 , 9996] = 1.0
[9996 , 9997] = 4.0
[9997 , 9997] = 9997.0
[9998 , 9997] = 1.0
[9997 , 9998] = 3.0
[9998 , 9998] = 9998.0
[9999 , 9998] = 1.0
[9998 , 9999] = 2.0
[9999 , 9999] = 9999.0
[10000, 9999] = 1.0
[9999 , 10000] = 1.0
[10000, 10000] = 10000.0
@time lu(A);
0.097763 seconds (274.35 k allocations: 23.452 MiB)