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 998 evaluations per sample.
Range (min … max): 16.216 ns … 36.457 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 16.779 ns ┊ GC (median): 0.00%
Time (mean ± σ): 16.876 ns ± 0.850 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▄▅▅▇█▇▂ ▂
▇███████▇▆▇▇█▇▆▅▅▅▄▆▅▆▆▃▁▃▃▁▄▄▁▄▁▁▁▃▄▃▁▁▁▃▁▃▃▁▁▁▃▁▁▃▁▄▅▄▅▇█ █
16.2 ns Histogram: log(frequency) by time 22.4 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 914 evaluations per sample.
Range (min … max): 114.602 ns … 180.586 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 118.700 ns ┊ GC (median): 0.00%
Time (mean ± σ): 120.324 ns ± 4.768 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▃▅█▆▃▁
▁▂▃▃▅▇███████▇▇▆▅▄▃▃▃▃▃▃▃▃▃▃▂▃▂▂▂▂▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▂▂▁▁▁▁ ▃
115 ns Histogram: frequency by time 137 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 9V45 96-Core Processor
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-15.0.7 (ORCJIT, generic)
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 `~/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 352 evaluations per sample.
Range (min … max): 255.952 ns … 331.148 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 262.210 ns ┊ GC (median): 0.00%
Time (mean ± σ): 263.490 ns ± 6.099 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▄▅▃▃▁ ▅█▆▅▅▄▃▁▁ ▁▁
▁▃▆█████▇▆█████████████▅▃▂▂▁▁▁▁▁▁▁▁▁▁▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▁ ▃
256 ns Histogram: frequency by time 283 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 9 evaluations per sample.
Range (min … max): 2.798 μs … 4.447 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 2.859 μs ┊ GC (median): 0.00%
Time (mean ± σ): 2.877 μs ± 110.077 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▁▅█████▇▆▄▂▁ ▁ ▁ ▃
████████████▇██▇█▇▆▆▅▆▅▃▃▃▁▁▄▃▄▁▁▁▃▄▁▃▁▃▁▁▁▃▁▁▁▁▁▃▁▃▃▇█████ █
2.8 μs Histogram: log(frequency) by time 3.48 μ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 6 evaluations per sample.
Range (min … max): 5.543 μs … 9.931 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 5.680 μs ┊ GC (median): 0.00%
Time (mean ± σ): 5.732 μs ± 231.362 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▄▇▇▇██▇▆▅▄▄▃▃▂▂▂▁▁▁ ▁▁▁▁▁ ▃
▄▆███████████████████████▇▇▆▆▅▆▇▃▅▅▅▄▅▃▄▃▃▅▃▁▃▁▃▃▄▇████████ █
5.54 μs Histogram: log(frequency) by time 6.67 μ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 9V45 96-Core Processor
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-15.0.7 (ORCJIT, generic)
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.12.0
[9f78cca6] SummationByPartsOperators v0.5.97 `~/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 340 evaluations per sample.
Range (min … max): 258.356 ns … 639.691 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 264.629 ns ┊ GC (median): 0.00%
Time (mean ± σ): 266.728 ns ± 8.647 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▅▄▁▁▃▇▇█▁
▁▃██████████▇▇▆▆▆▅▃▃▂▂▂▂▁▁▂▂▂▂▂▂▂▂▂▃▂▂▂▂▂▁▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▃
258 ns Histogram: frequency by time 297 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
Di_SBP:
BenchmarkTools.Trial: 10000 samples with 200 evaluations per sample.
Range (min … max): 402.355 ns … 700.695 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 410.360 ns ┊ GC (median): 0.00%
Time (mean ± σ): 413.962 ns ± 13.474 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▅█▃▂▃▆▅▂
▁▇████████▆▆▃▄▃▃▄▂▁▁▁▁▁▁▁▁▂▂▂▂▂▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
402 ns Histogram: frequency by time 468 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.317 μs … 7.469 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 3.371 μs ┊ GC (median): 0.00%
Time (mean ± σ): 3.406 μs ± 168.226 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▄██▇█▇▅▄▃▃▃▂▂▁▁▁ ▁ ▃
████████████████████▇▇▆▇▆▆▆▆▆▆▆▇▄▆▄▆▅▅▅▅▅▄▄▁▄▄▄▁▁▆▆██████▇▆ █
3.32 μs Histogram: log(frequency) by time 4.12 μs <
Memory estimate: 0 bytes, allocs estimate: 0.
Di_banded:
BenchmarkTools.Trial: 10000 samples with 7 evaluations per sample.
Range (min … max): 4.464 μs … 9.119 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 4.560 μs ┊ GC (median): 0.00%
Time (mean ± σ): 4.615 μs ± 222.981 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▃▆▇██▇▅▅▄▄▃▃▂▂▂▁▁▁ ▁▁ ▂
▆███████████████████████▇▇▆▆▆▆▅▄▆▅▃▅▅▄▅▁▄▁▁▁▅▇██████▇▆▇▅▆▆▆ █
4.46 μs Histogram: log(frequency) by time 5.54 μ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): 64.886 μs … 147.510 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 68.783 μs ┊ GC (median): 0.00%
Time (mean ± σ): 69.633 μs ± 2.927 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▂▆██▇▆▅▄▃▃▂▂▂▂▁ ▁▃▃▃▂▁ ▂
▄▆▆▆▆▅▁▅▄▄▄▃▄██████████████████▇▇▇▇▇▆▆▆▆▇▇▆▆▅████████▇▆▆▅▆▆▅ █
64.9 μs Histogram: log(frequency) by time 78.9 μ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 9V45 96-Core Processor
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-15.0.7 (ORCJIT, generic)
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.12.0
[9f78cca6] SummationByPartsOperators v0.5.97 `~/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 995 evaluations per sample.
Range (min … max): 28.938 ns … 401.221 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 29.582 ns ┊ GC (median): 0.00%
Time (mean ± σ): 31.237 ns ± 12.532 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
█▇▄▂▁ ▂▁ ▁
█████▇▇███▇▇▆▆▆▅▆▄▄▄▅▄▄▄▄▅▅▄▁▄▅▅▅▁▄▄▄▄▁▄▃▁▄▄▃▄▁▁▄▄▄▃▅▄▄▄▄▅▃▄ █
28.9 ns Histogram: log(frequency) by time 72.1 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 560 evaluations per sample.
Range (min … max): 202.320 ns … 813.159 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 208.650 ns ┊ GC (median): 0.00%
Time (mean ± σ): 209.825 ns ± 8.904 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▄█▃
▁▃▇██▅▅▄▄▄▄▇████▇▅▅▄▃▃▃▃▃▂▂▂▂▂▂▂▂▂▃▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
202 ns Histogram: frequency by time 231 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_full
BenchmarkTools.Trial: 10000 samples with 33 evaluations per sample.
Range (min … max): 932.303 ns … 3.084 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 954.152 ns ┊ GC (median): 0.00%
Time (mean ± σ): 963.287 ns ± 46.612 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▁▅▇██▇▇▅▄▃▃▂▂▁▁▁▁ ▁▁▁▁ ▂
▇████████████████████████▆▇▆▆▆▆▆▆▅▄▆▅▅▄▃▄▃▄▃▃▁▁▃▅▁▃▁▃▅██████ █
932 ns Histogram: log(frequency) by time 1.13 μ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 120 evaluations per sample.
Range (min … max): 752.950 ns … 3.534 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 765.392 ns ┊ GC (median): 0.00%
Time (mean ± σ): 773.691 ns ± 36.955 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▁██▇▆▂▂▁
▂▄████████▇▅▄▄▃▃▂▂▂▂▂▂▂▂▁▁▂▂▂▂▃▃▃▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
753 ns Histogram: frequency by time 858 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 10 evaluations per sample.
Range (min … max): 1.247 μs … 7.254 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 1.277 μs ┊ GC (median): 0.00%
Time (mean ± σ): 1.298 μs ± 106.903 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▄▇█▇▆▅▅▄▄▃▂▂▂▁▁▁ ▂
███████████████████▆▇▇▇▅▅▅▆▆▁▃▁▄▃▄▄▄▁▄▄▃▁▄▁▁▁▁▁▁▁▃▃▃▃▃▃▄▅▇▇ █
1.25 μs Histogram: log(frequency) by time 1.83 μs <
Memory estimate: 240 bytes, allocs estimate: 5.
D_full
BenchmarkTools.Trial: 10000 samples with 7 evaluations per sample.
Range (min … max): 4.771 μs … 9.486 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 4.975 μs ┊ GC (median): 0.00%
Time (mean ± σ): 5.027 μs ± 213.586 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▁▃▅▆▇▇███▇▆▆▅▄▃▃▂▁▁▁ ▁▁▁▁ ▁ ▃
▆▄▆██████████████████████▇▇██▇▇▇▇▆▇▅▇▆▆▅▆▅▅▇▅▄▆▆▇▇████████▆ █
4.77 μs Histogram: log(frequency) by time 5.87 μ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 944 evaluations per sample.
Range (min … max): 98.707 ns … 156.006 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 100.638 ns ┊ GC (median): 0.00%
Time (mean ± σ): 101.981 ns ± 4.519 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▅█▄▁▁
▃█████▇▅▃▃▃▃▃▃▃▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▂▂▂ ▃
98.7 ns Histogram: frequency by time 128 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 202 evaluations per sample.
Range (min … max): 380.965 ns … 11.249 μs ┊ GC (min … max): 0.00% … 93.88%
Time (median): 390.832 ns ┊ GC (median): 0.00%
Time (mean ± σ): 396.308 ns ± 111.196 ns ┊ GC (mean ± σ): 0.27% ± 0.94%
▂▅█▇▃▃▃▂
▁▃▅▇▇████████▇▅▄▃▃▂▂▂▂▂▂▂▁▁▁▂▂▂▂▂▂▂▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
381 ns Histogram: frequency by time 444 ns <
Memory estimate: 32 bytes, allocs estimate: 1.
D_full
BenchmarkTools.Trial: 10000 samples with 5 evaluations per sample.
Range (min … max): 6.854 μs … 11.693 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 7.027 μs ┊ GC (median): 0.00%
Time (mean ± σ): 7.076 μs ± 267.982 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▂▃▄▅▆▇██▇▆▄▃▂▁▁ ▁▁▁ ▂
█████████████████████▇▆▆▅▅▅▅▅▁▃▁▅▅▆▅▄▅▅▃▁▁▁▄▅▁▅▅▆▆▇████████ █
6.85 μs Histogram: log(frequency) by time 8.27 μ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 962 evaluations per sample.
Range (min … max): 85.012 ns … 155.720 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 86.886 ns ┊ GC (median): 0.00%
Time (mean ± σ): 87.479 ns ± 2.370 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▂▅███▆▅▇▇▆▅▄▁
▁▃▇█████████████▆▅▄▃▂▂▂▂▂▁▂▁▂▁▁▁▁▁▂▂▂▂▂▃▂▂▂▂▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁ ▃
85 ns Histogram: frequency by time 95.6 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 205 evaluations per sample.
Range (min … max): 372.800 ns … 705.346 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 381.493 ns ┊ GC (median): 0.00%
Time (mean ± σ): 386.956 ns ± 18.595 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▁▅█▄▁
▃█████▇▆▅▅▄▃▃▂▂▃▃▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▁▂▂▁▂▂▂▂▂▂▁▂▂▂▂▂ ▃
373 ns Histogram: frequency by time 489 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_full
BenchmarkTools.Trial: 10000 samples with 5 evaluations per sample.
Range (min … max): 6.990 μs … 16.603 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 7.129 μs ┊ GC (median): 0.00%
Time (mean ± σ): 7.185 μs ± 279.954 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▂▅▆▆▇███▆▅▄▃▂▁▁▁▁ ▁▁▁▁ ▃
████████████████████▇▇█▇▆▇▆▇▆▆▆▅▅▅▅▅▄▅▄▅▅▅▄▁▅▄▄▁▅▇▆████████ █
6.99 μs Histogram: log(frequency) by time 8.31 μ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 839 evaluations per sample.
Range (min … max): 146.799 ns … 191.502 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 150.570 ns ┊ GC (median): 0.00%
Time (mean ± σ): 151.773 ns ± 3.664 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▅▇▇█▄▄▄▃▁
▁▁▂▃▅▆▇█████████▇▅▄▃▃▃▂▂▂▃▃▃▃▃▄▃▃▃▃▃▃▂▂▂▂▁▂▂▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▃
147 ns Histogram: frequency by time 165 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
Range (min … max): 29.424 μs … 64.316 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 30.666 μs ┊ GC (median): 0.00%
Time (mean ± σ): 31.105 μs ± 1.632 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▁▅▆█▆▄▃▁
▁▂▃▆█████████▆▆▆▅▄▄▃▃▃▂▂▂▁▁▁▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▃
29.4 μs Histogram: frequency by time 37 μs <
Memory estimate: 0 bytes, allocs estimate: 0.
D_full
BenchmarkTools.Trial: 10000 samples with 8 evaluations per sample.
Range (min … max): 3.305 μs … 20.457 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 3.479 μs ┊ GC (median): 0.00%
Time (mean ± σ): 3.502 μs ± 258.358 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▄▇█▇▃
▂▂▂▃▅▆▆▇██████▇▄▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▁▁▂▂▂▁▂▂▁▁▁▁▁▂▁▂▂▂▂▂▂▂▂▂▂▂▂▂ ▃
3.3 μs Histogram: frequency by time 4.22 μ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 9V45 96-Core Processor
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-15.0.7 (ORCJIT, generic)
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.20
[09ab397b] StructArrays v0.7.3
[9f78cca6] SummationByPartsOperators v0.5.97 `~/work/SummationByPartsOperators.jl/SummationByPartsOperators.jl`