Unitful.jl attaches a physical unit to a number as part of its type: 2.0u"kg" is a Quantity{Float64, ...}, not a Float64 with a tag bolted on. That has a consequence for the compilation model from intro_julia.jl: multiple dispatch specialises on argument type, and the unit is part of the type, so a function called with unitful arguments compiles its own specialisation — once, the first time, same as any other type combination. No dimensional analysis happens at runtime.

The headline: after that one-off specialisation, a unitful call and its Float64 sibling lower, type-infer, and run near-identically. Then the catch — multiple dispatch is not automatic across every combination of types. Composing two abstractions (a sparse matrix, wrapped to say “symmetric”) can silently fall through to a generic, slower method the moment nobody has written the specific one for that combination — and the only fix is to write it.

using Unitful, SparseArrays, LinearAlgebra, InteractiveUtils, BenchmarkTools, Printf

kinetic_energy(m, v) = 0.5 * m * v^2

kinetic_energy(2.0, 3.0)                    # compile the Float64 specialisation
kinetic_energy(2.0u"kg", 3.0u"m/s")          # compile the Quantity specialisation

9.0 kg m^2 s^-2

Same function, two calls. m and v carry no unit annotation in the source — the units live entirely in the types of the values passed in.

println(kinetic_energy(2.0, 3.0))
println(kinetic_energy(2.0u"kg", 3.0u"m/s"))
println(typeof(kinetic_energy(2.0u"kg", 3.0u"m/s")))

9.0
9.0 kg m^2 s^-2
Quantity{Float64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(kg, m^2, s^-2), 𝐋^2 𝐌 𝐓^-2, nothing}}

9.0 kg m^2 s^-2 — correct energy dimensions, computed from a function that never mentions units. kg * (m/s)^2 combines to J-equivalent dimensions automatically because * and ^ on Quantity are ordinary multiple-dispatch methods that propagate units through the type.

Same lowered IR

If the unit only affects type, it shouldn’t affect lowering — lowering runs before any type is known.

lowered_f64 = @code_lowered kinetic_energy(2.0, 3.0)
lowered_u = @code_lowered kinetic_energy(2.0u"kg", 3.0u"m/s")
display(lowered_f64)
println("identical lowered IR: ", lowered_f64.code == lowered_u.code)

CodeInfo(
1 ─ %1 = Main.:*
│   %2 = Main.:^
│   %3 =   builtin Core.apply_type(Base.Val, 2)
│   %4 =   dynamic (%3)()
│   %5 =   dynamic Base.literal_pow(%2, v, %4)
│   %6 =   dynamic (%1)(0.5, m, %5)
└──      return %6
)

identical lowered IR: true

Not just similar — the exact same CodeInfo. kinetic_energy was parsed and lowered once; the Quantity call reuses that IR verbatim and only then forks into its own type-inferred, compiled specialisation. Contrast intro_julia.jl’s mul_naive vs mul_fused, where the dotted syntax produced genuinely different lowered IR — here nothing about the code changed, only the types flowing through it.

Nearly identical typed IR

Type inference is where the two calls finally diverge.

display(@code_typed kinetic_energy(2.0, 3.0))

CodeInfo(
1 ─ %1 = intrinsic Base.mul_float(v, v)::Float64
│   %2 = intrinsic Base.mul_float(0.5, m)::Float64
│   %3 = intrinsic Base.mul_float(%2, %1)::Float64
└──      return %3
) => Float64

display(@code_typed kinetic_energy(2.0u"kg", 3.0u"m/s"))

