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 995 evaluations per sample.
Range (min … max): 26.552 ns … 136.415 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 27.348 ns ┊ GC (median): 0.00%
Time (mean ± σ): 27.658 ns ± 1.949 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▂▆█▇▅▅▂
▄████████▇▅▃▃▃▂▂▂▂▂▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▁▂▂▂▁▁▂▂▂▁▂▂▂▂▂▂▂▂ ▃
26.6 ns Histogram: frequency by time 35.8 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 585 evaluations per sample.
Range (min … max): 203.458 ns … 345.552 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 213.716 ns ┊ GC (median): 0.00%
Time (mean ± σ): 215.515 ns ± 7.230 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▅▅▆▄█▇▆▂▄▂▃ ▂
▁▁▁▁▂▃▃▃▆▇█▇██████████████▇▄▅▃▄▃▂▃▃▂▂▂▃▃▃▃▄▃▃▃▃▃▂▂▂▂▂▂▁▂▁▁▁▁▁ ▃
203 ns Histogram: frequency by time 236 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 7763 64-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 `~/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 203 evaluations per sample.
Range (min … max): 382.438 ns … 527.335 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 387.128 ns ┊ GC (median): 0.00%
Time (mean ± σ): 390.005 ns ± 11.963 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▆▅ ▄█▃
▄██▇▆███▆▄▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▁▁▂▁▁▁▂▁▂▂▁▁▂▁▁▂▂▂▁▂▁▁▁▂▂▃▃▂▂▃▃▃▂▂▂▂ ▃
382 ns Histogram: frequency by time 432 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): 4.537 μs … 9.146 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 4.597 μs ┊ GC (median): 0.00%
Time (mean ± σ): 4.645 μs ± 252.349 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▃▇█▇▅▁ ▁▁ ▁▁▁ ▂
██████▇██▇▅▁▃▃▅▅▆▅▄▄▃▁▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▃▁▁▁▃▁▃▁▁▁▁▁▁▁▇████ █
4.54 μs Histogram: log(frequency) by time 5.81 μ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.497 μs … 35.657 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 9.527 μs ┊ GC (median): 0.00%
Time (mean ± σ): 9.663 μs ± 1.131 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
█ ▁
█▆▄▁▃▆▅▃▁▁▃▃▃▁▄▄▃▄▃▁▁▃▃▁▁▁▁▁▁▃▃▁▁▁▃▃▁▁▃▁▃▃▁▃▁▁▁▁▃▁▃▃▃▁▁▃▃▇ █
9.5 μs Histogram: log(frequency) by time 17.6 μ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 7763 64-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.14.0
[9f78cca6] SummationByPartsOperators v0.5.103 `~/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 204 evaluations per sample.
Range (min … max): 383.657 ns … 591.299 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 386.897 ns ┊ GC (median): 0.00%
Time (mean ± σ): 390.649 ns ± 12.566 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▄█▅
▃▆████▅▄▃▃▂▂▂▂▂▂▂▂▂▂▂▂▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▁▂▁▂▂▃▃▃▃▂▂▂▂▂ ▃
384 ns Histogram: frequency by time 432 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
Di_SBP:
BenchmarkTools.Trial: 10000 samples with 22 evaluations per sample.
Range (min … max): 954.955 ns … 2.437 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 965.000 ns ┊ GC (median): 0.00%
Time (mean ± σ): 977.868 ns ± 72.185 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▅█▆ ▁ ▁ ▁
███▇▅▄██▄▁▁▁▁▁▄▅▅▅▃▁▃▁▁▃▃▁▃▁▃▄▁▄▅▆▇▅▄▁▄▅▄▃▁▃▃▁▁▁▁▄▁▁▁▁▃▁▆██▇ █
955 ns Histogram: log(frequency) by time 1.34 μ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 7 evaluations per sample.
Range (min … max): 4.908 μs … 8.042 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 4.958 μs ┊ GC (median): 0.00%
Time (mean ± σ): 5.006 μs ± 243.362 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▂▇█▆ ▂ ▁▁▁ ▂
█████▄███▅▁▃▁▄▄▇▅▅▅▁▄▃▃▁▁▁▁▁▁▁▁▁▁▃▁▁▁▁▃▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▆████ █
4.91 μs Histogram: log(frequency) by time 6.16 μs <
Memory estimate: 0 bytes, allocs estimate: 0.
Di_banded:
BenchmarkTools.Trial: 10000 samples with 5 evaluations per sample.
Range (min … max): 6.911 μs … 12.173 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 6.961 μs ┊ GC (median): 0.00%
Time (mean ± σ): 7.029 μs ± 346.699 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▆█▆▁ ▁ ▁▁ ▁
████▃▆██▄▁▁▁▁▁▆▆▆▅▄▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▁▁▁▁▁▁▁▃▄▇████ █
6.91 μs Histogram: log(frequency) by time 8.69 μ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): 115.516 μs … 171.480 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 116.207 μs ┊ GC (median): 0.00%
Time (mean ± σ): 117.792 μs ± 4.110 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▂▆█▇▄▃▃ ▂▃▄▂▂ ▂
█████████▇████▇▆▅▅▄▄▄▄▄▁▁▄▃▃▄▃▃▃▃▁▅▆▆▅▅▄▇███████▇▇█▇▇▆▆▅▆▅▆▆▆ █
116 μs Histogram: log(frequency) by time 130 μ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 7763 64-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.14.0
[9f78cca6] SummationByPartsOperators v0.5.103 `~/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 988 evaluations per sample.
Range (min … max): 46.626 ns … 85.107 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 47.456 ns ┊ GC (median): 0.00%
Time (mean ± σ): 47.994 ns ± 2.123 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▃▅█▆▄▃▁
▄███████▆▅▅▃▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▂▂▁▁▁▂▂▂▁▁▁▁▂▂▁▂▂▂▂▃▂▂▂▂▂▂▂ ▃
46.6 ns Histogram: frequency by time 56.6 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 228 evaluations per sample.
Range (min … max): 324.728 ns … 1.734 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 377.329 ns ┊ GC (median): 0.00%
Time (mean ± σ): 402.222 ns ± 102.004 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▃▃▆▇██▇▆▅▄▃▂ ▂
██████████████▇█▇▆▇▅▇▇███▇███▇▇▇▅▆▆▆▆▅▆▆▆▆▆▅▅▆▆▅▅▅▄▅▂▆▅▆▄▆▆▇▆ █
325 ns Histogram: log(frequency) by time 912 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_full
BenchmarkTools.Trial: 10000 samples with 10 evaluations per sample.
Range (min … max): 1.174 μs … 2.506 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 1.188 μs ┊ GC (median): 0.00%
Time (mean ± σ): 1.200 μs ± 93.063 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
██ ▁ ▂
███▄▄▅█▅▁▁▁▁▁▁▄▅▄▃▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▇ █
1.17 μs Histogram: log(frequency) by time 2 μ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.277 μs … 3.048 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 1.286 μs ┊ GC (median): 0.00%
Time (mean ± σ): 1.300 μs ± 99.479 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
█▆ ▁
██▄▄▃▅█▃▁▁▁▁▁▅▄▄▄▃▁▁▁▁▁▁▁▁▁▃▃▁▃▃▁▄▄▃▁▃▁▁▁▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▅▇ █
1.28 μs Histogram: log(frequency) by time 2.08 μs <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 9 evaluations per sample.
Range (min … max): 2.436 μs … 5.340 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 2.550 μs ┊ GC (median): 0.00%
Time (mean ± σ): 2.577 μs ± 164.165 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▂▅▇██▇▆▅▃▂ ▁▁▁ ▂
▄▆███████████▇▆▄▅▅▅▆▅▃▄▆▄▃▃▃▁▁▁▁▁▁▃▄▁▁▁▁▁▃▁▁▁▁▁▁▃▁▁▁▁▄▅████ █
2.44 μs Histogram: log(frequency) by time 3.49 μs <
Memory estimate: 240 bytes, allocs estimate: 5.
D_full
BenchmarkTools.Trial: 10000 samples with 5 evaluations per sample.
Range (min … max): 6.250 μs … 13.631 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 6.288 μs ┊ GC (median): 0.00%
Time (mean ± σ): 6.358 μs ± 375.565 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▆█▄▁ ▂ ▁▁ ▁
█████▇█▆▃▄▄▁▁▁▅▆▆▅▄▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▁▁▁▃▁▁▁▁▃▇███ █
6.25 μs Histogram: log(frequency) by time 8.01 μ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 530 evaluations per sample.
Range (min … max): 212.851 ns … 328.651 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 216.177 ns ┊ GC (median): 0.00%
Time (mean ± σ): 218.253 ns ± 5.961 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▃▆█▇▆▄▄▃▁ ▂▃▄▃▂▁▁ ▂
▃▃▁▄▅███████████▆▇██▇▇▅▄▁▄▁▁▁▄▁▁▃▁▄▁▄▄▃▁▁▁▁▃▁▃▃▆█████████▆▆▇▆ █
213 ns Histogram: log(frequency) by time 235 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 185 evaluations per sample.
Range (min … max): 559.259 ns … 10.879 μs ┊ GC (min … max): 0.00% … 92.85%
Time (median): 588.124 ns ┊ GC (median): 0.00%
Time (mean ± σ): 595.060 ns ± 106.575 ns ┊ GC (mean ± σ): 0.17% ± 0.93%
▂▃▅▆▇▇▇██▇▇▇▅▃▂
▁▁▂▂▄▆█████████████████▇▆▆▅▅▄▃▄▃▃▃▃▃▃▃▃▃▃▃▃▃▃▃▃▂▂▁▂▂▁▁▁▁▁▁▁▁▁ ▄
559 ns Histogram: frequency by time 664 ns <
Memory estimate: 32 bytes, allocs estimate: 1.
D_full
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
Range (min … max): 13.895 μs … 32.881 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 13.966 μs ┊ GC (median): 0.00%
Time (mean ± σ): 14.107 μs ± 1.062 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
█▅ ▁
██▁▁▄▆▇▁▁▁▁▁▁▄▄▄▄▄▁▃▁▃▁▁▁▁▁▁▁▁▃▁▁▁▃▃▃▄▃▁▃▁▃▁▁▁▃▁▁▃▃▁▁▁▁▁▁▆▇ █
13.9 μs Histogram: log(frequency) by time 22 μ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 550 evaluations per sample.
Range (min … max): 208.736 ns … 295.805 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 212.105 ns ┊ GC (median): 0.00%
Time (mean ± σ): 214.075 ns ± 5.803 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▁▇█▄
▂▂▂▃▄████▇▅▃▃▂▂▂▂▂▂▂▂▂▂▂▂▁▁▁▁▂▁▁▂▂▂▁▂▂▂▃▄▃▃▃▂▂▂▂▂▂▂▂▁▂▁▂▁▂▁▂▂ ▃
209 ns Histogram: frequency by time 235 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 191 evaluations per sample.
Range (min … max): 524.016 ns … 1.146 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 547.251 ns ┊ GC (median): 0.00%
Time (mean ± σ): 553.114 ns ± 26.865 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▁▂▃▄▅▆██▆▅▄▄▁
▂▃▄▆█████████████▇▆▅▅▄▃▃▃▃▃▃▃▄▄▄▄▄▄▃▃▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂ ▄
524 ns Histogram: frequency by time 639 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_full
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
Range (min … max): 13.906 μs … 26.750 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 13.976 μs ┊ GC (median): 0.00%
Time (mean ± σ): 14.127 μs ± 1.048 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
█▅ ▁
██▅▅▄▇█▁▁▁▁▁▁▁▃▅▄▅▁▃▃▁▁▁▃▁▁▃▁▁▁▁▃▄▁▃▃▄▄▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▃▃▆▇ █
13.9 μs Histogram: log(frequency) by time 22.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_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 478 evaluations per sample.
Range (min … max): 222.843 ns … 315.987 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 226.533 ns ┊ GC (median): 0.00%
Time (mean ± σ): 228.575 ns ± 6.501 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▂▆█▇▇▅▂
▂▂▃▄▇████████▆▅▄▃▃▂▂▂▂▂▂▂▂▂▁▂▁▁▁▁▁▁▁▂▁▁▁▁▁▁▂▂▂▂▃▃▄▄▄▃▃▃▃▂▂▂▂▂ ▃
223 ns Histogram: frequency by time 247 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
Range (min … max): 49.201 μs … 152.906 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 49.382 μs ┊ GC (median): 0.00%
Time (mean ± σ): 49.920 μs ± 2.727 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▆█▃ ▂ ▁▂▁▁ ▁
███▅▅▆█▆▃▁▁▁▁▃▇▆▇▃▄▁▃▃▁▁▁▁▁▁▁▁▁▁▁▄▁▁▁▁▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▆█████▇ █
49.2 μs Histogram: log(frequency) by time 58.1 μs <
Memory estimate: 0 bytes, allocs estimate: 0.
D_full
BenchmarkTools.Trial: 10000 samples with 6 evaluations per sample.
Range (min … max): 5.363 μs … 11.179 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 5.415 μs ┊ GC (median): 0.00%
Time (mean ± σ): 5.473 μs ± 295.080 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▃██▄ ▁▂ ▁▁▁ ▂
████▇▄██▇▃▁▁▁▁▁▄▇▅▄▄▄▃▁▁▁▁▁▁▁▁▁▁▁▃▃▁▁▁▁▁▁▁▁▁▁▁▁▁▃▃▁▁▁▄▇████ █
5.36 μs Histogram: log(frequency) by time 6.86 μ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 7763 64-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.23
[09ab397b] StructArrays v0.7.4
[9f78cca6] SummationByPartsOperators v0.5.103 `~/work/SummationByPartsOperators.jl/SummationByPartsOperators.jl`