-
-
Notifications
You must be signed in to change notification settings - Fork 5.4k
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Fix (l/r)mul!
with Diagonal
/Bidiagonal
#55052
Conversation
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 |
The problem is that 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. |
I found it: #42321. |
This goes over the entire matrix and misses all the optimizations that we have in |
Yes, we would need to specialize |
(l/r)mul!
with Diagonal
(l/r)mul!
with Diagonal
/Bidiagonal
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. (cherry picked from commit 262b40a)
This seems to fail on CI on 1.11:
|
Backported PRs: - [x] #54201 <!-- Fix generic triangular solves with empty matrices --> - [x] #54358 <!-- Create `jl_clear_coverage_data` to dynamically reset coverage --> - [x] #54908 <!-- LazyString in interpolated error messages in threadingconstructs --> - [x] #54952 <!-- LAPACK: Avoid repr call in `chkvalidparam` --> - [x] #54898 <!-- fix concurrent module loading return value --> - [x] #55082 <!-- Add fast method for copyto!(::Memory, ::Memory) --> - [x] #55084 <!-- Use triple quotes in TOML.print when string contains newline --> - [x] #55141 <!-- Update the aarch64 devdocs to reflect the current state of its support --> - [x] #54955 <!-- don't throw EOFError from sleep --> - [x] #54871 <!-- Make warn missed transformations pass optional --> - [x] #55178 <!-- Compat for `Base.@nospecializeinfer` --> - [x] #55197 <!-- compat notice for a[begin] indexing --> - [x] #54917 <!-- Fix potential underrun with annotation merging --> - [x] #55209 <!-- correction to compat notice for a[begin] --> - [x] #55203 <!-- document mutable struct const fields --> - [x] #54791 <!-- Bump libblastrampoline to v5.10.1 --> - [x] #54950 <!-- SuiteSparse: Bump version --> - [x] #54956 <!-- Fix accidental early evaluation of imported `using` binding --> - [x] #54996 <!-- inference: add missing `MustAlias` widening in `_getfield_tfunc` --> - [x] #55070 <!-- LinearAlgebra: LazyString in error messages for Diagonal/Bidiagonal --> - [x] #54574 <!-- Make ScopedValues public --> - [x] #54739 <!-- finish implementation of upgradable stdlibs --> - [x] #54965 <!-- RFC: Make `include_dependency(path; track_content=true)` the default --> - [x] #53286 <!-- Raise an error when using `include_dependency` with non-existent file or directory --> - [x] #55066 <!-- fix loading of repeated/concurrent modules --> - [x] #52694 <!-- Reinstate similar for AbstractQ for backward compatibility --> - [x] #55218 <!-- Artifacts: use a different way of getting the UUID of a module --> - [x] #54891 <!-- #54739-related fixes for loading stdlibs --> - [x] #55072 <!-- trace-compile: don't generate `precompile` statements for OpaqueClosure methods --> - [x] #55188 <!-- Make Core.TypeofUnion use the type method table --> Need manual backport: - [ ] #55052 <!-- Fix `(l/r)mul!` with `Diagonal`/`Bidiagonal` --> Contains multiple commits, manual intervention needed: Non-merged PRs with backport label: - [ ] #55169 <!-- `propertynames` for SVD respects private argument --> - [ ] #55148 <!-- Random: Mark unexported public symbols as public --> - [ ] #55017 <!-- TOML: Make `Dates` a type parameter --> - [ ] #55013 <!-- [docs] change docstring to match code --> - [ ] #54919 <!-- Fix annotated join with non-concrete eltype iters --> - [ ] #54457 <!-- Make `String(::Memory)` copy --> - [ ] #53957 <!-- tweak how filtering is done for what packages should be precompiled --> - [ ] #51479 <!-- prevent code loading from lookin in the versioned environment when building Julia --> - [ ] #50813 <!-- More doctests for Sockets and capitalization fix --> - [ ] #50157 <!-- improve docs for `@inbounds` and `Base.@propagate_inbounds` --> - [ ] #41244 <!-- Fix shell `cd` error when working dir has been deleted -->
Currently,
rmul!(A::AbstractMatirx, D::Diagonal)
callsmul!(A, A, D)
, but this isn't a valid call, asmul!
assumes no aliasing between the destination and the matrices to be multiplied. As a consequence,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 fromreshape
ing anArray
.A similar problem also exists in
l/rmul!
withBidiaognal
, but that's a little harder to fix while remaining equally performant.