Rewrite the RS splitting algorithm

Author: ranjanan

src/classical.jl

@@ -42,12 +42,12 @@ function direct_interpolation{T,V}(A::T, S::T, splitting::Vector{V})
     fill!(S.nzval, 1.)
     S = A .* S
     Pp = rs_direct_interpolation_pass1(S, A, splitting)
-    Pp .= Pp .+ 1
+    Pp = Pp .+ 1
 
     Px, Pj = rs_direct_interpolation_pass2(A, S, splitting, Pp)
 
     # Px .= abs.(Px)
-    Pj .= Pj .+ 1
+    Pj = Pj .+ 1
 
     R = SparseMatrixCSC(maximum(Pj), size(A, 1), Pp, Pj, Px)
     P = R'
@@ -59,7 +59,7 @@ end
 function rs_direct_interpolation_pass1(S, A, splitting)
 
      Bp = zeros(Int, size(A.colptr))
-     Sp = S.colptr
+     #=Sp = S.colptr
      Sj = S.rowval
      n_nodes = size(A, 1)
      nnz = 0
@@ -75,6 +75,21 @@ function rs_direct_interpolation_pass1(S, A, splitting)
             end
         end
          Bp[i+1] = nnz
+     end=#
+     n = size(A, 1)
+     nnz = 0
+     for i = 1:n
+         if splitting[i] == C_NODE
+             nnz += 1
+         else
+             for j in nzrange(S, i)
+                 row = S.rowval[j]
+                 if splitting[row] == C_NODE && row != i
+                     nnz += 1
+                end
+             end
+        end
+        Bp[i+1] = nnz
      end
      Bp
  end

src/splitting.jl

@@ -12,12 +12,13 @@ function RS_CF_splitting(S::SparseMatrixCSC)
 
 	n_nodes = n
 	lambda = zeros(Int, n)
-	Tp = S.colptr
-	Tj = S.rowval
-	Sp = Tp
-	Sj = Tj
+	T = S'
+	Tp = T.colptr
+	Tj = T.rowval
+	Sp = S.colptr
+	Sj = S.rowval
 
-    # compute lambdas
+    # compute lambdas - number of neighbors
     for i = 1:n
         lambda[i] = Tp[i+1] - Tp[i]
     end
@@ -27,6 +28,7 @@ function RS_CF_splitting(S::SparseMatrixCSC)
 	index_to_node = zeros(Int,n)
 	node_to_index = zeros(Int,n)
 
+	# Number of nodes with a certain neighbor count
     for i = 1:n
         interval_count[lambda[i]+1] += 1
     end
@@ -52,9 +54,83 @@ function RS_CF_splitting(S::SparseMatrixCSC)
             splitting[i] = F_NODE
 		end
     end
+	@show lambda
+	@show node_to_index
+	@show index_to_node
+	@show interval_ptr
+
+	for top_index = n_nodes:-1:1
+		i = index_to_node[top_index]
+		lambda_i = lambda[i] + 1
+		interval_count[lambda_i] -= 1
+		if splitting[i] == F_NODE
+			continue
+		else
+			@assert splitting[i] == U_NODE
+			splitting[i] = C_NODE
+			for j in nzrange(T, i)
+				row = T.rowval[j]
+				if splitting[row] == U_NODE
+					splitting[row] = F_NODE
+
+					for k in nzrange(S, row)
+						rowk = S.rowval[k]
+						if splitting[rowk] == U_NODE
+							lambda[rowk] >= n_nodes - 1 && continue
+							lambda_k = lambda[rowk] + 1
+						  	old_pos  = node_to_index[rowk]
+						  	new_pos  = interval_ptr[lambda_k] + interval_count[lambda_k]# - 1
+
+						  	node_to_index[index_to_node[old_pos]] = new_pos
+						  	node_to_index[index_to_node[new_pos]] = old_pos
+						  	(index_to_node[old_pos], index_to_node[new_pos]) = (index_to_node[new_pos], index_to_node[old_pos])
+
+						  	# update intervals
+						  	interval_count[lambda_k]   -= 1
+						  	interval_count[lambda_k + 1] += 1 # invalid write!
+						  	interval_ptr[lambda_k + 1]    = new_pos - 1
+
+						 	# increment lambda_k
+						 	lambda[rowk] += 1
+						end
+					end
+				end
+			end
+			for j in nzrange(S, i)
+				row = S.rowval[j]
+				if splitting[row] == U_NODE
+
+					lambda[row] == 0 && continue
+
+					# assert(lambda[j] > 0);//this would cause problems!
+
+					# move j to the beginning of its current interval
+					lambda_j = lambda[row] + 1
+					old_pos  = node_to_index[row]
+					new_pos  = interval_ptr[lambda_j]
+
+					node_to_index[index_to_node[old_pos]] = new_pos
+					node_to_index[index_to_node[new_pos]] = old_pos
+					(index_to_node[old_pos],index_to_node[new_pos]) = (index_to_node[new_pos],index_to_node[old_pos])
+
+					# update intervals
+					interval_count[lambda_j]   -= 1
+					interval_count[lambda_j-1] += 1
+					interval_ptr[lambda_j]     += 1
+					interval_ptr[lambda_j-1]    = interval_ptr[lambda_j] - interval_count[lambda_j-1]
+
+					# decrement lambda_j
+					lambda[row] -= 1
+				end
+			end
+		end
+	end
+	splitting
+end
+
 
     # Now add elements to C and F, in descending order of lambda
