Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Make ishermitian and issym test for approx. symmetry for floats. fix #10298 #10369

Closed
wants to merge 1 commit into from
Closed
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
Add the issym/ishermitian based on norm(A - A', Inf)<tol*norm(A,Inf). f…
…ix #10298.
  • Loading branch information
dhoegh authored and dhoegh committed Oct 25, 2015
commit 80b7c207348944cc53cbaa81685dc00b4fdee4d0
8 changes: 4 additions & 4 deletions base/docs/helpdb.jl
Original file line number Diff line number Diff line change
Expand Up @@ -7211,9 +7211,9 @@ Given `@osx? a : b`, do `a` on OS X and `b` elsewhere. See documentation for Han
:@osx

doc"""
ishermitian(A) -> Bool
ishermitian(A; tol=eps(eltype(real(A)))*norm(A,1)) -> Bool
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think you mean real(eltype(A)). Doing real(A) makes a real copy of the matrix if it is complex.


Test whether a matrix is Hermitian.
Test whether a matrix is Hermitian. The `tol` keyword is only used for matrices of floating point or `Complex{T<:AbstractFloat}` numbers, where it tests `norm(A - A', Inf) < tol*norm(A,Inf)`.
"""
ishermitian

Expand Down Expand Up @@ -7345,9 +7345,9 @@ The arguments to a function or constructor are outside the valid domain.
DomainError

doc"""
issym(A) -> Bool
issym(A; tol=eps(eltype(real(A)))*norm(A,1)) -> Bool

Test whether a matrix is symmetric.
Test whether a matrix is symmetric. The `tol` keyword is only used for matrices of floating point or `Complex{T<:AbstractFloat}` numbers, where it tests `norm(A - A', Inf) < tol*norm(A,Inf)`.
"""
issym

Expand Down
35 changes: 33 additions & 2 deletions base/linalg/generic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -350,7 +350,7 @@ condskeel{T<:Integer}(A::AbstractMatrix{T}, p::Real=Inf) = norm(abs(inv(float(A)
condskeel(A::AbstractMatrix, x::AbstractVector, p::Real=Inf) = norm(abs(inv(A))*abs(A)*abs(x), p)
condskeel{T<:Integer}(A::AbstractMatrix{T}, x::AbstractVector, p::Real=Inf) = norm(abs(inv(float(A)))*abs(A)*abs(x), p)

function issym(A::AbstractMatrix)
function _issym(A::AbstractMatrix)
m, n = size(A)
if m != n
return false
Expand All @@ -365,7 +365,7 @@ end

issym(x::Number) = true

function ishermitian(A::AbstractMatrix)
function _ishermitian(A::AbstractMatrix)
m, n = size(A)
if m != n
return false
Expand All @@ -378,6 +378,37 @@ function ishermitian(A::AbstractMatrix)
return true
end

ishermitian(A::AbstractMatrix) = _ishermitian(A::AbstractMatrix)
issym(A::AbstractMatrix) = _issym(A::AbstractMatrix)

for (f,t) in ((:ishermitian, :ctranspose),(:issym, :transpose))
eval(quote
# ishermitian and issym functions test whether norm(A - A', Inf)<tol*norm(A,Inf) for FloatingPoint matrices
function $f{T<:AbstractFloat}(A::Union{AbstractMatrix{Complex{T}}, AbstractMatrix{T}}, tol=eps(T)*norm(A,1))
# This is performed due to if tol==0. _issym/_ishermitian is a lot faster for exact equality.
tol == 0. && return $(symbol("_", f))(A)
m,n = size(A)
Tnorm = typeof(float(real(zero(T))))
Tsum = promote_type(Float64,Tnorm)
nrm_diff::Tsum = 0
nrm_A::Tsum = 0
@inbounds begin
for i = 1:m
nrmi_diff::Tsum = 0
nrmi_A::Tsum = 0
for j = 1:n
nrmi_diff += abs(A[i,j]-$t(A[j,i]))
nrmi_A += abs(A[i,j])
end
nrm_diff = max(nrm_diff,nrmi_diff)
nrm_A = max(nrm_A,nrmi_A)
end
end
return nrm_diff < tol*nrm_A
end
end)
end

ishermitian(x::Number) = (x == conj(x))

function istriu(A::AbstractMatrix)
Expand Down
4 changes: 4 additions & 0 deletions base/linalg/symmetric.jl
Original file line number Diff line number Diff line change
Expand Up @@ -105,6 +105,10 @@ function triu(A::Symmetric, k::Integer=0)
end
end

#Test whether a matrix is positive-definite
isposdef!(A::HermOrSym) = isposdef!(A.data, symbol(A.uplo))
isposdef{T}(A::HermOrSym{T}) = (S = typeof(sqrt(one(T))); isposdef!(S == T ? copy(A.data) : convert(AbstractMatrix{S}, A.data), symbol(A.uplo)))

## Matvec
A_mul_B!{T<:BlasFloat,S<:StridedMatrix}(y::StridedVector{T}, A::Symmetric{T,S}, x::StridedVector{T}) = BLAS.symv!(A.uplo, one(T), A.data, x, zero(T), y)
A_mul_B!{T<:BlasComplex,S<:StridedMatrix}(y::StridedVector{T}, A::Hermitian{T,S}, x::StridedVector{T}) = BLAS.hemv!(A.uplo, one(T), A.data, x, zero(T), y)
Expand Down
22 changes: 22 additions & 0 deletions test/linalg/symmetric.jl
Original file line number Diff line number Diff line change
Expand Up @@ -209,3 +209,25 @@ let A = [1.0+im 2.0; 2.0 0.0]
@test !ishermitian(A)
@test_throws ArgumentError Hermitian(A)
end

# Tests for #10298
for elT in [Float32, Float64]
A1 = ones(elT,10,10)
A = convert(Array{Complex{elT}}, A1)
A[end,1] = Complex{elT}(nextfloat((A[end,1].re)), A[end,1].im)
@test issym(A)
A = A + triu(A1)im - tril(A1)im
@test ishermitian(A)

D=diagm(rand(100))*100
D[end,1] = nextfloat((D[end,1]))

for HOrS in [Hermitian, Symmetric]
@test isposdef(HOrS(D))
@test isposdef!(copy(HOrS(D)))
end
@test ishermitian(D)
@test issym(D)
@test isposdef(D)
end