This seems to cut down the runtime of `D .* D` by half for `5000x5000` diagonal matrices.
Using nightly `v"1.10.0-DEV.450"`
```julia
julia> using LinearAlgebra, BenchmarkTools

julia> D = Diagonal(rand(5000));

julia> @Btime $D * $D;
  3.648 μs (2 allocations: 39.11 KiB)

julia> @Btime $D .* $D;
  8.729 μs (2 allocations: 39.11 KiB)

julia> B = Broadcast.instantiate(Broadcast.Broadcasted(*, (D, D)));

julia> @Btime Base.copyto!($(copy(D)), $B);
  8.390 μs (0 allocations: 0 bytes)

julia> function copyto2!(dest::Diagonal, bc::Broadcast.Broadcasted{<:LinearAlgebra.StructuredMatrixStyle})
           !LinearAlgebra.isstructurepreserving(bc) && !LinearAlgebra.fzeropreserving(bc) && return copyto!(dest, convert(Broadcast.Broadcasted{Nothing}, bc))
           axs = axes(dest)
           axes(bc) == axs || Broadcast.throwdm(axes(bc), axs)
           @inbounds for i in axs[1]
               dest.diag[i] = Broadcast._broadcast_getindex(bc, CartesianIndex(i, i))
           end
           return dest
       end
copyto2! (generic function with 1 method)

julia> @Btime copyto2!($(copy(D)), $B);
  4.207 μs (0 allocations: 0 bytes)
```