-    for top_index = n_nodes:-1:1
+    #=for top_index = n_nodes:-1:1
         i        = index_to_node[top_index]
         lambda_i = lambda[i] + 1
 
@@ -74,6 +150,7 @@ function RS_CF_splitting(S::SparseMatrixCSC)
 
             # For each j in S^T_i /\ U
             for jj = Tp[i]:Tp[i+1]-1
+				#jj > length(Tp) && continue
 
                 j = Tj[jj]
 
@@ -82,6 +159,7 @@ function RS_CF_splitting(S::SparseMatrixCSC)
 
                     # For each k in S_j /\ U
                     for kk = Sp[j]: Sp[j+1]-1
+						# kk > length(Sj) && continue
                         k = Sj[kk]
 
                         if splitting[k] == U_NODE
@@ -109,7 +187,8 @@ function RS_CF_splitting(S::SparseMatrixCSC)
 			end
 
             # For each j in S_i /\ U
-            for jj = Sp[i]: Sp[i+1]-1
+            for jj = Sp[i]: Sp[i+1] - 1
+				# jj > length(Sj) && continue
 
                 j = Sj[jj]
 
@@ -141,4 +220,4 @@ function RS_CF_splitting(S::SparseMatrixCSC)
         end
     end
 	splitting
-end
+end=#

src/strength.jl

@@ -29,7 +29,9 @@ function strength_of_connection{T}(c::Classical{T}, A::SparseMatrixCSC)
     end
     S = sparse(I, J, V, m, n)
 
-    scale_cols_by_largest_entry(S)
+    scale_cols_by_largest_entry!(S)
+
+    S'
 end
 
 function find_max_off_diag(neighbors, col)
@@ -40,26 +42,15 @@ function find_max_off_diag(neighbors, col)
     return maxval
 end
 
-function scale_cols_by_largest_entry(A::SparseMatrixCSC)
-
-    m,n = size(A)
-
-    I = zeros(Int, size(A.nzval))
-    J = similar(I)
-    V = zeros(size(A.nzval))
+function scale_cols_by_largest_entry!(A::SparseMatrixCSC)
 
-    k = 1
+    n = size(A, 1)
     for i = 1:n
         _m = maximum(A[:,i])
         for j in nzrange(A, i)
-            row = A.rowval[j]
-            val = A.nzval[j]
-            I[k] = row
-            J[k] = i
-            V[k] = val / _m
-            k += 1
+            A.nzval[j] /= _m
         end
     end
 
-    sparse(I,J,V,m,n)
+    A
 end

test/runtests.jl

@@ -19,7 +19,7 @@ end
 # Ruge-Stuben splitting
 S = poisson(7)
 @test split_nodes(RS(), S) == [0, 1, 0, 1, 0, 1, 0]
-
+@show "buzz"
 srand(0)
 S = sprand(10,10,0.1); S = S + S'
 @test split_nodes(RS(), S) ==  [0, 1, 1, 0, 0, 0, 0, 0, 1, 1]