Symmetric(StaticMatrix) causes several problems

Hello,

Thank you for this package,

Starting to use it, I encountered several problems with Symmetric(StaticMatrix). This prevents me to have some generic code (working without modification for StaticMatrices and Julia regular matrix types).

Example 1: M\v strange behavior

M\v does not work with this invertible M matrix

julia> using LinearAlgebra, StaticArrays

julia> M = Symmetric(@SMatrix Float64[1 1;1 2])
2×2 Symmetric{Float64, SMatrix{2, 2, Float64, 4}}:
 1.0  1.0
 1.0  2.0

julia> v = @SVector Float64[1,2]
2-element SVector{2, Float64} with indices SOneTo(2):
 1.0
 2.0

julia> M\v
ERROR: setindex!(::SMatrix{2, 2, Float64, 4}, value, ::Int) is not defined.
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:33
  [2] setindex!(a::SMatrix{2, 2, Float64, 4}, value::Float64, i::Int64)
    @ StaticArrays ~/.julia/packages/StaticArrays/OWJK7/src/indexing.jl:3
  [3] macro expansion
    @ ~/.julia/packages/StaticArrays/OWJK7/src/indexing.jl:66 [inlined]
  [4] _setindex!_scalar
    @ ~/.julia/packages/StaticArrays/OWJK7/src/indexing.jl:46 [inlined]
  [5] setindex!
    @ ~/.julia/packages/StaticArrays/OWJK7/src/indexing.jl:42 [inlined]
  [6] copytri!
    @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/matmul.jl:514 [inlined]
  [7] copytri!
    @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/matmul.jl:510 [inlined]
  [8] lu!(A::Symmetric{Float64, SMatrix{2, 2, Float64, 4}}, pivot::Val{true}; check::Bool)
    @ LinearAlgebra /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lu.jl:89
  [9] lu(A::Symmetric{Float64, SMatrix{2, 2, Float64, 4}}, pivot::Val{true}; check::Bool)
    @ LinearAlgebra /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lu.jl:273
 [10] lu(A::Symmetric{Float64, SMatrix{2, 2, Float64, 4}}, pivot::Val{true}) (repeats 2 times)
    @ LinearAlgebra /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lu.jl:272
 [11] \(A::Symmetric{Float64, SMatrix{2, 2, Float64, 4}}, B::SVector{2, Float64})
    @ LinearAlgebra /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/generic.jl:1136
 [12] top-level scope
    @ REPL[119]:1

however, the same procedure works with the same matrix type but different component values:

julia> M = Symmetric(@SMatrix Float64[1 0;0 2])
2×2 Symmetric{Float64, SMatrix{2, 2, Float64, 4}}:
 1.0  0.0
 0.0  2.0

julia> M\v
2-element SVector{2, Float64} with indices SOneTo(2):
 1.0
 1.0

nb: same observation if one uses MMatrix instead of SMatrix.

Example 2: M + I

julia> using LinearAlgebra, StaticArrays

julia> M = Symmetric(@SMatrix Float64[1 0;0 2])
2×2 Symmetric{Float64, SMatrix{2, 2, Float64, 4}}:
 1.0  0.0
 0.0  2.0

julia> M + I
ERROR: setindex!(::SMatrix{2, 2, Float64, 4}, value, ::Int) is not defined.
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:33
 [2] setindex!(a::SMatrix{2, 2, Float64, 4}, value::Float64, i::Int64)
   @ StaticArrays ~/.julia/packages/StaticArrays/OWJK7/src/indexing.jl:3
 [3] macro expansion
   @ ~/.julia/packages/StaticArrays/OWJK7/src/indexing.jl:66 [inlined]
 [4] _setindex!_scalar
   @ ~/.julia/packages/StaticArrays/OWJK7/src/indexing.jl:46 [inlined]
 [5] setindex!
   @ ~/.julia/packages/StaticArrays/OWJK7/src/indexing.jl:42 [inlined]
 [6] setindex!(A::Symmetric{Float64, SMatrix{2, 2, Float64, 4}}, v::Float64, i::Int64, j::Int64)
   @ LinearAlgebra /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/symmetric.jl:226
 [7] +(A::Symmetric{Float64, SMatrix{2, 2, Float64, 4}}, J::UniformScaling{Bool})
   @ LinearAlgebra /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/uniformscaling.jl:219
 [8] top-level scope
   @ REPL[132]:1

Looking at the code I understand the origin of the problem:

function (+)(A::AbstractMatrix, J::UniformScaling)
    checksquare(A)
    B = copy_oftype(A, Base._return_type(+, Tuple{eltype(A), typeof(J)}))
    @inbounds for i in axes(A, 1)
        B[i,i] += J
    end
    return B
end

The function copy_oftype() create a copy of of M.data of type SMatrix then try to modify its diagonal.

As expected, this does not cause an issue with MMatrix:

Julia> M = Symmetric(@MMatrix Float64[1 1;1 2])
2×2 Symmetric{Float64, MMatrix{2, 2, Float64, 4}}:
 1.0  1.0
 1.0  2.0

julia> M + I
2×2 Symmetric{Float64, MMatrix{2, 2, Float64, 4}}:
 2.0  1.0
 1.0  3.0

Sorry to bring this to the table,
Thanks.
Vincent

Reproduce both examples with Julia, LinearAlgebra, and StaticArrays, then trace the mentioned paths through LinearAlgebra's lu.jl/matmul.jl and uniformscaling.jl, along with StaticArrays' indexing.jl. Done means Symmetric SMatrix and MMatrix values support the reported M\v and M + I operations without attempting unsupported mutation.