Structured array .= silently ignores unstructured non-zeros

[ararslan]

Consider the following:

julia> U = UpperTriangular(Float64[1 2 3; 0 1 2; 0 0 1])
3×3 UpperTriangular{Float64,Array{Float64,2}}:
 1.0  2.0  3.0
  ⋅   1.0  2.0
  ⋅    ⋅   1.0

julia> X = randn(3, 3)
3×3 Array{Float64,2}:
  0.9189    -0.942656  -0.825467
 -0.346463   0.829813   0.741945
 -0.239974  -0.232412   1.46933 

julia> U .+= X
3×3 UpperTriangular{Float64,Array{Float64,2}}:
 1.9189  1.05734  2.17453
  ⋅      1.82981  2.74194
  ⋅       ⋅       2.46933

Note that U was updated with the result of elementwise addition with X. However, U is still an UpperTriangular, even though the right hand side of the addition is not. Generally speaking, the result of adding two matrices is only upper triangular if both are.

It seems to me that the case of updating a structured matrix should not work except when the operation would produce a matrix of the same structure.

[mcognetta]

Seems like the most practical solution is to just disallow this with structured and non-structured arguments. Otherwise it would have to test for triangularity, diagonality, etc each time.

[ararslan]

Otherwise it would have to test for triangularity, diagonality, etc each time.

I'm not sure what you mean; we could simply whitelist types in a method signature so that dispatch knows which types it can use.

[mbauman]

Ah, dang. This is an issue for all structured array assignments. All those methods need to check for structure/zero-preserving-ness like the non-in-place broadcasts.

[mcognetta]

Sorry, we are on the same page. I misread this part:

It seems to me that the case of updating a structured matrix should not work except when the operation would produce a matrix of the same structure.

What I meant by this was that if we allowed any matrix type, then we would have to call istriu or whatever each time, which would be expensive in the non-structured cases.

A small problem is that we lose some cases where we have (for example) a LowerTriangular matrix that is actually diagonal as well.

[mbauman]

Note that it's not just argument types — it also depends on the function. For example, X .= sin.(X) works but X .= cos.(X) should not for all the structured matrices. We have this infrastructure, it just needed to be extended to the in-place methods.

JuliaLang/julia#31678