JuliaLang / LinearAlgebra.jl

Julia Linear Algebra standard library
Other
18 stars 5 forks source link

When not trying to set non-diagonal broadcast get: Cannot set a non-diagonal index in a Hermitian/Symmetric matrix #772

Open oxinabox opened 4 years ago

oxinabox commented 4 years ago
julia> A = [1 2; 3 4];

julia> AH = Symmetric(A);

julia> AH .+= Diagonal(A)
ERROR: ArgumentError: Cannot set a non-diagonal index in a symmetric matrix
Stacktrace:
  [1] setindex!
    @ /usr/local/src/julia/julia-master/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/symmetric.jl
:225 [inlined]
  [2] _setindex!
    @ ./abstractarray.jl:1247 [inlined]
  [3] setindex!
    @ ./abstractarray.jl:1217 [inlined]
  [4] macro expansion
    @ ./broadcast.jl:932 [inlined]
  [5] macro expansion
    @ ./simdloop.jl:77 [inlined]
  [6] copyto!
    @ ./broadcast.jl:931 [inlined]
  [7] copyto!
    @ ./broadcast.jl:886 [inlined]
  [8] materialize!
    @ ./broadcast.jl:848 [inlined]
  [9] materialize!(dest::Symmetric{Int64, Matrix{Int64}}, bc::Base.Broadcast.Broadcasted{LinearAlgebra.StructuredMatrixStyle{Diagonal{Int64, Vector{Int64}}}, Nothing, typeof(+), Tuple{Symmetric{Int64, Matrix{Int64}}, Diagonal{Int64, Vector{Int64}}}})
    @ Base.Broadcast ./broadcast.jl:845
 [10] top-level scope
    @ REPL[17]:1

Same happens to Hermian

I can do it it I am not broadcasting, but it alocates a new AH and just rebinds the name. Which is allocating.

julia> objectid(AH)
0x4b328c6c806ac9c2

julia> AH += Diagonal(A)
2×2 Symmetric{Int64, Matrix{Int64}}:
 4  2
 2  8

julia> objectid(AH)
0x67c658733e89b4a9

In contast UpperTrianguar is fine with this:

julia> AU = UpperTriangular(A)
2×2 UpperTriangular{Int64, Matrix{Int64}}:
 2  2
 ⋅  4

julia> AU .+= Diagonal(A)
2×2 UpperTriangular{Int64, Matrix{Int64}}:
 4  2
 ⋅  8

I can hack around it via:

julia> parent(AH) .+= Diagonal(A);

julia> AH
2×2 Symmetric{Int64, Matrix{Int64}}:
 8   2
 2  16
jishnub commented 3 years ago

What I find inconsistent about the present behavior is:

julia> using LinearAlgebra

julia> S = Symmetric(ones(2,2));

julia> S .= 4;

julia> S
2×2 Symmetric{Float64, Matrix{Float64}}:
 4.0  4.0
 4.0  4.0

julia> S2 = similar(S);

julia> S2 .= S
2×2 Symmetric{Float64, Matrix{Float64}}:
 4.0  4.0
 4.0  4.0

so .. it is possible to set off-diagonal elements after all? Not only is the error confusing, it is also misleading.