diff --git a/src/diagonal.jl b/src/diagonal.jl index e60fb009..63cff1c8 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -970,6 +970,10 @@ function inv(D::Diagonal{T}) where T Diagonal(Di) end +# Ensure doubly wrapped matrices use efficient diagonal methods and return a Symmetric/Hermitian type +inv(A::Symmetric{<:Number,<:Diagonal}) = Symmetric(inv(A.data), sym_uplo(A.uplo)) +inv(A::Hermitian{<:Number,<:Diagonal}) = Hermitian(inv(real(A.data)), sym_uplo(A.uplo)) + function pinv(D::Diagonal{T}) where T Di = similar(D.diag, typeof(inv(oneunit(T)))) for i = 1:length(D.diag) diff --git a/test/symmetric.jl b/test/symmetric.jl index 707b392d..ae10618f 100644 --- a/test/symmetric.jl +++ b/test/symmetric.jl @@ -874,6 +874,21 @@ end @test det(Hermitian(A))::Float64 == det(A) == 0.0 end +@testset "issue #1437: inverse of Symmetric|Hermitian{<:Any,<:Diagonal} returns of Symmetric|Hermitian{<:Any,<:Diagonal}" begin + Dreal = Diagonal(randn(3)) + Dcomplex = Diagonal(randn(ComplexF64, 3)) + # without wrapper + invDreal = inv(Dreal) + invDcomplex = inv(real(Dcomplex)) # because Hermitian implies a real diagonal + # with wrapper + SDreal = Symmetric(Dreal) + HDcomplex = Hermitian(Dcomplex) + @test inv(SDreal)::Symmetric{Float64,typeof(Dreal)} ≈ invDreal + @test inv(HDcomplex)::Hermitian{Float64,typeof(Dreal)} ≈ invDcomplex + Dcomplex[2,2] = 0 + @test_throws SingularException inv(HDcomplex) +end + @testset "symmetric()/hermitian() for Numbers" begin @test LinearAlgebra.symmetric(1) == LinearAlgebra.symmetric(1, :U) == 1 @test LinearAlgebra.symmetric_type(Int) == Int