JuliaLang / julia / #37835

26%

Build Type

push

local

Committed by web-flow

Commit Message

Fix `(l/r)mul!` with `Diagonal`/`Bidiagonal` (#55052)

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
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 `reshape`ing an `Array`.

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

Run Details

92 of 92 new or added lines in 2 files covered. (100.0%)

17 existing lines in 9 files now uncovered.

77496 of 88551 relevant lines covered (87.52%)

15114947.4 hits per line