DimensionError for determinant of generic, triangular, unitful matrices

[lxvm]

using Unitful
using LinearAlgebra
det([1 0; 0 1]*u"m")

leads to

ERROR: DimensionError: m and 1 m² are not dimensionally compatible.
Stacktrace:
 [1] convert(::Type{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}, x::Quantity{Int64, 𝐋², Unitful.FreeUnits{(m²,), 𝐋², nothing}})
   @ Unitful ~/.julia/packages/Unitful/sxiHA/src/conversion.jl:133
 [2] det(A::Matrix{Quantity{Int64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.11.6+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/generic.jl:1649
 [3] top-level scope
   @ REPL[24]:1

I suspect that it is because this type handling should be using oneunit instead of one

LinearAlgebra.jl/src/generic.jl

Line 1848 in 7e647cf

S = promote_type(T, typeof((one(T)*zero(T) + zero(T))/one(T)))

Although it isn't entirely clear to me why conversion is necessary, here is a potential fix that keeps the code similar but fixes this edge case

function mydet(A::AbstractMatrix{T}) where {T}
    if istriu(A) || istril(A)
        S = promote_type(T, typeof((oneunit(T)*zero(T) + zero(T)*zero(T))/one(T)))
        return convert(S, det(UpperTriangular(A)))
    end
    return det(lu(A; check = false))
end
mydet([1 0; 0 1]*u"m") # returns 1.0 m²

[stevengj]

Using oneunit doesn't seem enough here, because the correct units are oneunit(T)^m.

The reason for the convert is to ensure a type-stable result: det(UpperTriangular(A))) is division-free, and returns an integer result for integer matrices, whereas det(lu(A)) uses divisions and hence converts to a floating-point result. We want to return a floating-point result in both cases for an integer matrix even if it happens to be triangular.

[stevengj]

However, the current algorithm seems wrong for another reason: because it only converts to floating-point after computing the determinant in the triangular case, it is susceptible to spurious overflow. For example, it gives the wrong result in this case:

julia> A = diagm([typemax(Int), typemax(Int)])
2×2 Matrix{Int64}:
 9223372036854775807                    0
                   0  9223372036854775807

julia> det(A) # incorrect: spurious overflow
1.0

julia> det(float(A)) # correct
8.507059173023462e37

The following implementation should fix both the units issue and the overflow issue:

function det(A::AbstractMatrix{T}) where {T}
    if istriu(A) || istril(A)
        S = promote_type(T, typeof((one(T)*zero(T) + zero(T))/one(T)))
        return prod(S, @view A[diagind(A)])
    end
    return det(lu(A; check = false))
end

That is, you convert to S (which has the same units as T) before doing the product.

Feel like doing a PR with both the units test and the overflow test?

[lxvm]

Sure!