Skip to content
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

Specialize adding/subtracting mixed Upper/LowerTriangular #56149

Merged
merged 8 commits into from
Oct 18, 2024

Conversation

jishnub
Copy link
Contributor

@jishnub jishnub commented Oct 14, 2024

Fixes JuliaLang/LinearAlgebra.jl#1099
After this,

julia> using LinearAlgebra

julia> A = hermitianpart(rand(4, 4))
4×4 Hermitian{Float64, Matrix{Float64}}:
 0.387617  0.277226  0.67629   0.60678
 0.277226  0.894101  0.388416  0.489141
 0.67629   0.388416  0.100907  0.619955
 0.60678   0.489141  0.619955  0.452605

julia> B = UpperTriangular(A)
4×4 UpperTriangular{Float64, Hermitian{Float64, Matrix{Float64}}}:
 0.387617  0.277226  0.67629   0.60678
          0.894101  0.388416  0.489141
                   0.100907  0.619955
                            0.452605

julia> B - B'
4×4 Matrix{Float64}:
  0.0        0.277226   0.67629   0.60678
 -0.277226   0.0        0.388416  0.489141
 -0.67629   -0.388416   0.0       0.619955
 -0.60678   -0.489141  -0.619955  0.0

This preserves the band structure of the parent, if any:

julia> U = UpperTriangular(Diagonal(ones(4)))
4×4 UpperTriangular{Float64, Diagonal{Float64, Vector{Float64}}}:
 1.0  0.0  0.0  0.0
     1.0  0.0  0.0
         1.0  0.0
             1.0

julia> U - U'
4×4 Diagonal{Float64, Vector{Float64}}:
 0.0            
     0.0        
         0.0    
             0.0

This doesn't fully work with partly initialized matrices, and would need #55312 for that.

The abstract triangular methods now construct matrices using similar(parent(U), size(U)) so that the destinations are fully mutable.

julia> @invoke B::LinearAlgebra.AbstractTriangular - B'::LinearAlgebra.AbstractTriangular
4×4 Matrix{Float64}:
  0.0        0.277226   0.67629   0.60678
 -0.277226   0.0        0.388416  0.489141
 -0.67629   -0.388416   0.0       0.619955
 -0.60678   -0.489141  -0.619955  0.0

@jishnub jishnub added linear algebra Linear algebra arrays [a, r, r, a, y, s] bugfix This change fixes an existing bug labels Oct 14, 2024
@jishnub jishnub added the backport 1.11 Change should be backported to release-1.11 label Oct 14, 2024
@jishnub jishnub requested a review from dkarrasch October 17, 2024 06:31
@jishnub
Copy link
Contributor Author

jishnub commented Oct 17, 2024

The broken tests here would be fixed by #55312, so one of these should be rebased after the other is merged.

Copy link
Member

@dkarrasch dkarrasch left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure if these needs an update of the tests once the other PR is merged?

