Generalize 3-arg dot to HermOrSym

Author: dkarrasch

src/symmetric.jl

@@ -737,27 +737,28 @@ function mul(A::AdjOrTrans{<:BlasFloat,<:StridedMatrix}, B::HermOrSym{<:BlasFloa
             convert(AbstractMatrix{T}, B))
 end
 
-function dot(x::AbstractVector, A::RealHermSymComplexHerm, y::AbstractVector)
+function dot(x::AbstractVector, A::HermOrSym, y::AbstractVector)
     require_one_based_indexing(x, y)
     n = length(x)
     (n == length(y) == size(A, 1)) || throw(DimensionMismatch())
+    diagop = A isa Symmetric ? identity : real
     data = A.data
     r = dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
     iszero(n) && return r
     if A.uplo == 'U'
         @inbounds for j = 1:length(y)
-            r += dot(x[j], real(data[j,j]), y[j])
+            r += dot(x[j], diagop(data[j,j]), y[j])
             @simd for i = 1:j-1
                 Aij = data[i,j]
-                r += dot(x[i], Aij, y[j]) + dot(x[j], adjoint(Aij), y[i])
+                r += dot(x[i], Aij, y[j]) + dot(x[j], _conjugation(A)(Aij), y[i])
             end
         end
     else # A.uplo == 'L'
         @inbounds for j = 1:length(y)
-            r += dot(x[j], real(data[j,j]), y[j])
+            r += dot(x[j], diagop(data[j,j]), y[j])
             @simd for i = j+1:length(y)
                 Aij = data[i,j]
-                r += dot(x[i], Aij, y[j]) + dot(x[j], adjoint(Aij), y[i])
+                r += dot(x[i], Aij, y[j]) + dot(x[j], _conjugation(A)(Aij), y[i])
             end
         end
     end

test/symmetric.jl

@@ -582,13 +582,21 @@ end
 
 # bug identified in PR #52318: dot products of quaternionic Hermitian matrices,
 # or any number type where conj(a)*conj(b) ≠ conj(a*b):
-@testset "dot Hermitian quaternion #52318" begin
-    A, B = [Quaternion.(randn(3,3), randn(3, 3), randn(3, 3), randn(3,3)) |> t -> t + t' for i in 1:2]
+@testset "dot Hermitian quaternion" begin
+    A, B = [(randn(Quaternion{Float64},4,4)) |> t -> t + t' for i in 1:2]
     @test A == Hermitian(A) && B == Hermitian(B)
     @test dot(A, B) ≈ dot(Hermitian(A), Hermitian(B))
-    A, B = [Quaternion.(randn(3,3), randn(3, 3), randn(3, 3), randn(3,3)) |> t -> t + transpose(t) for i in 1:2]
+    A, B = [(randn(Quaternion{Float64},4,4)) |> t -> t + transpose(t) for i in 1:2]
     @test A == Symmetric(A) && B == Symmetric(B)
     @test dot(A, B) ≈ dot(Symmetric(A), Symmetric(B))
+    A = randn(Quaternion{Float64}, 4, 4)
+    x = randn(Quaternion{Float64}, 4)
+    y = randn(Quaternion{Float64}, 4)
+    for t in (Symmetric, Hermitian), uplo in (:U, :L)
+        M = t(A, uplo)
+        N = Matrix(M)
+        @test dot(x, M, y) ≈ dot(x, N, y)
+    end
 end
 
 # let's make sure the analogous bug will not show up with kronecker products