Fix (l/r)mul! with Diagonal/Bidiagonal

[jishnub]

Currently, rmul!(A::AbstractMatirx, D::Diagonal) calls mul!(A, A, D), but this isn't a valid call, as mul! assumes no aliasing between the destination and the matrices to be multiplied. As a consequence,

julia> B = Bidiagonal(rand(4), rand(3), :L)
4×4 Bidiagonal{Float64, Vector{Float64}}:
 0.476892   ⋅         ⋅         ⋅ 
 0.353756  0.139188   ⋅         ⋅ 
  ⋅        0.685839  0.309336   ⋅ 
  ⋅         ⋅        0.369038  0.304273

julia> D = Diagonal(rand(size(B,2)));

julia> rmul!(B, D)
4×4 Bidiagonal{Float64, Vector{Float64}}:
 0.0   ⋅    ⋅    ⋅ 
 0.0  0.0   ⋅    ⋅ 
  ⋅   0.0  0.0   ⋅ 
  ⋅    ⋅   0.0  0.0

julia> B
4×4 Bidiagonal{Float64, Vector{Float64}}:
 0.0   ⋅    ⋅    ⋅ 
 0.0  0.0   ⋅    ⋅ 
  ⋅   0.0  0.0   ⋅ 
  ⋅    ⋅   0.0  0.0

This is clearly nonsense, and happens because the internal _mul! function assumes that it can safely overwrite the destination with zeros before carrying out the multiplication. This is fixed in this PR by using broadcasting instead. The current implementation is generally equally performant, albeit occasionally with a minor allocation arising from reshapeing an Array.

A similar problem also exists in l/rmul! with Bidiaognal, but that's a little harder to fix while remaining equally performant.

[dkarrasch]

This is going back to where we came from, so it would be good to double check what the original concern was to make it call mul!(A, A, D). I guess this is a problem only with products of structured/banded matrices? One issue that I remember is that permuteddims allocates, and we didn't like that, so we wrote out the loops. Just to check: is there another option to fix this at the _mul! level?

[jishnub]

The problem is that mul! is documented to assume no aliasing, so a call like mul!(A, A, D) will invariably lead to problems downstream. We may fix this for the methods in LinearAlgebra, but libraries that define more specific methods may still encounter errors.

One option might be to write out the loop here instead of broadcasting. My idea was that array types that specialise broadcasting would be able to access optimised methods.

[dkarrasch]

I found it: #42321.

[dkarrasch]

This goes over the entire matrix and misses all the optimizations that we have in _mul! for the banded matrices. I guess that was (one of) the reason(s) to initially make [l/r]mul! run by mul!.

[jishnub]

Yes, we would need to specialize lmul! and rmul! for banded matrices. The fallback would probably need to go over the entire matrix, although we could introduce a bandwidth function to limit the loop bounds. I don't think that we can reliably call mul! here, although we may reuse some of the code from the specialized methods.

[KristofferC]

This seems to fail on CI on 1.11:

Error During Test at /cache/build/tester-amdci5-15/julialang/julia-release-1-dot-11/julia-f4b43765c5/share/julia/stdlib/v1.11/LinearAlgebra/test/tridiag.jl:846
--
  | Test threw exception
  | Expression: rmul!(copy(T), D) ≈ T * D
  | MethodError: no method matching zero(::Type{Main.OffsetArrays.OffsetMatrix{Float64, Matrix{Float64}}})
  | The function `zero` exists, but no method is defined for this combination of argument types.
  |  
  | Closest candidates are:
  | zero(!Matched::Type{Union{}}, Any...)
  | @ Base number.jl:310
  | zero(!Matched::Type{Missing})
  | @ Base missing.jl:104
  | zero(!Matched::Type{Main.Test26Main_LinearAlgebra_structuredbroadcast.TestStructuredBroadcast.Zero36193})
  | @ Main.Test26Main_LinearAlgebra_structuredbroadcast /cache/build/tester-amdci5-15/julialang/julia-release-1-dot-11/julia-f4b43765c5/share/julia/stdlib/v1.11/LinearAlgebra/test/structuredbroadcast.jl:255
  | ...
  |  
  | Stacktrace:
  | [1] getindex
  | @ /cache/build/tester-amdci5-15/julialang/julia-release-1-dot-11/julia-f4b43765c5/share/julia/stdlib/v1.11/LinearAlgebra/src/tridiag.jl:689 [inlined]
  | [2] _getindex
  | @ ./abstractarray.jl:1361 [inlined]
  | [3] getindex
  | @ ./abstractarray.jl:1315 [inlined]
  | [4] iterate
  | @ ./abstractarray.jl:1212 [inlined]
  | [5] _foldl_impl(op::Base.MappingRF{typeof(LinearAlgebra.norm), Base.BottomRF{typeof(max)}}, init::Base._InitialValue, itr::LinearAlgebra.Tridiagonal{Main.OffsetArrays.OffsetMatrix{Float64, Matrix{Float64}}, Vector{Main.OffsetArrays.OffsetMatrix{Float64, Matrix{Float64}}}})