make Diagonal+Symmetric return Symmetric (#35333)

* make Diagonal+Symmetric return Symmetric

* introduce generic +/- AbstractMatrix with Diagonal

* fix typos in tests

* remove generic addition/subtraction with Diagonals

stdlib/LinearAlgebra/src/diagonal.jl

@@ -156,6 +156,21 @@ end
 (+)(Da::Diagonal, Db::Diagonal) = Diagonal(Da.diag + Db.diag)
 (-)(Da::Diagonal, Db::Diagonal) = Diagonal(Da.diag - Db.diag)
 
+for f in (:+, :-)
+    @eval function $f(D::Diagonal, S::Symmetric)
+        return Symmetric($f(D, S.data), sym_uplo(S.uplo))
+    end
+    @eval function $f(S::Symmetric, D::Diagonal)
+        return Symmetric($f(S.data, D), sym_uplo(S.uplo))
+    end
+    @eval function $f(D::Diagonal{<:Real}, H::Hermitian)
+        return Hermitian($f(D, H.data), sym_uplo(H.uplo))
+    end
+    @eval function $f(H::Hermitian, D::Diagonal{<:Real})
+        return Hermitian($f(H.data, D), sym_uplo(H.uplo))
+    end
+end
+
 (*)(x::Number, D::Diagonal) = Diagonal(x * D.diag)
 (*)(D::Diagonal, x::Number) = Diagonal(D.diag * x)
 (/)(D::Diagonal, x::Number) = Diagonal(D.diag / x)

stdlib/LinearAlgebra/test/diagonal.jl

@@ -200,6 +200,19 @@ Random.seed!(1)
         A     = rand(elty, n, n)
         Asym  = Symmetric(A + transpose(A), :U)
         Aherm = Hermitian(A + adjoint(A), :U)
+        for op in (+, -)
+            @test op(Asym, D) isa Symmetric
+            @test Array(op(Asym, D)) ≈ Array(Symmetric(op(Array(Asym), Array(D))))
+            @test op(D, Asym) isa Symmetric
+            @test Array(op(D, Asym)) ≈ Array(Symmetric(op(Array(D), Array(Asym))))
+            if !(elty <: Real)
+                Dr = real(D)
+                @test op(Aherm, Dr) isa Hermitian
+                @test Array(op(Aherm, Dr)) ≈ Array(Hermitian(op(Array(Aherm), Array(Dr))))
+                @test op(Dr, Aherm) isa Hermitian
+                @test Array(op(Dr, Aherm)) ≈ Array(Hermitian(op(Array(Dr), Array(Aherm))))
+            end
+        end
         @test Array(D*Transpose(Asym)) ≈ Array(D) * Array(transpose(Asym))
         @test Array(D*Adjoint(Asym)) ≈ Array(D) * Array(adjoint(Asym))
         @test Array(D*Transpose(Aherm)) ≈ Array(D) * Array(transpose(Aherm))