Do not use zero(...) in generic_lu! algorithm to support Unitful

[haampie]

This is a minimal effort attempt to work around the bigger problem that factorizations only define a single number type.

Suppose eltype(A) == typeof(1.0u"kg"), then by design of the algorithm the LU factorization should end up with A = L * U where eltype(L) = Float64 and eltype(U) == eltype(A) == typeof(1.0u"kg").

In julia the algorithm assumes it can work in-place, so we have to promote/union number types of L and U:

julia> promote_type(Float64, typeof(1.0u"kg"))
Quantity{Float64, D, U} where U where D

and in this case zero(promote_type(Float64, typeof(1.0u"kg"))) will error. This PR tries to avoid using zero then.

A better solution IMHO is to allocate L and U separately (or at the very least parametrize the eltype of L and U separately in the factorization) to avoid having to work with non-concrete number types. It's nice to get a unitful LU-factorization, but pretty much any linalg operation \, *, +, etc will run oneunit(eltype(factorization)) and throw:

julia> F = lu(rand(typeof(1.0u"kg"), 3, 3))
LU{Quantity{Float64, D, U} where U where D, Matrix{Quantity{Float64, D, U} where U where D}}
L factor:
3×3 Matrix{Quantity{Float64, D, U} where U where D}:
      1.0     0.0 kg  0.0 kg
 0.268046        1.0  0.0 kg
 0.512366  -0.566135     1.0
U factor:
3×3 Matrix{Quantity{Float64, D, U} where U where D}:
 0.654136 kg  0.762665 kg  0.637124 kg
         0.0  0.640381 kg  0.682195 kg
         0.0          0.0  0.102111 kg

julia> F \ rand(3)
ERROR: MethodError: no method matching (Quantity{Float64, D, U} where U where D)(::Float64)
Closest candidates are:

[Keno]

Could use a test. We have Furlong in the test code for this purpose.

[Keno]

This hit during CI downtime, but otherwise LGTM. Should have a test though, the addition of which will also re-trigger CI.

[haampie]

I think we can't use Furlong here, since * is not a closed operation. It would mean the eltype would depend on the indices, meaning that the matrix would have some abstract type, in which case simple stuff like one(T) won't work.

Edit: I take that back, Furlongs work fine:

julia> lu(Furlongs.Furlong.(rand(6, 6))).factors
6×6 Matrix{Main.Furlongs.Furlong{p, Float64} where p}:
 Furlong{1, Float64}(0.789404)  Furlong{1, Float64}(0.0483697)   Furlong{1, Float64}(0.156466)  …   Furlong{1, Float64}(0.146989)   Furlong{1, Float64}(0.305877)
 Furlong{0, Float64}(0.107174)   Furlong{1, Float64}(0.567589)   Furlong{1, Float64}(0.479531)      Furlong{1, Float64}(0.710281)   Furlong{1, Float64}(0.892259)
 Furlong{0, Float64}(0.686722)   Furlong{0, Float64}(0.813328)   Furlong{1, Float64}(0.217378)      Furlong{1, Float64}(-0.40261)  Furlong{1, Float64}(-0.816591)
 Furlong{0, Float64}(0.714551)   Furlong{0, Float64}(0.964134)  Furlong{0, Float64}(-0.435512)     Furlong{1, Float64}(-0.518144)  Furlong{1, Float64}(-0.824331)
 Furlong{0, Float64}(0.587981)   Furlong{0, Float64}(0.342699)  Furlong{0, Float64}(-0.216329)     Furlong{1, Float64}(-0.662855)  Furlong{1, Float64}(-0.629475)
 Furlong{0, Float64}(0.799651)   Furlong{0, Float64}(0.217718)  Furlong{0, Float64}(-0.147971)  …  Furlong{0, Float64}(-0.130078)  Furlong{1, Float64}(-0.826502)

but you end up with a matrix of some abstract type.

[Keno]

I'm confused what this change does then. Isn't that the same issue with Unitful also?

[haampie]

Sorry for the confusion, it's exactly the same as Unitful yes. There was an error displaying the LU-factorization of Furlongs in the REPL where one(T) and zero(T) were used, which made me think the computation failed, but in fact it works fine. Somehow this printing error doesn't occur for Unitful.

[dkarrasch]

The only failure is a timeout, but the LU tests passed. So this is ready to go, right?