diff --git a/src/symmetric.jl b/src/symmetric.jl index ed5fe3b5..6d928dd1 100644 --- a/src/symmetric.jl +++ b/src/symmetric.jl @@ -224,11 +224,15 @@ const RealHermSymSymTri{T<:Real} = Union{RealHermSym{T}, SymTridiagonal{T}} const RealHermSymComplexHerm{T<:Real,S} = Union{Hermitian{T,S}, Symmetric{T,S}, Hermitian{Complex{T},S}} const RealHermSymComplexSym{T<:Real,S} = Union{Hermitian{T,S}, Symmetric{T,S}, Symmetric{Complex{T},S}} const RealHermSymSymTriComplexHerm{T<:Real} = Union{RealHermSymComplexSym{T}, SymTridiagonal{T}} -const SelfAdjoint = Union{Symmetric{<:Real}, Hermitian{<:Number}} +const SelfAdjoint = Union{SymTridiagonal{<:Real}, Symmetric{<:Real}, Hermitian} wrappertype(::Union{Symmetric, SymTridiagonal}) = Symmetric wrappertype(::Hermitian) = Hermitian +nonhermitianwrappertype(::SymSymTri{<:Real}) = Symmetric +nonhermitianwrappertype(::Hermitian{<:Real}) = Symmetric +nonhermitianwrappertype(::Hermitian) = identity + size(A::HermOrSym) = size(A.data) axes(A::HermOrSym) = axes(A.data) @inline function Base.isassigned(A::HermOrSym, i::Int, j::Int) @@ -834,119 +838,74 @@ end ^(A::Symmetric{<:Complex}, p::Integer) = sympow(A, p) ^(A::SymTridiagonal{<:Real}, p::Integer) = sympow(A, p) ^(A::SymTridiagonal{<:Complex}, p::Integer) = sympow(A, p) -function sympow(A::SymSymTri, p::Integer) - if p < 0 - return Symmetric(Base.power_by_squaring(inv(A), -p)) - else - return Symmetric(Base.power_by_squaring(A, p)) - end -end -for hermtype in (:Symmetric, :SymTridiagonal) - @eval begin - function ^(A::$hermtype{<:Real}, p::Real) - isinteger(p) && return integerpow(A, p) - F = eigen(A) - if all(λ -> λ ≥ 0, F.values) - return Symmetric((F.vectors * Diagonal((F.values).^p)) * F.vectors') - else - return Symmetric((F.vectors * Diagonal(complex.(F.values).^p)) * F.vectors') - end - end - function ^(A::$hermtype{<:Complex}, p::Real) - isinteger(p) && return integerpow(A, p) - return Symmetric(schurpow(A, p)) - end - end -end -function ^(A::Hermitian, p::Integer) +^(A::Hermitian, p::Integer) = sympow(A, p) +function sympow(A, p::Integer) if p < 0 retmat = Base.power_by_squaring(inv(A), -p) else retmat = Base.power_by_squaring(A, p) end - return Hermitian(retmat) + return wrappertype(A)(retmat) end -function ^(A::Hermitian{T}, p::Real) where T +function ^(A::SelfAdjoint, p::Real) isinteger(p) && return integerpow(A, p) F = eigen(A) if all(λ -> λ ≥ 0, F.values) retmat = (F.vectors * Diagonal((F.values).^p)) * F.vectors' - return Hermitian(retmat) + return wrappertype(A)(retmat) else - retmat = (F.vectors * Diagonal((complex.(F.values).^p))) * F.vectors' - if T <: Real - return Symmetric(retmat) - else - return retmat - end + retmat = (F.vectors * Diagonal(complex.(F.values).^p)) * F.vectors' + return nonhermitianwrappertype(A)(retmat) end end +function ^(A::SymSymTri{<:Complex}, p::Real) + isinteger(p) && return integerpow(A, p) + return Symmetric(schurpow(A, p)) +end for func in (:exp, :cos, :sin, :tan, :cosh, :sinh, :tanh, :atan, :asinh, :atanh, :cbrt) @eval begin - function ($func)(A::RealHermSymSymTri) - F = eigen(A) - return wrappertype(A)((F.vectors * Diagonal(($func).(F.values))) * F.vectors') - end - function ($func)(A::Hermitian{<:Complex}) + function ($func)(A::SelfAdjoint) F = eigen(A) retmat = (F.vectors * Diagonal(($func).(F.values))) * F.vectors' - return Hermitian(retmat) + return wrappertype(A)(retmat) end end end -function cis(A::RealHermSymSymTri) +function cis(A::SelfAdjoint) F = eigen(A) - return Symmetric(F.vectors .* cis.(F.values') * F.vectors') + retmat = F.vectors .* cis.(F.values') * F.vectors' + return nonhermitianwrappertype(A)(retmat) end -function cis(A::Hermitian{<:Complex}) - F = eigen(A) - return F.vectors .* cis.(F.values') * F.vectors' -end - for func in (:acos, :asin) @eval begin - function ($func)(A::RealHermSymSymTri) - F = eigen(A) - if all(λ -> -1 ≤ λ ≤ 1, F.values) - return wrappertype(A)((F.vectors * Diagonal(($func).(F.values))) * F.vectors') - else - return Symmetric((F.vectors * Diagonal(($func).(complex.(F.values)))) * F.vectors') - end - end - function ($func)(A::Hermitian{<:Complex}) + function ($func)(A::SelfAdjoint) F = eigen(A) if all(λ -> -1 ≤ λ ≤ 1, F.values) retmat = (F.vectors * Diagonal(($func).(F.values))) * F.vectors' - return Hermitian(retmat) + return wrappertype(A)(retmat) else - return (F.vectors * Diagonal(($func).(complex.(F.values)))) * F.vectors' + retmat = (F.vectors * Diagonal(($func).(complex.(F.values)))) * F.vectors' + return nonhermitianwrappertype(A)(retmat) end end end end -function acosh(A::RealHermSymSymTri) - F = eigen(A) - if all(λ -> λ ≥ 1, F.values) - return wrappertype(A)((F.vectors * Diagonal(acosh.(F.values))) * F.vectors') - else - return Symmetric((F.vectors * Diagonal(acosh.(complex.(F.values)))) * F.vectors') - end -end -function acosh(A::Hermitian{<:Complex}) +function acosh(A::SelfAdjoint) F = eigen(A) if all(λ -> λ ≥ 1, F.values) retmat = (F.vectors * Diagonal(acosh.(F.values))) * F.vectors' - return Hermitian(retmat) + return wrappertype(A)(retmat) else - return (F.vectors * Diagonal(acosh.(complex.(F.values)))) * F.vectors' + retmat = (F.vectors * Diagonal(acosh.(complex.(F.values)))) * F.vectors' + return nonhermitianwrappertype(A)(retmat) end end -function sincos(A::RealHermSymSymTri) +function sincos(A::SelfAdjoint) n = checksquare(A) F = eigen(A) T = float(eltype(F.values)) @@ -956,49 +915,28 @@ function sincos(A::RealHermSymSymTri) end return wrappertype(A)((F.vectors * S) * F.vectors'), wrappertype(A)((F.vectors * C) * F.vectors') end -function sincos(A::Hermitian{<:Complex}) - n = checksquare(A) + +function log(A::SelfAdjoint) F = eigen(A) - T = float(eltype(F.values)) - S, C = Diagonal(similar(A, T, (n,))), Diagonal(similar(A, T, (n,))) - for i in eachindex(S.diag, C.diag, F.values) - S.diag[i], C.diag[i] = sincos(F.values[i]) - end - retmatS, retmatC = (F.vectors * S) * F.vectors', (F.vectors * C) * F.vectors' - for i in diagind(retmatS, IndexStyle(retmatS)) - retmatS[i] = real(retmatS[i]) - retmatC[i] = real(retmatC[i]) + if all(λ -> λ ≥ 0, F.values) + retmat = (F.vectors * Diagonal(log.(F.values))) * F.vectors' + return wrappertype(A)(retmat) + else + retmat = (F.vectors * Diagonal(log.(complex.(F.values)))) * F.vectors' + return nonhermitianwrappertype(A)(retmat) end - return Hermitian(retmatS), Hermitian(retmatC) end - -for func in (:log, :sqrt) - # sqrt has rtol arg to handle matrices that are semidefinite up to roundoff errors - rtolarg = func === :sqrt ? Any[Expr(:kw, :(rtol::Real), :(eps(real(float(one(T))))*size(A,1)))] : Any[] - rtolval = func === :sqrt ? :(-maximum(abs, F.values) * rtol) : 0 - @eval begin - function ($func)(A::RealHermSymSymTri{T}; $(rtolarg...)) where {T<:Real} - F = eigen(A) - λ₀ = $rtolval # treat λ ≥ λ₀ as "zero" eigenvalues up to roundoff - if all(λ -> λ ≥ λ₀, F.values) - return wrappertype(A)((F.vectors * Diagonal(($func).(max.(0, F.values)))) * F.vectors') - else - return Symmetric((F.vectors * Diagonal(($func).(complex.(F.values)))) * F.vectors') - end - end - function ($func)(A::Hermitian{T}; $(rtolarg...)) where {T<:Complex} - n = checksquare(A) - F = eigen(A) - λ₀ = $rtolval # treat λ ≥ λ₀ as "zero" eigenvalues up to roundoff - if all(λ -> λ ≥ λ₀, F.values) - retmat = (F.vectors * Diagonal(($func).(max.(0, F.values)))) * F.vectors' - return Hermitian(retmat) - else - retmat = (F.vectors * Diagonal(($func).(complex.(F.values)))) * F.vectors' - return retmat - end - end +# sqrt has rtol kwarg to handle matrices that are semidefinite up to roundoff errors +function sqrt(A::SelfAdjoint; rtol = eps(real(float(eltype(A)))) * size(A, 1)) + F = eigen(A) + λ₀ = -maximum(abs, F.values) * rtol # treat λ ≥ λ₀ as "zero" eigenvalues up to roundoff + if all(λ -> λ ≥ λ₀, F.values) + retmat = (F.vectors * Diagonal(sqrt.(max.(0, F.values)))) * F.vectors' + return wrappertype(A)(retmat) + else + retmat = (F.vectors * Diagonal(sqrt.(complex.(F.values)))) * F.vectors' + return nonhermitianwrappertype(A)(retmat) end end diff --git a/test/symmetric.jl b/test/symmetric.jl index 9f727b8c..3f1397c0 100644 --- a/test/symmetric.jl +++ b/test/symmetric.jl @@ -1199,4 +1199,14 @@ end end end +@testset "asin/acos/acosh for matrix outside the real domain" begin + M = [0 2;2 0] #eigenvalues are ±2 + for T ∈ (Float32, Float64, ComplexF32, ComplexF64) + M2 = Hermitian(T.(M)) + @test sin(asin(M2)) ≈ M2 + @test cos(acos(M2)) ≈ M2 + @test cosh(acosh(M2)) ≈ M2 + end +end + end # module TestSymmetric