Compute the number of self-loops

Author: amontoison

src/graph.jl

@@ -179,6 +179,7 @@ function build_edge_to_index(S::SparsityPatternCSC{T}) where {T}
     # edge_to_index gives an index for each edge
     edge_to_index = Vector{T}(undef, nnz(S))
     offsets = zeros(T, S.n)
+    nb_self_loops = 0
     counter = 0
     rvS = rowvals(S)
     for j in axes(S, 2)
@@ -193,10 +194,11 @@ function build_edge_to_index(S::SparsityPatternCSC{T}) where {T}
             elseif i == j
                 # this should never be used, make sure it errors
                 edge_to_index[k] = 0
+                nb_self_loops += 1
             end
         end
     end
-    return edge_to_index
+    return edge_to_index, nb_self_loops
 end
 
 ## Adjacency graph
@@ -227,16 +229,23 @@ The adjacency graph of a symmetric matrix `A ∈ ℝ^{n × n}` is `G(A) = (V, E)
 struct AdjacencyGraph{T<:Integer,augmented_graph}
     S::SparsityPatternCSC{T}
     edge_to_index::Vector{T}
+    nb_self_loops::Int
 end
 
 Base.eltype(::AdjacencyGraph{T}) where {T} = T
 
 function AdjacencyGraph(
     S::SparsityPatternCSC{T},
-    edge_to_index::Vector{T}=build_edge_to_index(S);
+    edge_to_index::Vector{T},
+    nb_self_loops::Int;
     augmented_graph::Bool=false,
 ) where {T}
-    return AdjacencyGraph{T,augmented_graph}(S, edge_to_index)
+    return AdjacencyGraph{T,augmented_graph}(S, edge_to_index, nb_self_loops)
+end
+
+function AdjacencyGraph(S::SparsityPatternCSC; augmented_graph::Bool=false)
+    edge_to_index, nb_self_loops = build_edge_to_index(S)
+    return AdjacencyGraph(S, edge_to_index, nb_self_loops; augmented_graph)
 end
 
 function AdjacencyGraph(A::SparseMatrixCSC; augmented_graph::Bool=false)

src/interface.jl

@@ -323,7 +323,7 @@ function _coloring(
     forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing,
 ) where {R}
     A_and_Aᵀ, edge_to_index = bidirectional_pattern(A; symmetric_pattern)
-    ag = AdjacencyGraph(A_and_Aᵀ, edge_to_index; augmented_graph=true)
+    ag = AdjacencyGraph(A_and_Aᵀ, edge_to_index, 0; augmented_graph=true)
     outputs_by_order = map(algo.orders) do order
         vertices_in_order = vertices(ag, order)
         _color, _star_set = star_coloring(
@@ -372,7 +372,7 @@ function _coloring(
     symmetric_pattern::Bool,
 ) where {R}
     A_and_Aᵀ, edge_to_index = bidirectional_pattern(A; symmetric_pattern)
-    ag = AdjacencyGraph(A_and_Aᵀ, edge_to_index; augmented_graph=true)
+    ag = AdjacencyGraph(A_and_Aᵀ, edge_to_index, 0; augmented_graph=true)
     outputs_by_order = map(algo.orders) do order
         vertices_in_order = vertices(ag, order)
         _color, _tree_set = acyclic_coloring(ag, vertices_in_order, algo.postprocessing)

src/result.jl

@@ -334,16 +334,15 @@ function star_csc_indices(
     decompression_uplo::Symbol,
 ) where {T}
     (; star, hub) = star_set
-    S = pattern(ag)
-    edge_to_index = edge_indices(ag)
+    (; S, edge_to_index, nb_self_loops) = ag
     n = S.n
     rvS = rowvals(S)
-    l = 0
-    if augmented_graph(ag) && (decompression_uplo != :F)
-        compressed_indices = zeros(T, nnz(S) ÷ 2)
-    else
-        compressed_indices = zeros(T, nnz(S))  # needs to be independent from the storage in the graph, in case the graph gets reused
+    nb_indices = nnz(S)
+    if decompression_uplo != :F
+        nb_indices = (nb_indices - nb_self_loops) ÷ 2 + nb_self_loops
     end
+    compressed_indices = zeros(T, nb_indices)  # needs to be independent from the storage in the graph, in case the graph gets reused
+    l = 0
     for j in axes(S, 2)
         for k in nzrange(S, j)
             i = rvS[k]
@@ -374,9 +373,6 @@ function star_csc_indices(
             end
         end
     end
-    if !augmented_graph(ag) && (decompression_uplo != :F)
-        resize!(compressed_indices, l)
-    end
     return compressed_indices
 end
 
@@ -422,26 +418,14 @@ function TreeSetColoringResult(
     decompression_uplo::Symbol=:F,
 ) where {T<:Integer,R}
     (; reverse_bfs_orders, tree_edge_indices, nt) = tree_set
-    (; S) = ag
+    (; S, nb_self_loops) = ag
     nvertices = length(color)
     group = group_by_color(T, color)
     rv = rowvals(S)
 
-    ndiag = 0
-    if !augmented_graph(ag)
-        for j in axes(S, 2)
-            for k in nzrange(S, j)
-                i = rv[k]
-                if i == j
-                    ndiag += 1
-                end
-            end
-        end
-    end
-
     # Vector for the decompression of the diagonal coefficients
-    diagonal_indices = Vector{T}(undef, ndiag)
-    diagonal_nzind = (decompression_uplo == :F) ? Vector{T}(undef, ndiag) : T[]
+    diagonal_indices = Vector{T}(undef, nb_self_loops)
+    diagonal_nzind = (decompression_uplo == :F) ? Vector{T}(undef, nb_self_loops) : T[]
 
     if !augmented_graph(ag)
         l = 0
@@ -460,7 +444,7 @@ function TreeSetColoringResult(
     end
 
     # Vectors for the decompression of the off-diagonal coefficients
-    nedges = (nnz(S) - ndiag) ÷ 2
+    nedges = (nnz(S) - nb_self_loops) ÷ 2
     lower_triangle_offsets = decompression_uplo == :U ? T[] : Vector{T}(undef, nedges)
     upper_triangle_offsets = decompression_uplo == :L ? T[] : Vector{T}(undef, nedges)