JuliaLang
Can't sum partially initialized Hermitian matrices that are not isbits with Hermitian Diagonal matrices
Edge case of JuliaLang/julia#52895, missed by the bugfix JuliaLang/julia#52942. @jishnub I think this will be easy for you.
id = Hermitian(I(2))
m = Hermitian(Matrix{BigFloat}(undef,2,2))
m.data[1,1] = 1
m.data[1,2] = 2
m.data[2,2] = 3
m+id
The issue seems to be with summing a partially uninitialised Matrix and a Diagonal, which lacks a specialised method. I'll try to come up with a fix, which may also help with performance.
I've looked into it. The problem is that when you call UpperTriangular(id.data) the result is a UpperTriangular{Bool, Diagonal{Bool, Vector{Bool}}}, that is, it is not a strided matrix, and therefore doesn't get caught by the specialized method for summing triangular matrices, but instead goes to the default sum, which cannot handle the undefs.
Argh, this means that whenever Hermitian is wrapping an abstract matrix the same problem will appear. There must be a better solution than writing a specialized method for every single type of abstract matrix. Maybe it's better to do the sum over the upper (or lower) triangular for everything, and let the sparse arrays define their own specialized methods?
Yes, that makes sense. Perhaps we need to expose some function that falls back to summing over the triangular parts, and one that may be specialized by packages. The reason the current implementation was limited to StridedArrays was to avoid breakages by using scalar indexing for custom arrays, so we need to be a bit careful about that.
