[ANN] UnitfulTensors.jl: Efficient arrays with units

Unitful.jl sometimes works well for arrays with homogeneous units, but as soon as you have heterogeneous arrays, you are doomed.
Even with homogeneous arrays there can be slowdowns by orders of magnitude, see the matrix*matrix benchmark below.
\ doesn’t work with Unitful.jl.

using Unitful, LinearAlgebra, SparseArrays, BenchmarkTools

Nx = 10
N = Nx^2

Δ1d = Tridiagonal(fill(1., Nx - 1), fill(-2., Nx), fill(1., Nx - 1))
𝟙 = I(Nx)

Δ = kron(Δ1d, 𝟙) + kron(𝟙, Δ1d)
Δs = sparse(Δ)
ρ = randn(N)

uΔ = Δ * u"m^-2"
uΔs = Δs * u"m^-2"
uρ = ρ * u"C/m^3"
uΔhet = hcat(uΔ[:, 1:N÷2], uΔ[:, N÷2+1:end] * u"m^2")
uΔshet = hcat(uΔs[:, 1:N÷2], uΔs[:, N÷2+1:end] * u"m^2")
uρhet = [uρ[1:N÷2] * u"m^2"; uρ[N÷2+1:end]]

println("Matrix * vector:")
@btime $Δ * $ρ
@btime $uΔ * $uρ
@btime $Δs * $ρ
@btime $uΔs * $uρ
println("Matrix * matrix:")
@btime $Δ * $Δ
@btime $uΔ * $uΔ
@btime $Δs * $Δs
@btime $uΔs * $uΔs
println("Matrix * vector, heterogeneous units:")
@btime $Δ * $ρ
@btime $uΔhet * $uρhet
@btime $Δs * $ρ
@btime $uΔshet * $uρhet

# output

Matrix * vector:
  1.625 μs (1 allocation: 896 bytes)
  2.599 μs (1 allocation: 896 bytes)
  755.117 ns (1 allocation: 896 bytes)
  754.833 ns (1 allocation: 896 bytes)
Matrix * matrix:
  57.442 μs (2 allocations: 78.17 KiB)
  1.156 ms (6 allocations: 108.94 KiB)
  16.610 μs (6 allocations: 35.50 KiB)
  16.015 μs (6 allocations: 35.50 KiB)
Matrix * vector, heterogeneous units:
  1.622 μs (1 allocation: 896 bytes)
  655.666 μs (20301 allocations: 318.06 KiB)
  756.780 ns (1 allocation: 896 bytes)

MethodError: no method matching zero(::Type{Any})
1 Like