Disambiguate structured and abstract matrix multiplication

[jishnub]

The various methods of the form *(::AbstractMatrix, ::Diagonal) for structured matrices lead to method ambiguities in packages that define *(::MyMatrix, ::AbstractMatrix). However, these methods only seem to be defined here to use the first argument to construct the destination, as opposed to the fallback method that uses the second argument. Instead of specializing *, we might as well specialize the call to similar that chooses which matrix to use to construct the destination. This reduces code duplication, as well as avoids these ambiguities altogether. This also makes the destination types of A * B and B * A identical if one of the two is a StructuredMatrix.

After this, for example, both of these produce sparse matrices:

julia> S = sprand(4,4,0.2); B = Bidiagonal(ones(4), ones(3), :U);

julia> S * B
4×4 SparseMatrixCSC{Float64, Int64} with 2 stored entries:
  ⋅    ⋅   0.706461  0.706461
  ⋅    ⋅    ⋅         ⋅ 
  ⋅    ⋅    ⋅         ⋅ 
  ⋅    ⋅    ⋅         ⋅ 

julia> B * S
4×4 SparseMatrixCSC{Float64, Int64} with 1 stored entry:
  ⋅    ⋅   0.706461   ⋅ 
  ⋅    ⋅    ⋅         ⋅ 
  ⋅    ⋅    ⋅         ⋅ 
  ⋅    ⋅    ⋅         ⋅ 

The former is a Matrix on master.

Close #51354
Address JuliaStats/PDMats.jl#176 on nightly

This PR changes the following behavior:

julia> T = UpperTriangular(Number[i+j for i in 1:4, j in 1:4]);

julia> v = rand(4);

julia> T * v |> eltype
Float64

julia> Array(T) * v |> eltype
Any

to

julia> T * v |> eltype
Any

I'm unsure if the choice of _init_eltype was deliberate here (there is a test that specifically checks that _init_eltype isn't being used for general matrices).

[dkarrasch]

That is nice. I believe we have similar * methods in SparseArrays.jl that tell it to call similar on the first argument. We should use this mechanism over there as well.

[dkarrasch]

I have pushed a commit that removes the clumsy _init* stuff, since unitful arithmetic seems to work just fine with promote_op(matprod, eltype(A), eltype(B)). The inconsistency between triangular and plain matrices was unintended. The promotion mechanism was introduced in a big triangular overhaul, and back then, I didn't extrapolate to the whole LinearAlgebra library.

[dkarrasch]

Hm, in fact, matprod_dest could be equally well used in the solves, seems like...?

[jishnub]

Looks possible. I'll try to get to this in a couple of days. Thanks for making this simpler!

[dkarrasch]

BTW, note #35166.

[jishnub]

Is there anything left to be done here?

[dkarrasch]

I have some local stuff for the determination of the eltype of the destination array, but that could go into another PR.

[singularitti]

Nice work @jishnub @dkarrasch, thank you!