diff --git a/src/symmetric.jl b/src/symmetric.jl index 6d928dd1..ed5fe3b5 100644 --- a/src/symmetric.jl +++ b/src/symmetric.jl @@ -224,15 +224,11 @@ 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{SymTridiagonal{<:Real}, Symmetric{<:Real}, Hermitian} +const SelfAdjoint = Union{Symmetric{<:Real}, Hermitian{<:Number}} 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) @@ -838,74 +834,119 @@ 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) -^(A::Hermitian, p::Integer) = sympow(A, p) -function sympow(A, p::Integer) +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) if p < 0 retmat = Base.power_by_squaring(inv(A), -p) else retmat = Base.power_by_squaring(A, p) end - return wrappertype(A)(retmat) + return Hermitian(retmat) end -function ^(A::SelfAdjoint, p::Real) +function ^(A::Hermitian{T}, p::Real) where T isinteger(p) && return integerpow(A, p) F = eigen(A) if all(λ -> λ ≥ 0, F.values) retmat = (F.vectors * Diagonal((F.values).^p)) * F.vectors' - return wrappertype(A)(retmat) + return Hermitian(retmat) else - retmat = (F.vectors * Diagonal(complex.(F.values).^p)) * F.vectors' - return nonhermitianwrappertype(A)(retmat) + retmat = (F.vectors * Diagonal((complex.(F.values).^p))) * F.vectors' + if T <: Real + return Symmetric(retmat) + else + return retmat + end 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::SelfAdjoint) + function ($func)(A::RealHermSymSymTri) + F = eigen(A) + return wrappertype(A)((F.vectors * Diagonal(($func).(F.values))) * F.vectors') + end + function ($func)(A::Hermitian{<:Complex}) F = eigen(A) retmat = (F.vectors * Diagonal(($func).(F.values))) * F.vectors' - return wrappertype(A)(retmat) + return Hermitian(retmat) end end end -function cis(A::SelfAdjoint) +function cis(A::RealHermSymSymTri) F = eigen(A) - retmat = F.vectors .* cis.(F.values') * F.vectors' - return nonhermitianwrappertype(A)(retmat) + return Symmetric(F.vectors .* cis.(F.values') * F.vectors') 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::SelfAdjoint) + 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}) F = eigen(A) if all(λ -> -1 ≤ λ ≤ 1, F.values) retmat = (F.vectors * Diagonal(($func).(F.values))) * F.vectors' - return wrappertype(A)(retmat) + return Hermitian(retmat) else - retmat = (F.vectors * Diagonal(($func).(complex.(F.values)))) * F.vectors' - return nonhermitianwrappertype(A)(retmat) + return (F.vectors * Diagonal(($func).(complex.(F.values)))) * F.vectors' end end end end -function acosh(A::SelfAdjoint) +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}) F = eigen(A) if all(λ -> λ ≥ 1, F.values) retmat = (F.vectors * Diagonal(acosh.(F.values))) * F.vectors' - return wrappertype(A)(retmat) + return Hermitian(retmat) else - retmat = (F.vectors * Diagonal(acosh.(complex.(F.values)))) * F.vectors' - return nonhermitianwrappertype(A)(retmat) + return (F.vectors * Diagonal(acosh.(complex.(F.values)))) * F.vectors' end end -function sincos(A::SelfAdjoint) +function sincos(A::RealHermSymSymTri) n = checksquare(A) F = eigen(A) T = float(eltype(F.values)) @@ -915,28 +956,49 @@ function sincos(A::SelfAdjoint) end return wrappertype(A)((F.vectors * S) * F.vectors'), wrappertype(A)((F.vectors * C) * F.vectors') end - -function log(A::SelfAdjoint) +function sincos(A::Hermitian{<:Complex}) + n = checksquare(A) F = eigen(A) - 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) + 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]) end + return Hermitian(retmatS), Hermitian(retmatC) 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) + +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 end end diff --git a/test/symmetric.jl b/test/symmetric.jl index 3f1397c0..9f727b8c 100644 --- a/test/symmetric.jl +++ b/test/symmetric.jl @@ -1199,14 +1199,4 @@ 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