CodeInfo(
1 ─ %1 =   builtin Base.getfield(v, :val)::Float64
│   %2 = intrinsic Base.mul_float(%1, %1)::Float64
│   %3 =   builtin Base.getfield(m, :val)::Float64
│   %4 = intrinsic Base.mul_float(0.5, %3)::Float64
│   %5 = intrinsic Base.mul_float(%4, %2)::Float64
│   %6 = %new(Quantity{Float64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(kg, m^2, s^-2), 𝐋^2 𝐌 𝐓^-2, nothing}}, %5)::Quantity{Float64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(kg, m^2, s^-2), 𝐋^2 𝐌 𝐓^-2, nothing}}
└──      return %6
) => Quantity{Float64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(kg, m^2, s^-2), 𝐋^2 𝐌 𝐓^-2, nothing}}

Both bodies are the same two mul_float intrinsics. The Quantity version adds two getfields to unwrap .val from m and v and one %new to rewrap the product back into a Quantity — and all three are compiled away to nothing more than moving bits into and out of an unboxed struct. No branch, no dictionary lookup, no runtime unit check: kg * (m/s)^2 was resolved to a concrete result type during type inference, so the “unit arithmetic” is over before the function body starts running.

kinetic_energy(2.0, 3.0)
kinetic_energy(2.0u"kg", 3.0u"m/s")
t_f64 = @benchmark kinetic_energy($2.0, $3.0)
t_u = @benchmark kinetic_energy($(2.0u"kg"), $(3.0u"m/s"))
@printf("Float64:  %s\n", BenchmarkTools.prettytime(median(t_f64).time))
@printf("Unitful:  %s\n", BenchmarkTools.prettytime(median(t_u).time))

Float64:  1.443 ns
Unitful:  1.443 ns

Sub-2ns, no daylight between them within noise. The type/dispatch/JIT machinery paid for the unit exactly once, at first call, and charges nothing per call after that.

Where the automagic ends

Build a small sparse matrix, attach a unit, and strip it back off with ustrip. (broadcast — the idiomatic Unitful spelling):

S = spzeros(5, 5)
S[1, 1] = 1.0
S[2, 2] = 2.0
S[1, 2] = S[2, 1] = 0.5
Su = S .* 1.0u"MHz"

5×5 SparseMatrixCSC{Quantity{Float64, 𝐓^-1, Unitful.FreeUnits{(MHz,), 𝐓^-1, nothing}}, Int64} with 4 stored entries:
 1.0 MHz  0.5 MHz     ⋅        ⋅        ⋅   
 0.5 MHz  2.0 MHz     ⋅        ⋅        ⋅   
    ⋅        ⋅        ⋅        ⋅        ⋅   
    ⋅        ⋅        ⋅        ⋅        ⋅   
    ⋅        ⋅        ⋅        ⋅        ⋅   

ustrip.(u"MHz", Su)

5×5 SparseMatrixCSC{Float64, Int64} with 4 stored entries:
 1.0  0.5   ⋅    ⋅    ⋅ 
 0.5  2.0   ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅ 

println(typeof(ustrip.(u"MHz", Su)))

SparseMatrixCSC{Float64, Int64}

SparseMatrixCSC{Float64} — sparsity preserved, as expected: SparseArrays taught broadcast how to only touch stored entries. Now wrap the same matrix to mark it symmetric — a completely ordinary thing to do with a Hamiltonian or a covariance matrix — and strip the unit again:

Hsym = Symmetric(Su)

5×5 Symmetric{Quantity{Float64, 𝐓^-1, Unitful.FreeUnits{(MHz,), 𝐓^-1, nothing}}, SparseMatrixCSC{Quantity{Float64, 𝐓^-1, Unitful.FreeUnits{(MHz,), 𝐓^-1, nothing}}, Int64}}:
 1.0 MHz  0.5 MHz     ⋅        ⋅        ⋅   
 0.5 MHz  2.0 MHz     ⋅        ⋅        ⋅   
    ⋅        ⋅        ⋅        ⋅        ⋅   
    ⋅        ⋅        ⋅        ⋅        ⋅   
    ⋅        ⋅        ⋅        ⋅        ⋅   

ustrip.(u"MHz", Hsym)

5×5 Matrix{Float64}:
 1.0  0.5  0.0  0.0  0.0
 0.5  2.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

