Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 14 additions & 0 deletions src/LinearAlgebra.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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'
Expand All @@ -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)"))

Expand Down
2 changes: 1 addition & 1 deletion src/bidiag.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
4 changes: 3 additions & 1 deletion src/bunchkaufman.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
12 changes: 6 additions & 6 deletions src/special.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
72 changes: 36 additions & 36 deletions src/symmetric.jl
Original file line number Diff line number Diff line change
Expand Up @@ -206,15 +206,15 @@ 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) &&
any(!isreal, A.data[i] for i in diagind(A.data, IndexStyle(A.data)))

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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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'
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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))
Expand All @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down
12 changes: 8 additions & 4 deletions src/symmetriceigen.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)

"""
Expand Down Expand Up @@ -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}
Expand Down
6 changes: 6 additions & 0 deletions test/symmetric.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down