From 230c8f142e7d242dc3702057cae37cfd22d7d736 Mon Sep 17 00:00:00 2001
From: Kristoffer Carlsson <kristoffer.carlsson@chalmers.se>
Date: Wed, 20 May 2015 13:41:50 +0200
Subject: [PATCH] improve performance of issym and ishermitian

---
 base/functors.jl            |  3 +++
 base/sparse.jl              |  2 +-
 base/sparse/sparsematrix.jl | 38 ++++++++++++++++++++++++++++++-------
 test/sparsedir/sparse.jl    | 38 +++++++++++++++++++++++++++++++++++++
 4 files changed, 73 insertions(+), 8 deletions(-)

diff --git a/base/functors.jl b/base/functors.jl
index 24c10dfd43f69..646facee9d749 100644
--- a/base/functors.jl
+++ b/base/functors.jl
@@ -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
 
diff --git a/base/sparse.jl b/base/sparse.jl
index bf82e81d3f712..d931e43082e4e 100644
--- a/base/sparse.jl
+++ b/base/sparse.jl
@@ -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
 
diff --git a/base/sparse/sparsematrix.jl b/base/sparse/sparsematrix.jl
index 8cc9f676a3525..ab9092a4562fd 100644
--- a/base/sparse/sparsematrix.jl
+++ b/base/sparse/sparsematrix.jl
@@ -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]
+          val == 0 && continue # In case of explicit zeros
+          row = rowval[p]
+          row < col && return false
+          if row == col
+              if val != check(val)
+                  return false
+              end
+          else
+              row2 = rowval[tracker[row]]
+              row2 > col && return false
+              if row2 == col
+                 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})
diff --git a/test/sparsedir/sparse.jl b/test/sparsedir/sparse.jl
index 0c2e3d1fb41e2..db22999119389 100644
--- a/test/sparsedir/sparse.jl
+++ b/test/sparsedir/sparse.jl
@@ -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