Make LU factorization work for more types

[blegat]

When determining the element type need, lufact does the following

julia/stdlib/LinearAlgebra/src/lu.jl

Line 180 in 1a870b4

S = typeof(zero(T)/one(T))

However LU factorization also multiplies elements of type T together and sum them together.
If the type T is not invariant under multiplication and/or addition, lufact does work, see JuliaAlgebra/MultivariatePolynomials.jl#69.
Note the strong similary between this PR and #18218 which did a similar change but to matrix multiplication.

[andreasnoack]

👍 Minor detail but maybe we should call it something like rationalop since it would probably be the right operation for all triangular-triangular factorizations.

[fredrikekre]

Better as something like

lutype(T) = (zero(T)*zero(T) + zero(T)*zero(T)) / (one(T)*one(T) + one(T)*one(T))

or is it OK to rely on promote_op?

[blegat]

we should call it something like rationalop since it would probably be the right operation for all triangular-triangular factorizations

Agreed, I have renamed it

is it OK to rely on promote_op ?

I have assumed that it was ok since matrix product relies on it

julia/stdlib/LinearAlgebra/src/matmul.jl

Lines 123 to 124 in 4ed4637

	function (*)(A::AbstractMatrix, B::AbstractMatrix)
	TS = promote_op(matprod, eltype(A), eltype(B))

but you think otherwise, I can change, let me know :)

[martinholters]

Use of promote_op should be avoided if possible. Ceterum censeo promote_op esse delendam.

[andreasnoack]

User of promote_op should be avoided if possible.

The price of using instances through zero(T) and one(T) here is that it might fail for types where these don't exist. That was the issue for matrix multiplication so it could also happen here although it probably is less likely.

@blegat After thinking a little more about this, I think we might be able to compute the type in an even better way here. This was discussed in JuliaLang/LinearAlgebra.jl#128 and actually, the issue isn't completely fixed. The problem is that L is unitless while U has units but they'd have to be stored in the same array. I think something like

UT = typeof(zero(T) - oneunit(T)/oneunit(T)*oneunit(T))
LT = typeof(oneunit(T)/oneunit(T))
S = promote_type(T, UT, LT)

should be correct. In most cases it should just collapse to a single type and for unitful elements, it should provide the right union.

[blegat]

@andreasnoack I have updated the PR based on your idea. The code now also works for quantities with units (tested with Unitful)