More modifications

Author: amontoison

src/decompression.jl

@@ -537,6 +537,44 @@ function decompress!(
     return A
 end
 
+function decompress_single_color!(
+    A::SparseMatrixCSC,
+    b::AbstractVector,
+    c::Integer,
+    result::StarSetColoringResult,
+    uplo::Symbol=:F,
+)
+    (; ag, compressed_indices) = result
+    (; S) = ag
+    lower_index = (c - 1) * S.n + 1
+    upper_index = c * S.n
+    nzA = nonzeros(A)
+    if result.decompression_uplo == uplo
+        uplo == :F && check_same_pattern(A, S)
+        for k in eachindex(nzA, compressed_indices)
+            if lower_index <= compressed_indices[k] <= upper_index
+                nzA[k] = b[compressed_indices[k] - lower_index + 1]
+            end
+        end
+    else
+        @assert result.decompression_uplo == :F
+        rvS = rowvals(S)
+        l = 0  # assume A has the same pattern as the triangle
+        for j in axes(S, 2)
+            for k in nzrange(S, j)
+                i = rvS[k]
+                if in_triangle(i, j, uplo)
+                    l += 1
+                    if lower_index <= compressed_indices[k] <= upper_index
+                        nzA[l] = b[i]
+                    end
+                end
+            end
+        end
+    end
+    return A
+end
+
 ## TreeSetColoringResult
 
 function decompress!(

src/graph.jl

@@ -284,15 +284,7 @@ function degree(g::AdjacencyGraph{T,false}, v::Integer) where {T}
     return g.S.colptr[v + 1] - g.S.colptr[v] - has_selfloop
 end
 
-nb_edges(g::AdjacencyGraph{T,true}) where {T} = nnz(g.S) ÷ 2
-
-function nb_edges(g::AdjacencyGraph{T,false}) where {T}
-    ne = 0
-    for v in vertices(g)
-        ne += degree(g, v)
-    end
-    return ne ÷ 2
-end
+nb_edges(g::AdjacencyGraph)  = (nnz(g.S) - g.nb_self_loops) ÷ 2
 
 maximum_degree(g::AdjacencyGraph) = maximum(Base.Fix1(degree, g), vertices(g))
 minimum_degree(g::AdjacencyGraph) = minimum(Base.Fix1(degree, g), vertices(g))

src/interface.jl

@@ -190,10 +190,11 @@ function coloring(
     A::AbstractMatrix,
     problem::ColoringProblem,
     algo::GreedyColoringAlgorithm;
-    decompression_eltype::Type{R}=Float64,
     symmetric_pattern::Bool=false,
+    decompression_eltype::Type{R}=Float64,
+    decompression_uplo::Symbol=:F,
 ) where {R}
-    return _coloring(WithResult(), A, problem, algo, R, symmetric_pattern)
+    return _coloring(WithResult(), A, problem, algo, symmetric_pattern, R, decompression_uplo)
 end
 
 """
@@ -229,8 +230,9 @@ function _coloring(
     A::AbstractMatrix,
     ::ColoringProblem{:nonsymmetric,:column},
     algo::GreedyColoringAlgorithm,
+    symmetric_pattern::Bool,
     decompression_eltype::Type,
-    symmetric_pattern::Bool;
+    decompression_uplo::Symbol;
     forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing,
 )
     symmetric_pattern = symmetric_pattern || A isa Union{Symmetric,Hermitian}
@@ -252,8 +254,9 @@ function _coloring(
     A::AbstractMatrix,
     ::ColoringProblem{:nonsymmetric,:row},
     algo::GreedyColoringAlgorithm,
+    symmetric_pattern::Bool,
     decompression_eltype::Type,
-    symmetric_pattern::Bool;
+    decompression_uplo::Symbol;
     forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing,
 )
     symmetric_pattern = symmetric_pattern || A isa Union{Symmetric,Hermitian}
@@ -275,8 +278,9 @@ function _coloring(
     A::AbstractMatrix,
     ::ColoringProblem{:symmetric,:column},
     algo::GreedyColoringAlgorithm{:direct},
+    symmetric_pattern::Bool,
     decompression_eltype::Type,
-    symmetric_pattern::Bool;
+    decompression_uplo::Symbol;
     forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing,
 )
     ag = AdjacencyGraph(A; augmented_graph=false)
@@ -286,7 +290,7 @@ function _coloring(
     end
     color, star_set = argmin(maximum ∘ first, color_and_star_set_by_order)
     if speed_setting isa WithResult
-        return StarSetColoringResult(A, ag, color, star_set)
+        return StarSetColoringResult(A, ag, color, star_set; decompression_uplo)
     else
         return color
     end
@@ -297,8 +301,9 @@ function _coloring(
     A::AbstractMatrix,
     ::ColoringProblem{:symmetric,:column},
     algo::GreedyColoringAlgorithm{:substitution},
-    decompression_eltype::Type{R},
     symmetric_pattern::Bool,
+    decompression_eltype::Type{R},
+    decompression_uplo::Symbol,
 ) where {R}
     ag = AdjacencyGraph(A; augmented_graph=false)
     color_and_tree_set_by_order = map(algo.orders) do order
@@ -307,7 +312,7 @@ function _coloring(
     end
     color, tree_set = argmin(maximum ∘ first, color_and_tree_set_by_order)
     if speed_setting isa WithResult
-        return TreeSetColoringResult(A, ag, color, tree_set, R)
+        return TreeSetColoringResult(A, ag, color, tree_set, R; decompression_uplo)
     else
         return color
     end
@@ -318,8 +323,9 @@ function _coloring(
     A::AbstractMatrix,
     ::ColoringProblem{:nonsymmetric,:bidirectional},
     algo::GreedyColoringAlgorithm{:direct},
+    symmetric_pattern::Bool,
     decompression_eltype::Type{R},
-    symmetric_pattern::Bool;
+    decompression_uplo::Symbol;
     forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing,
 ) where {R}
     A_and_Aᵀ, edge_to_index = bidirectional_pattern(A; symmetric_pattern)
@@ -368,8 +374,9 @@ function _coloring(
     A::AbstractMatrix,
     ::ColoringProblem{:nonsymmetric,:bidirectional},
     algo::GreedyColoringAlgorithm{:substitution},
-    decompression_eltype::Type{R},
     symmetric_pattern::Bool,
+    decompression_eltype::Type{R},
+    decompression_uplo::Symbol,
 ) where {R}
     A_and_Aᵀ, edge_to_index = bidirectional_pattern(A; symmetric_pattern)
     ag = AdjacencyGraph(A_and_Aᵀ, edge_to_index, 0; augmented_graph=true)

src/result.jl

@@ -339,7 +339,7 @@ function star_csc_indices(
     rvS = rowvals(S)
     nb_indices = nnz(S)
     if decompression_uplo != :F
-        nb_indices = (nb_indices - nb_self_loops) ÷ 2 + nb_self_loops
+        nb_indices = nb_edges(ag) + 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