# Correct way to benchmark StaticVectors

**URL:** <https://discourse.julialang.org/t/correct-way-to-benchmark-staticvectors/119357>\
**Category:** General Usage\
**Tags:** benchmarktools, staticarrays\
**Created:** [September 13, 2024, 12:26am UTC](https://discourse.julialang.org/t/correct-way-to-benchmark-staticvectors/119357 "2024-09-13T00:26:24Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![bthorne93](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bthorne93/32/19570_2.png) [@bthorne93](https://discourse.julialang.org/u/bthorne93)\
**Post date:** [September 13, 2024, 12:26am UTC](https://discourse.julialang.org/t/correct-way-to-benchmark-staticvectors/119357/1 "2024-09-13T00:26:24Z")

</div>

I am writing some code to do geometric intersection tests using StaticArrays.jl. I would like to keep some benchmarks as I am developing and would like to know if I am benchmarking static arrays correctly.

The documentation for BenchmarkTools mentions that operations that can be simplified at compile time will not be measured properly, and that you therefore have to reference and dereference the interpolated variable. I think I am doing that properly below, however the measured output is still very fast for this operation (around 6 ns), so I am suspicious that I am not measuring the right thing.

Would the below benchmark correctly measure the execution speed of the BVH.intersection function, and would anyone be able to explain / point to a source on how this interpolation and referencing / dereferencing works?

```julia
using BenchmarkTools
using Revise, StaticArrays, LinearAlgebra

function benchmark_intersection()
    for T in [Float64, Float32]

        v₁ = SVector{3, T}(1.0, 0.0, 1.0)
        v₂ = SVector{3, T}(-1.0, 1.0, 1.0)
        v₃ = SVector{3, T}(-1.0, -1.0, 1.0)

        O = SVector{3, T}(0.0, 0.0, 0.0)
        direction = SVector{3, T}(0.0, 0.0, 1.0)

        A = [-direction (v₂-v₁) (v₃-v₁)]
        b = O - v₁

        ray = Ray(O, direction)
        triangle = Triangle(v₁, v₂, v₃)
          
        println("Benchmarking $T")
        display(@benchmark intersection($(Ref(ray))[], $(Ref(triangle))[]))
    end
end

struct Ray{N, T}
    origin::SVector{N, T}
    direction::SVector{N, T}
end

struct Triangle{N, T}
    v₁::SVector{N, T}
    v₂::SVector{N, T}
    v₃::SVector{N, T}
end

function intersection(ray::Ray{N, T}, triangle::Triangle{N, T}) :: IntersectionResult{N, T} where {N, T}
    (; origin, direction) = ray
    (; v₁, v₂, v₃) = triangle

    A = [-direction (v₂-v₁) (v₃-v₁)]
    b = origin - v₁

    if det(A) ≈ 0
        return IntersectionResult(false, nothing)
    end

    (t, _, _) = A \ b
    point = origin + direction * t

    return IntersectionResult(true, point)
end

```

---

<div class="post-metadata">

**Author:** ![matthias314](https://avatars.discourse-cdn.com/v4/letter/m/a88e4f/32.png) [@matthias314](https://discourse.julialang.org/u/matthias314)\
**Post date:** [September 13, 2024, 1:18am UTC](https://discourse.julialang.org/t/correct-way-to-benchmark-staticvectors/119357/2 "2024-09-13T01:18:15Z")

</div>

What is `BVH`?

---

<div class="post-metadata">

**Author:** ![bthorne93](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bthorne93/32/19570_2.png) [@bthorne93](https://discourse.julialang.org/u/bthorne93)\
**Post date:** [September 13, 2024, 3:35am UTC](https://discourse.julialang.org/t/correct-way-to-benchmark-staticvectors/119357/3 "2024-09-13T03:35:50Z")

</div>

The module I am writing - sorry forgot to add that code. Edited above to include.

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [September 13, 2024, 6:59am UTC](https://discourse.julialang.org/t/correct-way-to-benchmark-staticvectors/119357/4 "2024-09-13T06:59:49Z")

</div>

Benchmarking looks fine to me. Generally, timings `<1ns` are suspicious and indicate a complete constant-folding of your function. But here you are using StaticArrays.jl with vectors of size 3 so each function are just a handful of FLOPS and modern CPUs are quite fast 🙂

However this

> [@bthorne93](#):
>
> `det(A) ≈ 0`

is nonsense because it will always evaluate to `false` as just being discussed [in another thread](https://discourse.julialang.org/t/history-of-isapprox-in-languages/119347).  
You’ll need to think about an absolute scale, e.g. given by the largest element of `A` and then set `abstol` accordingly. Though I am not sure what “accordingly” should mean concretely when considering `det`. @stevengj probably has an idea what an appropriate scale should be. I would propose `abstol=sqrt(eps(eltype(A)))*maximum(A)`

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [September 13, 2024, 5:01pm UTC](https://discourse.julialang.org/t/correct-way-to-benchmark-staticvectors/119357/5 "2024-09-13T17:01:22Z")

</div>

> [@abraemer](#):
>
> I would propose `abstol=sqrt(eps(eltype(A)))*maximum(A)`

That won’t work, because it has the wrong units. (The units of the determinant of an m \times m matrix are the matrix elements to the m'th power.) An absolute tolerance has to have the same units as the quantity you are comparing it to.

Also `maximum(A)` returns the largest element _algebraically_, whereas I think you want the largest _magnitude_.

The “right” way to determine whether a matrix is close to singular is to compute its condition number, which usually is computed via the SVD, though you could in principle also use `norm(A)*norm(inv(A))` (except that `inv` might throw). A nearly as good way is to use rank-revealing QR.

But for a 3 \times 3 matrix, checking that `det(A) ≤ rtol * maximum(abs ∘ float, A)^size(A,1)` should be fairly reasonable, I guess (without thinking about it too hard), for some decent tolerance, e.g. `rtol = sqrt(eps(float(eltype(A))))`.

(Determinants can also spuriously overflow/underflow because they can grow exponentially with matrix size. But for 3x3 matrices with reasonable-magnitude entries it shouldn’t be an issue.)