diff --git a/src/LinearAlgebra.jl b/src/LinearAlgebra.jl index 661672cd..8b190e39 100644 --- a/src/LinearAlgebra.jl +++ b/src/LinearAlgebra.jl @@ -364,7 +364,14 @@ function char_uplo(uplo::Symbol) end end +""" + sym_uplo(uplo::Char) + +Return the `Symbol` corresponding the `uplo` by checking for validity. +""" function sym_uplo(uplo::Char) + # This method is called by other packages, and isn't used within LinearAlgebra + # It's retained here for backward compatibility. if uplo == 'U' return :U elseif uplo == 'L' @@ -373,6 +380,13 @@ function sym_uplo(uplo::Char) throw_uplo() end end +""" + _sym_uplo(uplo::Char) + +Return the `Symbol` corresponding to `uplo` without checking for validity. +See also `sym_uplo`, which checks for validity. +""" +_sym_uplo(uplo::Char) = uplo == 'U' ? (:U) : (:L) @noinline throw_uplo() = throw(ArgumentError("uplo argument must be either :U (upper) or :L (lower)")) diff --git a/src/bidiag.jl b/src/bidiag.jl index ed50a2dd..88b5a763 100644 --- a/src/bidiag.jl +++ b/src/bidiag.jl @@ -302,7 +302,7 @@ function show(io::IO, M::Bidiagonal) print(io, ", ") show(io, M.ev) print(io, ", ") - show(io, sym_uplo(M.uplo)) + show(io, _sym_uplo(M.uplo)) print(io, ")") end diff --git a/src/bunchkaufman.jl b/src/bunchkaufman.jl index 7efe57eb..69811fde 100644 --- a/src/bunchkaufman.jl +++ b/src/bunchkaufman.jl @@ -130,7 +130,9 @@ function bunchkaufman!(A::StridedMatrix{<:BlasFloat}, rook::Bool = false; check: end bkcopy_oftype(A, S) = eigencopy_oftype(A, S) -bkcopy_oftype(A::Symmetric{<:Complex}, S) = Symmetric(copytrito!(similar(parent(A), S, size(A)), A.data, A.uplo), sym_uplo(A.uplo)) +function bkcopy_oftype(A::Symmetric{<:Complex}, S) + Symmetric(copytrito!(similar(parent(A), S, size(A)), A.data, A.uplo), _sym_uplo(A.uplo)) +end """ bunchkaufman(A, rook::Bool=false; check = true) -> S::BunchKaufman diff --git a/src/special.jl b/src/special.jl index 83880ca5..9dff10a4 100644 --- a/src/special.jl +++ b/src/special.jl @@ -292,19 +292,19 @@ end for f in (:+, :-) @eval function $f(D::Diagonal{<:Number}, S::Symmetric) - uplo = sym_uplo(S.uplo) + uplo = _sym_uplo(S.uplo) return Symmetric(parentof_applytri($f, Symmetric(D, uplo), S), uplo) end @eval function $f(S::Symmetric, D::Diagonal{<:Number}) - uplo = sym_uplo(S.uplo) + uplo = _sym_uplo(S.uplo) return Symmetric(parentof_applytri($f, S, Symmetric(D, uplo)), uplo) end @eval function $f(D::Diagonal{<:Real}, H::Hermitian) - uplo = sym_uplo(H.uplo) + uplo = _sym_uplo(H.uplo) return Hermitian(parentof_applytri($f, Hermitian(D, uplo), H), uplo) end @eval function $f(H::Hermitian, D::Diagonal{<:Real}) - uplo = sym_uplo(H.uplo) + uplo = _sym_uplo(H.uplo) return Hermitian(parentof_applytri($f, H, Hermitian(D, uplo)), uplo) end end @@ -608,8 +608,8 @@ end # tridiagonal cholesky factorization function cholesky(S::RealSymHermitian{<:BiTriSym}, ::NoPivot = NoPivot(); check::Bool = true) T = choltype(S) - B = Bidiagonal{T}(diag(S, 0), diag(S, S.uplo == 'U' ? 1 : -1), sym_uplo(S.uplo)) - cholesky!(Hermitian(B, sym_uplo(S.uplo)), NoPivot(); check = check) + B = Bidiagonal{T}(diag(S, 0), diag(S, S.uplo == 'U' ? 1 : -1), _sym_uplo(S.uplo)) + cholesky!(Hermitian(B, _sym_uplo(S.uplo)), NoPivot(); check = check) end # istriu/istril for triangular wrappers of structured matrices diff --git a/src/symmetric.jl b/src/symmetric.jl index 089eefa2..c04c54ff 100644 --- a/src/symmetric.jl +++ b/src/symmetric.jl @@ -206,7 +206,7 @@ for (S, H) in ((:Symmetric, :Hermitian), (:Hermitian, :Symmetric)) throw(ArgumentError("Cannot construct $($S); uplo doesn't match")) end end - $S(A::$H) = $S(A, sym_uplo(A.uplo)) + $S(A::$H) = $S(A, _sym_uplo(A.uplo)) function $S(A::$H, uplo::Symbol) if A.uplo == char_uplo(uplo) if $H === Hermitian && !(eltype(A) <: Real) && @@ -214,7 +214,7 @@ for (S, H) in ((:Symmetric, :Hermitian), (:Hermitian, :Symmetric)) throw(ArgumentError("Cannot construct $($S)($($H))); diagonal contains complex values")) end - return $S(A.data, sym_uplo(A.uplo)) + return $S(A.data, _sym_uplo(A.uplo)) else throw(ArgumentError("Cannot construct $($S); uplo doesn't match")) end @@ -286,7 +286,7 @@ end @inline function getindex(A::Symmetric, i::Int, j::Int) @boundscheck checkbounds(A, i, j) @inbounds if i == j - return symmetric(A.data[i, j], sym_uplo(A.uplo))::symmetric_type(eltype(A.data)) + return symmetric(A.data[i, j], _sym_uplo(A.uplo))::symmetric_type(eltype(A.data)) elseif (A.uplo == 'U') == (i < j) return A.data[i, j] else @@ -296,7 +296,7 @@ end @inline function getindex(A::Hermitian, i::Int, j::Int) @boundscheck checkbounds(A, i, j) @inbounds if i == j - return hermitian(A.data[i, j], sym_uplo(A.uplo))::hermitian_type(eltype(A.data)) + return hermitian(A.data[i, j], _sym_uplo(A.uplo))::hermitian_type(eltype(A.data)) elseif (A.uplo == 'U') == (i < j) return A.data[i, j] else @@ -329,14 +329,14 @@ Base._reverse(A::Hermitian, ::Colon) = Hermitian(reverse(A.data), A.uplo == 'U' end Base.dataids(A::HermOrSym) = Base.dataids(parent(A)) -Base.unaliascopy(A::Hermitian) = Hermitian(Base.unaliascopy(parent(A)), sym_uplo(A.uplo)) -Base.unaliascopy(A::Symmetric) = Symmetric(Base.unaliascopy(parent(A)), sym_uplo(A.uplo)) +Base.unaliascopy(A::Hermitian) = Hermitian(Base.unaliascopy(parent(A)), _sym_uplo(A.uplo)) +Base.unaliascopy(A::Symmetric) = Symmetric(Base.unaliascopy(parent(A)), _sym_uplo(A.uplo)) _conjugation(::Union{Symmetric, Hermitian{<:Real}}) = transpose _conjugation(::Hermitian) = adjoint -diag(A::Symmetric) = symmetric.(diag(parent(A)), sym_uplo(A.uplo)) -diag(A::Hermitian) = hermitian.(diag(parent(A)), sym_uplo(A.uplo)) +diag(A::Symmetric) = symmetric.(diag(parent(A)), _sym_uplo(A.uplo)) +diag(A::Hermitian) = hermitian.(diag(parent(A)), _sym_uplo(A.uplo)) function applytri(f, A::HermOrSym) if A.uplo == 'U' @@ -374,15 +374,15 @@ similar(A::Union{Symmetric,Hermitian}, ::Type{T}, dims::Dims{N}) where {T,N} = s parent(A::HermOrSym) = A.data Symmetric{T,S}(A::Symmetric{T,S}) where {T,S<:AbstractMatrix{T}} = A Symmetric{T,S}(A::Symmetric) where {T,S<:AbstractMatrix{T}} = Symmetric{T,S}(convert(S,A.data),A.uplo) -AbstractMatrix{T}(A::Symmetric) where {T} = Symmetric(convert(AbstractMatrix{T}, A.data), sym_uplo(A.uplo)) +AbstractMatrix{T}(A::Symmetric) where {T} = Symmetric(convert(AbstractMatrix{T}, A.data), _sym_uplo(A.uplo)) AbstractMatrix{T}(A::Symmetric{T}) where {T} = copy(A) Hermitian{T,S}(A::Hermitian{T,S}) where {T,S<:AbstractMatrix{T}} = A Hermitian{T,S}(A::Hermitian) where {T,S<:AbstractMatrix{T}} = Hermitian{T,S}(convert(S,A.data),A.uplo) -AbstractMatrix{T}(A::Hermitian) where {T} = Hermitian(convert(AbstractMatrix{T}, A.data), sym_uplo(A.uplo)) +AbstractMatrix{T}(A::Hermitian) where {T} = Hermitian(convert(AbstractMatrix{T}, A.data), _sym_uplo(A.uplo)) AbstractMatrix{T}(A::Hermitian{T}) where {T} = copy(A) -copy(A::Symmetric) = (Symmetric(parentof_applytri(copy, A), sym_uplo(A.uplo))) -copy(A::Hermitian) = (Hermitian(parentof_applytri(copy, A), sym_uplo(A.uplo))) +copy(A::Symmetric) = (Symmetric(parentof_applytri(copy, A), _sym_uplo(A.uplo))) +copy(A::Hermitian) = (Hermitian(parentof_applytri(copy, A), _sym_uplo(A.uplo))) function copyto!(dest::Symmetric, src::Symmetric) if axes(dest) != axes(src) @@ -423,13 +423,13 @@ end end @inline function _symmetrize_diagonal!(B, A::Symmetric) for i = 1:size(A, 1) - B[i,i] = symmetric(A[i,i], sym_uplo(A.uplo))::symmetric_type(eltype(A.data)) + B[i,i] = symmetric(A[i,i], _sym_uplo(A.uplo))::symmetric_type(eltype(A.data)) end return B end @inline function _symmetrize_diagonal!(B, A::Hermitian) for i = 1:size(A, 1) - B[i,i] = hermitian(A[i,i], sym_uplo(A.uplo))::hermitian_type(eltype(A.data)) + B[i,i] = hermitian(A[i,i], _sym_uplo(A.uplo))::hermitian_type(eltype(A.data)) end return B end @@ -501,9 +501,9 @@ transpose(A::Hermitian{<:Real}) = A real(A::Symmetric{<:Real}) = A real(A::Hermitian{<:Real}) = A -real(A::Symmetric) = Symmetric(parentof_applytri(real, A), sym_uplo(A.uplo)) -real(A::Hermitian) = Hermitian(parentof_applytri(real, A), sym_uplo(A.uplo)) -imag(A::Symmetric) = Symmetric(parentof_applytri(imag, A), sym_uplo(A.uplo)) +real(A::Symmetric) = Symmetric(parentof_applytri(real, A), _sym_uplo(A.uplo)) +real(A::Hermitian) = Hermitian(parentof_applytri(real, A), _sym_uplo(A.uplo)) +imag(A::Symmetric) = Symmetric(parentof_applytri(imag, A), _sym_uplo(A.uplo)) Base.copy(A::Adjoint{<:Any,<:Symmetric}) = Symmetric(copy(adjoint(A.parent.data)), ifelse(A.parent.uplo == 'U', :L, :U)) @@ -513,8 +513,8 @@ Base.copy(A::Transpose{<:Any,<:Hermitian}) = tr(A::Symmetric{<:Number}) = tr(A.data) # to avoid AbstractMatrix fallback (incl. allocations) tr(A::Hermitian{<:Number}) = real(tr(A.data)) -Base.conj(A::Symmetric) = Symmetric(parentof_applytri(conj, A), sym_uplo(A.uplo)) -Base.conj(A::Hermitian) = Hermitian(parentof_applytri(conj, A), sym_uplo(A.uplo)) +Base.conj(A::Symmetric) = Symmetric(parentof_applytri(conj, A), _sym_uplo(A.uplo)) +Base.conj(A::Hermitian) = Hermitian(parentof_applytri(conj, A), _sym_uplo(A.uplo)) Base.conj!(A::HermOrSym) = typeof(A)(parentof_applytri(conj!, A), A.uplo) # tril/triu @@ -731,25 +731,25 @@ function _hermkron!(C, A, B, conj, real, Auplo, Buplo) end end -(-)(A::Symmetric) = Symmetric(parentof_applytri(-, A), sym_uplo(A.uplo)) -(-)(A::Hermitian) = Hermitian(parentof_applytri(-, A), sym_uplo(A.uplo)) +(-)(A::Symmetric) = Symmetric(parentof_applytri(-, A), _sym_uplo(A.uplo)) +(-)(A::Hermitian) = Hermitian(parentof_applytri(-, A), _sym_uplo(A.uplo)) ## Addition/subtraction for f ∈ (:+, :-), Wrapper ∈ (:Hermitian, :Symmetric) @eval function $f(A::$Wrapper, B::$Wrapper) - uplo = A.uplo == B.uplo ? sym_uplo(A.uplo) : (:U) + uplo = A.uplo == B.uplo ? _sym_uplo(A.uplo) : (:U) $Wrapper(parentof_applytri($f, A, B), uplo) end end for f in (:+, :-) @eval begin - $f(A::Hermitian, B::Symmetric{<:Real}) = $f(A, Hermitian(parent(B), sym_uplo(B.uplo))) - $f(A::Symmetric{<:Real}, B::Hermitian) = $f(Hermitian(parent(A), sym_uplo(A.uplo)), B) - $f(A::SymTridiagonal, B::Symmetric) = $f(Symmetric(A, sym_uplo(B.uplo)), B) - $f(A::Symmetric, B::SymTridiagonal) = $f(A, Symmetric(B, sym_uplo(A.uplo))) - $f(A::SymTridiagonal{<:Real}, B::Hermitian) = $f(Hermitian(A, sym_uplo(B.uplo)), B) - $f(A::Hermitian, B::SymTridiagonal{<:Real}) = $f(A, Hermitian(B, sym_uplo(A.uplo))) + $f(A::Hermitian, B::Symmetric{<:Real}) = $f(A, Hermitian(parent(B), _sym_uplo(B.uplo))) + $f(A::Symmetric{<:Real}, B::Hermitian) = $f(Hermitian(parent(A), _sym_uplo(A.uplo)), B) + $f(A::SymTridiagonal, B::Symmetric) = $f(Symmetric(A, _sym_uplo(B.uplo)), B) + $f(A::Symmetric, B::SymTridiagonal) = $f(A, Symmetric(B, _sym_uplo(A.uplo))) + $f(A::SymTridiagonal{<:Real}, B::Hermitian) = $f(Hermitian(A, _sym_uplo(B.uplo)), B) + $f(A::Hermitian, B::SymTridiagonal{<:Real}) = $f(A, Hermitian(B, _sym_uplo(A.uplo))) end end @@ -799,12 +799,12 @@ function dot(x::AbstractVector, A::HermOrSym, y::AbstractVector) end # Scaling with Number -*(A::Symmetric, x::Number) = Symmetric(parentof_applytri(y -> y * x, A), sym_uplo(A.uplo)) -*(x::Number, A::Symmetric) = Symmetric(parentof_applytri(y -> x * y, A), sym_uplo(A.uplo)) -*(A::Hermitian, x::Real) = Hermitian(parentof_applytri(y -> y * x, A), sym_uplo(A.uplo)) -*(x::Real, A::Hermitian) = Hermitian(parentof_applytri(y -> x * y, A), sym_uplo(A.uplo)) -/(A::Symmetric, x::Number) = Symmetric(parentof_applytri(y -> y/x, A), sym_uplo(A.uplo)) -/(A::Hermitian, x::Real) = Hermitian(parentof_applytri(y -> y/x, A), sym_uplo(A.uplo)) +*(A::Symmetric, x::Number) = Symmetric(parentof_applytri(y -> y * x, A), _sym_uplo(A.uplo)) +*(x::Number, A::Symmetric) = Symmetric(parentof_applytri(y -> x * y, A), _sym_uplo(A.uplo)) +*(A::Hermitian, x::Real) = Hermitian(parentof_applytri(y -> y * x, A), _sym_uplo(A.uplo)) +*(x::Real, A::Hermitian) = Hermitian(parentof_applytri(y -> x * y, A), _sym_uplo(A.uplo)) +/(A::Symmetric, x::Number) = Symmetric(parentof_applytri(y -> y/x, A), _sym_uplo(A.uplo)) +/(A::Hermitian, x::Real) = Hermitian(parentof_applytri(y -> y/x, A), _sym_uplo(A.uplo)) factorize(A::HermOrSym) = _factorize(A) function _factorize(A::HermOrSym{T}; check::Bool=true) where T @@ -850,8 +850,8 @@ function _inv(A::HermOrSym) B end # StridedMatrix restriction seems necessary due to inv! call in _inv above -inv(A::Hermitian{<:Any,<:StridedMatrix}) = Hermitian(_inv(A), sym_uplo(A.uplo)) -inv(A::Symmetric{<:Any,<:StridedMatrix}) = Symmetric(_inv(A), sym_uplo(A.uplo)) +inv(A::Hermitian{<:Any,<:StridedMatrix}) = Hermitian(_inv(A), _sym_uplo(A.uplo)) +inv(A::Symmetric{<:Any,<:StridedMatrix}) = Symmetric(_inv(A), _sym_uplo(A.uplo)) function svd(A::RealHermSymComplexHerm; full::Bool=false) vals, vecs = eigen(A) diff --git a/src/symmetriceigen.jl b/src/symmetriceigen.jl index f6640bc1..a8a54037 100644 --- a/src/symmetriceigen.jl +++ b/src/symmetriceigen.jl @@ -2,8 +2,12 @@ # preserve HermOrSym wrapper # Call `copytrito!` instead of `copy_similar` to only copy the matching triangular half -eigencopy_oftype(A::Hermitian, ::Type{S}) where S = Hermitian(copytrito!(similar(parent(A), S, size(A)), A.data, A.uplo), sym_uplo(A.uplo)) -eigencopy_oftype(A::Symmetric, ::Type{S}) where S = Symmetric(copytrito!(similar(parent(A), S, size(A)), A.data, A.uplo), sym_uplo(A.uplo)) +function eigencopy_oftype(A::Hermitian, ::Type{S}) where S + Hermitian(copytrito!(similar(parent(A), S, size(A)), A.data, A.uplo), _sym_uplo(A.uplo)) +end +function eigencopy_oftype(A::Symmetric, ::Type{S}) where S + Symmetric(copytrito!(similar(parent(A), S, size(A)), A.data, A.uplo), _sym_uplo(A.uplo)) +end eigencopy_oftype(A::Symmetric{<:Complex}, ::Type{S}) where S = copyto!(similar(parent(A), S), A) """ @@ -314,8 +318,8 @@ end # Perform U' \ A / U in-place, where U::Union{UpperTriangular,Diagonal} UtiAUi!(A, U) = _UtiAUi!(A, U) -UtiAUi!(A::Symmetric, U) = Symmetric(_UtiAUi!(copytri!(parent(A), A.uplo), U), sym_uplo(A.uplo)) -UtiAUi!(A::Hermitian, U) = Hermitian(_UtiAUi!(copytri!(parent(A), A.uplo, true), U), sym_uplo(A.uplo)) +UtiAUi!(A::Symmetric, U) = Symmetric(_UtiAUi!(copytri!(parent(A), A.uplo), U), _sym_uplo(A.uplo)) +UtiAUi!(A::Hermitian, U) = Hermitian(_UtiAUi!(copytri!(parent(A), A.uplo, true), U), _sym_uplo(A.uplo)) _UtiAUi!(A, U) = rdiv!(ldiv!(U', A), U) function eigvals(A::HermOrSym{TA}, B::HermOrSym{TB}; kws...) where {TA,TB} diff --git a/test/symmetric.jl b/test/symmetric.jl index 45125591..031dac27 100644 --- a/test/symmetric.jl +++ b/test/symmetric.jl @@ -1343,6 +1343,12 @@ end @test_throws msg LinearAlgebra.fillband!(Symmetric(A), 2, 0, 1) end +@testset "sym_uplo" begin + @test LinearAlgebra.sym_uplo('U') == :U + @test LinearAlgebra.sym_uplo('L') == :L + @test_throws ArgumentError LinearAlgebra.sym_uplo('N') +end + @testset "uplo" begin S = Symmetric([1 2; 3 4], :U) @test LinearAlgebra.uplo(S) == :U