diff --git a/src/symmetric.jl b/src/symmetric.jl index 51a1d94b..14f4ab49 100644 --- a/src/symmetric.jl +++ b/src/symmetric.jl @@ -229,6 +229,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 @@ -737,27 +740,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 diff --git a/test/symmetric.jl b/test/symmetric.jl index 50c72e04..707b392d 100644 --- a/test/symmetric.jl +++ b/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)),