Inverse Triangular methods missing

[btmit]

The inverse of a Triangular SMatrix does not preserve type. The resulting triangular is based on Matrix not SMatrix because the fallback method is hard-coded to Matrix. It seems like this should be a new method in StaticArrays.jl rather than something in new Base since it'll probably need to use inv() instead of ldiv!().

using LinearAlgebra, StaticArrays

a = rand(SMatrix{3,3})
ua = UpperTriangular(a) 
la = LowerTriangular(a)

typeof(inv(ua)) == typeof(ua)  # false UpperTriangular{Float64, Matrix{Float64}}
typeof(inv(la)) == typeof(la)  # false LowerTriangular{Float64, Matrix{Float64}}

[btmit]

In addition to unexpected typing there is a noteworthy performance penalty.

ua_raw = SMatrix{3,3}(ua)
@btime UpperTriangular(inv($ua_raw))  # 4.666 ns (0 allocations: 0 bytes)

@btime inv($ua)  # 162.800 ns (1 allocation: 128 bytes)

[RomeoV]

Here's an implementation

"""
    inv(U::UpperTriangular{T, <:SMatrix}) where T

Custom inverse for upper triangular static matrices using back-substitution.
Preserves SMatrix type instead of converting to Matrix.
"""
function LinearAlgebra.inv(U::UpperTriangular{T, <:SMatrix{N,N}}) where {T,N}
    A = parent(U)

    # Build columns as SVectors, then construct SMatrix from tuple
    cols = ntuple(N) do j
        # Standard basis vector for column j
        b = SVector{N}(i == j ? one(T) : zero(T) for i in 1:N)

        # Back-substitution for column j
        x = MVector{N,T}(undef)
        for i in N:-1:1
            s = b[i]
            for k in i+1:N
                s -= A[i,k] * x[k]
            end
            x[i] = s / A[i,i]
        end

        SVector{N}(x)
    end

    # Construct matrix from column vectors
    return hcat(cols...)
end

"""
    inv(L::LowerTriangular{T, <:SMatrix}) where T

Custom inverse for lower triangular static matrices using forward-substitution.
Preserves SMatrix type instead of converting to Matrix.
"""
function LinearAlgebra.inv(L::LowerTriangular{T, <:SMatrix{N,N}}) where {T,N}
    A = parent(L)

    # Build columns as SVectors, then construct SMatrix from tuple
    cols = ntuple(N) do j
        # Standard basis vector for column j
        b = SVector{N}(i == j ? one(T) : zero(T) for i in 1:N)

        # Forward-substitution for column j
        x = MVector{N,T}(undef)
        for i in 1:N
            s = b[i]
            for k in 1:i-1
                s -= A[i,k] * x[k]
            end
            x[i] = s / A[i,i]
        end

        SVector{N}(x)
    end

    # Construct matrix from column vectors
    return hcat(cols...)
end