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

Improve performance of issym and ishermitian for sparse matrices #11371

Merged
merged 1 commit into from
May 27, 2015
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
3 changes: 3 additions & 0 deletions base/functors.jl
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,9 @@ call(::ExpFun, x) = exp(x)
immutable LogFun <: Func{1} end
call(::LogFun, x) = log(x)

immutable ConjFun <: Func{1} end
call(::ConjFun, x) = conj(x)

immutable AndFun <: Func{2} end
call(::AndFun, x, y) = x & y

Expand Down
2 changes: 1 addition & 1 deletion base/sparse.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@

module SparseMatrix

using Base: Func, AddFun, OrFun
using Base: Func, AddFun, OrFun, ConjFun, IdFun
using Base.Sort: Forward
using Base.LinAlg: AbstractTriangular

Expand Down
38 changes: 31 additions & 7 deletions base/sparse/sparsematrix.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2376,17 +2376,41 @@ end
squeeze(S::SparseMatrixCSC, dims::Dims) = throw(ArgumentError("squeeze is not available for sparse matrices"))

## Structure query functions
issym(A::SparseMatrixCSC) = is_hermsym(A, IdFun())

function issym(A::SparseMatrixCSC)
m, n = size(A)
if m != n; return false; end
return countnz(A - A.') == 0
end
ishermitian(A::SparseMatrixCSC) = is_hermsym(A, ConjFun())

function ishermitian(A::SparseMatrixCSC)
function is_hermsym(A::SparseMatrixCSC, check::Func)
m, n = size(A)
if m != n; return false; end
return countnz(A - A') == 0

colptr = A.colptr
rowval = A.rowval
nzval = A.nzval
tracker = copy(A.colptr)
@inbounds for col = 1:A.n
for p = tracker[col]:colptr[col+1]-1
val = nzval[p]
if val == 0; continue; end # In case of explicit zeros
row = rowval[p]
if row < col; return false; end
if row == col # Diagonal element
if val != check(val)
return false
end
else
row2 = rowval[tracker[row]]
if row2 > col; return false; end
if row2 == col # A[i,j] and A[j,i] exists
if val != check(nzval[tracker[row]])
return false
end
tracker[row] += 1
end
end
end
end
return true
end

function istriu{Tv}(A::SparseMatrixCSC{Tv})
Expand Down
38 changes: 38 additions & 0 deletions test/sparsedir/sparse.jl
Original file line number Diff line number Diff line change
Expand Up @@ -930,3 +930,41 @@ A = sparse([1.0])
@test norm(A) == 1.0
@test_throws ArgumentError norm(sprand(5,5,0.2),3)
@test_throws ArgumentError norm(sprand(5,5,0.2),2)

# test ishermitian and issym real matrices
A = speye(5,5)
@test ishermitian(A) == true
@test issym(A) == true
A[1,3] = 1.0
@test ishermitian(A) == false
@test issym(A) == false
A[3,1] = 1.0
@test ishermitian(A) == true
@test issym(A) == true

# test ishermitian and issym complex matrices
A = speye(5,5) + im*speye(5,5)
@test ishermitian(A) == false
@test issym(A) == true
A[1,4] = 1.0 + im
@test ishermitian(A) == false
@test issym(A) == false

A = speye(Complex128, 5,5)
A[3,2] = 1.0 + im
@test ishermitian(A) == false
@test issym(A) == false
A[2,3] = 1.0 - im
@test ishermitian(A) == true
@test issym(A) == false

A = sparse(zeros(5,5))
@test ishermitian(A) == true
@test issym(A) == true

# Test with explicit zeros
A = speye(Complex128, 5,5)
A[3,1] = 2
A.nzval[2] = 0.0
@test ishermitian(A) == true
@test issym(A) == true