@KristofferC KristofferC mentioned this pull request Oct 18, 2024
43 tasks
@jishnub jishnub merged commit e33c6a8 into master Oct 18, 2024
6 of 7 checks passed
@jishnub jishnub deleted the jishnub/abstri_similar branch October 18, 2024 11:34
@jishnub jishnub removed the backport 1.11 Change should be backported to release-1.11 label Oct 21, 2024
KristofferC pushed a commit that referenced this pull request Oct 21, 2024
Fixes https://github.com/JuliaLang/julia/issues/56134
After this,
```julia
julia> using LinearAlgebra

julia> A = hermitianpart(rand(4, 4))
4×4 Hermitian{Float64, Matrix{Float64}}:
 0.387617  0.277226  0.67629   0.60678
 0.277226  0.894101  0.388416  0.489141
 0.67629   0.388416  0.100907  0.619955
 0.60678   0.489141  0.619955  0.452605

julia> B = UpperTriangular(A)
4×4 UpperTriangular{Float64, Hermitian{Float64, Matrix{Float64}}}:
 0.387617  0.277226  0.67629   0.60678
  ⋅        0.894101  0.388416  0.489141
  ⋅         ⋅        0.100907  0.619955
  ⋅         ⋅         ⋅        0.452605

julia> B - B'
4×4 Matrix{Float64}:
  0.0        0.277226   0.67629   0.60678
 -0.277226   0.0        0.388416  0.489141
 -0.67629   -0.388416   0.0       0.619955
 -0.60678   -0.489141  -0.619955  0.0
```
This preserves the band structure of the parent, if any:
```julia
julia> U = UpperTriangular(Diagonal(ones(4)))
4×4 UpperTriangular{Float64, Diagonal{Float64, Vector{Float64}}}:
 1.0  0.0  0.0  0.0
  ⋅   1.0  0.0  0.0
  ⋅    ⋅   1.0  0.0
  ⋅    ⋅    ⋅   1.0

julia> U - U'
4×4 Diagonal{Float64, Vector{Float64}}:
 0.0   ⋅    ⋅    ⋅ 
  ⋅   0.0   ⋅    ⋅ 
  ⋅    ⋅   0.0   ⋅ 
  ⋅    ⋅    ⋅   0.0
```
This doesn't fully work with partly initialized matrices, and would need
#55312 for that.

The abstract triangular methods now construct matrices using
`similar(parent(U), size(U))` so that the destinations are fully
mutable.
```julia
julia> @invoke B::LinearAlgebra.AbstractTriangular - B'::LinearAlgebra.AbstractTriangular
4×4 Matrix{Float64}:
  0.0        0.277226   0.67629   0.60678
 -0.277226   0.0        0.388416  0.489141
 -0.67629   -0.388416   0.0       0.619955
 -0.60678   -0.489141  -0.619955  0.0
```

---------

Co-authored-by: Daniel Karrasch <[email protected]>
KristofferC added a commit that referenced this pull request Nov 21, 2024
Backported PRs:
- [x] #55886 <!-- irrationals: restrict assume effects annotations to
known types -->
- [x] #55867 <!-- update `hash` doc string: `widen` not required any
more -->
- [x] #56084 <!-- slightly improve inference in precompilation code -->
- [x] #56088 <!-- make `Base.ANSIIterator` have a concrete field -->
- [x] #54093 <!-- Fix `JULIA_CPU_TARGET` being propagated to workers
precompiling stdlib pkgimages -->
- [x] #56165 <!-- Fix markdown list in installation.md -->
- [x] #56148 <!-- Make loading work when stdlib deps are missing in the
manifest -->
- [x] #56174 <!-- Fix implicit `convert(String, ...)` in several places
-->
- [x] #56159 <!-- Add invalidation barriers for `displaysize` and
`implicit_typeinfo` -->
- [x] #56089 <!-- Call `MulAddMul` instead of multiplication in
_generic_matmatmul! -->
- [x] #56195 <!-- Include default user depot when JULIA_DEPOT_PATH has
leading empty entry -->
- [x] #56215 <!-- [REPL] fix lock ordering mistake in load_pkg -->
- [x] #56251 <!-- REPL: run repl hint generation for modeswitch chars
when not switching -->
- [x] #56092 <!-- stream: fix reading LibuvStream into array -->
- [x] #55870 <!-- fix infinite recursion in `promote_type` for
`Irrational` -->
- [x] #56227 <!-- Do not call `rand` during sysimage precompilation -->
- [x] #55741 <!-- Change annotations to use a NamedTuple -->
- [x] #56149 <!-- Specialize adding/subtracting mixed
Upper/LowerTriangular -->
- [x] #56214 <!-- fix precompile process flags -->
- [x] #54471
- [x] #55622
- [x] #55704
- [x] #55764
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
arrays [a, r, r, a, y, s] bugfix This change fixes an existing bug linear algebra Linear algebra
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Bug Report: Error when Subtracting AbstractTriangular wrapping Symmetric/Hermitian Matrices in LinearAlgebra
2 participants