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): 20.029 ns … 53.520 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 20.903 ns ┊ GC (median): 0.00%
Time (mean ± σ): 21.196 ns ± 1.312 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▂▄█▄▃ ▂
▃▇▇██████▇█▇▄▅▄▂▃▃▃▃▃▃▂▂▂▂▂▂▂▂▁▁▂▂▁▁▂▁▂▁▁▂▂▁▂▁▂▂▂▂▂▂▂▂▂▂▂▂▂ ▃
20 ns Histogram: frequency by time 27.7 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 750 evaluations per sample.
Range (min … max): 168.076 ns … 235.763 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 169.932 ns ┊ GC (median): 0.00%
Time (mean ± σ): 171.228 ns ± 3.935 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▄▆█▆▄▂
▂▃▅▇██████▇▅▄▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▁▂▂▂▂▃▃▃▃▃▃▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ ▃
168 ns Histogram: frequency by time 183 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.12
Commit d93beab124c (2026-08-15 10:29 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.103-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 222 evaluations per sample.
Range (min … max): 334.005 ns … 546.486 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 336.986 ns ┊ GC (median): 0.00%
Time (mean ± σ): 340.207 ns ± 11.355 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▄▇██▇▆▆▅▄▃▂ ▁▂▂▂▂▁▁▁ ▂
████████████▇█▇▇▆▆▆▅▅▃▅▅▃▁▃▅▃▄▄▄▃▃▆████████████▇▇▇▆▆▆▆▆▅▅▅▅▅▆ █
334 ns Histogram: log(frequency) by time 383 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 7 evaluations per sample.
Range (min … max): 3.992 μs … 6.962 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 4.123 μs ┊ GC (median): 0.00%
Time (mean ± σ): 4.157 μs ± 187.683 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▂▅▆███▇▆▄▂ ▂
▆▇██████████▇▆▆▅▄▅▃▄▄▆▆▃▆▄▁▁▁▁▅▁▄▁▁▁▃▁▁▅▁▃▁▁▁▃▃▅▆▇██████▇▇▇ █
3.99 μs Histogram: log(frequency) by time 5.13 μ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.764 μs … 15.929 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 7.769 μs ┊ GC (median): 0.00%
Time (mean ± σ): 7.834 μs ± 370.211 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
█▂ ▁
███▇▆▄▆▄▄▄▁▃▁▁▄▅▆▅▅▄▄▃▃▁▁▃▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▃▆██▇▇▇▆▆▆ █
7.76 μs Histogram: log(frequency) by time 9.53 μ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.12
Commit d93beab124c (2026-08-15 10:29 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.13.0
[9f78cca6] SummationByPartsOperators v0.5.103-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 227 evaluations per sample.
Range (min … max): 327.890 ns … 735.890 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 330.709 ns ┊ GC (median): 0.00%
Time (mean ± σ): 333.462 ns ± 10.396 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▃▇███▇▆▅▄▃▂▁ ▁▂▂▂▂▁▁▁▁ ▃
██████████████▇███▇▇▇▆▄▅▅▃▁▃▄▃▁▁▁▄▁▁▅▁▃▁▁▆██████████▇█▇█▇▇▆▆▇ █
328 ns Histogram: log(frequency) by time 368 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
Di_SBP:
BenchmarkTools.Trial: 10000 samples with 170 evaluations per sample.
Range (min … max): 630.935 ns … 891.082 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 650.788 ns ┊ GC (median): 0.00%
Time (mean ± σ): 655.473 ns ± 16.709 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▂▄▆▆▇█▇▄▄▃▁
▁▁▁▁▂▃▄▆█████████████▇▅▅▄▃▃▂▂▂▂▁▁▁▁▁▁▁▂▂▂▂▃▃▃▃▃▃▂▂▂▂▂▂▁▁▁▁▁▁▁ ▃
631 ns Histogram: frequency by time 709 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 7 evaluations per sample.
Range (min … max): 4.056 μs … 8.950 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 4.102 μs ┊ GC (median): 0.00%
Time (mean ± σ): 4.136 μs ± 188.531 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▂▆█▇▄▁ ▂
██████▇▇▆▅▅▃▄▃▄▄▄▅▅▅▄▅▄▄▁▃▁▁▁▁▃▁▁▃▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▅▆███▇▇▇▇▆ █
4.06 μs Histogram: log(frequency) by time 5.09 μs <
Memory estimate: 0 bytes, allocs estimate: 0.
Di_banded:
BenchmarkTools.Trial: 10000 samples with 6 evaluations per sample.
Range (min … max): 5.625 μs … 11.113 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 5.667 μs ┊ GC (median): 0.00%
Time (mean ± σ): 5.717 μs ± 278.540 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▄██▄▁ ▂
█████▇▆▅▄▄▃▅▄▁▄▃▆▆▆▅▅▃▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▃▁▃▁▁▁▃▃▁▁▃▆█████▇▆▆▆▆ █
5.62 μs Histogram: log(frequency) by time 6.89 μ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): 97.374 μs … 138.255 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 98.295 μs ┊ GC (median): 0.00%
Time (mean ± σ): 99.946 μs ± 3.700 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▅██▇▆▅▄▃▃▃▃▂▃▃▃▃▂▂▁▁ ▂▃▃▂▂▁▁▁▁ ▂
█████████████████████▇▇▇▅▆▆▅▆▇▅███████████▇███▇▇▇▇▆▇▆▆▆▆▆▆▄▅ █
97.4 μs Histogram: log(frequency) by time 113 μ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.12
Commit d93beab124c (2026-08-15 10:29 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.13.0
[9f78cca6] SummationByPartsOperators v0.5.103-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 992 evaluations per sample.
Range (min … max): 38.181 ns … 69.438 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 39.292 ns ┊ GC (median): 0.00%
Time (mean ± σ): 39.698 ns ± 1.742 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▅▅█▆█▄▇▅█
▇███████████▇▆▆▆▄▄▄▃▃▃▃▂▃▃▂▂▂▂▂▂▂▂▂▂▁▁▁▂▂▂▂▂▂▂▂▃▂▂▂▂▂▂▂▂▂▂▂ ▄
38.2 ns Histogram: frequency by time 47.1 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 243 evaluations per sample.
Range (min … max): 309.058 ns … 477.658 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 337.535 ns ┊ GC (median): 0.00%
Time (mean ± σ): 339.702 ns ± 11.051 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▁▂▄▆▇▇██▇▇▄▄▂▂
▂▁▁▁▁▁▁▂▂▂▂▂▂▃▄▅▆▇███████████████▇▅▄▄▃▃▃▃▃▃▃▃▃▃▃▃▃▃▃▃▃▃▃▃▃▂▂▂ ▄
309 ns Histogram: frequency by time 376 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_full
BenchmarkTools.Trial: 10000 samples with 10 evaluations per sample.
Range (min … max): 1.036 μs … 3.219 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 1.049 μs ┊ GC (median): 0.00%
Time (mean ± σ): 1.057 μs ± 75.432 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
█
██▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▂ ▂
1.04 μs Histogram: frequency by time 1.64 μ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.221 μs … 3.636 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 1.231 μs ┊ GC (median): 0.00%
Time (mean ± σ): 1.242 μs ± 86.928 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
█▆
██▃▂▂▂▂▂▂▁▁▂▂▁▁▁▂▂▂▂▂▂▂▂▂▁▁▁▁▂▂▂▁▁▁▁▁▁▂▂▁▁▁▁▂▂▂▂▁▁▂▂▂▂▁▂▂▂ ▂
1.22 μs Histogram: frequency by time 1.83 μs <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 9 evaluations per sample.
Range (min … max): 2.022 μs … 7.266 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 2.062 μs ┊ GC (median): 0.00%
Time (mean ± σ): 2.083 μs ± 126.539 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▁█▆
▃███▇▅▄▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▁▁▂▁▁▂▂▂▂▂▂▂▁▂▁▂▁▁▁▁▁▁▂▂▂▂▂▂▂ ▂
2.02 μs Histogram: frequency by time 2.78 μs <
Memory estimate: 240 bytes, allocs estimate: 5.
D_full
BenchmarkTools.Trial: 10000 samples with 6 evaluations per sample.
Range (min … max): 6.037 μs … 11.305 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 6.096 μs ┊ GC (median): 0.00%
Time (mean ± σ): 6.145 μs ± 242.004 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▁▆█▇▅▂ ▁ ▂
███████▇▆▇▇▆▆▆▇▆▇▅▆▆▅▄▃▄▁▃▄▃▁▁▁▁▁▁▃▁▁▁▃▁▁▃▃▃▁▃▁▃▃▅█████▇▇▇▆ █
6.04 μs Histogram: log(frequency) by time 7.29 μ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 709 evaluations per sample.
Range (min … max): 177.598 ns … 275.261 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 179.801 ns ┊ GC (median): 0.00%
Time (mean ± σ): 181.246 ns ± 4.545 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▁▅█▅
▁▂▃▅█████▆▄▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▂▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
178 ns Histogram: frequency by time 197 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 193 evaluations per sample.
Range (min … max): 487.663 ns … 11.456 μs ┊ GC (min … max): 0.00% … 92.99%
Time (median): 562.228 ns ┊ GC (median): 0.00%
Time (mean ± σ): 564.956 ns ± 111.920 ns ┊ GC (mean ± σ): 0.19% ± 0.93%
▁▂▄▅▇██▇▇▅▅▃▂
▁▁▁▁▁▁▁▁▁▁▂▂▂▂▃▃▃▃▃▃▄▅▆▆▇███████████████▇▆▆▅▅▄▅▄▄▄▃▃▂▂▂▂▂▂▂▁▁ ▄
488 ns Histogram: frequency by time 625 ns <
Memory estimate: 32 bytes, allocs estimate: 1.
D_full
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
Range (min … max): 12.638 μs … 27.280 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 12.759 μs ┊ GC (median): 0.00%
Time (mean ± σ): 12.865 μs ± 810.440 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▇█▄ ▁
███▇▄▅▃▄▃▅▄▄▁▄▁▃▃▄▄▅▅▁▃▁▁▁▅▆▃▃▃▃▃▃▁▁▁▁▁▁▁▁▁▁▃▁▁▃▁▁▁▁▁▁▁▁▁▄▆▇ █
12.6 μs Histogram: log(frequency) by time 18.9 μ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 709 evaluations per sample.
Range (min … max): 176.509 ns … 222.262 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 178.460 ns ┊ GC (median): 0.00%
Time (mean ± σ): 179.758 ns ± 3.613 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▃▇█▆▃
▂▂▃▄▆█████▇▅▄▃▃▂▃▂▂▂▂▂▂▂▂▂▂▂▂▁▁▁▂▂▁▁▂▂▂▂▃▃▃▃▃▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂ ▃
177 ns Histogram: frequency by time 191 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 197 evaluations per sample.
Range (min … max): 454.381 ns … 673.690 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 520.315 ns ┊ GC (median): 0.00%
Time (mean ± σ): 521.729 ns ± 19.267 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▂▃▄▆██▇▃▂
▁▁▁▁▁▁▁▁▁▁▂▂▂▃▂▂▂▂▂▃▂▄▅▅█▆▇██████████▇▇▇▅▅▄▃▄▃▃▃▃▃▃▃▃▂▂▂▂▁▁▁▁ ▃
454 ns Histogram: frequency by time 576 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_full
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
Range (min … max): 12.568 μs … 34.481 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 12.699 μs ┊ GC (median): 0.00%
Time (mean ± σ): 12.822 μs ± 883.733 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▇█▄ ▁
███▇▆▅▃▁▄▄▄▃▁▄▃▁▄▄▄▄▃▄▄▅▅▅▄▄▆▆▅▅▁▄▅▅▅▁▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▃▁▁▁▇▆ █
12.6 μs Histogram: log(frequency) by time 18.8 μ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 656 evaluations per sample.
Range (min … max): 187.213 ns … 272.584 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 190.848 ns ┊ GC (median): 0.00%
Time (mean ± σ): 192.315 ns ± 5.117 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▁▄▃█▆▇▅▃▄▂
▂▂▃▅▇███████████▇▅▅▄▄▃▂▂▂▁▁▁▁▁▁▁▁▁▂▂▂▂▃▂▃▂▃▂▃▂▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁ ▃
187 ns Histogram: frequency by time 206 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
Range (min … max): 39.599 μs … 77.164 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 39.819 μs ┊ GC (median): 0.00%
Time (mean ± σ): 40.189 μs ± 1.648 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▆█▆▆▅▂ ▁ ▂
███████▆▆▆▅▅▅▅▁▅▆▆▆▅▄▃▁▅▄▄▃▁▃▁▄▁▁▄▃▁▁▁▁▁▁▅▁▁▁▃▁▄▆▇████▇█▇▇▇ █
39.6 μs Histogram: log(frequency) by time 46.9 μs <
Memory estimate: 0 bytes, allocs estimate: 0.
D_full
BenchmarkTools.Trial: 10000 samples with 6 evaluations per sample.
Range (min … max): 5.209 μs … 10.666 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 5.256 μs ┊ GC (median): 0.00%
Time (mean ± σ): 5.299 μs ± 234.366 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▃██▆▂ ▂
██████▇▇▆▄▅▅▅▄▄▄▆▆▅▅▄▃▃▁▃▁▁▁▁▁▁▃▁▁▃▁▁▃▁▃▁▁▁▁▁▁▁▄▆████▇▆▇▆▅▆ █
5.21 μs Histogram: log(frequency) by time 6.44 μ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.12
Commit d93beab124c (2026-08-15 10:29 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.22
[09ab397b] StructArrays v0.7.3
[9f78cca6] SummationByPartsOperators v0.5.103-DEV `~/work/SummationByPartsOperators.jl/SummationByPartsOperators.jl`