Skip to content

Commit

Permalink
Merge pull request JuliaLang#11371 from KristofferC/kc/herm_sym
Browse files Browse the repository at this point in the history
Improve performance of issym and ishermitian for sparse matrices
  • Loading branch information
andreasnoack committed May 27, 2015
2 parents 4c7bb25 + 48d5e5f commit c71cf48
Show file tree
Hide file tree
Showing 4 changed files with 73 additions and 8 deletions.
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

0 comments on commit c71cf48

Please sign in to comment.