Add dot method for symmetric matrices

stdlib/LinearAlgebra/src/symmetric.jl

@@ -414,6 +414,84 @@ function triu(A::Symmetric, k::Integer=0)
     end
 end
 
+function dot(A::Symmetric, B::Symmetric)
+    n = size(A, 2)
+    if n != size(B, 2)
+        throw(DimensionMismatch("A has dimensions $(size(A)) but B has dimensions $(size(B))"))
+    end
+
+    dotprod = zero(dot(first(A), first(B)))
+    @inbounds if A.uplo == 'U' && B.uplo == 'U'
+        for j in 1:n
+            for i in 1:(j - 1)
+                dotprod += 2 * dot(A.data[i, j], B.data[i, j])
+            end
+            dotprod += dot(A[j, j], B[j, j])
+        end
+    elseif A.uplo == 'L' && B.uplo == 'L'
+        for j in 1:n
+            dotprod += dot(A[j, j], B[j, j])
+            for i in (j + 1):n
+                dotprod += 2 * dot(A.data[i, j], B.data[i, j])
+            end
+        end
+    elseif A.uplo == 'U' && B.uplo == 'L'
+        for j in 1:n
+            for i in 1:(j - 1)
+                dotprod += 2 * dot(A.data[i, j], transpose(B.data[j, i]))
+            end
+            dotprod += dot(A[j, j], B[j, j])
+        end
+    else
+        for j in 1:n
+            dotprod += dot(A[j, j], B[j, j])
+            for i in (j + 1):n
+                dotprod += 2 * dot(A.data[i, j], transpose(B.data[j, i]))
+            end
+        end
+    end
+    return dotprod
+end
+
+function dot(A::Hermitian, B::Hermitian)
+    n = size(A, 2)
+    if n != size(B, 2)
+        throw(DimensionMismatch("A has dimensions $(size(A)) but B has dimensions $(size(B))"))
+    end
+
+    dotprod = zero(dot(first(A), first(B)))
+    @inbounds if A.uplo == 'U' && B.uplo == 'U'
+        for j in 1:n
+            for i in 1:(j - 1)
+                dotprod += 2 * real(dot(A.data[i, j], B.data[i, j]))
+            end
+            dotprod += dot(A[j, j], B[j, j])
+        end
+    elseif A.uplo == 'L' && B.uplo == 'L'
+        for j in 1:n
+            dotprod += dot(A[j, j], B[j, j])
+            for i in (j + 1):n
+                dotprod += 2 * real(dot(A.data[i, j], B.data[i, j]))
+            end
+        end
+    elseif A.uplo == 'U' && B.uplo == 'L'
+        for j in 1:n
+            for i in 1:(j - 1)
+                dotprod += 2 * real(dot(A.data[i, j], adjoint(B.data[j, i])))
+            end
+            dotprod += dot(A[j, j], B[j, j])
+        end
+    else
+        for j in 1:n
+            dotprod += dot(A[j, j], B[j, j])
+            for i in (j + 1):n
+                dotprod += 2 * real(dot(A.data[i, j], adjoint(B.data[j, i])))
+            end
+        end
+    end
+    return dotprod
+end
+
 (-)(A::Symmetric) = Symmetric(-A.data, sym_uplo(A.uplo))
 (-)(A::Hermitian) = Hermitian(-A.data, sym_uplo(A.uplo))
 

stdlib/LinearAlgebra/test/symmetric.jl

@@ -377,6 +377,40 @@ end
                 end
             end
         end
+
+        @testset "dot product of symmetric and Hermitian matrices" begin
+            for mtype in (Symmetric, Hermitian)
+                symau = mtype(a, :U)
+                symal = mtype(a, :L)
+                msymau = Matrix(symau)
+                msymal = Matrix(symal)
+                @test_throws DimensionMismatch dot(symau, mtype(zeros(eltya, n-1, n-1)))
+                for eltyc in (Float32, Float64, ComplexF32, ComplexF64, BigFloat, Int)
+                    creal = randn(n, n)/2
+                    cimag = randn(n, n)/2
+                    c = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(creal, cimag) : creal)
+                    symcu = mtype(c, :U)
+                    symcl = mtype(c, :L)
+                    msymcu = Matrix(symcu)
+                    msymcl = Matrix(symcl)
+                    @test dot(symau, symcu) ≈ dot(msymau, msymcu)
+                    @test dot(symau, symcl) ≈ dot(msymau, msymcl)
+                    @test dot(symal, symcu) ≈ dot(msymal, msymcu)
+                    @test dot(symal, symcl) ≈ dot(msymal, msymcl)
+                end
+
+                # block matrices
+                blockm = [eltya == Int ? rand(1:7, 3, 3) : convert(Matrix{eltya}, eltya <: Complex ? complex.(randn(3, 3)/2, randn(3, 3)/2) : randn(3, 3)/2) for _ in 1:3, _ in 1:3]
+                symblockmu = mtype(blockm, :U)
+                symblockml = mtype(blockm, :L)
+                msymblockmu = Matrix(symblockmu)
+                msymblockml = Matrix(symblockml)
+                @test dot(symblockmu, symblockmu) ≈ dot(msymblockmu, msymblockmu)
+                @test dot(symblockmu, symblockml) ≈ dot(msymblockmu, msymblockml)
+                @test dot(symblockml, symblockmu) ≈ dot(msymblockml, msymblockmu)
+                @test dot(symblockml, symblockml) ≈ dot(msymblockml, msymblockml)
+            end
+        end
     end
 end