Avoid leaking memory between structs

Author: gdalle

src/coloring.jl

@@ -81,8 +81,7 @@ function star_coloring(g::AdjacencyGraph, order::AbstractOrder, postprocessing::
     nv = nb_vertices(g)
     ne = nb_edges(g)
     color = zeros(Int, nv)
-    forbidden_colors = g.vertex_buffer  # Vector of length nv
-    fill!(forbidden_colors, 0)
+    forbidden_colors = zeros(Int, nv)
     first_neighbor = fill((0, 0, 0), nv)  # at first no neighbors have been encountered
     treated = zeros(Int, nv)
     star = Vector{Int}(undef, ne)
@@ -123,7 +122,9 @@ function star_coloring(g::AdjacencyGraph, order::AbstractOrder, postprocessing::
     end
     star_set = StarSet(star, hub)
     if postprocessing
-        postprocess!(color, star_set, g)
+        # Reuse the vector forbidden_colors to compute offsets during post-processing
+        offsets = forbidden_colors
+        postprocess!(color, star_set, g, offsets)
     end
     return color, star_set
 end
@@ -230,8 +231,7 @@ function acyclic_coloring(g::AdjacencyGraph, order::AbstractOrder, postprocessin
     nv = nb_vertices(g)
     ne = nb_edges(g)
     color = zeros(Int, nv)
-    forbidden_colors = g.vertex_buffer  # Vector of length nv
-    fill!(forbidden_colors, 0)
+    forbidden_colors = zeros(Int, nv)
     first_neighbor = fill((0, 0, 0), nv)  # at first no neighbors have been encountered
     first_visit_to_tree = fill((0, 0), ne)
     forest = Forest{Int}(ne)
@@ -289,7 +289,9 @@ function acyclic_coloring(g::AdjacencyGraph, order::AbstractOrder, postprocessin
 
     tree_set = TreeSet(g, forest)
     if postprocessing
-        postprocess!(color, tree_set, g)
+        # Reuse the vector forbidden_colors to compute offsets during post-processing
+        offsets = forbidden_colors
+        postprocess!(color, tree_set, g, offsets)
     end
     return color, tree_set
 end
@@ -517,6 +519,7 @@ function postprocess!(
     color::AbstractVector{<:Integer},
     star_or_tree_set::Union{StarSet,TreeSet},
     g::AdjacencyGraph,
+    offsets::Vector{Int},
 )
     S = pattern(g)
     edge_to_index = edge_indices(g)
@@ -621,7 +624,6 @@ function postprocess!(
 
     # if at least one of the colors is useless, modify the color assignments of vertices
     if any(!, color_used)
-        offsets = g.vertex_buffer
         num_colors_useless = 0
 
         # determine what are the useless colors and compute the offsets

src/graph.jl

@@ -36,7 +36,7 @@ SparseArrays.nzrange(S::SparsityPatternCSC, j::Integer) = S.colptr[j]:(S.colptr[
 
 Return a [`SparsityPatternCSC`](@ref) corresponding to the transpose of `S`.
 """
-function Base.transpose(S::SparsityPatternCSC{T}) where {T<:Integer}
+function Base.transpose(S::SparsityPatternCSC{T}) where {T}
     m, n = size(S)
     nnzA = nnz(S)
     A_colptr = S.colptr
@@ -94,16 +94,14 @@ function Base.getindex(S::SparsityPatternCSC, i0::Integer, i1::Integer)
 end
 
 """
-    bidirectional_pattern(A::AbstractMatrix; symmetric_pattern::Bool=false)
+    bidirectional_pattern(A::AbstractMatrix; symmetric_pattern::Bool)
 
 Return a [`SparsityPatternCSC`](@ref) corresponding to the matrix `[0 Aᵀ; A 0]`, with a minimum of allocations.
 """
-bidirectional_pattern(A::AbstractMatrix; symmetric_pattern::Bool=false) =
+bidirectional_pattern(A::AbstractMatrix; symmetric_pattern::Bool) =
     bidirectional_pattern(SparsityPatternCSC(SparseMatrixCSC(A)); symmetric_pattern)
 
-function bidirectional_pattern(
-    S::SparsityPatternCSC{T}; symmetric_pattern::Bool=false
-) where {T<:Integer}
+function bidirectional_pattern(S::SparsityPatternCSC{T}; symmetric_pattern::Bool) where {T}
     m, n = size(S)
     p = m + n
     nnzS = nnz(S)
@@ -175,7 +173,7 @@ function bidirectional_pattern(
     return S_and_Sᵀ, edge_to_index
 end
 
-function build_edge_to_index(S::SparsityPatternCSC{T}) where {T<:Integer}
+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)
@@ -196,7 +194,7 @@ function build_edge_to_index(S::SparsityPatternCSC{T}) where {T<:Integer}
             end
         end
     end
-    return edge_to_index, offsets
+    return edge_to_index
 end
 
 ## Adjacency graph
@@ -219,7 +217,6 @@ The adjacency graph of a symmetric matrix `A ∈ ℝ^{n × n}` is `G(A) = (V, E)
 
 - `S::SparsityPatternCSC{T}`: Underlying sparsity pattern, whose diagonal is empty whenever `has_diagonal` is `false`
 - `edge_to_index::Vector{T}`: A vector mapping each nonzero of `S` to a unique edge index (ignoring diagonal and accounting for symmetry, so that `(i, j)` and `(j, i)` get the same index)
-- `vertex_buffer::Vector{T}`: A buffer of length `S.n`, which is the number of vertices in the graph
 
 # References
 
@@ -228,28 +225,14 @@ The adjacency graph of a symmetric matrix `A ∈ ℝ^{n × n}` is `G(A) = (V, E)
 struct AdjacencyGraph{T,has_diagonal}
     S::SparsityPatternCSC{T}
     edge_to_index::Vector{T}
-    vertex_buffer::Vector{T}
 end
 
 function AdjacencyGraph(
     S::SparsityPatternCSC{T},
-    edge_to_index::Vector{T},
-    vertex_buffer::Vector{T};
+    edge_to_index::Vector{T}=build_edge_to_index(S);
     has_diagonal::Bool=true,
 ) where {T}
-    return AdjacencyGraph{T,has_diagonal}(S, edge_to_index, vertex_buffer)
-end
-
-function AdjacencyGraph(
-    S::SparsityPatternCSC{T}, edge_to_index::Vector{T}; has_diagonal::Bool=true
-) where {T}
-    vertex_buffer = Vector{Int}(undef, S.n)
-    return AdjacencyGraph(S, edge_to_index, vertex_buffer; has_diagonal)
-end
-
-function AdjacencyGraph(S::SparsityPatternCSC; has_diagonal::Bool=true)
-    edge_to_index, vertex_buffer = build_edge_to_index(S)
-    return AdjacencyGraph(S, edge_to_index, vertex_buffer; has_diagonal)
+    return AdjacencyGraph{T,has_diagonal}(S, edge_to_index)
 end
 
 function AdjacencyGraph(A::SparseMatrixCSC; has_diagonal::Bool=true)

src/result.jl

@@ -241,15 +241,12 @@ end
 function StarSetColoringResult(
     A::AbstractMatrix, ag::AdjacencyGraph, color::Vector{Int}, star_set::StarSet
 )
-    # Reuse edge_to_index to store the compressed indices for decompression
-    edge_to_index = edge_indices(ag)
-    compressed_indices = edge_to_index
-
     (; star, hub) = star_set
     S = pattern(ag)
     n = S.n
     group = group_by_color(color)
     rvS = rowvals(S)
+    compressed_indices = zeros(Int, nnz(S))  # needs to be independent from the storage in the graph, in case the graph gets reused
     for j in axes(S, 2)
         for k in nzrange(S, j)
             i = rvS[k]

test/graph.jl

@@ -16,42 +16,47 @@ using Test
 
 @testset "SparsityPatternCSC" begin
     @testset "Transpose" begin
-        @test all(1:1000) do _
+        for _ in 1:1000
             m, n = rand(100:1000), rand(100:1000)
             p = 0.05 * rand()
             A = sprand(m, n, p)
             S = SparsityPatternCSC(A)
             Sᵀ = transpose(S)
             Sᵀ_true = SparsityPatternCSC(sparse(transpose(A)))
-            Sᵀ.colptr == Sᵀ_true.colptr && Sᵀ.rowval == Sᵀ_true.rowval
+            @test Sᵀ.colptr == Sᵀ_true.colptr
+            @test Sᵀ.rowval == Sᵀ_true.rowval
         end
     end
     @testset "Bidirectional" begin
         @testset "symmetric_pattern = false" begin
-            m, n = rand(100:1000), rand(100:1000)
-            p = 0.05 * rand()
-            A = sprand(Bool, m, n, p)
-            A_and_Aᵀ = [spzeros(Bool, n, n) transpose(A); A spzeros(Bool, m, m)]
-            S_and_Sᵀ, edge_to_index = bidirectional_pattern(A; symmetric_pattern=false)
-            @test S_and_Sᵀ.colptr == A_and_Aᵀ.colptr
-            @test S_and_Sᵀ.rowval == A_and_Aᵀ.rowval
-            M = SparseMatrixCSC(
-                m + n, m + n, S_and_Sᵀ.colptr, S_and_Sᵀ.rowval, edge_to_index
-            )
-            @test issymmetric(M)
+            for _ in 1:1000
+                m, n = rand(100:1000), rand(100:1000)
+                p = 0.05 * rand()
+                A = sprand(Bool, m, n, p)
+                A_and_Aᵀ = [spzeros(Bool, n, n) transpose(A); A spzeros(Bool, m, m)]
+                S_and_Sᵀ, edge_to_index = bidirectional_pattern(A; symmetric_pattern=false)
+                @test S_and_Sᵀ.colptr == A_and_Aᵀ.colptr
+                @test S_and_Sᵀ.rowval == A_and_Aᵀ.rowval
+                M = SparseMatrixCSC(
+                    m + n, m + n, S_and_Sᵀ.colptr, S_and_Sᵀ.rowval, edge_to_index
+                )
+                @test issymmetric(M)
+            end
         end
         @testset "symmetric_pattern = true" begin
-            m = rand(100:1000)
-            p = 0.05 * rand()
-            A = sparse(Symmetric(sprand(Bool, m, m, p)))
-            A_and_Aᵀ = [spzeros(Bool, m, m) transpose(A); A spzeros(Bool, m, m)]
-            S_and_Sᵀ, edge_to_index = bidirectional_pattern(A; symmetric_pattern=true)
-            @test S_and_Sᵀ.colptr == A_and_Aᵀ.colptr
-            @test S_and_Sᵀ.rowval == A_and_Aᵀ.rowval
-            M = SparseMatrixCSC(
-                2 * m, 2 * m, S_and_Sᵀ.colptr, S_and_Sᵀ.rowval, edge_to_index
-            )
-            @test issymmetric(M)
+            for _ in 1:1000
+                m = rand(100:1000)
+                p = 0.05 * rand()
+                A = sparse(Symmetric(sprand(Bool, m, m, p)))
+                A_and_Aᵀ = [spzeros(Bool, m, m) transpose(A); A spzeros(Bool, m, m)]
+                S_and_Sᵀ, edge_to_index = bidirectional_pattern(A; symmetric_pattern=true)
+                @test S_and_Sᵀ.colptr == A_and_Aᵀ.colptr
+                @test S_and_Sᵀ.rowval == A_and_Aᵀ.rowval
+                M = SparseMatrixCSC(
+                    2 * m, 2 * m, S_and_Sᵀ.colptr, S_and_Sᵀ.rowval, edge_to_index
+                )
+                @test issymmetric(M)
+            end
         end
     end
     @testset "size" begin
@@ -137,15 +142,35 @@ end;
         0 0 0 1 1 1 1 0
     ])
 
-    B = transpose(A) * A
-    g = AdjacencyGraph(B - Diagonal(B))
+    g = AdjacencyGraph(transpose(A) * A)
     @test nb_vertices(g) == 8
-    @test collect(neighbors(g, 1)) == [6, 7, 8]
-    @test collect(neighbors(g, 2)) == [5, 7, 8]
-    @test collect(neighbors(g, 3)) == [5, 6, 8]
-    @test collect(neighbors(g, 4)) == [5, 6, 7]
-    @test collect(neighbors(g, 5)) == [2, 3, 4, 6, 7, 8]
-    @test collect(neighbors(g, 6)) == [1, 3, 4, 5, 7, 8]
-    @test collect(neighbors(g, 7)) == [1, 2, 4, 5, 6, 8]
-    @test collect(neighbors(g, 8)) == [1, 2, 3, 5, 6, 7]
+    # wrong neighbors, it's okay they are filtered after
+    @test collect(neighbors(g, 1)) == sort(vcat(1, [6, 7, 8]))
+    @test collect(neighbors(g, 2)) == sort(vcat(2, [5, 7, 8]))
+    @test collect(neighbors(g, 3)) == sort(vcat(3, [5, 6, 8]))
+    @test collect(neighbors(g, 4)) == sort(vcat(4, [5, 6, 7]))
+    @test collect(neighbors(g, 5)) == sort(vcat(5, [2, 3, 4, 6, 7, 8]))
+    @test collect(neighbors(g, 6)) == sort(vcat(6, [1, 3, 4, 5, 7, 8]))
+    @test collect(neighbors(g, 7)) == sort(vcat(7, [1, 2, 4, 5, 6, 8]))
+    @test collect(neighbors(g, 8)) == sort(vcat(8, [1, 2, 3, 5, 6, 7]))
+    # right degree
+    @test degree(g, 1) == 3
+    @test degree(g, 2) == 3
+    @test degree(g, 3) == 3
+    @test degree(g, 4) == 3
+    @test degree(g, 5) == 6
+    @test degree(g, 6) == 6
+    @test degree(g, 7) == 6
+    @test degree(g, 8) == 6
+
+    g = AdjacencyGraph(transpose(A) * A; has_diagonal=false)
+    # wrong degree
+    @test degree(g, 1) == 4
+    @test degree(g, 2) == 4
+    @test degree(g, 3) == 4
+    @test degree(g, 4) == 4
+    @test degree(g, 5) == 7
+    @test degree(g, 6) == 7
+    @test degree(g, 7) == 7
+    @test degree(g, 8) == 7
 end