println(typeof(ustrip.(u"MHz", Hsym)))

Matrix{Float64}

Matrix, not SparseMatrixCSC. Silently — no error, no warning, just a container swap. Nothing about Symmetric or SparseMatrixCSC is individually broken; broadcast just doesn’t know how to combine them. Broadcast.combine_styles shows why: Julia picks a BroadcastStyle for the outermost wrapper — SparseArrays registered a style that keeps . operations sparse for a bare SparseMatrixCSC, but nobody registered one for “Symmetric around a sparse parent”, so it falls back to the same dense style a Symmetric{Float64,Matrix} would get:

println("bare sparse  broadcast style: ", Broadcast.combine_styles(Su))
println("Symmetric(.) broadcast style: ", Broadcast.combine_styles(Hsym))

bare sparse  broadcast style: SparseArrays.HigherOrderFns.SparseMatStyle()
Symmetric(.) broadcast style: Base.Broadcast.DefaultArrayStyle{2}()

The consequence isn’t cosmetic. A dense copy feeds a dense mul!, which is O(N²) instead of O(nnz):

N = 4000
Sbig = spdiagm(0 => fill(2.0, N), 1 => fill(0.3, N - 1), -1 => fill(0.3, N - 1))
x = rand(ComplexF64, N)

H_dense = ustrip.(u"MHz", Symmetric(Sbig .* 1.0u"MHz"))  # densified by the broadcast above
H_sparse = Sbig                                            # what it should have stayed as

mul!(similar(x), H_dense, x)
mul!(similar(x), H_sparse, x)
t_dense = @benchmark mul!(y, $H_dense, $x) setup = (y = similar($x))
t_sparse = @benchmark mul!(y, $H_sparse, $x) setup = (y = similar($x))
@printf("dense:  %s\n", BenchmarkTools.prettytime(median(t_dense).time))
@printf("sparse: %s\n", BenchmarkTools.prettytime(median(t_sparse).time))
@printf("speedup: %.0fx\n", median(t_dense).time / median(t_sparse).time)

dense:  3.614 ms
sparse: 7.946 μs
speedup: 455x

N = 4000 with a 3-wide band is already several hundred times slower; a matrix large enough to need sparsity in the first place wouldn’t fit in memory as a dense copy at all. This is the same scale argument as mul_native! in intro_julia.jl, just for algorithmic complexity instead of temporaries — except here the wrong dispatch doesn’t just cost a constant factor, it changes the exponent, and nothing in the program said so.

The only fix: write the missing method

There’s no way to make broadcast “just know” how to handle every wrapper combination in advance — multiple dispatch gives you a method table, not an inference engine. The fix is to add the specific method for the specific combination that matters here, bypassing broadcast entirely by overloading the non-broadcast, whole-array form of ustrip (called without the dot):

function Unitful.ustrip(u::Unitful.Units, H::Symmetric{<:Unitful.Quantity,<:SparseMatrixCSC})
    P = parent(H)
    bare = SparseMatrixCSC(P.m, P.n, P.colptr, P.rowval, ustrip.(u, P.nzval))
    return Symmetric(bare, Symbol(H.uplo))
end

fixed = ustrip(u"MHz", Hsym)                  # no dot: calls the new method directly
println(typeof(fixed))
println(fixed == ustrip.(u"MHz", Matrix(Hsym)))  # same values, right storage

Symmetric{Float64, SparseMatrixCSC{Float64, Int64}}
true

One method, scoped exactly to the type combination that broke, restores both the sparsity and the Symmetric wrapper — colptr/rowval are shared with the original, not copied; only nzval is touched. Nothing upstream needs to change: every caller that already wrote ustrip(u, H) picks this up automatically, because that’s the one spelling that was never routed through broadcast in the first place. That’s the trade multiple dispatch makes: it never guesses a fast path for a combination you haven’t told it about, but the fix, once you’ve found where it broke, is always this local — one more method, nothing else touched.