diff --git a/SubgraphIsomorphism/ktruss/code/julia/ktruss.jl b/SubgraphIsomorphism/ktruss/code/julia/ktruss.jl index 9bfff12..3195a63 100644 --- a/SubgraphIsomorphism/ktruss/code/julia/ktruss.jl +++ b/SubgraphIsomorphism/ktruss/code/julia/ktruss.jl @@ -1,34 +1,71 @@ +const Chol = if isdefined(Base, :SparseArrays) + Base.SparseArrays.CHOLMOD +else + Base.SparseMatrix.CHOLMOD +end + +colptr0(A::Chol.Sparse, col) = unsafe_load(unsafe_load(A.p).p, col) +rowval(A::Chol.Sparse, j) = unsafe_load(unsafe_load(A.p).i, j) + 1 +nzval(A::Chol.Sparse, j) = unsafe_load(unsafe_load(A.p).x, j) + function calcx(E, m, n, k) - tmp = E.'*E; - R = E * ( tmp - spdiagm( diag(tmp) ) ); - r,c,v = findnz(R); - id = v.==2; - A = sparse( r[id], c[id], 1, m, n); - s = sum(A, 2); - x = s .< (k-2); - return(x, !x); + CSE = Chol.Sparse(E) + tmp = CSE'*CSE + + # subtract spdiagm( diag(tmp) ) from tmp in-place by setting diagonals to 0 + tmp_ptr = unsafe_load(tmp.p) + for col in 1:n + for j in colptr0(tmp, col) + 1 : colptr0(tmp, col+1) + if rowval(tmp, j) == col + unsafe_store!(tmp_ptr.x, 0, j) + end + end + end + + R = CSE * tmp + s = zeros(Int, m) + + @inbounds for col in 1:n + for j in colptr0(R, col) + 1 : colptr0(R, col+1) + if nzval(R, j) == 2 + s[rowval(R, j)] += 1 + end + end + end + + rowflags = falses(m) # rows of E with at least one nonzero in them + @inbounds for col in 1:n + for j in colptr0(CSE, col) + 1 : colptr0(CSE, col+1) + if nzval(CSE, j) != 0 + rowflags[rowval(CSE, j)] = true + end + end + end + + return find(s .>= (k-2)), sum(rowflags) end function ktruss(inc_mtx_file, k) - if ~isfile( inc_mtx_file ) - println("unable to open input file"); - return (-1); + if !isfile( inc_mtx_file ) + println("unable to open input file") + return (-1) end # load input data - t_read_inc=@elapsed ii = readdlm( inc_mtx_file, '\t', Int64); - println("incidence matrix read time : ", t_read_inc); - - t_create_inc=@elapsed E = sparse( ii[:,1], ii[:,2], ii[:,3] ); - println("sparse adj. matrix creation time : ", t_create_inc); - - # - tic(); - m,n = size(E); - x, xc = calcx(E, m, n, k); - while sum(xc) != sum( any(E,2) ) - E[find(x), :] = 0 - x, xc = calcx(E, m, n, k); + t_read_inc=@elapsed ii = readdlm( inc_mtx_file, '\t', Int64) + println("incidence matrix read time : ", t_read_inc) + + t_create_inc=@elapsed E = sparse( ii[:,1], ii[:,2], ii[:,3] ) + println("sparse adj. matrix creation time : ", t_create_inc) + + tic() + m,n = size(E) + xcinds, rowcount = calcx(E, m, n, k) + while length(xcinds) != rowcount + # keep rows of E where x is false, equivalent to E[find(x), :] = 0 + Ekeep = E[xcinds, :] + E = SparseMatrixCSC(m, n, Ekeep.colptr, xcinds[Ekeep.rowval], Ekeep.nzval) + xcinds, rowcount = calcx(E, m, n, k) end return E