Support structured array types

[putianyi889]

Since #222 , fill +/- loses support for most structured arrays.

julia> s=Symmetric(rand(3,3))
3×3 Symmetric{Float64, Matrix{Float64}}:
 0.110199  0.379641  0.536482
 0.379641  0.518226  0.384931
 0.536482  0.384931  0.955319

julia> s+Zeros(3,3)
3×3 Matrix{Float64}:
 0.110199  0.379641  0.536482
 0.379641  0.518226  0.384931
 0.536482  0.384931  0.955319

It's not practical to list the types since ToeplitzMatrices, LazyArrays, StaticArrays, etc., are not dependencies. This could be solved in Julia 1.9 using package extensions.

[jishnub]

See JuliaLang/julia#49033 for a different approach.

In general, perhaps this should be done through broadcast styles?

[putianyi889]

In general, perhaps this should be done through broadcast styles?

That's true, but I don't know how broadcast works. Will broadcasting not be affected by ambiguity? StaticArrays uses map so will be fine, but what about general cases?

[putianyi889]

Prior to 1.9, can we have mini packages that resolve the ambiguities? For example, FillArraysXStaticArrays.jl

[putianyi889]

one could use fill_add to work around for now.

[dlfivefifty]

I think broadcasting is the right approach. I believe if there's an ambiguity in broadcasting it falls back to the default behaviour.

That is: we should overload BroadcastStyle and do intelligent broadcasting there.

[dlfivefifty]

Note Symmetric addition is handled at the +:

julia> S .+ S
5×5 Matrix{Float64}:
 -1.90614    2.27016   -0.412719   0.331858   3.62299
  2.27016   -1.48817    0.456745   2.18015    0.599337
 -0.412719   0.456745   4.30436   -2.37771    0.431586
  0.331858   2.18015   -2.37771    0.966284   0.971221
  3.62299    0.599337   0.431586   0.971221  -1.9591

julia> S + S
5×5 Symmetric{Float64, Matrix{Float64}}:
 -1.90614    2.27016   -0.412719   0.331858   3.62299
  2.27016   -1.48817    0.456745   2.18015    0.599337
 -0.412719   0.456745   4.30436   -2.37771    0.431586
  0.331858   2.18015   -2.37771    0.966284   0.971221
  3.62299    0.599337   0.431586   0.971221  -1.9591

So I'm inclined not to bother atm

[jishnub]

It appears that the fzeros approach isn't sufficient here, as

julia> typeof(s)
Symmetric{Float64, Matrix{Float64}}

julia> Broadcast.BroadcastStyle(typeof(s))
Base.Broadcast.DefaultArrayStyle{2}()

Perhaps Symmetric should define its own broadcast style in LinearAlgebra