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 BenchmarkToolsusing LinearAlgebra, SparseArraysusing SummationByPartsOperatorsBLAS.set_num_threads(1) # make sure that BLAS is serial to be fairT = Float64xmin, 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 * ufunction 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 997 evaluations per sample.
Range (min … max): 18.008 ns … 72.853 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 18.176 ns ┊ GC (median): 0.00%
Time (mean ± σ): 18.447 ns ± 1.548 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▇█▃ ▁
███▆▅▄▆▄▃▅▄▅▆▆▆▅▅▄▃▆▅▅▄▄▅▇▄▅▃▄▄▄▅▅▆▁▁▁▃▄▁▁▃▁▄▁▃▁▁█▇▁▁▁▃▅▇▇▇ █
18 ns Histogram: log(frequency) by time 26.9 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 800 evaluations per sample.
Range (min … max): 157.064 ns … 300.830 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 162.683 ns ┊ GC (median): 0.00%
Time (mean ± σ): 164.046 ns ± 3.820 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▁ ▁▆█▅▂▁ ▁▁▁▁▂▄▃▁ ▂▂▁▁▁▁▃▂▁ ▁
▆▅▄▅▇▆▇▇▇▆▇▇█████████████████████████████████▇▇█▇▇▇▆▆▅▆▅▅▅▅▅▅ █
157 ns Histogram: log(frequency) by time 175 ns <
Memory estimate: 0 bytes, allocs estimate: 0.These results were obtained using the following versions.
using InteractiveUtilsversioninfo()using PkgPkg.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 × Intel(R) Xeon(R) 6973P-C
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-15.0.7 (ORCJIT, icelake-client)
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.97-DEV `~/work/SummationByPartsOperators.jl/SummationByPartsOperators.jl`Bounded domains
We start again by setting up some benchmark code.
using BenchmarkToolsusing LinearAlgebra, SparseArraysusing SummationByPartsOperators, BandedMatricesBLAS.set_num_threads(1) # make sure that BLAS is serial to be fairT = Float64xmin, 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 * ufunction 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 287 evaluations per sample.
Range (min … max): 282.373 ns … 408.237 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 283.955 ns ┊ GC (median): 0.00%
Time (mean ± σ): 287.713 ns ± 10.808 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▂██ ▁ ▂▂▁▁▁▁ ▁
███▇▄▄▅▄▆▇▇▆▄▄▄▅▄▄▄▃▃▅▅▅▄▄▄▃▆▇▇▇▆▅▆█▆▂▅████████▇▆▅▅▅▄▃▅▄▃▄▄▃▄ █
282 ns Histogram: log(frequency) by time 329 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 8 evaluations per sample.
Range (min … max): 3.141 μs … 6.322 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 3.193 μs ┊ GC (median): 0.00%
Time (mean ± σ): 3.227 μs ± 166.245 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▄▇█▇▅▄▃▂▂ ▁ ▂
▇██████████▇▄▅▁▁▆▆▆▄▄▃▃▁▁▁▄▁▃▁▃▃▃▄▁▁▄▄▃▃▄▃▁▃▅▇████▇█▇▇▆▆▆▆▆ █
3.14 μs Histogram: log(frequency) by time 3.94 μ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 4 evaluations per sample.
Range (min … max): 7.058 μs … 11.945 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 7.126 μs ┊ GC (median): 0.00%
Time (mean ± σ): 7.174 μs ± 279.070 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
█▂█▃▅▃▂▄ ▄▁ ▂
███████████▄▄▄▅▄▄▄▄▄▅▄▅▃▅▃▄▃▃▁▁▁▁▁▃▃▁▃▁▁▁▁▃▄▄▁▅▅▅▆▇▆▇█▇█▇▇▇ █
7.06 μs Histogram: log(frequency) by time 8.63 μs <
Memory estimate: 0 bytes, allocs estimate: 0.These results were obtained using the following versions.
using InteractiveUtilsversioninfo()using PkgPkg.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 × Intel(R) Xeon(R) 6973P-C
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-15.0.7 (ORCJIT, icelake-client)
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.97-DEV `~/work/SummationByPartsOperators.jl/SummationByPartsOperators.jl`Dissipation operators
We follow the same structure as before. At first, we set up some benchmark code.
using BenchmarkToolsusing LinearAlgebra, SparseArraysusing SummationByPartsOperators, BandedMatricesBLAS.set_num_threads(1) # make sure that BLAS is serial to be fairT = Float64xmin, 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 * ufunction 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 262 evaluations per sample.
Range (min … max): 296.756 ns … 657.302 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 298.214 ns ┊ GC (median): 0.00%
Time (mean ± σ): 301.963 ns ± 11.800 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▁██▁ ▂▂▁▁ ▁
████▅▄▄▄▅▅▆▇▇▆▅▄▃▄▅▅▆▄▄▄▄▄▅▆▇▇▆▅▅▄▄▅▆▅▄▄█▇▆▅▅███████▇▇▆▅▄▅▃▃▃ █
297 ns Histogram: log(frequency) by time 340 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
Di_SBP:
BenchmarkTools.Trial: 10000 samples with 129 evaluations per sample.
Range (min … max): 733.450 ns … 1.659 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 756.818 ns ┊ GC (median): 0.00%
Time (mean ± σ): 768.888 ns ± 46.227 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
█▆▅▆▆▅▅▅▅█▃▁ ▁ ▃▃▃▃▂▂▂▂▁▁ ▂
███████████████████████████▇▆▇▇▇▇▇▇▇▅▆▅▅▅▅▅▅▃▅▄▁▄▅▅▅▄▆▅▇▇▇▇▆ █
733 ns Histogram: log(frequency) by time 995 ns <
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 8 evaluations per sample.
Range (min … max): 3.636 μs … 5.569 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 3.712 μs ┊ GC (median): 0.00%
Time (mean ± σ): 3.745 μs ± 155.958 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▃▆▇██▇▆▅▄▃▂▁ ▃
▇████████████▇▆▃▅▁▅▅▆▅▅▅▆▃▅▅▅▃▁▄▁▃▁▄▃▄▅▆▆█▇▇▇▆▆▇▇███████▇▇▇ █
3.64 μs Histogram: log(frequency) by time 4.47 μs <
Memory estimate: 0 bytes, allocs estimate: 0.
Di_banded:
BenchmarkTools.Trial: 10000 samples with 6 evaluations per sample.
Range (min … max): 5.147 μs … 7.361 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 5.192 μs ┊ GC (median): 0.00%
Time (mean ± σ): 5.235 μs ± 218.510 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
█▃█▃▄▁▂▁ ▂ ▁ ▂ ▁ ▁ ▂
██████████▇█▇█▇▁▃▃▄▃▁▄▃▄▅▃▅▅▃▃▄▃▄▃▁▃▃▁▁▁▁▆█▆█▅▆▆▅▆▇▇▆▇▇██▇█ █
5.15 μs Histogram: log(frequency) by time 6.24 μ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): 260.518 μs … 367.391 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 267.664 μs ┊ GC (median): 0.00%
Time (mean ± σ): 271.064 μs ± 6.384 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▁▂▂▂▂▂▆██▇▄▁ ▁▁ ▃▅▅▅▄▄▃▃▂▂▁▁ ▂
▃▁▅▅▅▄▆████████████▆▅▄▃▅▅█▇████▇▆▇▆▆███████████████▇▆▆▇▇▇▇▇▇▆ █
261 μs Histogram: log(frequency) by time 288 μs <
Memory estimate: 0 bytes, allocs estimate: 0.These results were obtained using the following versions.
using InteractiveUtilsversioninfo()using PkgPkg.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 × Intel(R) Xeon(R) 6973P-C
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-15.0.7 (ORCJIT, icelake-client)
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.97-DEV `~/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 BenchmarkToolsusing StaticArrays, StructArraysusing LinearAlgebra, SparseArraysusing SummationByPartsOperatorsBLAS.set_num_threads(1) # make sure that BLAS is serial to be fairstruct Vec5{T} <: FieldVector{5,T} x1::T x2::T x3::T x4::T x5::Tend# Apply `mul!` to each component of a plain array of structures one after anotherfunction mul_aos!(du, D, u, args...) for i in 1:size(du, 1) mul!(view(du, i, :), D, view(u, i, :), args...) endendT = Float64xmin, 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 994 evaluations per sample.
Range (min … max): 31.694 ns … 93.209 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 33.089 ns ┊ GC (median): 0.00%
Time (mean ± σ): 33.413 ns ± 1.971 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▂ ▃▃ █▄ ▁
█▃▁██▁▄██▆▅▅▆▅▅▆▇▆▆▆▆▆▆▅▆▆▅▄▅▄▅▅▅▁▃▄▁▁▁▃▁▃▄▃▃▄▁▄▄▄▆▅▆███▇▇▇ █
31.7 ns Histogram: log(frequency) by time 42.4 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 401 evaluations per sample.
Range (min … max): 241.411 ns … 313.015 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 243.439 ns ┊ GC (median): 0.00%
Time (mean ± σ): 244.694 ns ± 4.666 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▃▆▅▄▃▅█▇▄▁ ▂▃▂▂▁▁▁▁ ▂
██████████▇▆▄▁▄▄▆█▇▇▅▅▃▄▁▁▃▁▃▃▃▃▄▁▄▁▁▁▆▇▇██████████████▇▇▆▅▅▆ █
241 ns Histogram: log(frequency) by time 260 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_full
BenchmarkTools.Trial: 10000 samples with 10 evaluations per sample.
Range (min … max): 1.104 μs … 5.497 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 1.114 μs ┊ GC (median): 0.00%
Time (mean ± σ): 1.125 μs ± 96.926 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
█
█▇▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▁▁▁▂▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▁▂▁▁▁▁▂▁▂▂▂ ▂
1.1 μs Histogram: frequency by time 1.82 μ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 194 evaluations per sample.
Range (min … max): 493.108 ns … 891.361 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 498.613 ns ┊ GC (median): 0.00%
Time (mean ± σ): 504.074 ns ± 16.368 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▂▅▇██▇▆▅▅▄▄▃▂ ▂▂▂▃▂▃▂▂▂▁ ▃
▇███████████████▇▆▆▅▅▅▃▁▃▁▁▃▁▁▁▁▃▁▁▃▅▅▄▄▅▅▆▆▆▇████████████▇▇▆ █
493 ns Histogram: log(frequency) by time 554 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 10 evaluations per sample.
Range (min … max): 1.603 μs … 6.062 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 1.626 μs ┊ GC (median): 0.00%
Time (mean ± σ): 1.648 μs ± 96.836 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
█▂
▃██▅▅▅▄▃▃▃▂▂▂▂▂▂▂▁▂▁▂▂▂▂▂▂▂▂▂▂▂▂▂▁▁▁▁▁▁▁▂▁▁▁▁▂▁▂▂▁▂▂▂▂▂▂▂▂ ▂
1.6 μs Histogram: frequency by time 2.16 μs <
Memory estimate: 240 bytes, allocs estimate: 5.
D_full
BenchmarkTools.Trial: 10000 samples with 5 evaluations per sample.
Range (min … max): 6.355 μs … 10.674 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 6.410 μs ┊ GC (median): 0.00%
Time (mean ± σ): 6.466 μs ± 297.511 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▃██▂ ▁ ▂
████▇▅▄▁▃▁▁▄▁▃▆▆▅▅▁▃▄▁▄▃▁▃▃▁▃▃▄▄▁▁▃▄▅▁▁▃▁▁▁▁▁▁▁▁▃▃▄▄▅█████▇ █
6.36 μs Histogram: log(frequency) by time 7.99 μ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_rmul!(du_aos_r, D_SBP, u_aos_r)@show reinterpret(reshape, T, du_aos_r) ≈ du_aos_plainprintln("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 907 evaluations per sample.
Range (min … max): 119.821 ns … 886.234 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 120.393 ns ┊ GC (median): 0.00%
Time (mean ± σ): 122.170 ns ± 10.743 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▃█▇▁ ▃▂▂▁▁ ▁
████▄▅▄▅▇█▇▅▅▅▄▄▃▅▇▇█▇▆▅▃▂▂▄▆▅▅▅▅▆███████▇▇▆▆▅▅▃▄▄▃▃▂▅▅▅▅▃▅▄▄ █
120 ns Histogram: log(frequency) by time 138 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 207 evaluations per sample.
Range (min … max): 364.812 ns … 17.290 μs ┊ GC (min … max): 0.00% … 96.26%
Time (median): 373.932 ns ┊ GC (median): 0.00%
Time (mean ± σ): 378.837 ns ± 169.496 ns ┊ GC (mean ± σ): 0.44% ± 0.96%
▁▄▅▆███▇▆▅▄▄▃▂ ▂▂▃▂▂▂▂▂▁▁ ▃
▃▅▅▇▇███████████████▇███▆▅▆▅▃▁▃▆▇▇▇██████████████▇▇▆▆▆▆▆▄▃▄▃▅ █
365 ns Histogram: log(frequency) by time 415 ns <
Memory estimate: 32 bytes, allocs estimate: 1.
D_full
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
Range (min … max): 9.312 μs … 25.455 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 9.354 μs ┊ GC (median): 0.00%
Time (mean ± σ): 9.430 μs ± 641.586 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
█▃ ▁
██▄▅▆▄▄▄▄▃▄▄▁▁▃▃▁▁▃▃▄▁▁▃▁▄▃▃▃▄▁▁▃▁▁▃▁▁▁▁▁▁▁▁▁▁▃▁▁▃▁▁▃▃▄▅▅▄▆ █
9.31 μs Histogram: log(frequency) by time 13.4 μ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_aosmul!(du_aos, D_SBP, u_aos)@show du_aos ≈ du_aos_rprintln("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 907 evaluations per sample.
Range (min … max): 120.180 ns … 173.050 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 120.465 ns ┊ GC (median): 0.00%
Time (mean ± σ): 121.945 ns ± 3.904 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▇█▃ ▁ ▃▃▂▂▂▁ ▂
███▅▅▄▄▆▆██▇▅▅▅▅▅▅▅▅▇██▇█▆▅▅▅▃▄▅▆▅▄▆▅▅▅████████▇▇▇▅▆▄▅▅▅▄▄▄▁▅ █
120 ns Histogram: log(frequency) by time 135 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 211 evaluations per sample.
Range (min … max): 351.815 ns … 727.137 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 421.777 ns ┊ GC (median): 0.00%
Time (mean ± σ): 419.500 ns ± 35.132 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▆▇▂▁▂ ▁▂ ▆█▄▁▃▆▆▅▄▄▄▄▄▃▃▂▂▂▂▁ ▂
▆▄███████▇█████▇▇▇█▆▇█▇▇▇██████████████████████████▆▆▆▆▅▅▆▅▅▅ █
352 ns Histogram: log(frequency) by time 510 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_full
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
Range (min … max): 9.762 μs … 34.654 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 10.310 μs ┊ GC (median): 0.00%
Time (mean ± σ): 10.283 μs ± 807.419 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▇ ██ ▁▂ ▂
█▅▅▆▃██▆██▅▅▇▃▅▅▁▅▄▅▆▆▃▁▄▄▄▁▁▃▁▁▁▁▄▃▁▃▃▁▁▃▁▄▄▄▅▄▁▅▆▃▄▇▅▁▁▁▃▄ █
9.76 μs Histogram: log(frequency) by time 15.3 μ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_soamul!(du_soa, D_SBP, u_soa)@show du_soa ≈ du_aosprintln("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 676 evaluations per sample.
Range (min … max): 177.975 ns … 274.596 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 185.938 ns ┊ GC (median): 0.00%
Time (mean ± σ): 188.001 ns ± 6.090 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▁▃▂ ▆███▄▂▂▂▃▂▂▂▁ ▂▃▄▄▃▂▁ ▂
▃▁████▆▅▆▆▇▅▇███████████████▇▇▆▆▅███▆▇▅▆█████████▇▇▇▆▇▆▅▄▅▄▆▆ █
178 ns Histogram: log(frequency) by time 209 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
Range (min … max): 34.229 μs … 91.749 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 34.643 μs ┊ GC (median): 0.00%
Time (mean ± σ): 36.554 μs ± 3.092 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▄█▆ ▂ ▂▄▆ ▁
████▆██▅▆███▇██▇▆▆▇▇▇▆▆██▇██████▇▆▆▅▅▆▇█▇▆▇▇▇▆▅▆▇▅▅▅▅▅▃▄▅▅▆ █
34.2 μs Histogram: log(frequency) by time 46.2 μs <
Memory estimate: 0 bytes, allocs estimate: 0.
D_full
BenchmarkTools.Trial: 10000 samples with 7 evaluations per sample.
Range (min … max): 4.666 μs … 11.230 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 4.911 μs ┊ GC (median): 0.00%
Time (mean ± σ): 5.077 μs ± 548.969 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▇█ ▂▃▁▂▇▅ ▃▁ ▄▅ ▂
██▅▅▅██████▇▇█▇▇██▇▆▇█▇▇▇█▇▆▇▆▆▅▇█████▇▆▆▆▇██▇▆▅▄▄▅▃▄▄▄▄▄▃▅ █
4.67 μs Histogram: log(frequency) by time 7.02 μs <
Memory estimate: 0 bytes, allocs estimate: 0.These results were obtained using the following versions.
using InteractiveUtilsversioninfo()using PkgPkg.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 × Intel(R) Xeon(R) 6973P-C
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-15.0.7 (ORCJIT, icelake-client)
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.97-DEV `~/work/SummationByPartsOperators.jl/SummationByPartsOperators.jl`