Symmetric matrices do not permit editing off-diagonal elements. I have two matrices S1 and S2 and I’m trying to evaluate S1 + S2 and save the result into S2. Typically S2 .+= S1 would achieve this but this doesn’t work for Symmetric matrices.
julia> S1 = Symmetric(ones(3,3), :U);
julia> S2 = similar(S1);
julia> S2 .+= S1
ERROR: ArgumentError: Cannot set a non-diagonal index in a symmetric matrix
Stacktrace:
[1] setindex!
@ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.7/LinearAlgebra/src/symmetric.jl:225 [inlined]
[2] _setindex!
@ ./abstractarray.jl:1297 [inlined]
[3] setindex!
@ ./abstractarray.jl:1267 [inlined]
[4] macro expansion
@ ./broadcast.jl:983 [inlined]
[5] macro expansion
@ ./simdloop.jl:77 [inlined]
[6] copyto!
@ ./broadcast.jl:982 [inlined]
[7] copyto!
@ ./broadcast.jl:935 [inlined]
[8] materialize!
@ ./broadcast.jl:893 [inlined]
[9] materialize!(dest::Symmetric{Float64, Matrix{Float64}}, bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2}, Nothing, typeof(+), Tuple{Symmetric{Float64, Matrix{Float64}}, Symmetric{Float64, Matrix{Float64}}}})
@ Base.Broadcast ./broadcast.jl:890
[10] top-level scope
@ REPL[24]:1
Strangely assignment works:
julia> S2 .= S1
3×3 Symmetric{Float64, Matrix{Float64}}:
1.0 1.0 1.0
1.0 1.0 1.0
1.0 1.0 1.0
julia> all(isone, S2)
true
so perhaps one may set off-diagonal elements after all?
I realize that there’s a way around this by using the parent arrays:
julia> parent(S2) .= 0
julia> parent(S2) .+= parent(S1);
julia> S2
3×3 Symmetric{Float64, Matrix{Float64}}:
1.0 1.0 1.0
1.0 1.0 1.0
1.0 1.0 1.0
but this seems to be an uncessary foray into the internals.
Also
julia> S2 += S1
3×3 Symmetric{Float64, Matrix{Float64}}:
3.0 2.0 2.0
2.0 2.0 2.0
2.0 2.0 2.0
works as expected, but this creates a new binding. This is not what I want.