Skip to content

Conversation

jishnub
Copy link
Member

@jishnub jishnub commented Jan 17, 2024

This ensures that only the triangular indices are accessed for strided parent matrices. Fix #52895

julia> M = Matrix{Complex{BigFloat}}(undef, 2, 2);

julia> M[1,1] = M[2,2] = M[1,2] = 2;

julia> H = Hermitian(M)
2×2 Hermitian{Complex{BigFloat}, Matrix{Complex{BigFloat}}}:
 2.0+0.0im  2.0+0.0im
 2.0-0.0im  2.0+0.0im

julia> H + H # works after this
2×2 Hermitian{Complex{BigFloat}, Matrix{Complex{BigFloat}}}:
 4.0+0.0im  4.0+0.0im
 4.0-0.0im  4.0+0.0im

This also provides a speed-up in several common cases (allocations mentioned only when they differ):

julia> H = Hermitian(rand(ComplexF64,1000,1000));

julia> H2 = Hermitian(rand(ComplexF64,1000,1000),:L);
Operation master PR
-H 2.247 ms 1.384 ms
real(H) 1.544 ms 1.175 ms
H + H 2.288 ms 1.978 ms
H + H2 5.139 ms 3.287 ms
isdiag(H) 23.042 ns (1 allocation: 16 bytes) 16.778 ns (0 allocations: 0 bytes)

I'm not entirely certain why isdiag(H) allocates on master, as union splitting should handle this automatically, but manually splitting the union appears to help.

@jishnub jishnub added linear algebra Linear algebra arrays [a, r, r, a, y, s] labels Jan 17, 2024
@jishnub jishnub changed the title Reroute symmetric algebraic functions through triangular Reroute algebraic functions for Symmetric/Hermitian through triangular Jan 17, 2024
@araujoms
Copy link
Contributor

This PR exceeds my most optimistic expectations. Note that it doesn't work for mismatched uplo, though, because that will be wrapped by Adjoint, and that falls back to the non-symmetric +. Example:

function matrix_sum(T)
   d = 2
   a = Hermitian(Matrix{T}(undef,d,d),:U)
   b = Hermitian(Matrix{T}(undef,d,d),:L)
   for i=1:d
       for j=i:d
           a.data[i,j] = T(1)
           b.data[j,i] = T(1)
       end
   end
   a+b
end    

@jishnub
Copy link
Member Author

jishnub commented Jan 18, 2024

Yes, fixing the mismatched uplo case would require #52530

@jishnub
Copy link
Member Author

jishnub commented Jan 29, 2024

The addition/subtraction with mismatched uplo is fixed now:

julia> matrix_sum(BigFloat)
2×2 Hermitian{BigFloat, Matrix{BigFloat}}:
 2.0  2.0
 2.0  2.0

@jishnub jishnub force-pushed the jishnub/hermsymtri branch from 195cea2 to 2c5bd94 Compare February 1, 2024 18:21
@jishnub
Copy link
Member Author

jishnub commented Feb 2, 2024

This should be good to go

@jishnub
Copy link
Member Author

jishnub commented Feb 5, 2024

Gentle bump

Copy link
Member

@oscardssmith oscardssmith left a comment

Choose a reason for hiding this comment

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

This looks reasonable to me (although i do kind of hate that we are doing all of this with the uplo flag. Using characters is ugly. None of that is the fault of this PR though.)

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.

Nice work!

@dkarrasch dkarrasch added merge me PR is reviewed. Merge when all tests are passing and removed status: waiting for PR reviewer labels Feb 8, 2024
@jishnub jishnub merged commit 27b31d1 into master Feb 9, 2024
@jishnub jishnub deleted the jishnub/hermsymtri branch February 9, 2024 10:56
@oscardssmith oscardssmith removed the merge me PR is reviewed. Merge when all tests are passing label Feb 9, 2024
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] linear algebra Linear algebra

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Cannot sum partially initialized Hermitian matrices that are not an isbits type

4 participants