Benchmarks
Here are some simple benchmarks. Take them with a grain of salt since they run on virtual machines in the cloud to generate the documentation automatically.
First-derivative operators
Periodic domains
Let's set up some benchmark code.
using BenchmarkTools
using LinearAlgebra, SparseArrays
using SummationByPartsOperators
BLAS.set_num_threads(1) # make sure that BLAS is serial to be fair
T = Float64
xmin, xmax = T(0), T(1)
D_SBP = periodic_derivative_operator(derivative_order=1, accuracy_order=2,
xmin=xmin, xmax=xmax, N=100)
x = grid(D_SBP)
D_sparse = sparse(D_SBP)
u = randn(eltype(D_SBP), length(x)); du = similar(u);
@show D_SBP * u ≈ D_sparse * u
function doit(D, text, du, u)
println(text)
sleep(0.1)
show(stdout, MIME"text/plain"(), @benchmark mul!($du, $D, $u))
println()
enddoit (generic function with 1 method)First, we benchmark the implementation from SummationByPartsOperators.jl.
doit(D_SBP, "D_SBP:", du, u)D_SBP:
BenchmarkTools.Trial: 10000 samples with 995 evaluations per sample.
Range (min … max): 27.478 ns … 52.723 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 28.354 ns ┊ GC (median): 0.00%
Time (mean ± σ): 28.709 ns ± 1.673 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▃▄▅██▂▄▂
▅████████▃▄▅▄▃▂▂▃▃▂▂▂▁▂▂▁▂▂▁▂▂▂▁▁▁▁▂▁▁▁▁▂▂▂▁▁▁▂▂▂▂▂▂▂▂▂▂▂▂▂ ▃
27.5 ns Histogram: frequency by time 37.3 ns <
Memory estimate: 0 bytes, allocs estimate: 0.Next, we compare this to the runtime obtained using a sparse matrix representation of the derivative operator. Depending on the hardware etc., this can be an order of magnitude slower than the optimized implementation from SummationByPartsOperators.jl.
doit(D_sparse, "D_sparse:", du, u)D_sparse:
BenchmarkTools.Trial: 10000 samples with 482 evaluations per sample.
Range (min … max): 223.033 ns … 386.494 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 228.060 ns ┊ GC (median): 0.00%
Time (mean ± σ): 231.004 ns ± 9.021 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▆█▅▂▂▂
▁▃████████▇▅▆▆▅▄▄▃▃▂▂▂▂▂▁▂▂▂▁▁▂▂▃▂▂▂▂▂▂▂▂▂▁▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
223 ns Histogram: frequency by time 258 ns <
Memory estimate: 0 bytes, allocs estimate: 0.These results were obtained using the following versions.
using InteractiveUtils
versioninfo()
using Pkg
Pkg.status(["SummationByPartsOperators"],
mode=PKGMODE_MANIFEST)Julia Version 1.10.11
Commit a2b11907d7b (2026-03-09 14:59 UTC)
Build Info:
Official https://julialang.org/ release
Platform Info:
OS: Linux (x86_64-linux-gnu)
CPU: 4 × AMD EPYC 9V74 80-Core Processor
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-15.0.7 (ORCJIT, znver3)
Threads: 2 default, 0 interactive, 1 GC (on 4 virtual cores)
Environment:
JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
Status `~/work/SummationByPartsOperators.jl/SummationByPartsOperators.jl/docs/Manifest.toml`
[9f78cca6] SummationByPartsOperators v0.5.96 `~/work/SummationByPartsOperators.jl/SummationByPartsOperators.jl`Bounded domains
We start again by setting up some benchmark code.
using BenchmarkTools
using LinearAlgebra, SparseArrays
using SummationByPartsOperators, BandedMatrices
BLAS.set_num_threads(1) # make sure that BLAS is serial to be fair
T = Float64
xmin, xmax = T(0), T(1)
D_SBP = derivative_operator(MattssonNordström2004(), derivative_order=1,
accuracy_order=6, xmin=xmin, xmax=xmax, N=10^3)
D_sparse = sparse(D_SBP)
D_banded = BandedMatrix(D_SBP)
u = randn(eltype(D_SBP), size(D_SBP, 1)); du = similar(u);
@show D_SBP * u ≈ D_sparse * u ≈ D_banded * u
function doit(D, text, du, u)
println(text)
sleep(0.1)
show(stdout, MIME"text/plain"(), @benchmark mul!($du, $D, $u))
println()
enddoit (generic function with 1 method)First, we benchmark the implementation from SummationByPartsOperators.jl.
doit(D_SBP, "D_SBP:", du, u)D_SBP:
BenchmarkTools.Trial: 10000 samples with 199 evaluations per sample.
Range (min … max): 424.809 ns … 1.116 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 429.392 ns ┊ GC (median): 0.00%
Time (mean ± σ): 433.730 ns ± 15.480 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▂▅▇███▇▆▅▃▄▃▂ ▁ ▁▁ ▂▂▃▃▂▂▂▁▁ ▃
█████████████▆█████▇▆▄▄▅▄▃▃▁▁▁▁▁▁▄▃▃▃▁▁▁▁▄▆████████████▇▇▆▇▇ █
425 ns Histogram: log(frequency) by time 481 ns <
Memory estimate: 0 bytes, allocs estimate: 0.Again, we compare this to a representation of the derivative operator as a sparse matrix. No surprise - it is again much slower, as in periodic domains.
doit(D_sparse, "D_sparse:", du, u)D_sparse:
BenchmarkTools.Trial: 10000 samples with 6 evaluations per sample.
Range (min … max): 5.189 μs … 10.968 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 5.395 μs ┊ GC (median): 0.00%
Time (mean ± σ): 5.445 μs ± 282.155 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▃▇█▄
▂▂▂▃▆█████▄▃▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▁▁▂▂▂▂▁▁▁▁▁▁▁▂▁▁▁▁▂▂▂▂▂▂▂▂▂▂▂▂▂ ▃
5.19 μs Histogram: frequency by time 6.85 μs <
Memory estimate: 0 bytes, allocs estimate: 0.Finally, we compare it to a representation as a banded matrix. Disappointingly, this is still much slower than the optimized implementation from SummationByPartsOperators.jl.
doit(D_banded, "D_banded:", du, u)D_banded:
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
Range (min … max): 9.995 μs … 37.426 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 10.025 μs ┊ GC (median): 0.00%
Time (mean ± σ): 10.133 μs ± 983.450 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
█ ▁
██▅▄▁▃▃▄▄▁▁▁▁▁▃▁▁▄▄▃▄▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▁▁▁▁▁▃▁▁▃▁▁▁▁▁▁▁▃▄▅ █
10 μs Histogram: log(frequency) by time 18 μs <
Memory estimate: 0 bytes, allocs estimate: 0.These results were obtained using the following versions.
using InteractiveUtils
versioninfo()
using Pkg
Pkg.status(["SummationByPartsOperators", "BandedMatrices"],
mode=PKGMODE_MANIFEST)Julia Version 1.10.11
Commit a2b11907d7b (2026-03-09 14:59 UTC)
Build Info:
Official https://julialang.org/ release
Platform Info:
OS: Linux (x86_64-linux-gnu)
CPU: 4 × AMD EPYC 9V74 80-Core Processor
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-15.0.7 (ORCJIT, znver3)
Threads: 2 default, 0 interactive, 1 GC (on 4 virtual cores)
Environment:
JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
Status `~/work/SummationByPartsOperators.jl/SummationByPartsOperators.jl/docs/Manifest.toml`
[aae01518] BandedMatrices v1.11.0
[9f78cca6] SummationByPartsOperators v0.5.96 `~/work/SummationByPartsOperators.jl/SummationByPartsOperators.jl`Dissipation operators
We follow the same structure as before. At first, we set up some benchmark code.
using BenchmarkTools
using LinearAlgebra, SparseArrays
using SummationByPartsOperators, BandedMatrices
BLAS.set_num_threads(1) # make sure that BLAS is serial to be fair
T = Float64
xmin, xmax = T(0), T(1)
D_SBP = derivative_operator(MattssonNordström2004(), derivative_order=1,
accuracy_order=6, xmin=xmin, xmax=xmax, N=10^3)
Di_SBP = dissipation_operator(MattssonSvärdNordström2004(), D_SBP)
Di_sparse = sparse(Di_SBP)
Di_banded = BandedMatrix(Di_SBP)
Di_full = Matrix(Di_SBP)
u = randn(eltype(D_SBP), size(D_SBP, 1)); du = similar(u);
@show Di_SBP * u ≈ Di_sparse * u ≈ Di_banded * u ≈ Di_full * u
function doit(D, text, du, u)
println(text)
sleep(0.1)
show(stdout, MIME"text/plain"(), @benchmark mul!($du, $D, $u))
println()
enddoit (generic function with 1 method)At first, let us benchmark the derivative and dissipation operators implemented in SummationByPartsOperators.jl.
doit(D_SBP, "D_SBP:", du, u)
doit(Di_SBP, "Di_SBP:", du, u)D_SBP:
BenchmarkTools.Trial: 10000 samples with 199 evaluations per sample.
Range (min … max): 421.588 ns … 678.658 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 426.523 ns ┊ GC (median): 0.00%
Time (mean ± σ): 430.492 ns ± 13.823 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▃▇█▅█▆▂
▄███████▆▅▄▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▁▂▁▂▁▁▁▂▁▁▁▂▂▂▂▃▃▃▃▃▃▃▂▂▂▂▂▂▂▂▂▂ ▃
422 ns Histogram: frequency by time 479 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
Di_SBP:
BenchmarkTools.Trial: 10000 samples with 17 evaluations per sample.
Range (min … max): 975.588 ns … 2.017 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 985.000 ns ┊ GC (median): 0.00%
Time (mean ± σ): 994.447 ns ± 67.611 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▇█▂ ▁
███▆▅▅▅▃▄▅▃▃▃▁▁▃▃▅▄▁▃▄▄▃▃▄▃▁▄▁▁▁▃▃▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▁▄▇▇▇ █
976 ns Histogram: log(frequency) by time 1.45 μs <
Memory estimate: 0 bytes, allocs estimate: 0.Next, we compare the results to sparse matrix representations. It will not come as a surprise that these are again much (around an order of magnitude) slower.
doit(Di_sparse, "Di_sparse:", du, u)
doit(Di_banded, "Di_banded:", du, u)Di_sparse:
BenchmarkTools.Trial: 10000 samples with 5 evaluations per sample.
Range (min … max): 6.313 μs … 13.438 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 6.456 μs ┊ GC (median): 0.00%
Time (mean ± σ): 6.515 μs ± 331.596 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▃▆██▇▅▃ ▂
▇███████▇▆▅▄▃▃▁▃▄▄▅▆▄▅▅▄▃▁▃▁▁▃▃▁▃▃▃▃▃▁▃▁▄▄▁▁▃▁▃▁▁▅▇█████▇▇▇ █
6.31 μs Histogram: log(frequency) by time 8.21 μs <
Memory estimate: 0 bytes, allocs estimate: 0.
Di_banded:
BenchmarkTools.Trial: 10000 samples with 4 evaluations per sample.
Range (min … max): 7.296 μs … 14.812 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 7.373 μs ┊ GC (median): 0.00%
Time (mean ± σ): 7.453 μs ± 415.913 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▄█▇▄▄▃ ▂
███████▅▅▁▃▃▄▄▃▄▅▄▅▄▁▅▁▄▄▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▁▁▁▁▅▇███▇▇▇ █
7.3 μs Histogram: log(frequency) by time 9.6 μs <
Memory estimate: 0 bytes, allocs estimate: 0.Finally, let's benchmark the same computation if a full (dense) matrix is used to represent the derivative operator. This is obviously a bad idea but 🤷
doit(Di_full, "Di_full:", du, u)Di_full:
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
Range (min … max): 128.253 μs … 428.004 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 132.719 μs ┊ GC (median): 0.00%
Time (mean ± σ): 136.124 μs ± 10.805 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▁▇█▇▅▃▃▃▁ ▃▄▄▃▂▁ ▂
▆▇▆██████████▇▆▇▅▇████████████▇▇▇▇▆▅▆▆▆▆▆▆▆▅▆▄▅▄▅▅▄▅▄▅▅▅▅▂▄▄▄ █
128 μs Histogram: log(frequency) by time 173 μs <
Memory estimate: 0 bytes, allocs estimate: 0.These results were obtained using the following versions.
using InteractiveUtils
versioninfo()
using Pkg
Pkg.status(["SummationByPartsOperators", "BandedMatrices"],
mode=PKGMODE_MANIFEST)Julia Version 1.10.11
Commit a2b11907d7b (2026-03-09 14:59 UTC)
Build Info:
Official https://julialang.org/ release
Platform Info:
OS: Linux (x86_64-linux-gnu)
CPU: 4 × AMD EPYC 9V74 80-Core Processor
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-15.0.7 (ORCJIT, znver3)
Threads: 2 default, 0 interactive, 1 GC (on 4 virtual cores)
Environment:
JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
Status `~/work/SummationByPartsOperators.jl/SummationByPartsOperators.jl/docs/Manifest.toml`
[aae01518] BandedMatrices v1.11.0
[9f78cca6] SummationByPartsOperators v0.5.96 `~/work/SummationByPartsOperators.jl/SummationByPartsOperators.jl`Structure-of-Arrays (SoA) and Array-of-Structures (AoS)
SummationByPartsOperators.jl tries to provide efficient support of
StaticVectors from StaticArrays.jl- StructArrays.jl
To demonstrate this, let us set up some benchmark code.
using BenchmarkTools
using StaticArrays, StructArrays
using LinearAlgebra, SparseArrays
using SummationByPartsOperators
BLAS.set_num_threads(1) # make sure that BLAS is serial to be fair
struct Vec5{T} <: FieldVector{5,T}
x1::T
x2::T
x3::T
x4::T
x5::T
end
# Apply `mul!` to each component of a plain array of structures one after another
function mul_aos!(du, D, u, args...)
for i in 1:size(du, 1)
mul!(view(du, i, :), D, view(u, i, :), args...)
end
end
T = Float64
xmin, xmax = T(0), T(1)
D_SBP = derivative_operator(MattssonNordström2004(), derivative_order=1,
accuracy_order=4, xmin=xmin, xmax=xmax, N=101)
D_sparse = sparse(D_SBP)
D_full = Matrix(D_SBP)101×101 Matrix{Float64}:
-141.176 173.529 -23.5294 … 0.0 0.0 0.0
-50.0 0.0 50.0 0.0 0.0 0.0
9.30233 -68.6047 0.0 0.0 0.0 0.0
3.06122 0.0 -60.2041 0.0 0.0 0.0
0.0 0.0 8.33333 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 -8.33333 0.0 0.0
0.0 0.0 0.0 60.2041 0.0 -3.06122
0.0 0.0 0.0 0.0 68.6047 -9.30233
0.0 0.0 0.0 -50.0 0.0 50.0
0.0 0.0 0.0 … 23.5294 -173.529 141.176At first, we benchmark the application of the operators implemented in SummationByPartsOperators.jl and their representations as sparse and dense matrices in the scalar case. As before, the sparse matrix representation is around an order of magnitude slower and the dense matrix representation is far off.
println("Scalar case")
u = randn(T, size(D_SBP, 1)); du = similar(u)
println("D_SBP")
show(stdout, MIME"text/plain"(), @benchmark mul!($du, $D_SBP, $u))
println("\nD_sparse")
show(stdout, MIME"text/plain"(), @benchmark mul!($du, $D_sparse, $u))
println("\nD_full")
show(stdout, MIME"text/plain"(), @benchmark mul!($du, $D_full, $u))Scalar case
D_SBP
BenchmarkTools.Trial: 10000 samples with 987 evaluations per sample.
Range (min … max): 49.832 ns … 92.733 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 51.922 ns ┊ GC (median): 0.00%
Time (mean ± σ): 52.588 ns ± 2.709 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▂▄▄▆█▇▆▅▅▅▄▄▁
▂█████████████▇▇▅▅▅▅▆▅▄▄▃▃▃▂▃▃▂▂▂▁▂▂▂▂▂▃▃▂▃▃▃▃▃▃▃▂▂▂▂▂▂▂▂▂▂ ▄
49.8 ns Histogram: frequency by time 62.7 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 204 evaluations per sample.
Range (min … max): 389.358 ns … 1.043 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 428.044 ns ┊ GC (median): 0.00%
Time (mean ± σ): 432.504 ns ± 22.702 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▃▄▆███▇▅▄▁
▂▂▁▂▁▂▂▂▂▂▂▂▃▃▄▅▇████████████▆▅▄▃▃▃▂▃▂▂▂▂▃▂▃▃▃▃▃▃▃▃▃▃▃▂▂▂▂▂▂ ▄
389 ns Histogram: frequency by time 488 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_full
BenchmarkTools.Trial: 10000 samples with 10 evaluations per sample.
Range (min … max): 1.222 μs … 10.087 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 1.236 μs ┊ GC (median): 0.00%
Time (mean ± σ): 1.253 μs ± 162.241 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
██▂ ▂
███▇▆▆▆▆▅▅▃▄▄▄▁▄▃▄▄▅▃▄▃▄▄▁▁▃▁▃▆▆▃▁▃▃▁▁▁▁▁▁▁▁▁▃▃▁▁▁▁▁▄▁▁▁▆▆▆ █
1.22 μs Histogram: log(frequency) by time 2.05 μs <
Memory estimate: 0 bytes, allocs estimate: 0.Next, we use a plain array of structures (AoS) in the form of a two-dimensional array and our custom mul_aos! implementation that loops over each component, using mul! on views. Here, the differences between the timings are less pronounced.
println("Plain Array of Structures")
u_aos_plain = randn(T, 5, size(D_SBP, 1)); du_aos_plain = similar(u_aos_plain)
println("D_SBP")
show(stdout, MIME"text/plain"(), @benchmark mul_aos!($du_aos_plain, $D_SBP, $u_aos_plain))
println("\nD_sparse")
show(stdout, MIME"text/plain"(), @benchmark mul_aos!($du_aos_plain, $D_sparse, $u_aos_plain))
println("\nD_full")
show(stdout, MIME"text/plain"(), @benchmark mul_aos!($du_aos_plain, $D_full, $u_aos_plain))Plain Array of Structures
D_SBP
BenchmarkTools.Trial: 10000 samples with 10 evaluations per sample.
Range (min … max): 1.305 μs … 4.137 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 1.315 μs ┊ GC (median): 0.00%
Time (mean ± σ): 1.329 μs ± 108.091 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
█
█▄▂▂▂▂▁▁▂▂▂▂▁▁▂▁▂▂▂▁▂▂▂▁▂▂▂▂▂▂▂▂▂▂▂▁▂▂▁▂▂▁▁▁▁▂▁▂▁▂▁▁▁▁▁▁▂▂▂ ▂
1.3 μs Histogram: frequency by time 2.1 μs <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 9 evaluations per sample.
Range (min … max): 2.822 μs … 5.831 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 2.903 μs ┊ GC (median): 0.00%
Time (mean ± σ): 2.934 μs ± 167.216 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▁▄▆███▆▅▄▃▁ ▂
███████████▆▄▅▄▄▄▁▄▅▄▄▅▅▅▄▃▃▅▃▁▁▃▃▁▃▁▃▁▃▃▃▃▁▁▁▃▁▁▁▄▄▆██████ █
2.82 μs Histogram: log(frequency) by time 3.85 μs <
Memory estimate: 240 bytes, allocs estimate: 5.
D_full
BenchmarkTools.Trial: 10000 samples with 4 evaluations per sample.
Range (min … max): 6.905 μs … 14.011 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 6.953 μs ┊ GC (median): 0.00%
Time (mean ± σ): 7.026 μs ± 403.387 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▇█▃ ▁
████▇▇▇▇▇▆▅▅▅▁▁▃▁▆▆▃▄▁▃▃▃▁▄▃▁▁▁▄▁▃▁▁▁▃▁▁▁▁▁▁▄▁▃▃▁▁▁▁▅▇███▇▇ █
6.91 μs Histogram: log(frequency) by time 9.16 μs <
Memory estimate: 0 bytes, allocs estimate: 0.Now, we use an array of structures (AoS) based on reinterpret and standard mul!. This is much more efficient for the implementation in SummationByPartsOperators.jl. In Julia v1.6, this is also more efficient for sparse matrices but less efficient for dense matrices (compared to the plain AoS approach with mul_aos! above).
println("Array of Structures (reinterpreted array)")
u_aos_r = reinterpret(reshape, Vec5{T}, u_aos_plain); du_aos_r = similar(u_aos_r)
@show D_SBP * u_aos_r ≈ D_sparse * u_aos_r ≈ D_full * u_aos_r
mul!(du_aos_r, D_SBP, u_aos_r)
@show reinterpret(reshape, T, du_aos_r) ≈ du_aos_plain
println("D_SBP")
show(stdout, MIME"text/plain"(), @benchmark mul!($du_aos_r, $D_SBP, $u_aos_r))
println("\nD_sparse")
show(stdout, MIME"text/plain"(), @benchmark mul!($du_aos_r, $D_sparse, $u_aos_r))
println("\nD_full")
show(stdout, MIME"text/plain"(), @benchmark mul!($du_aos_r, $D_full, $u_aos_r))Array of Structures (reinterpreted array)
D_SBP * u_aos_r ≈ D_sparse * u_aos_r ≈ D_full * u_aos_r = true
reinterpret(reshape, T, du_aos_r) ≈ du_aos_plain = true
D_SBP
BenchmarkTools.Trial: 10000 samples with 409 evaluations per sample.
Range (min … max): 232.868 ns … 362.037 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 234.411 ns ┊ GC (median): 0.00%
Time (mean ± σ): 236.634 ns ± 6.800 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▁▇█▆▄▂ ▂▃▃▂▁▁ ▂
▃███████▇█▄▄▄▅▆▇▇▇▆▆▄▄▃▃▃▁▁▁▁▁▁▃▁▁▄▁▁▁▃▁▁▃▃▅██████████▇▇▇▇▇▆▆ █
233 ns Histogram: log(frequency) by time 259 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 179 evaluations per sample.
Range (min … max): 588.318 ns … 16.667 μs ┊ GC (min … max): 0.00% … 93.82%
Time (median): 652.209 ns ┊ GC (median): 0.00%
Time (mean ± σ): 662.461 ns ± 169.358 ns ┊ GC (mean ± σ): 0.24% ± 0.94%
▁▁▄▆█▆▃▁
▁▂▃▂▂▂▃▅▇██████████▇▆▅▄▄▅▄▃▃▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▃
588 ns Histogram: frequency by time 849 ns <
Memory estimate: 32 bytes, allocs estimate: 1.
D_full
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
Range (min … max): 13.300 μs … 31.858 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 13.370 μs ┊ GC (median): 0.00%
Time (mean ± σ): 13.511 μs ± 1.059 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
█▁ ▁
██▆▄▄▃▄▃▄▁▃▁▃▁▁▁▄▄▄▄▄▁▁▃▃▁▃▁▃▁▁▃▁▁▄▃▄▁▁▄▁▁▃▃▃▁▃▃▁▁▃▃▁▃▁▃▄▇▇ █
13.3 μs Histogram: log(frequency) by time 21.2 μs <
Memory estimate: 0 bytes, allocs estimate: 0.Next, we still use an array of structures (AoS), but copy the data into a plain Array instead of using the reinterpreted versions. There is no significant difference to the previous version in this case.
println("Array of Structures")
u_aos = Array(u_aos_r); du_aos = similar(u_aos)
@show D_SBP * u_aos ≈ D_sparse * u_aos ≈ D_full * u_aos
mul!(du_aos, D_SBP, u_aos)
@show du_aos ≈ du_aos_r
println("D_SBP")
show(stdout, MIME"text/plain"(), @benchmark mul!($du_aos, $D_SBP, $u_aos))
println("\nD_sparse")
show(stdout, MIME"text/plain"(), @benchmark mul!($du_aos, $D_sparse, $u_aos))
println("\nD_full")
show(stdout, MIME"text/plain"(), @benchmark mul!($du_aos, $D_full, $u_aos))Array of Structures
D_SBP * u_aos ≈ D_sparse * u_aos ≈ D_full * u_aos = true
du_aos ≈ du_aos_r = true
D_SBP
BenchmarkTools.Trial: 10000 samples with 413 evaluations per sample.
Range (min … max): 232.869 ns … 334.329 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 238.421 ns ┊ GC (median): 0.00%
Time (mean ± σ): 241.363 ns ± 7.675 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▄▇█▄▃▁
▁▁▁▂▇███████▇▆▅▃▄▃▃▃▃▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▂▂▂▂▂▂▂▂▂▁▂▂▁▁▁▁▁▁▁▁▁▁ ▂
233 ns Histogram: frequency by time 267 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 171 evaluations per sample.
Range (min … max): 617.655 ns … 1.314 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 690.222 ns ┊ GC (median): 0.00%
Time (mean ± σ): 693.249 ns ± 30.397 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▁▂▂▄▄▇▇▇█▅▃▃▃▁ ▁▁
▁▁▂▂▂▃▂▂▂▂▃▃▃▄▄▆▆▇█▇█████████████████▇▅▅▆▅▄▄▄▄▃▄▄▃▃▃▃▃▂▂▂▂▂▁ ▄
618 ns Histogram: frequency by time 772 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_full
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
Range (min … max): 13.320 μs … 38.347 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 13.390 μs ┊ GC (median): 0.00%
Time (mean ± σ): 13.524 μs ± 1.072 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
█ ▁
██▆▅▄▄▄▃▁▁▃▄▁▁▁▄▄▃▄▄▄▃▁▁▁▁▁▁▁▃▃▃▁▁▃▁▃▁▃▃▃▁▁▁▃▃▄▁▁▃▁▄▃▄▁▁▃▆▇ █
13.3 μs Histogram: log(frequency) by time 21.2 μs <
Memory estimate: 0 bytes, allocs estimate: 0.Finally, let's look at a structure of arrays (SoA). Interestingly, this is slower than the array of structures we used above. On Julia v1.6, the sparse matrix representation performs particularly bad in this case.
println("Structure of Arrays")
u_soa = StructArray(u_aos); du_soa = similar(u_soa)
@show D_SBP * u_soa ≈ D_sparse * u_soa ≈ D_full * u_soa
mul!(du_soa, D_SBP, u_soa)
@show du_soa ≈ du_aos
println("D_SBP")
show(stdout, MIME"text/plain"(), @benchmark mul!($du_soa, $D_SBP, $u_soa))
println("\nD_sparse")
show(stdout, MIME"text/plain"(), @benchmark mul!($du_soa, $D_sparse, $u_soa))
println("\nD_full")
show(stdout, MIME"text/plain"(), @benchmark mul!($du_soa, $D_full, $u_soa))Structure of Arrays
D_SBP * u_soa ≈ D_sparse * u_soa ≈ D_full * u_soa = true
du_soa ≈ du_aos = true
D_SBP
BenchmarkTools.Trial: 10000 samples with 376 evaluations per sample.
Range (min … max): 248.085 ns … 391.439 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 253.149 ns ┊ GC (median): 0.00%
Time (mean ± σ): 255.352 ns ± 7.835 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▁▅▇▅▆█▅▂▁
▃▄▅▆█████████▇▅▄▃▂▃▂▂▂▂▂▂▂▂▂▂▂▂▁▁▁▁▂▁▁▂▂▂▂▃▂▃▃▃▃▃▃▃▃▂▂▂▂▂▂▂▂▂ ▃
248 ns Histogram: frequency by time 282 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
Range (min … max): 51.087 μs … 110.767 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 51.377 μs ┊ GC (median): 0.00%
Time (mean ± σ): 51.911 μs ± 2.241 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▅█▇▅▂ ▁▁▁ ▂
██████▆▆▆▅▆▄▅▄▅▇▆▆▅▅▁▄▄▃▁▃▃▁▁▁▁▃▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▆▇██████▇▇▇▆▆ █
51.1 μs Histogram: log(frequency) by time 60.9 μs <
Memory estimate: 0 bytes, allocs estimate: 0.
D_full
BenchmarkTools.Trial: 10000 samples with 6 evaluations per sample.
Range (min … max): 5.627 μs … 10.628 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 5.654 μs ┊ GC (median): 0.00%
Time (mean ± σ): 5.711 μs ± 285.667 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▇█▄ ▁ ▂
█████▆▅▄▅▅▄▃▁▄▃▅▆▅▅▄▃▃▄▄▁▁▃▃▁▁▁▁▁▁▄▃▃▃▁▃▁▄▁▃▁▄▁▁▁▁▃▃█████▇▇ █
5.63 μs Histogram: log(frequency) by time 7.08 μs <
Memory estimate: 0 bytes, allocs estimate: 0.These results were obtained using the following versions.
using InteractiveUtils
versioninfo()
using Pkg
Pkg.status(["SummationByPartsOperators", "StaticArrays", "StructArrays"],
mode=PKGMODE_MANIFEST)Julia Version 1.10.11
Commit a2b11907d7b (2026-03-09 14:59 UTC)
Build Info:
Official https://julialang.org/ release
Platform Info:
OS: Linux (x86_64-linux-gnu)
CPU: 4 × AMD EPYC 9V74 80-Core Processor
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-15.0.7 (ORCJIT, znver3)
Threads: 2 default, 0 interactive, 1 GC (on 4 virtual cores)
Environment:
JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
Status `~/work/SummationByPartsOperators.jl/SummationByPartsOperators.jl/docs/Manifest.toml`
[90137ffa] StaticArrays v1.9.18
[09ab397b] StructArrays v0.7.3
[9f78cca6] SummationByPartsOperators v0.5.96 `~/work/SummationByPartsOperators.jl/SummationByPartsOperators.jl`