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

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

Closed
araujoms opened this issue Jan 13, 2024 · 1 comment · Fixed by #52942
Closed

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

araujoms opened this issue Jan 13, 2024 · 1 comment · Fixed by #52942

Comments

@araujoms
Copy link
Contributor

If T is an isbitstype, creating an uninitialized matrix a = Matrix{T}(undef,d,d) will initialize it with whatever is in memory, so the resulting matrix will be just a regular numerical matrix. Otherwise the matrix will be filled with #undef, and I'll get an error if I try to access them. So far so good.

I wanted to create, however, a Hermitian matrix, and initialize only the upper triangular. If I try to sum it with another Hermitian matrix, however, I get the error, even though I don't actually want to access the uninitialized elements. The code below shows what I mean, if you call it with matrix_sum(BigFloat).

Looking at the code of +, I see that it just adds the underlying parent and wraps it again with Hermitian, which by necessity accesses all elements. If one would sum instead only the upper (or lower) triangular this bug would be fixed, and presumably we would get a speedup as well.

A bit unrelated, but I noticed that GenericLinearAlgebra doesn't choke on such matrices, it can calculate the eigenvalues, for example.

using LinearAlgebra

function matrix_sum(T)
   d = 2
   a = Hermitian(Matrix{T}(undef,d,d))
   for i=1:d
       for j=i:d
           a.data[i,j] = T(1)
       end
   end
   b = a + a
end        
@jishnub
Copy link
Contributor

jishnub commented Jan 14, 2024

This should benefit from the recently incorporated changes for AbstractTriangular. We may wrap the parents in an appropriate triangular wrapper before adding them, and finally unwrap the result.

jishnub added a commit that referenced this issue Feb 9, 2024
…gular (#52942)

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

```julia
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
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.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants