Make *Triangular handle units

[dkarrasch]

This is an attempt to make triangular linear algebra handle unitful objects well, like (most of) UpperHessenberg and Bidiagonal. The first commit rewrites the generic functions in terms of 3-arg mul!, ldiv! and the internal _rdiv! and doesn't change any tests to convince ourselves that this doesn't change behavior. 2-arg versions for BLAS-able argument types remain fully functional. Some of the loops had to be reorder to keep the variables consistent with units.

Closes JuliaLang/LinearAlgebra.jl#875.

EDIT (copied from below):

The features of the PR are:

• multiplication and solving with triangular matrices involving units (the triangular matrix, the other factor, or both) works
• inversion of unitful triangular matrices works
• matrix square root of unitful UnitTriangular matrices works
• right in-place multiplication by a triangular matrix does no longer fall back to generic code (see mul!(::Matrix{T}, ::Matrix{T}, ::AbstractTriangular{T, <:Matrix{T}}) where T<:BlasFloat does not use BLAS LinearAlgebra.jl#875)
• 3-arg triangular multiplication/division of triangular matrices is directed towards BLAS if possible
• reduces mul! methods down to 55 (we have 138 in v1.9)
• there is a "_ustrip" mechanism that packages can hook into to have generic substitution codes run without overhead by division by some oneunit

[dkarrasch]

@nanosoldier runbenchmarks("linalg", vs = ":master")

[dkarrasch]

@nanosoldier runtests(ALL, vs = ":master")

[nanosoldier]

Your benchmark job has completed - possible performance regressions were detected. A full report can be found here.

[nanosoldier]

Your package evaluation job has completed - possible new issues were detected. A full report can be found here.

[KlausC]

Isn't there a more efficient way to move the correct physical dimensions to bi?

When the element types of A and b are unitful integers, the result could also be unitful integers, without any floating point operations performed (no division by 1).

[dkarrasch]

I don't know of any. convert and alike fail because you can't convert, for instance, square meters to meters or a unitless number. Also, I'm not sure if there is a generic way of identifying whether a unitful quantity is integer-based. For unitless numbers you can see that I made an effort to preserve the integer return type.

[dkarrasch]

I worked around the issue. First, the target c is initialized with an integer unitful eltype, and second I replaced the division by a multiplication. I'll add a test for that also.

[dkarrasch]

Actually, no. For unitful quantities, we need to "perform" the division to get the right unit: assume both A and b have integer units of kilometers, then the result is "unit-free". In terms of types, this means (Int*km) \ (Int*km) = Float64, otherwise the kilometers don't cancel out.

[dkarrasch]

@nanosoldier runtests(ALL, vs = ":master")

[dkarrasch]

what do you need to make progress here?

Feedback. 😜

This is a pretty big refactoring, and it increases CI runtime in the triangular test file (due to compilation of more methods) and in the bidiag test file due to quite some more tests involving units. I'm wondering if this is all worth it or if we have other options to handle unitful linear algebra (and whether this is all of sufficient interest).

[dkarrasch]

@nanosoldier runtests()

[nanosoldier]

Your package evaluation job has completed - possible new issues were detected.
A full report can be found here.

[dkarrasch]

@nanosoldier runtests(["ReusableFunctions", "Pyehtim", "CommutativeRings", "HomotopyOpt", "GradientRobustMultiPhysics", "HTTP", "JuMP", "LazyBandedMatrices", "SignalTablesInterface_CairoMakie", "Pickle", "SignalTablesInterface_WGLMakie", "CompressiveLearning", "SummationByPartsOperators", "UnfoldMakie", "Decapodes", "Estapir", "LiteQTL", "OptimizationMetaheuristics", "Kronecker", "ClimaCore"])

[nanosoldier]

Your package evaluation job has completed - possible new issues were detected.
A full report can be found here.

[dkarrasch]

pkgeval says we're good. So let's go?

[dkarrasch]

@nanosoldier runbenchmarks()

[dkarrasch]

@nanosoldier runbenchmarks("linalg")

[nanosoldier]

Your benchmark job has completed - successfully executed benchmarks. A full report can be found here.

[dkarrasch]

@nanosoldier runbenchmarks("linalg")

[nanosoldier]

Your benchmark job has completed - successfully executed benchmarks. A full report can be found here.

[dkarrasch]

I must have missed some update for nanosoldier. We used to have green checkmarks and red crosses in the benchmarks output. How do I read the output now? Does the absence of these signs mean "everything within expectation"? @maleadt @KristofferC

[maleadt]

That's because you didn't do a comparison benchmark. Only the PkgEval Nanosoldier instance supports leaving off the vs keyword.
@vtjnash This has happened a couple of times now. Can you rebase your Nanosoldier instance on top of the latest Nanosoldier.jl?

[dkarrasch]

I wasn't reading carefully, but just looked at the examples at https://github.com/JuliaCI/Nanosoldier.jl, and copy-pasted. Now I see way below that there are more details to be specified. Let's try again.

@nanosoldier runbenchmarks("linalg", vs = ":master")

[nanosoldier]

Your benchmark job has completed - possible performance regressions were detected. A full report can be found here.

[dkarrasch]

The only reported regression must be fake. The code path (mul! for 2x2 matrices) is not touched here, and I also can't reproduce locally.