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 996 evaluations per sample.
Range (min … max): 25.791 ns … 49.279 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 26.545 ns ┊ GC (median): 0.00%
Time (mean ± σ): 26.850 ns ± 1.592 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▄█▁▂▃▂
▃███████▆▃▄▃▃▃▃▂▂▂▂▂▂▂▂▁▂▁▂▁▂▁▂▂▁▂▁▁▂▂▂▁▂▂▂▁▁▁▁▂▂▁▁▂▂▂▂▂▂▂▂ ▃
25.8 ns Histogram: frequency by time 34.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 585 evaluations per sample.
Range (min … max): 201.369 ns … 475.164 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 212.415 ns ┊ GC (median): 0.00%
Time (mean ± σ): 214.598 ns ± 10.600 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▁▂▅▆▆██▇▆▅▃▂▁
▁▁▁▁▁▂▃▄▅▆█████████████▇▆▅▄▄▃▃▃▃▂▃▃▃▄▃▃▄▃▃▃▃▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁ ▃
201 ns Histogram: frequency by time 239 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.101 `~/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): 383.133 ns … 761.177 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 388.611 ns ┊ GC (median): 0.00%
Time (mean ± σ): 392.313 ns ± 13.194 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▄█▆
▂▃▃▅███▇▄▃▃▂▂▂▂▂▂▂▂▂▂▂▂▁▂▁▁▁▁▂▁▂▁▁▁▁▂▂▁▁▁▂▂▂▂▂▂▃▃▃▂▂▂▂▂▂▂▂▂▂▂ ▃
383 ns Histogram: frequency by time 439 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.474 μs … 17.630 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 4.609 μs ┊ GC (median): 0.00%
Time (mean ± σ): 4.701 μs ± 581.499 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▃▆██▆▃▁ ▁▂▂ ▂
▆████████▄▄▆▆▅▅▅▄▁▁▁▁▅▃▄▅▁▄▅▄▄▅▄▄▃▃▃▆████▇▆▅▁▄▄▄▅▆▅▄▅▄▄▁▁▃▃ █
4.47 μs Histogram: log(frequency) by time 6.47 μ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 … 32.040 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 9.518 μs ┊ GC (median): 0.00%
Time (mean ± σ): 9.643 μs ± 1.016 μ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.12.0
[9f78cca6] SummationByPartsOperators v0.5.101 `~/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 203 evaluations per sample.
Range (min … max): 383.182 ns … 530.355 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 386.143 ns ┊ GC (median): 0.00%
Time (mean ± σ): 389.835 ns ± 12.281 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▃▇██▇▆▄▄▃▂ ▂▃▃▂▁ ▂
███████████▆▅▆▇▇▇▆▅▃▄▁▃▁▁▁▁▁▃▃▄▁▁▁▁▃▁▃▃▁▁▁▃▃▃▁▃▁▇███████▇▇▆▅▅ █
383 ns Histogram: log(frequency) by time 434 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
Di_SBP:
BenchmarkTools.Trial: 10000 samples with 16 evaluations per sample.
Range (min … max): 984.938 ns … 2.253 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 993.125 ns ┊ GC (median): 0.00%
Time (mean ± σ): 1.004 μs ± 71.342 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
██ ▁ ▂
███▄▁▃█▆▁▃▁▁▁▃▁▆▄▄▄▁▄▄▃▄▅▄▄▃▄▅▁▄▁▃▁▁▃▁▄▄▁▁▁▁▁▃▁▄▁▁▁▁▁▃▁▁▃▁▅█ █
985 ns Histogram: log(frequency) by time 1.5 μ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 6 evaluations per sample.
Range (min … max): 5.016 μs … 10.159 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 5.091 μs ┊ GC (median): 0.00%
Time (mean ± σ): 5.141 μs ± 273.065 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▂▆██▆▂ ▁▁ ▁▁▁ ▂
██████▇███▅▁▁▁▁▅▆▅▆▃▄▁▁▁▁▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▁▁▁▁▁▁▁▁▁▇███ █
5.02 μs Histogram: log(frequency) by time 6.5 μs <
Memory estimate: 0 bytes, allocs estimate: 0.
Di_banded:
BenchmarkTools.Trial: 10000 samples with 5 evaluations per sample.
Range (min … max): 6.891 μs … 13.974 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 6.949 μs ┊ GC (median): 0.00%
Time (mean ± σ): 7.024 μs ± 401.688 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▃█▇▁ ▁▁ ▁▁ ▂
████▄▅██▄▁▁▁▁▁▃▅▆▅▄▅▃▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄████▇ █
6.89 μ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): 117.730 μs … 199.134 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 148.288 μs ┊ GC (median): 0.00%
Time (mean ± σ): 150.016 μs ± 5.911 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▁▁▁▂▂▄▄▆▇██▇▅▄▃▂▂▁▁▂▃▄▄▅▅▃▂▁ ▃
▇▆▁▁▄▃▁▁▁▄▁▃▁▃▄▄▃▃▁▁▁▁▁▁▃▇███████████████████████████████▇█▇█ █
118 μs Histogram: log(frequency) by time 169 μ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.12.0
[9f78cca6] SummationByPartsOperators v0.5.101 `~/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.656 ns … 79.623 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 48.056 ns ┊ GC (median): 0.00%
Time (mean ± σ): 48.579 ns ± 2.175 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▃█▇▅▃▃▁
▂▂▃▅███████▇▅▄▃▃▂▂▃▂▂▂▂▂▂▂▁▂▁▁▁▁▁▁▂▁▂▂▂▂▁▁▁▁▁▂▂▂▂▂▃▂▂▂▂▂▂▂▂ ▃
46.7 ns Histogram: frequency by time 57.4 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 224 evaluations per sample.
Range (min … max): 326.951 ns … 595.804 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 390.326 ns ┊ GC (median): 0.00%
Time (mean ± σ): 391.388 ns ± 22.592 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▂▃▃▂▅▆▆█▇▇▆▅▃▁▁
▁▂▂▂▂▁▁▁▁▂▃▃▄▅▄▄▄▄▄▅▅████████████████▇▆▆▅▅▅▅▄▄▃▃▃▃▃▃▂▂▂▂▂▂▁▂▁ ▄
327 ns Histogram: frequency by time 455 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_full
BenchmarkTools.Trial: 10000 samples with 10 evaluations per sample.
Range (min … max): 1.100 μs … 2.751 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 1.112 μs ┊ GC (median): 0.00%
Time (mean ± σ): 1.124 μs ± 99.502 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
█▇▁ ▁
███▅▄▆▇▃▃▁▃▁▁▁▁▁▄▁▃▁▁▁▁▁▁▁▁▁▁▃▃▁▃▁▁▁▃▁▁▁▃▁▃▃▁▄▃▃▁▄▁▁▁▁▁▁▃▆ █
1.1 μs Histogram: log(frequency) by time 1.95 μ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.278 μs … 3.937 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 1.285 μs ┊ GC (median): 0.00%
Time (mean ± σ): 1.302 μs ± 107.286 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.471 μs … 6.490 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 2.574 μs ┊ GC (median): 0.00%
Time (mean ± σ): 2.604 μs ± 174.986 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▁▇█▅▂
▂▂▄█████▇▄▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▁▁▁▂▁▂▁▂▂▁▂▁▂▁▁▁▁▁▁▁▁▁▂▁▂▂▂▂▂▂▂ ▃
2.47 μs Histogram: frequency by time 3.52 μ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 … 12.233 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 6.290 μs ┊ GC (median): 0.00%
Time (mean ± σ): 6.358 μs ± 346.597 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 501 evaluations per sample.
Range (min … max): 210.295 ns … 297.443 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 213.735 ns ┊ GC (median): 0.00%
Time (mean ± σ): 215.676 ns ± 5.798 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▅██▆▁
▂▂▂▃▄██████▇▆▄▃▃▂▂▂▂▂▂▂▂▂▁▂▂▁▁▁▂▁▁▁▁▁▁▁▁▁▂▁▂▂▂▂▃▄▄▃▃▃▂▂▂▂▂▂▂▂ ▃
210 ns Histogram: frequency by time 234 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 187 evaluations per sample.
Range (min … max): 549.856 ns … 11.819 μs ┊ GC (min … max): 0.00% … 93.17%
Time (median): 568.556 ns ┊ GC (median): 0.00%
Time (mean ± σ): 576.331 ns ± 115.243 ns ┊ GC (mean ± σ): 0.19% ± 0.93%
▃▆██▇▄▂
▁▁▁▂▂▄▅█████████▆▅▄▃▃▃▃▂▂▂▁▁▁▁▁▁▁▁▁▁▂▂▂▃▃▃▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▃
550 ns Histogram: frequency by time 642 ns <
Memory estimate: 32 bytes, allocs estimate: 1.
D_full
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
Range (min … max): 13.936 μs … 60.043 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 14.016 μs ┊ GC (median): 0.00%
Time (mean ± σ): 14.239 μs ± 1.374 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
█
█▄▂▂▂▂▂▂▁▂▂▂▂▁▂▂▂▂▂▂▂▂▂▂▂▁▁▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▁▂▂▂▂▂▂▁▂▂ ▂
13.9 μs Histogram: frequency by time 22.2 μs <
Memory estimate: 0 bytes, allocs estimate: 0.Next, we still use an array of structures (AoS), but copy the data into a plain Array instead of using the reinterpreted versions. There is no significant difference to the previous version in this case.
println("Array of Structures")u_aos = Array(u_aos_r); du_aos = similar(u_aos)@show D_SBP * u_aos ≈ D_sparse * u_aos ≈ D_full * u_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 506 evaluations per sample.
Range (min … max): 212.711 ns … 398.393 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 215.879 ns ┊ GC (median): 0.00%
Time (mean ± σ): 218.055 ns ± 6.711 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▁▄▇███▆▅▄▃▂ ▂▃▄▄▃▂▁ ▃
▄▅█████████████▇▇███▇▇▅▄▁▁▃▆▆▁▄▁▃▁▁▃▄▃▃▃▅▃▆▇█████████▇▇▇▇▇▅▆▅ █
213 ns Histogram: log(frequency) by time 237 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 141 evaluations per sample.
Range (min … max): 674.454 ns … 1.043 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 732.504 ns ┊ GC (median): 0.00%
Time (mean ± σ): 738.607 ns ± 23.175 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▁▃▇█▇▆▅▂
▂▁▁▂▁▂▁▁▂▁▂▂▁▂▂▂▂▃▄▅▇████████▇▅▅▄▃▃▂▂▂▂▂▂▂▂▂▂▂▂▃▃▄▄▄▃▃▃▃▃▂▂▂ ▃
674 ns Histogram: frequency by time 810 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_full
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
Range (min … max): 13.895 μs … 42.580 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 13.966 μs ┊ GC (median): 0.00%
Time (mean ± σ): 14.350 μs ± 1.793 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
█▄ ▁
██▄▅▄█▆▁▄▄▃▁▁▁▄▅▄▄▄▄▄▄▄▄▄▅▃▄▄▄▆▅▃▅▅▅▅▅▄▃▄▇█▇▆▆▆▄▆▇▆▃▄▅▇▇▇▆▆ █
13.9 μs Histogram: log(frequency) by time 22.7 μ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 487 evaluations per sample.
Range (min … max): 222.945 ns … 312.002 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 226.236 ns ┊ GC (median): 0.00%
Time (mean ± σ): 228.339 ns ± 6.569 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▂▆█▆▅▃
▂▃▅▇███████▇▅▅▄▃▃▂▂▂▂▂▂▂▂▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▂▂▂▂▃▃▄▄▃▃▃▃▂▃▂▂▂▂▂▂▂ ▃
223 ns Histogram: frequency by time 248 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
D_sparse
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
Range (min … max): 48.901 μs … 108.503 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 49.092 μs ┊ GC (median): 0.00%
Time (mean ± σ): 49.644 μs ± 2.641 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▆█▅▁ ▂ ▂▂ ▁
████▅▇█▇▄▁▁▁▃▅▆▆▆▄▅▃▁▄▁▃▃▁▁▁▁▁▁▁▁▃▁▁▁▁▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▇███▇▆▆ █
48.9 μs Histogram: log(frequency) by time 57.8 μs <
Memory estimate: 0 bytes, allocs estimate: 0.
D_full
BenchmarkTools.Trial: 10000 samples with 6 evaluations per sample.
Range (min … max): 5.570 μs … 11.797 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 5.612 μs ┊ GC (median): 0.00%
Time (mean ± σ): 5.677 μs ± 331.523 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▅█▅ ▁ ▁▁ ▁
███▄▃██▅▃▁▁▁▄▅▅▄▅▁▄▁▁▁▁▁▁▃▃▁▁▄▁▄▁▁▃▃▁▃▃▃▁▃▁▁▃▃▃▁▁▁▇███▇▅▅▄▄ █
5.57 μs Histogram: log(frequency) by time 7.17 μ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.22
[09ab397b] StructArrays v0.7.3
[9f78cca6] SummationByPartsOperators v0.5.101 `~/work/SummationByPartsOperators.jl/SummationByPartsOperators.jl`