Add *(::Diagonal, ::Diagonal, ::Diagonal) (#49005)

[dlfivefifty]

This adds a special case for triple multiplication of diagonals to ensure that the result is a Diagonal, see JuliaLang/LinearAlgebra.jl#992.

[jishnub]

Currently, this uses broadcasting, and I wonder if it might not be better to make the result of the broadcast operation return a Diagonal instead of special-casing this 3-term multiplication?

julia> (1:4) .* Diagonal(1:4) .* (1:4)' # may return a Diagonal
4×4 Matrix{Int64}:
 1  0   0   0
 0  8   0   0
 0  0  27   0
 0  0   0  64

In fact, each component may be evaluated to a Diagonal here:

julia> (1:4) .* Diagonal(1:4)
4×4 Matrix{Int64}:
 1  0  0   0
 0  4  0   0
 0  0  9   0
 0  0  0  16

julia> Diagonal(1:4) .* (1:4)'
4×4 Matrix{Int64}:
 1  0  0   0
 0  4  0   0
 0  0  9   0
 0  0  0  16

[dlfivefifty]

That would require special case-ing broadcast for * which would be a bigger project. I'm just trying to fix the regression introduced by a special case *(::Diagonal, ::AbstractMatrix, ::Diagonal) added.

[jishnub]

This is tricky with floating-point values, as NaN and Inf don't preserve the Diagonal nature:

julia> Diagonal([1,NaN]) * Diagonal(1:2) * Diagonal([1,NaN])
2×2 Matrix{Float64}:
   1.0  NaN
 NaN    NaN

[KristofferC]

This is tricky with floating-point values, as NaN and Inf don't preserve the Diagonal nature:

This is already the case with e.g. sparse arrays. And it was like that on 1.8 as well.

[martinholters]

How does * with even more Diagonal arguments fare?

[jishnub]

How does * with even more Diagonal arguments fare?

julia> Diagonal(1:2) * Diagonal(1:2) * Diagonal(1:2) * Diagonal(1:2)
2×2 Diagonal{Int64, Vector{Int64}}:
 1   ⋅
 ⋅  16

The 4-term multiplication, as implemented, uses combinations of 2-term multiplications, which returns a Diagonal. Higher number of terms fall back to afoldl, which uses the 2-term multiplication as well. It's only the 3-term multiplication that uses broadcasting and materializes the matrix.

[dlfivefifty]

I added a 4-term test so any future regressions will be picked up