Generalize 3-arg dot to HermOrSym (#1410)

Author: dkarrasch

src/symmetric.jl

@@ -230,6 +230,9 @@ const SelfAdjoint = Union{SymTridiagonal{<:Real}, Symmetric{<:Real}, Hermitian}
 wrappertype(::Union{Symmetric, SymTridiagonal}) = Symmetric
 wrappertype(::Hermitian) = Hermitian
 
+hswrapperop(::Symmetric) = symmetric
+hswrapperop(::Hermitian) = hermitian
+
 nonhermitianwrappertype(::SymSymTri{<:Real}) = Symmetric
 nonhermitianwrappertype(::Hermitian{<:Real}) = Symmetric
 nonhermitianwrappertype(::Hermitian) = identity
@@ -738,27 +741,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())
     data = A.data
-    r = dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
+    s = dot(first(x), first(A), first(y))
+    r = zero(s+s)
     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], hswrapperop(A)(data[j,j], :U), 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], hswrapperop(A)(data[j,j], :L), 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, M*y) ≈ dot(x, N, y)
+    end
 end
 
 # let's make sure the analogous bug will not show up with kronecker products
@@ -610,6 +618,18 @@ end
     end
 end
 
+@testset "3-arg dot with Symmetric/Hermitian matrix of matrices" begin
+    for m in (Symmetric([randn(ComplexF64, 2, 2) for i in 1:2, j in 1:2]),
+                Symmetric([randn(ComplexF64, 2, 2) for i in 1:2, j in 1:2], :L),
+                Hermitian([randn(ComplexF64, 2, 2) for i in 1:2, j in 1:2]),
+                Hermitian([randn(ComplexF64, 2, 2) for i in 1:2, j in 1:2], :L)
+            )
+        x = [randn(ComplexF64, 2) for i in 1:2]
+        y = [randn(ComplexF64, 2) for i in 1:2]
+        @test dot(x, m, y) ≈ dot(x, m*y) ≈ dot(x, Matrix(m), y)
+    end
+end
+
 #Issue #7647: test xsyevr, xheevr, xstevr drivers.
 @testset "Eigenvalues in interval for $(typeof(Mi7647))" for Mi7647 in
         (Symmetric(diagm(0 => 1.0:3.0)),