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 #56134
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
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
41 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 #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]>
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