From 7bfb9f528dfafdf7c75d77b9f75b908630bb1ab0 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Tue, 5 Aug 2025 09:30:35 +0200 Subject: [PATCH 01/39] Bump julia-actions/julia-downgrade-compat from 1 to 2 (#258) Bumps [julia-actions/julia-downgrade-compat](https://github.com/julia-actions/julia-downgrade-compat) from 1 to 2. - [Release notes](https://github.com/julia-actions/julia-downgrade-compat/releases) - [Commits](https://github.com/julia-actions/julia-downgrade-compat/compare/v1...v2) --- updated-dependencies: - dependency-name: julia-actions/julia-downgrade-compat dependency-version: '2' dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/Test-GPU.yml | 2 +- .github/workflows/Test.yml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/Test-GPU.yml b/.github/workflows/Test-GPU.yml index cd67c6e1..3557a520 100644 --- a/.github/workflows/Test-GPU.yml +++ b/.github/workflows/Test-GPU.yml @@ -31,7 +31,7 @@ jobs: with: version: ${{ matrix.julia-version }} arch: x64 - - uses: julia-actions/julia-downgrade-compat@v1 + - uses: julia-actions/julia-downgrade-compat@v2 if: ${{ matrix.version == '1.10' }} with: skip: LinearAlgebra, Random, SparseArrays diff --git a/.github/workflows/Test.yml b/.github/workflows/Test.yml index 1f0046a1..bed25d9d 100644 --- a/.github/workflows/Test.yml +++ b/.github/workflows/Test.yml @@ -28,7 +28,7 @@ jobs: with: version: ${{ matrix.julia-version }} arch: x64 - - uses: julia-actions/julia-downgrade-compat@v1 + - uses: julia-actions/julia-downgrade-compat@v2 if: ${{ matrix.version == '1.10' }} with: skip: LinearAlgebra, Random, SparseArrays From cbed1d58fa3d25eacb35ff6ef1f282d5b3d06f6d Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 18 Aug 2025 07:58:09 +0200 Subject: [PATCH 02/39] Bump actions/checkout from 4 to 5 (#259) --- .github/workflows/Benchmark.yml | 2 +- .github/workflows/Documentation.yml | 2 +- .github/workflows/Test-GPU.yml | 2 +- .github/workflows/Test.yml | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/.github/workflows/Benchmark.yml b/.github/workflows/Benchmark.yml index fda51ed6..12de3ffd 100644 --- a/.github/workflows/Benchmark.yml +++ b/.github/workflows/Benchmark.yml @@ -12,7 +12,7 @@ jobs: runs-on: ubuntu-latest if: contains(github.event.pull_request.labels.*.name, 'benchmark') steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 - uses: julia-actions/setup-julia@v2 with: version: "1" diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml index a8d1110e..c5527ae9 100644 --- a/.github/workflows/Documentation.yml +++ b/.github/workflows/Documentation.yml @@ -19,7 +19,7 @@ jobs: statuses: write runs-on: ubuntu-latest steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 - uses: julia-actions/setup-julia@v2 with: version: '1' diff --git a/.github/workflows/Test-GPU.yml b/.github/workflows/Test-GPU.yml index 3557a520..05628de0 100644 --- a/.github/workflows/Test-GPU.yml +++ b/.github/workflows/Test-GPU.yml @@ -26,7 +26,7 @@ jobs: matrix: julia-version: ['1.10', '1'] steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.julia-version }} diff --git a/.github/workflows/Test.yml b/.github/workflows/Test.yml index bed25d9d..ef161b41 100644 --- a/.github/workflows/Test.yml +++ b/.github/workflows/Test.yml @@ -23,7 +23,7 @@ jobs: julia-version: ['1.10', '1'] steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.julia-version }} From ccfcd351c26152683a028a75619521ad7914dab4 Mon Sep 17 00:00:00 2001 From: Neven Sajko <4944410+nsajko@users.noreply.github.com> Date: Thu, 11 Sep 2025 09:09:31 +0200 Subject: [PATCH 03/39] ensure `UnsupportedDecompressionError <: Exception` (#260) Not sure if there is any benefit, but I suppose the intention is for this type, as an exception type, to subtype `Exception`. --- src/decompression.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/decompression.jl b/src/decompression.jl index a3dd471a..a48c7e34 100644 --- a/src/decompression.jl +++ b/src/decompression.jl @@ -106,7 +106,7 @@ function compress( return Br, Bc end -struct UnsupportedDecompressionError +struct UnsupportedDecompressionError <: Exception msg::String end From 56a0c4f8025743bb4933c4f3075bda14525463c7 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Sun, 5 Oct 2025 23:12:19 -0500 Subject: [PATCH 04/39] Bump peter-evans/create-or-update-comment from 4 to 5 (#261) --- .github/workflows/Benchmark.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/Benchmark.yml b/.github/workflows/Benchmark.yml index 12de3ffd..3c7d70b8 100644 --- a/.github/workflows/Benchmark.yml +++ b/.github/workflows/Benchmark.yml @@ -53,7 +53,7 @@ jobs: comment-author: 'github-actions[bot]' body-includes: Benchmark Results - name: Comment on PR - uses: peter-evans/create-or-update-comment@v4 + uses: peter-evans/create-or-update-comment@v5 with: comment-id: ${{ steps.fcbenchmark.outputs.comment-id }} issue-number: ${{ github.event.pull_request.number }} From 7efc6bc7650c442c326aa40fe3063d27f6bd8ea5 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 6 Oct 2025 08:30:15 +0200 Subject: [PATCH 05/39] Bump peter-evans/find-comment from 3 to 4 (#262) Bumps [peter-evans/find-comment](https://github.com/peter-evans/find-comment) from 3 to 4. - [Release notes](https://github.com/peter-evans/find-comment/releases) - [Commits](https://github.com/peter-evans/find-comment/compare/v3...v4) --- updated-dependencies: - dependency-name: peter-evans/find-comment dependency-version: '4' dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/Benchmark.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/Benchmark.yml b/.github/workflows/Benchmark.yml index 3c7d70b8..23dd4973 100644 --- a/.github/workflows/Benchmark.yml +++ b/.github/workflows/Benchmark.yml @@ -46,7 +46,7 @@ jobs: echo '' >> body.md echo '' >> body.md - name: Find Comment - uses: peter-evans/find-comment@v3 + uses: peter-evans/find-comment@v4 id: fcbenchmark with: issue-number: ${{ github.event.pull_request.number }} From b3cf23ded1060e48921babd8d4529a01ccd36048 Mon Sep 17 00:00:00 2001 From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com> Date: Tue, 14 Oct 2025 07:17:24 +0200 Subject: [PATCH 06/39] Allow picking the best of several orders for GreedyColoringAlgorithm (#265) * Allow picking the best of several orders for GreedyColoringAlgorithm * Fix type params * Tuple * Better tests * Better test * Fix foc * Fix type inference inside closure * Fix seed * Test on 1.11 * Avoid duplicate remap_colors --- .github/workflows/Test.yml | 4 +- docs/src/tutorial.md | 2 +- src/interface.jl | 125 +++++++++++++++++++++++++++---------- src/result.jl | 7 ++- test/order.jl | 96 ++++++++++++++++++++++++++++ test/random.jl | 8 +-- test/type_stability.jl | 10 +++ test/utils.jl | 36 +++++++++++ 8 files changed, 246 insertions(+), 42 deletions(-) diff --git a/.github/workflows/Test.yml b/.github/workflows/Test.yml index ef161b41..3e17ae97 100644 --- a/.github/workflows/Test.yml +++ b/.github/workflows/Test.yml @@ -20,7 +20,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - julia-version: ['1.10', '1'] + julia-version: ['1.10', '1.11'] steps: - uses: actions/checkout@v5 @@ -40,4 +40,4 @@ jobs: with: files: lcov.info token: ${{ secrets.CODECOV_TOKEN }} - fail_ci_if_error: false \ No newline at end of file + fail_ci_if_error: false diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md index fec54f13..6e803511 100644 --- a/docs/src/tutorial.md +++ b/docs/src/tutorial.md @@ -31,7 +31,7 @@ problem = ColoringProblem() The algorithm defines how you want to solve it. It can be either a [`GreedyColoringAlgorithm`](@ref) or a [`ConstantColoringAlgorithm`](@ref). For `GreedyColoringAlgorithm`, you can select options such as -- the order in which vertices are processed (a subtype of [`AbstractOrder`](@ref SparseMatrixColorings.AbstractOrder)) +- the order in which vertices are processed (a subtype of [`AbstractOrder`](@ref SparseMatrixColorings.AbstractOrder) , or a tuple of such objects) - the type of decompression you want (`:direct` or `:substitution`) ```@example tutorial diff --git a/src/interface.jl b/src/interface.jl index 79fe71e9..c1a56c29 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -72,7 +72,7 @@ It is passed as an argument to the main function [`coloring`](@ref). GreedyColoringAlgorithm{decompression}(order=NaturalOrder(); postprocessing=false) GreedyColoringAlgorithm(order=NaturalOrder(); postprocessing=false, decompression=:direct) -- `order::AbstractOrder`: the order in which the columns or rows are colored, which can impact the number of colors. +- `order::Union{AbstractOrder,Tuple}`: the order in which the columns or rows are colored, which can impact the number of colors. Can also be a tuple of different orders to try out, from which the best order (the one with the lowest total number of colors) will be used. - `postprocessing::Bool`: whether or not the coloring will be refined by assigning the neutral color `0` to some vertices. - `decompression::Symbol`: either `:direct` or `:substitution`. Usually `:substitution` leads to fewer colors, at the cost of a more expensive coloring (and decompression). When `:substitution` is not applicable, it falls back on `:direct` decompression. @@ -94,26 +94,31 @@ See their respective docstrings for details. - [`AbstractOrder`](@ref) - [`decompress`](@ref) """ -struct GreedyColoringAlgorithm{decompression,O<:AbstractOrder} <: +struct GreedyColoringAlgorithm{decompression,N,O<:NTuple{N,AbstractOrder}} <: ADTypes.AbstractColoringAlgorithm - order::O + orders::O postprocessing::Bool -end -function GreedyColoringAlgorithm{decompression}( - order::AbstractOrder=NaturalOrder(); postprocessing::Bool=false -) where {decompression} - check_valid_algorithm(decompression) - return GreedyColoringAlgorithm{decompression,typeof(order)}(order, postprocessing) + function GreedyColoringAlgorithm{decompression}( + order_or_orders::Union{AbstractOrder,Tuple}=NaturalOrder(); + postprocessing::Bool=false, + ) where {decompression} + check_valid_algorithm(decompression) + if order_or_orders isa AbstractOrder + orders = (order_or_orders,) + else + orders = order_or_orders + end + return new{decompression,length(orders),typeof(orders)}(orders, postprocessing) + end end function GreedyColoringAlgorithm( - order::AbstractOrder=NaturalOrder(); + order_or_orders::Union{AbstractOrder,Tuple}=NaturalOrder(); postprocessing::Bool=false, decompression::Symbol=:direct, ) - check_valid_algorithm(decompression) - return GreedyColoringAlgorithm{decompression,typeof(order)}(order, postprocessing) + return GreedyColoringAlgorithm{decompression}(order_or_orders; postprocessing) end ## Coloring @@ -229,8 +234,11 @@ function _coloring( ) symmetric_pattern = symmetric_pattern || A isa Union{Symmetric,Hermitian} bg = BipartiteGraph(A; symmetric_pattern) - vertices_in_order = vertices(bg, Val(2), algo.order) - color = partial_distance2_coloring(bg, Val(2), vertices_in_order) + color_by_order = map(algo.orders) do order + vertices_in_order = vertices(bg, Val(2), order) + return partial_distance2_coloring(bg, Val(2), vertices_in_order) + end + color = argmin(maximum, color_by_order) if speed_setting isa WithResult return ColumnColoringResult(A, bg, color) else @@ -248,8 +256,11 @@ function _coloring( ) symmetric_pattern = symmetric_pattern || A isa Union{Symmetric,Hermitian} bg = BipartiteGraph(A; symmetric_pattern) - vertices_in_order = vertices(bg, Val(1), algo.order) - color = partial_distance2_coloring(bg, Val(1), vertices_in_order) + color_by_order = map(algo.orders) do order + vertices_in_order = vertices(bg, Val(1), order) + return partial_distance2_coloring(bg, Val(1), vertices_in_order) + end + color = argmin(maximum, color_by_order) if speed_setting isa WithResult return RowColoringResult(A, bg, color) else @@ -266,8 +277,11 @@ function _coloring( symmetric_pattern::Bool, ) ag = AdjacencyGraph(A; has_diagonal=true) - vertices_in_order = vertices(ag, algo.order) - color, star_set = star_coloring(ag, vertices_in_order, algo.postprocessing) + color_and_star_set_by_order = map(algo.orders) do order + vertices_in_order = vertices(ag, order) + return star_coloring(ag, vertices_in_order, algo.postprocessing) + 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) else @@ -284,8 +298,11 @@ function _coloring( symmetric_pattern::Bool, ) where {R} ag = AdjacencyGraph(A; has_diagonal=true) - vertices_in_order = vertices(ag, algo.order) - color, tree_set = acyclic_coloring(ag, vertices_in_order, algo.postprocessing) + color_and_tree_set_by_order = map(algo.orders) do order + vertices_in_order = vertices(ag, order) + return acyclic_coloring(ag, vertices_in_order, algo.postprocessing) + 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) else @@ -303,15 +320,37 @@ function _coloring( ) where {R} A_and_Aᵀ, edge_to_index = bidirectional_pattern(A; symmetric_pattern) ag = AdjacencyGraph(A_and_Aᵀ, edge_to_index; has_diagonal=false) - vertices_in_order = vertices(ag, algo.order) - color, star_set = star_coloring(ag, vertices_in_order, algo.postprocessing) + outputs_by_order = map(algo.orders) do order + vertices_in_order = vertices(ag, order) + _color, _star_set = star_coloring(ag, vertices_in_order, algo.postprocessing) + (_row_color, _column_color, _symmetric_to_row, _symmetric_to_column) = remap_colors( + eltype(ag), _color, maximum(_color), size(A)... + ) + return ( + _color, + _star_set, + _row_color, + _column_color, + _symmetric_to_row, + _symmetric_to_column, + ) + end + (color, star_set, row_color, column_color, symmetric_to_row, symmetric_to_column) = argmin( + t -> maximum(t[3]) + maximum(t[4]), outputs_by_order + ) # can't use ncolors without computing the full result if speed_setting isa WithResult symmetric_result = StarSetColoringResult(A_and_Aᵀ, ag, color, star_set) - return BicoloringResult(A, ag, symmetric_result, R) - else - row_color, column_color, _ = remap_colors( - eltype(ag), color, maximum(color), size(A)... + return BicoloringResult( + A, + ag, + symmetric_result, + row_color, + column_color, + symmetric_to_row, + symmetric_to_column, + R, ) + else return row_color, column_color end end @@ -326,15 +365,37 @@ function _coloring( ) where {R} A_and_Aᵀ, edge_to_index = bidirectional_pattern(A; symmetric_pattern) ag = AdjacencyGraph(A_and_Aᵀ, edge_to_index; has_diagonal=false) - vertices_in_order = vertices(ag, algo.order) - color, tree_set = acyclic_coloring(ag, vertices_in_order, algo.postprocessing) + 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) + (_row_color, _column_color, _symmetric_to_row, _symmetric_to_column) = remap_colors( + eltype(ag), _color, maximum(_color), size(A)... + ) + return ( + _color, + _tree_set, + _row_color, + _column_color, + _symmetric_to_row, + _symmetric_to_column, + ) + end + (color, tree_set, row_color, column_color, symmetric_to_row, symmetric_to_column) = argmin( + t -> maximum(t[3]) + maximum(t[4]), outputs_by_order + ) # can't use ncolors without computing the full result if speed_setting isa WithResult symmetric_result = TreeSetColoringResult(A_and_Aᵀ, ag, color, tree_set, R) - return BicoloringResult(A, ag, symmetric_result, R) - else - row_color, column_color, _ = remap_colors( - eltype(ag), color, maximum(color), size(A)... + return BicoloringResult( + A, + ag, + symmetric_result, + row_color, + column_color, + symmetric_to_row, + symmetric_to_column, + R, ) + else return row_color, column_color end end diff --git a/src/result.jl b/src/result.jl index b5b9cdd2..7a3e3b68 100644 --- a/src/result.jl +++ b/src/result.jl @@ -686,14 +686,15 @@ function BicoloringResult( A::AbstractMatrix, ag::AdjacencyGraph{T}, symmetric_result::AbstractColoringResult{:symmetric,:column}, + row_color::Vector{T}, + column_color::Vector{T}, + symmetric_to_row::Vector{T}, + symmetric_to_column::Vector{T}, decompression_eltype::Type{R}, ) where {T,R} m, n = size(A) symmetric_color = column_colors(symmetric_result) num_sym_colors = maximum(symmetric_color) - row_color, column_color, symmetric_to_row, symmetric_to_column = remap_colors( - T, symmetric_color, num_sym_colors, m, n - ) column_group = group_by_color(T, column_color) row_group = group_by_color(T, row_color) Br_and_Bc = Matrix{R}(undef, n + m, num_sym_colors) diff --git a/test/order.jl b/test/order.jl index 9bb9e6dc..838c32a9 100644 --- a/test/order.jl +++ b/test/order.jl @@ -146,3 +146,99 @@ end; @test isperm(π) end end + +@testset "Multiple orders" begin + # I used brute force to find examples where LargestFirst is *strictly* better than NaturalOrder, just to check that the best order is indeed selected when multiple orders are provided + @testset "Column coloring" begin + A = [ + 0 0 1 1 + 0 1 0 1 + 0 0 1 1 + 1 1 0 0 + ] + problem = ColoringProblem{:nonsymmetric,:column}() + algo = GreedyColoringAlgorithm(NaturalOrder()) + better_algo = GreedyColoringAlgorithm((NaturalOrder(), LargestFirst())) + @test ncolors(coloring(A, problem, better_algo)) < + ncolors(coloring(A, problem, algo)) + end + @testset "Row coloring" begin + A = [ + 1 0 0 0 + 0 0 1 0 + 0 1 1 1 + 1 0 0 1 + ] + problem = ColoringProblem{:nonsymmetric,:row}() + algo = GreedyColoringAlgorithm(NaturalOrder()) + better_algo = GreedyColoringAlgorithm((NaturalOrder(), LargestFirst())) + @test ncolors(coloring(A, problem, better_algo)) < + ncolors(coloring(A, problem, algo)) + end + @testset "Star coloring" begin + A = [ + 0 1 0 1 1 + 1 1 0 1 0 + 0 0 1 0 1 + 1 1 0 1 0 + 1 0 1 0 0 + ] + problem = ColoringProblem{:symmetric,:column}() + algo = GreedyColoringAlgorithm(NaturalOrder()) + better_algo = GreedyColoringAlgorithm((NaturalOrder(), LargestFirst())) + @test ncolors(coloring(A, problem, better_algo)) < + ncolors(coloring(A, problem, algo)) + end + @testset "Acyclic coloring" begin + A = [ + 1 0 0 0 0 1 0 + 0 0 0 1 0 0 0 + 0 0 0 1 0 0 0 + 0 1 1 1 0 1 1 + 0 0 0 0 0 0 1 + 1 0 0 1 0 0 1 + 0 0 0 1 1 1 1 + ] + problem = ColoringProblem{:symmetric,:column}() + algo = GreedyColoringAlgorithm{:substitution}(NaturalOrder()) + better_algo = GreedyColoringAlgorithm{:substitution}(( + NaturalOrder(), LargestFirst() + )) + @test ncolors(coloring(A, problem, better_algo)) < + ncolors(coloring(A, problem, algo)) + end + @testset "Star bicoloring" begin + A = [ + 0 1 0 0 0 + 1 0 1 0 0 + 0 1 0 0 1 + 0 0 0 0 0 + 0 0 1 0 1 + ] + problem = ColoringProblem{:nonsymmetric,:bidirectional}() + algo = GreedyColoringAlgorithm(NaturalOrder()) + better_algo = GreedyColoringAlgorithm((NaturalOrder(), LargestFirst())) + @test ncolors(coloring(A, problem, better_algo)) < + ncolors(coloring(A, problem, algo)) + end + @testset "Acyclic bicoloring" begin + A = [ + 0 1 0 1 1 0 1 0 1 + 1 0 0 0 0 0 0 0 1 + 0 0 0 0 0 0 0 0 0 + 1 0 0 1 1 0 1 0 0 + 1 0 0 1 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 + 1 0 0 1 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 + 1 1 0 0 0 0 0 0 0 + ] + problem = ColoringProblem{:nonsymmetric,:bidirectional}() + algo = GreedyColoringAlgorithm{:substitution}(NaturalOrder()) + better_algo = GreedyColoringAlgorithm{:substitution}(( + NaturalOrder(), LargestFirst() + )) + @test ncolors(coloring(A, problem, better_algo)) < + ncolors(coloring(A, problem, algo)) + end +end diff --git a/test/random.jl b/test/random.jl index 02c1a3a1..406dfea9 100644 --- a/test/random.jl +++ b/test/random.jl @@ -81,10 +81,10 @@ end; problem = ColoringProblem(; structure=:nonsymmetric, partition=:bidirectional) @testset for algo in ( GreedyColoringAlgorithm( - RandomOrder(rng); postprocessing=false, decompression=:direct + RandomOrder(StableRNG(0), 0); postprocessing=false, decompression=:direct ), GreedyColoringAlgorithm( - RandomOrder(rng); postprocessing=true, decompression=:direct + RandomOrder(StableRNG(0), 0); postprocessing=true, decompression=:direct ), ) @testset "$((; m, n, p))" for (m, n, p) in asymmetric_params @@ -102,10 +102,10 @@ end; problem = ColoringProblem(; structure=:nonsymmetric, partition=:bidirectional) @testset for algo in ( GreedyColoringAlgorithm( - RandomOrder(rng); postprocessing=false, decompression=:substitution + RandomOrder(StableRNG(0), 0); postprocessing=false, decompression=:substitution ), GreedyColoringAlgorithm( - RandomOrder(rng); postprocessing=true, decompression=:substitution + RandomOrder(StableRNG(0), 0); postprocessing=true, decompression=:substitution ), ) @testset "$((; m, n, p))" for (m, n, p) in asymmetric_params diff --git a/test/type_stability.jl b/test/type_stability.jl index 5a488596..3c6f9f02 100644 --- a/test/type_stability.jl +++ b/test/type_stability.jl @@ -40,11 +40,21 @@ rng = StableRNG(63) ColoringProblem(; structure, partition), GreedyColoringAlgorithm(order; decompression), ) + @test_opt coloring( + A, + ColoringProblem(; structure, partition), + GreedyColoringAlgorithm((NaturalOrder(), order); decompression), + ) @inferred coloring( A, ColoringProblem(; structure, partition), GreedyColoringAlgorithm(order; decompression), ) + @inferred coloring( + A, + ColoringProblem(; structure, partition), + GreedyColoringAlgorithm((NaturalOrder(), order); decompression), + ) end end end; diff --git a/test/utils.jl b/test/utils.jl index bb80f95f..2462f8c1 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -15,6 +15,10 @@ using SparseMatrixColorings: structurally_biorthogonal using Test +const _ALL_ORDERS = ( + NaturalOrder(), LargestFirst(), SmallestLast(), IncidenceDegree(), DynamicLargestFirst() +) + function test_coloring_decompression( A0::AbstractMatrix, problem::ColoringProblem{structure,partition}, @@ -169,6 +173,22 @@ function test_coloring_decompression( @show color_vec end end + + @testset "More orders is better" begin + more_orders = (algo.orders..., _ALL_ORDERS...) + better_algo = GreedyColoringAlgorithm{decompression}( + more_orders; algo.postprocessing + ) + all_algos = [ + GreedyColoringAlgorithm{decompression}(order; algo.postprocessing) for + order in more_orders + ] + result = coloring(A0, problem, algo) + better_result = coloring(A0, problem, better_algo) + all_results = [coloring(A0, problem, _algo) for _algo in all_algos] + @test ncolors(better_result) <= ncolors(result) + @test ncolors(better_result) == minimum(ncolors, all_results) + end end function test_bicoloring_decompression( @@ -216,6 +236,22 @@ function test_bicoloring_decompression( end end end + + @testset "More orders is better" begin + more_orders = (algo.orders..., _ALL_ORDERS...) + better_algo = GreedyColoringAlgorithm{decompression}( + more_orders; algo.postprocessing + ) + all_algos = [ + GreedyColoringAlgorithm{decompression}(order; algo.postprocessing) for + order in more_orders + ] + result = coloring(A0, problem, algo) + better_result = coloring(A0, problem, better_algo) + all_results = [coloring(A0, problem, _algo) for _algo in all_algos] + @test ncolors(better_result) <= ncolors(result) + @test ncolors(better_result) == minimum(ncolors, all_results) + end end function test_structured_coloring_decompression(A::AbstractMatrix) From 734e953d1dbc47449e2f76b36340f0beb0dca5c2 Mon Sep 17 00:00:00 2001 From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com> Date: Tue, 14 Oct 2025 09:26:45 +0200 Subject: [PATCH 07/39] chore: bump version to 0.4.22 (#268) --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 6656e125..003eab87 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SparseMatrixColorings" uuid = "0a514795-09f3-496d-8182-132a7b665d35" authors = ["Guillaume Dalle", "Alexis Montoison"] -version = "0.4.21" +version = "0.4.22" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" @@ -14,7 +14,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [weakdeps] CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8" -Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" +Colors = "5ae59095-9ag9b-59fe-a467-6f913c188581" [extensions] SparseMatrixColoringsCUDAExt = "CUDA" From e77226af7dcc1e180395b3dd2eef6397172349ab Mon Sep 17 00:00:00 2001 From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com> Date: Tue, 14 Oct 2025 09:29:10 +0200 Subject: [PATCH 08/39] chore: fix typo in package UUID (#269) --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 003eab87..3eb2711d 100644 --- a/Project.toml +++ b/Project.toml @@ -14,7 +14,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [weakdeps] CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8" -Colors = "5ae59095-9ag9b-59fe-a467-6f913c188581" +Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" [extensions] SparseMatrixColoringsCUDAExt = "CUDA" From b53dbf9849f33b7830d6f9e219f42916d7536f21 Mon Sep 17 00:00:00 2001 From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com> Date: Wed, 22 Oct 2025 18:48:41 +0200 Subject: [PATCH 09/39] Make any coloring algorithm compatible with SMC (#263) * Make any coloring algorithm compatible with SMC * Fix StackOverflow * Run tests on 1.11 * Start working on optimal coloring * Fix tests * Improve feasibility check * Format * Fix * Fix * Typo * Use HiGHS * Remove OptimalColoringAlgorithm * Remove test deps * Don't handle forced colors in acyclic * Run GPU CI on 1.11 too * Actually use forced colors * Format * Deactivate JuliaFormatter (can't figure out why it fails) --- .github/workflows/Test-GPU.yml | 2 +- src/adtypes.jl | 46 ++++++++++------- src/coloring.jl | 70 ++++++++++++++++++++------ src/constant.jl | 91 +++++++++++++++++----------------- src/interface.jl | 22 +++++--- test/adtypes.jl | 61 ++++++++++++++++------- test/constant.jl | 39 +++++++++++---- test/result.jl | 3 +- test/runtests.jl | 11 ++-- 9 files changed, 226 insertions(+), 119 deletions(-) diff --git a/.github/workflows/Test-GPU.yml b/.github/workflows/Test-GPU.yml index 05628de0..f5a4fbc8 100644 --- a/.github/workflows/Test-GPU.yml +++ b/.github/workflows/Test-GPU.yml @@ -24,7 +24,7 @@ jobs: JULIA_SMC_TEST_GROUP: "GPU" strategy: matrix: - julia-version: ['1.10', '1'] + julia-version: ['1.10', '1.11'] steps: - uses: actions/checkout@v5 - uses: julia-actions/setup-julia@v2 diff --git a/src/adtypes.jl b/src/adtypes.jl index 7500b1e2..7c464433 100644 --- a/src/adtypes.jl +++ b/src/adtypes.jl @@ -1,21 +1,33 @@ -function coloring( - A::AbstractMatrix, - ::ColoringProblem{:nonsymmetric,:column}, - algo::ADTypes.NoColoringAlgorithm; - kwargs..., -) - bg = BipartiteGraph(A) - color = convert(Vector{eltype(bg)}, ADTypes.column_coloring(A, algo)) - return ColumnColoringResult(A, bg, color) -end +## From ADTypes to SMC function coloring( A::AbstractMatrix, - ::ColoringProblem{:nonsymmetric,:row}, - algo::ADTypes.NoColoringAlgorithm; - kwargs..., -) - bg = BipartiteGraph(A) - color = convert(Vector{eltype(bg)}, ADTypes.row_coloring(A, algo)) - return RowColoringResult(A, bg, color) + problem::ColoringProblem{structure,partition}, + algo::ADTypes.AbstractColoringAlgorithm; + decompression_eltype::Type{R}=Float64, + symmetric_pattern::Bool=false, +) where {structure,partition,R} + symmetric_pattern = symmetric_pattern || A isa Union{Symmetric,Hermitian} + if structure == :nonsymmetric + if partition == :column + forced_colors = ADTypes.column_coloring(A, algo) + elseif partition == :row + forced_colors = ADTypes.row_coloring(A, algo) + else + # TODO: improve once https://github.com/SciML/ADTypes.jl/issues/69 is done + A_and_Aᵀ, _ = bidirectional_pattern(A; symmetric_pattern) + forced_colors = ADTypes.symmetric_coloring(A_and_Aᵀ, algo) + end + else + forced_colors = ADTypes.symmetric_coloring(A, algo) + end + return _coloring( + WithResult(), + A, + problem, + GreedyColoringAlgorithm(), + R, + symmetric_pattern; + forced_colors, + ) end diff --git a/src/coloring.jl b/src/coloring.jl index 7eec55c0..7288b773 100644 --- a/src/coloring.jl +++ b/src/coloring.jl @@ -1,5 +1,10 @@ +struct InvalidColoringError <: Exception end + """ - partial_distance2_coloring(bg::BipartiteGraph, ::Val{side}, vertices_in_order::AbstractVector) + partial_distance2_coloring( + bg::BipartiteGraph, ::Val{side}, vertices_in_order::AbstractVector; + forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing + ) Compute a distance-2 coloring of the given `side` (`1` or `2`) in the bipartite graph `bg` and return a vector of integer colors. @@ -7,6 +12,8 @@ A _distance-2 coloring_ is such that two vertices have different colors if they The vertices are colored in a greedy fashion, following the order supplied. +The optional `forced_colors` keyword argument is used to enforce predefined vertex colors (e.g. coming from another optimization algorithm) but still run the distance-2 coloring procedure to verify correctness. + # See also - [`BipartiteGraph`](@ref) @@ -17,11 +24,16 @@ The vertices are colored in a greedy fashion, following the order supplied. > [_What Color Is Your Jacobian? Graph Coloring for Computing Derivatives_](https://epubs.siam.org/doi/10.1137/S0036144504444711), Gebremedhin et al. (2005), Algorithm 3.2 """ function partial_distance2_coloring( - bg::BipartiteGraph{T}, ::Val{side}, vertices_in_order::AbstractVector{<:Integer} + bg::BipartiteGraph{T}, + ::Val{side}, + vertices_in_order::AbstractVector{<:Integer}; + forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing, ) where {T,side} color = Vector{T}(undef, nb_vertices(bg, Val(side))) forbidden_colors = Vector{T}(undef, nb_vertices(bg, Val(side))) - partial_distance2_coloring!(color, forbidden_colors, bg, Val(side), vertices_in_order) + partial_distance2_coloring!( + color, forbidden_colors, bg, Val(side), vertices_in_order; forced_colors + ) return color end @@ -30,7 +42,8 @@ function partial_distance2_coloring!( forbidden_colors::AbstractVector{<:Integer}, bg::BipartiteGraph, ::Val{side}, - vertices_in_order::AbstractVector{<:Integer}, + vertices_in_order::AbstractVector{<:Integer}; + forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing, ) where {side} color .= 0 forbidden_colors .= 0 @@ -44,17 +57,32 @@ function partial_distance2_coloring!( end end end - for i in eachindex(forbidden_colors) - if forbidden_colors[i] != v - color[v] = i - break + if isnothing(forced_colors) + for i in eachindex(forbidden_colors) + if forbidden_colors[i] != v + color[v] = i + break + end + end + else + f = forced_colors[v] + if ( + (f == 0 && length(neighbors(bg, Val(side), v)) > 0) || + (f > 0 && forbidden_colors[f] == v) + ) + throw(InvalidColoringError()) + else + color[v] = f end end end end """ - star_coloring(g::AdjacencyGraph, vertices_in_order::AbstractVector, postprocessing::Bool) + star_coloring( + g::AdjacencyGraph, vertices_in_order::AbstractVector, postprocessing::Bool; + forced_colors::Union{AbstractVector,Nothing}=nothing + ) Compute a star coloring of all vertices in the adjacency graph `g` and return a tuple `(color, star_set)`, where @@ -67,6 +95,8 @@ The vertices are colored in a greedy fashion, following the order supplied. If `postprocessing=true`, some colors might be replaced with `0` (the "neutral" color) as long as they are not needed during decompression. +The optional `forced_colors` keyword argument is used to enforce predefined vertex colors (e.g. coming from another optimization algorithm) but still run the star coloring procedure to verify correctness and build auxiliary data structures, useful during decompression. + # See also - [`AdjacencyGraph`](@ref) @@ -77,7 +107,10 @@ If `postprocessing=true`, some colors might be replaced with `0` (the "neutral" > [_New Acyclic and Star Coloring Algorithms with Application to Computing Hessians_](https://epubs.siam.org/doi/abs/10.1137/050639879), Gebremedhin et al. (2007), Algorithm 4.1 """ function star_coloring( - g::AdjacencyGraph{T}, vertices_in_order::AbstractVector{<:Integer}, postprocessing::Bool + g::AdjacencyGraph{T}, + vertices_in_order::AbstractVector{<:Integer}, + postprocessing::Bool; + forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing, ) where {T<:Integer} # Initialize data structures nv = nb_vertices(g) @@ -115,10 +148,18 @@ function star_coloring( end end end - for i in eachindex(forbidden_colors) - if forbidden_colors[i] != v - color[v] = i - break + if isnothing(forced_colors) + for i in eachindex(forbidden_colors) + if forbidden_colors[i] != v + color[v] = i + break + end + end + else + if forbidden_colors[forced_colors[v]] == v # TODO: handle forced_colors[v] == 0 + throw(InvalidColoringError()) + else + color[v] = forced_colors[v] end end _update_stars!(star, hub, g, v, color, first_neighbor) @@ -271,6 +312,7 @@ function acyclic_coloring( end end end + # TODO: handle forced colors for i in eachindex(forbidden_colors) if forbidden_colors[i] != v color[v] = i diff --git a/src/constant.jl b/src/constant.jl index d79caab0..56bf4101 100644 --- a/src/constant.jl +++ b/src/constant.jl @@ -4,20 +4,24 @@ Coloring algorithm which always returns the same precomputed vector of colors. Useful when the optimal coloring of a matrix can be determined a priori due to its specific structure (e.g. banded). -It is passed as an argument to the main function [`coloring`](@ref), but will only work if the associated `problem` has `:nonsymmetric` structure. -Indeed, for symmetric coloring problems, we need more than just the vector of colors to allow fast decompression. +It is passed as an argument to the main function [`coloring`](@ref), but will only work if the associated `problem` has a `:column` or `:row` partition. # Constructors ConstantColoringAlgorithm{partition}(matrix_template, color) - ConstantColoringAlgorithm(matrix_template, color; partition=:column) + ConstantColoringAlgorithm{partition,structure}(matrix_template, color) + ConstantColoringAlgorithm( + matrix_template, color; + structure=:nonsymmetric, partition=:column + ) - `partition::Symbol`: either `:row` or `:column`. +- `structure::Symbol`: either `:nonsymmetric` or `:symmetric`. - `matrix_template::AbstractMatrix`: matrix for which the vector of colors was precomputed (the algorithm will only accept matrices of the exact same size). - `color::Vector{<:Integer}`: vector of integer colors, one for each row or column (depending on `partition`). !!! warning - The second constructor (based on keyword arguments) is type-unstable. + The constructor based on keyword arguments is type-unstable if these arguments are not compile-time constants. We do not necessarily verify consistency between the matrix template and the vector of colors, this is the responsibility of the user. @@ -63,71 +67,68 @@ julia> column_colors(result) - [`ADTypes.column_coloring`](@extref ADTypes.column_coloring) - [`ADTypes.row_coloring`](@extref ADTypes.row_coloring) +- [`ADTypes.symmetric_coloring`](@extref ADTypes.symmetric_coloring) """ -struct ConstantColoringAlgorithm{ - partition, - M<:AbstractMatrix, - T<:Integer, - R<:AbstractColoringResult{:nonsymmetric,partition,:direct}, -} <: ADTypes.AbstractColoringAlgorithm +struct ConstantColoringAlgorithm{partition,structure,M<:AbstractMatrix,T<:Integer} <: + ADTypes.AbstractColoringAlgorithm matrix_template::M color::Vector{T} - result::R -end -function ConstantColoringAlgorithm{:column}( - matrix_template::AbstractMatrix, color::Vector{<:Integer} -) - bg = BipartiteGraph(matrix_template) - result = ColumnColoringResult(matrix_template, bg, color) - T, M, R = eltype(bg), typeof(matrix_template), typeof(result) - return ConstantColoringAlgorithm{:column,M,T,R}(matrix_template, color, result) + function ConstantColoringAlgorithm{partition,structure}( + matrix_template::AbstractMatrix, color::Vector{<:Integer} + ) where {partition,structure} + check_valid_problem(structure, partition) + return new{partition,structure,typeof(matrix_template),eltype(color)}( + matrix_template, color + ) + end end -function ConstantColoringAlgorithm{:row}( +function ConstantColoringAlgorithm{partition}( matrix_template::AbstractMatrix, color::Vector{<:Integer} -) - bg = BipartiteGraph(matrix_template) - result = RowColoringResult(matrix_template, bg, color) - T, M, R = eltype(bg), typeof(matrix_template), typeof(result) - return ConstantColoringAlgorithm{:row,M,T,R}(matrix_template, color, result) +) where {partition} + return ConstantColoringAlgorithm{partition,:nonsymmetric}(matrix_template, color) end function ConstantColoringAlgorithm( - matrix_template::AbstractMatrix, color::Vector{<:Integer}; partition::Symbol=:column + matrix_template::AbstractMatrix, + color::Vector{<:Integer}; + structure::Symbol=:nonsymmetric, + partition::Symbol=:column, ) - return ConstantColoringAlgorithm{partition}(matrix_template, color) + return ConstantColoringAlgorithm{partition,structure}(matrix_template, color) end -function coloring( - A::AbstractMatrix, - ::ColoringProblem{:nonsymmetric,partition}, - algo::ConstantColoringAlgorithm{partition}; - decompression_eltype::Type=Float64, - symmetric_pattern::Bool=false, -) where {partition} - (; matrix_template, result) = algo +function check_template(algo::ConstantColoringAlgorithm, A::AbstractMatrix) + (; matrix_template) = algo if size(A) != size(matrix_template) throw( DimensionMismatch( "`ConstantColoringAlgorithm` expected matrix of size $(size(matrix_template)) but got matrix of size $(size(A))", ), ) - else - return result end end function ADTypes.column_coloring( - A::AbstractMatrix, algo::ConstantColoringAlgorithm{:column} + A::AbstractMatrix, algo::ConstantColoringAlgorithm{:column,:nonsymmetric} +) + check_template(algo, A) + return algo.color +end + +function ADTypes.row_coloring( + A::AbstractMatrix, algo::ConstantColoringAlgorithm{:row,:nonsymmetric} ) - problem = ColoringProblem{:nonsymmetric,:column}() - result = coloring(A, problem, algo) - return column_colors(result) + check_template(algo, A) + return algo.color end -function ADTypes.row_coloring(A::AbstractMatrix, algo::ConstantColoringAlgorithm) - problem = ColoringProblem{:nonsymmetric,:row}() - result = coloring(A, problem, algo) - return row_colors(result) +function ADTypes.symmetric_coloring( + A::AbstractMatrix, algo::ConstantColoringAlgorithm{:column,:symmetric} +) + check_template(algo, A) + return algo.color end + +# TODO: handle bidirectional once https://github.com/SciML/ADTypes.jl/issues/69 is done diff --git a/src/interface.jl b/src/interface.jl index c1a56c29..b8d4ce54 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -230,13 +230,14 @@ function _coloring( ::ColoringProblem{:nonsymmetric,:column}, algo::GreedyColoringAlgorithm, decompression_eltype::Type, - symmetric_pattern::Bool, + symmetric_pattern::Bool; + forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing, ) symmetric_pattern = symmetric_pattern || A isa Union{Symmetric,Hermitian} bg = BipartiteGraph(A; symmetric_pattern) color_by_order = map(algo.orders) do order vertices_in_order = vertices(bg, Val(2), order) - return partial_distance2_coloring(bg, Val(2), vertices_in_order) + return partial_distance2_coloring(bg, Val(2), vertices_in_order; forced_colors) end color = argmin(maximum, color_by_order) if speed_setting isa WithResult @@ -252,13 +253,14 @@ function _coloring( ::ColoringProblem{:nonsymmetric,:row}, algo::GreedyColoringAlgorithm, decompression_eltype::Type, - symmetric_pattern::Bool, + symmetric_pattern::Bool; + forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing, ) symmetric_pattern = symmetric_pattern || A isa Union{Symmetric,Hermitian} bg = BipartiteGraph(A; symmetric_pattern) color_by_order = map(algo.orders) do order vertices_in_order = vertices(bg, Val(1), order) - return partial_distance2_coloring(bg, Val(1), vertices_in_order) + return partial_distance2_coloring(bg, Val(1), vertices_in_order; forced_colors) end color = argmin(maximum, color_by_order) if speed_setting isa WithResult @@ -274,12 +276,13 @@ function _coloring( ::ColoringProblem{:symmetric,:column}, algo::GreedyColoringAlgorithm{:direct}, decompression_eltype::Type, - symmetric_pattern::Bool, + symmetric_pattern::Bool; + forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing, ) ag = AdjacencyGraph(A; has_diagonal=true) color_and_star_set_by_order = map(algo.orders) do order vertices_in_order = vertices(ag, order) - return star_coloring(ag, vertices_in_order, algo.postprocessing) + return star_coloring(ag, vertices_in_order, algo.postprocessing; forced_colors) end color, star_set = argmin(maximum ∘ first, color_and_star_set_by_order) if speed_setting isa WithResult @@ -316,13 +319,16 @@ function _coloring( ::ColoringProblem{:nonsymmetric,:bidirectional}, algo::GreedyColoringAlgorithm{:direct}, decompression_eltype::Type{R}, - symmetric_pattern::Bool, + symmetric_pattern::Bool; + 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; has_diagonal=false) outputs_by_order = map(algo.orders) do order vertices_in_order = vertices(ag, order) - _color, _star_set = star_coloring(ag, vertices_in_order, algo.postprocessing) + _color, _star_set = star_coloring( + ag, vertices_in_order, algo.postprocessing; forced_colors + ) (_row_color, _column_color, _symmetric_to_row, _symmetric_to_column) = remap_colors( eltype(ag), _color, maximum(_color), size(A)... ) diff --git a/test/adtypes.jl b/test/adtypes.jl index bbb124a6..da0807d3 100644 --- a/test/adtypes.jl +++ b/test/adtypes.jl @@ -1,26 +1,49 @@ using ADTypes: ADTypes using SparseArrays +using LinearAlgebra using SparseMatrixColorings using Test -@testset "Column coloring" begin - problem = ColoringProblem(; structure=:nonsymmetric, partition=:column) - algo = ADTypes.NoColoringAlgorithm() - A = sprand(10, 20, 0.1) - result = coloring(A, problem, algo) - B = compress(A, result) - @test size(B) == size(A) - @test column_colors(result) == ADTypes.column_coloring(A, algo) - @test decompress(B, result) == A -end +@testset "NoColoringAlgorithm" begin + @testset "Column coloring" begin + problem = ColoringProblem(; structure=:nonsymmetric, partition=:column) + algo = ADTypes.NoColoringAlgorithm() + A = sprand(10, 20, 0.3) + result = coloring(A, problem, algo) + B = compress(A, result) + @test size(B) == size(A) + @test column_colors(result) == ADTypes.column_coloring(A, algo) + @test decompress(B, result) == A + end + + @testset "Row coloring" begin + problem = ColoringProblem(; structure=:nonsymmetric, partition=:row) + algo = ADTypes.NoColoringAlgorithm() + A = sprand(10, 20, 0.3) + result = coloring(A, problem, algo) + B = compress(A, result) + @test size(B) == size(A) + @test row_colors(result) == ADTypes.row_coloring(A, algo) + @test decompress(B, result) == A + end + + @testset "Symmetric coloring" begin + problem = ColoringProblem(; structure=:symmetric, partition=:column) + algo = ADTypes.NoColoringAlgorithm() + A = Symmetric(sprand(20, 20, 0.3)) + result = coloring(A, problem, algo) + B = compress(A, result) + @test size(B) == size(A) + @test column_colors(result) == ADTypes.column_coloring(A, algo) + @test decompress(B, result) == A + end -@testset "Row coloring" begin - problem = ColoringProblem(; structure=:nonsymmetric, partition=:row) - algo = ADTypes.NoColoringAlgorithm() - A = sprand(10, 20, 0.1) - result = coloring(A, problem, algo) - B = compress(A, result) - @test size(B) == size(A) - @test row_colors(result) == ADTypes.row_coloring(A, algo) - @test decompress(B, result) == A + @testset "Bicoloring" begin + problem = ColoringProblem(; structure=:nonsymmetric, partition=:bidirectional) + algo = ADTypes.NoColoringAlgorithm() + A = sprand(10, 20, 0.3) + result = coloring(A, problem, algo) + Br, Bc = compress(A, result) + @test decompress(Br, Bc, result) == A + end end diff --git a/test/constant.jl b/test/constant.jl index 5de881b8..97c7f238 100644 --- a/test/constant.jl +++ b/test/constant.jl @@ -1,16 +1,20 @@ using ADTypes: ADTypes using SparseMatrixColorings +using SparseMatrixColorings: InvalidColoringError using Test -matrix_template = ones(100, 200) +matrix_template = ones(Bool, 10, 20) +sym_matrix_template = ones(Bool, 10, 10) @testset "Column coloring" begin problem = ColoringProblem(; structure=:nonsymmetric, partition=:column) - color = rand(1:5, size(matrix_template, 2)) + color = collect(1:20) algo = ConstantColoringAlgorithm(matrix_template, color; partition=:column) - wrong_algo = ConstantColoringAlgorithm(matrix_template, color; partition=:row) + wrong_algo = ConstantColoringAlgorithm{:row}(matrix_template, color) + wrong_color = ConstantColoringAlgorithm{:column}(matrix_template, ones(Int, 20)) @test_throws DimensionMismatch coloring(transpose(matrix_template), problem, algo) @test_throws MethodError coloring(matrix_template, problem, wrong_algo) + @test_throws InvalidColoringError coloring(matrix_template, problem, wrong_color) result = coloring(matrix_template, problem, algo) @test column_colors(result) == color @test ADTypes.column_coloring(matrix_template, algo) == color @@ -19,9 +23,13 @@ end @testset "Row coloring" begin problem = ColoringProblem(; structure=:nonsymmetric, partition=:row) - color = rand(1:5, size(matrix_template, 1)) + color = collect(1:10) algo = ConstantColoringAlgorithm(matrix_template, color; partition=:row) + wrong_algo = ConstantColoringAlgorithm{:column}(matrix_template, color) + wrong_color = ConstantColoringAlgorithm{:row}(matrix_template, ones(Int, 10)) @test_throws DimensionMismatch coloring(transpose(matrix_template), problem, algo) + @test_throws MethodError coloring(matrix_template, problem, wrong_algo) + @test_throws InvalidColoringError coloring(matrix_template, problem, wrong_color) result = coloring(matrix_template, problem, algo) @test row_colors(result) == color @test ADTypes.row_coloring(matrix_template, algo) == color @@ -29,9 +37,22 @@ end end @testset "Symmetric coloring" begin - wrong_problem = ColoringProblem(; structure=:symmetric, partition=:column) - color = rand(1:5, size(matrix_template, 2)) - algo = ConstantColoringAlgorithm(matrix_template, color; partition=:column) - @test_throws MethodError coloring(matrix_template, wrong_problem, algo) - @test_throws MethodError ADTypes.symmetric_coloring(matrix_template, algo) + problem = ColoringProblem(; structure=:symmetric, partition=:column) + color = collect(1:10) + algo = ConstantColoringAlgorithm( + sym_matrix_template, color; partition=:column, structure=:symmetric + ) + wrong_algo = ConstantColoringAlgorithm{:column,:nonsymmetric}( + sym_matrix_template, color + ) + wrong_color = ConstantColoringAlgorithm{:column,:symmetric}( + sym_matrix_template, ones(Int, 20) + ) + @test_throws DimensionMismatch coloring(matrix_template, problem, algo) + @test_throws MethodError coloring(sym_matrix_template, problem, wrong_algo) + @test_throws InvalidColoringError coloring(sym_matrix_template, problem, wrong_color) + result = coloring(sym_matrix_template, problem, algo) + @test column_colors(result) == color + @test ADTypes.symmetric_coloring(sym_matrix_template, algo) == color + @test_throws MethodError ADTypes.column_coloring(sym_matrix_template, algo) end diff --git a/test/result.jl b/test/result.jl index 0bea3353..49344aab 100644 --- a/test/result.jl +++ b/test/result.jl @@ -1,3 +1,4 @@ +using SparseMatrixColorings using SparseMatrixColorings: group_by_color, UnsupportedDecompressionError using Test @@ -20,7 +21,7 @@ using Test end @testset "Empty compression" begin - A = rand(10, 10) + A = zeros(Bool, 10, 10) color = zeros(Int, 10) problem = ColoringProblem{:nonsymmetric,:column}() algo = ConstantColoringAlgorithm(A, color; partition=:column) diff --git a/test/runtests.jl b/test/runtests.jl index 81c4ecc9..cfd51dce 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -24,11 +24,12 @@ include("utils.jl") @testset "JET" begin JET.test_package(SparseMatrixColorings; target_defined_modules=true) end - @testset "JuliaFormatter" begin - @test JuliaFormatter.format( - SparseMatrixColorings; verbose=false, overwrite=false - ) - end + # @testset "JuliaFormatter" begin + # TODO: switch to Runic (temporarily deactivated) + # @test JuliaFormatter.format( + # SparseMatrixColorings; verbose=false, overwrite=false + # ) + # end @testset "Doctests" begin Documenter.doctest(SparseMatrixColorings) end From 17bc04461cf4c730b0fc91797db1effae25d07b4 Mon Sep 17 00:00:00 2001 From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com> Date: Fri, 31 Oct 2025 22:44:22 +0100 Subject: [PATCH 10/39] Optimal coloring algorithm with JuMP formulation (#271) --- Project.toml | 5 ++ docs/src/api.md | 6 +++ ext/SparseMatrixColoringsJuMPExt.jl | 80 +++++++++++++++++++++++++++++ src/SparseMatrixColorings.jl | 2 + src/optimal.jl | 34 ++++++++++++ test/Project.toml | 3 ++ test/optimal.jl | 46 +++++++++++++++++ test/runtests.jl | 3 ++ 8 files changed, 179 insertions(+) create mode 100644 ext/SparseMatrixColoringsJuMPExt.jl create mode 100644 src/optimal.jl create mode 100644 test/optimal.jl diff --git a/Project.toml b/Project.toml index 3eb2711d..864fae83 100644 --- a/Project.toml +++ b/Project.toml @@ -15,11 +15,14 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8" Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" +MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" [extensions] SparseMatrixColoringsCUDAExt = "CUDA" SparseMatrixColoringsCliqueTreesExt = "CliqueTrees" SparseMatrixColoringsColorsExt = "Colors" +SparseMatrixColoringsJuMPExt = ["JuMP", "MathOptInterface"] [compat] ADTypes = "1.2.1" @@ -27,7 +30,9 @@ CUDA = "5.8.2" CliqueTrees = "1" Colors = "0.12.11, 0.13" DocStringExtensions = "0.8,0.9" +JuMP = "1.29.1" LinearAlgebra = "<0.0.1, 1" +MathOptInterface = "1.45.0" PrecompileTools = "1.2.1" Random = "<0.0.1, 1" SparseArrays = "<0.0.1, 1" diff --git a/docs/src/api.md b/docs/src/api.md index c121d1e3..f90c888f 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -17,8 +17,14 @@ SparseMatrixColorings coloring fast_coloring ColoringProblem +``` + +## Coloring algorithms + +```@docs GreedyColoringAlgorithm ConstantColoringAlgorithm +OptimalColoringAlgorithm ``` ## Result analysis diff --git a/ext/SparseMatrixColoringsJuMPExt.jl b/ext/SparseMatrixColoringsJuMPExt.jl new file mode 100644 index 00000000..ca4f5e7b --- /dev/null +++ b/ext/SparseMatrixColoringsJuMPExt.jl @@ -0,0 +1,80 @@ +module SparseMatrixColoringsJuMPExt + +using ADTypes: ADTypes +using JuMP: + Model, + is_solved_and_feasible, + optimize!, + primal_status, + set_silent, + set_start_value, + value, + @variable, + @constraint, + @objective +using JuMP +import MathOptInterface as MOI +using SparseMatrixColorings: + BipartiteGraph, OptimalColoringAlgorithm, nb_vertices, neighbors, pattern, vertices + +function optimal_distance2_coloring( + bg::BipartiteGraph, + ::Val{side}, + optimizer::O; + silent::Bool=true, + assert_solved::Bool=true, +) where {side,O} + other_side = 3 - side + n = nb_vertices(bg, Val(side)) + model = Model(optimizer) + silent && set_silent(model) + # one variable per vertex to color, removing some renumbering symmetries + @variable(model, 1 <= color[i=1:n] <= i, Int) + # one variable to count the number of distinct colors + @variable(model, ncolors, Int) + @constraint(model, [ncolors; color] in MOI.CountDistinct(n + 1)) + # distance-2 coloring: neighbors of the same vertex must have distinct colors + for i in vertices(bg, Val(other_side)) + neigh = neighbors(bg, Val(other_side), i) + @constraint(model, color[neigh] in MOI.AllDifferent(length(neigh))) + end + # minimize the number of distinct colors (can't use maximum because they are not necessarily numbered contiguously) + @objective(model, Min, ncolors) + # actual solving step where time is spent + optimize!(model) + if assert_solved + # assert feasibility and optimality + @assert is_solved_and_feasible(model) + else + # only assert feasibility + @assert primal_status(model) == MOI.FEASIBLE_POINT + end + # native solver solutions are floating point numbers + color_int = round.(Int, value.(color)) + # remap to 1:cmax in case they are not contiguous + true_ncolors = 0 + remap = fill(0, maximum(color_int)) + for c in color_int + if remap[c] == 0 + true_ncolors += 1 + remap[c] = true_ncolors + end + end + return remap[color_int] +end + +function ADTypes.column_coloring(A::AbstractMatrix, algo::OptimalColoringAlgorithm) + bg = BipartiteGraph(A) + return optimal_distance2_coloring( + bg, Val(2), algo.optimizer; algo.silent, algo.assert_solved + ) +end + +function ADTypes.row_coloring(A::AbstractMatrix, algo::OptimalColoringAlgorithm) + bg = BipartiteGraph(A) + return optimal_distance2_coloring( + bg, Val(1), algo.optimizer; algo.silent, algo.assert_solved + ) +end + +end diff --git a/src/SparseMatrixColorings.jl b/src/SparseMatrixColorings.jl index 7c25bac4..fca70229 100644 --- a/src/SparseMatrixColorings.jl +++ b/src/SparseMatrixColorings.jl @@ -56,6 +56,7 @@ include("decompression.jl") include("check.jl") include("examples.jl") include("show_colors.jl") +include("optimal.jl") include("precompile.jl") @@ -64,6 +65,7 @@ export DynamicDegreeBasedOrder, SmallestLast, IncidenceDegree, DynamicLargestFir export PerfectEliminationOrder export ColoringProblem, GreedyColoringAlgorithm, AbstractColoringResult export ConstantColoringAlgorithm +export OptimalColoringAlgorithm export coloring, fast_coloring export column_colors, row_colors, ncolors export column_groups, row_groups diff --git a/src/optimal.jl b/src/optimal.jl new file mode 100644 index 00000000..2a70a59d --- /dev/null +++ b/src/optimal.jl @@ -0,0 +1,34 @@ +""" + OptimalColoringAlgorithm + +Coloring algorithm that relies on mathematical programming with [JuMP](https://jump.dev/) to find an optimal coloring. + +!!! warning + This algorithm is only available when JuMP is loaded. If you encounter a method error, run `import JuMP` in your REPL and try again. + It only works for nonsymmetric, unidirectional colorings problems. + +!!! danger + The coloring problem is NP-hard, so it is unreasonable to expect an optimal solution in reasonable time for large instances. + +# Constructor + + OptimalColoringAlgorithm(optimizer; silent::Bool=true, assert_solved::Bool=true) + +The `optimizer` argument can be any JuMP-compatible optimizer. +However, the problem formulation is best suited to CP-SAT optimizers like [MiniZinc](https://github.com/jump-dev/MiniZinc.jl). +You can use [`optimizer_with_attributes`](https://jump.dev/JuMP.jl/stable/api/JuMP/#optimizer_with_attributes) to set solver-specific parameters. + +# Keyword arguments + +- `silent`: whether to suppress solver output +- `assert_solved`: whether to check that the solver found an optimal solution (as opposed to running out of time for example) +""" +struct OptimalColoringAlgorithm{O} <: ADTypes.AbstractColoringAlgorithm + optimizer::O + silent::Bool + assert_solved::Bool +end + +function OptimalColoringAlgorithm(optimizer; silent::Bool=true, assert_solved::Bool=true) + return OptimalColoringAlgorithm(optimizer, silent, assert_solved) +end diff --git a/test/Project.toml b/test/Project.toml index ab595aba..217d62c9 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -12,10 +12,13 @@ CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8" Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MatrixDepot = "b51810bb-c9f3-55da-ae3c-350fc1fbce05" +MiniZinc = "a7f392d2-6c35-496e-b8cc-0974fbfcbf91" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35" diff --git a/test/optimal.jl b/test/optimal.jl new file mode 100644 index 00000000..10fc32ef --- /dev/null +++ b/test/optimal.jl @@ -0,0 +1,46 @@ +using SparseArrays +using SparseMatrixColorings +using StableRNGs +using Test +using JuMP +using MiniZinc +using HiGHS + +rng = StableRNG(0) + +asymmetric_params = vcat( + [(10, 20, p) for p in (0.0:0.1:0.5)], [(20, 10, p) for p in (0.0:0.1:0.5)] +) + +algo = GreedyColoringAlgorithm() +optalgo = OptimalColoringAlgorithm(() -> MiniZinc.Optimizer{Float64}("highs"); silent=false) + +@testset "Column coloring" begin + problem = ColoringProblem(; structure=:nonsymmetric, partition=:column) + for (m, n, p) in asymmetric_params + A = sprand(rng, m, n, p) + result = coloring(A, problem, algo) + optresult = coloring(A, problem, optalgo) + @test ncolors(result) >= ncolors(optresult) + end +end + +@testset "Row coloring" begin + problem = ColoringProblem(; structure=:nonsymmetric, partition=:row) + for (m, n, p) in asymmetric_params + A = sprand(rng, m, n, p) + result = coloring(A, problem, algo) + optresult = coloring(A, problem, optalgo) + @test ncolors(result) >= ncolors(optresult) + end +end + +@testset "Too big" begin + A = sprand(rng, Bool, 100, 100, 0.1) + optalgo_timelimit = OptimalColoringAlgorithm( + optimizer_with_attributes(HiGHS.Optimizer, "time_limit" => 10.0); # 1 second + silent=false, + assert_solved=false, + ) + @test_throws AssertionError coloring(A, ColoringProblem(), optalgo_timelimit) +end diff --git a/test/runtests.jl b/test/runtests.jl index cfd51dce..fa5cc7ff 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -59,6 +59,9 @@ include("utils.jl") @testset "Constant coloring" begin include("constant.jl") end + @testset "Optimal coloring" begin + include("optimal.jl") + end @testset "ADTypes coloring algorithms" begin include("adtypes.jl") end From 33a1ef1a6b7bab5b5dcdf2a5929e6e55bebfd887 Mon Sep 17 00:00:00 2001 From: Alexis Montoison <35051714+amontoison@users.noreply.github.com> Date: Sat, 1 Nov 2025 01:06:55 -0500 Subject: [PATCH 11/39] Add a file postprocessing.jl (#275) --- src/SparseMatrixColorings.jl | 1 + src/coloring.jl | 147 ----------------------------------- src/postprocessing.jl | 146 ++++++++++++++++++++++++++++++++++ 3 files changed, 147 insertions(+), 147 deletions(-) create mode 100644 src/postprocessing.jl diff --git a/src/SparseMatrixColorings.jl b/src/SparseMatrixColorings.jl index fca70229..481c8584 100644 --- a/src/SparseMatrixColorings.jl +++ b/src/SparseMatrixColorings.jl @@ -47,6 +47,7 @@ include("graph.jl") include("forest.jl") include("order.jl") include("coloring.jl") +include("postprocessing.jl") include("result.jl") include("matrices.jl") include("interface.jl") diff --git a/src/coloring.jl b/src/coloring.jl index 7288b773..2746b30b 100644 --- a/src/coloring.jl +++ b/src/coloring.jl @@ -643,150 +643,3 @@ function TreeSet( return TreeSet(reverse_bfs_orders, is_star, tree_edge_indices, nt) end - -## Postprocessing, mirrors decompression code - -function postprocess!( - color::AbstractVector{<:Integer}, - star_or_tree_set::Union{StarSet,TreeSet}, - g::AdjacencyGraph, - offsets::AbstractVector{<:Integer}, -) - S = pattern(g) - edge_to_index = edge_indices(g) - # flag which colors are actually used during decompression - nb_colors = maximum(color) - color_used = zeros(Bool, nb_colors) - - # nonzero diagonal coefficients force the use of their respective color (there can be no neutral colors if the diagonal is fully nonzero) - if has_diagonal(g) - for i in axes(S, 1) - if !iszero(S[i, i]) - color_used[color[i]] = true - end - end - end - - if star_or_tree_set isa StarSet - # only the colors of the hubs are used - (; star, hub) = star_or_tree_set - nb_trivial_stars = 0 - - # Iterate through all non-trivial stars - for s in eachindex(hub) - h = hub[s] - if h > 0 - color_used[color[h]] = true - else - nb_trivial_stars += 1 - end - end - - # Process the trivial stars (if any) - if nb_trivial_stars > 0 - rvS = rowvals(S) - for j in axes(S, 2) - for k in nzrange(S, j) - i = rvS[k] - if i > j - index_ij = edge_to_index[k] - s = star[index_ij] - h = hub[s] - if h < 0 - h = abs(h) - spoke = h == j ? i : j - if color_used[color[spoke]] - # Switch the hub and the spoke to possibly avoid adding one more used color - hub[s] = spoke - else - # Keep the current hub - color_used[color[h]] = true - end - end - end - end - end - end - else - # only the colors of non-leaf vertices are used - (; reverse_bfs_orders, is_star, tree_edge_indices, nt) = star_or_tree_set - nb_trivial_trees = 0 - - # Iterate through all non-trivial trees - for k in 1:nt - # Position of the first edge in the tree - first = tree_edge_indices[k] - - # Total number of edges in the tree - ne_tree = tree_edge_indices[k + 1] - first - - # Check if we have more than one edge in the tree (non-trivial tree) - if ne_tree > 1 - # Determine if the tree is a star - if is_star[k] - # It is a non-trivial star and only the color of the hub is needed - (_, hub) = reverse_bfs_orders[first] - color_used[color[hub]] = true - else - # It is not a star and both colors are needed during the decompression - (i, j) = reverse_bfs_orders[first] - color_used[color[i]] = true - color_used[color[j]] = true - end - else - nb_trivial_trees += 1 - end - end - - # Process the trivial trees (if any) - if nb_trivial_trees > 0 - for k in 1:nt - # Position of the first edge in the tree - first = tree_edge_indices[k] - - # Total number of edges in the tree - ne_tree = tree_edge_indices[k + 1] - first - - # Check if we have exactly one edge in the tree - if ne_tree == 1 - (i, j) = reverse_bfs_orders[first] - if color_used[color[i]] - # Make i the root to avoid possibly adding one more used color - # Switch it with the (only) leaf - reverse_bfs_orders[first] = (j, i) - else - # Keep j as the root - color_used[color[j]] = true - end - end - end - end - end - - # if at least one of the colors is useless, modify the color assignments of vertices - if any(!, color_used) - num_colors_useless = 0 - - # determine what are the useless colors and compute the offsets - for ci in 1:nb_colors - if color_used[ci] - offsets[ci] = num_colors_useless - else - num_colors_useless += 1 - end - end - - # assign the neutral color to every vertex with a useless color and remap the colors - for i in eachindex(color) - ci = color[i] - if !color_used[ci] - # assign the neutral color - color[i] = 0 - else - # remap the color to not have any gap - color[i] -= offsets[ci] - end - end - end - return color -end diff --git a/src/postprocessing.jl b/src/postprocessing.jl new file mode 100644 index 00000000..217b663b --- /dev/null +++ b/src/postprocessing.jl @@ -0,0 +1,146 @@ +## Postprocessing + +function postprocess!( + color::AbstractVector{<:Integer}, + star_or_tree_set::Union{StarSet,TreeSet}, + g::AdjacencyGraph, + offsets::AbstractVector{<:Integer}, +) + S = pattern(g) + edge_to_index = edge_indices(g) + # flag which colors are actually used during decompression + nb_colors = maximum(color) + color_used = zeros(Bool, nb_colors) + + # nonzero diagonal coefficients force the use of their respective color (there can be no neutral colors if the diagonal is fully nonzero) + if has_diagonal(g) + for i in axes(S, 1) + if !iszero(S[i, i]) + color_used[color[i]] = true + end + end + end + + if star_or_tree_set isa StarSet + # only the colors of the hubs are used + (; star, hub) = star_or_tree_set + nb_trivial_stars = 0 + + # Iterate through all non-trivial stars + for s in eachindex(hub) + h = hub[s] + if h > 0 + color_used[color[h]] = true + else + nb_trivial_stars += 1 + end + end + + # Process the trivial stars (if any) + if nb_trivial_stars > 0 + rvS = rowvals(S) + for j in axes(S, 2) + for k in nzrange(S, j) + i = rvS[k] + if i > j + index_ij = edge_to_index[k] + s = star[index_ij] + h = hub[s] + if h < 0 + h = abs(h) + spoke = h == j ? i : j + if color_used[color[spoke]] + # Switch the hub and the spoke to possibly avoid adding one more used color + hub[s] = spoke + else + # Keep the current hub + color_used[color[h]] = true + end + end + end + end + end + end + else + # only the colors of non-leaf vertices are used + (; reverse_bfs_orders, is_star, tree_edge_indices, nt) = star_or_tree_set + nb_trivial_trees = 0 + + # Iterate through all non-trivial trees + for k in 1:nt + # Position of the first edge in the tree + first = tree_edge_indices[k] + + # Total number of edges in the tree + ne_tree = tree_edge_indices[k + 1] - first + + # Check if we have more than one edge in the tree (non-trivial tree) + if ne_tree > 1 + # Determine if the tree is a star + if is_star[k] + # It is a non-trivial star and only the color of the hub is needed + (_, hub) = reverse_bfs_orders[first] + color_used[color[hub]] = true + else + # It is not a star and both colors are needed during the decompression + (i, j) = reverse_bfs_orders[first] + color_used[color[i]] = true + color_used[color[j]] = true + end + else + nb_trivial_trees += 1 + end + end + + # Process the trivial trees (if any) + if nb_trivial_trees > 0 + for k in 1:nt + # Position of the first edge in the tree + first = tree_edge_indices[k] + + # Total number of edges in the tree + ne_tree = tree_edge_indices[k + 1] - first + + # Check if we have exactly one edge in the tree + if ne_tree == 1 + (i, j) = reverse_bfs_orders[first] + if color_used[color[i]] + # Make i the root to avoid possibly adding one more used color + # Switch it with the (only) leaf + reverse_bfs_orders[first] = (j, i) + else + # Keep j as the root + color_used[color[j]] = true + end + end + end + end + end + + # if at least one of the colors is useless, modify the color assignments of vertices + if any(!, color_used) + num_colors_useless = 0 + + # determine what are the useless colors and compute the offsets + for ci in 1:nb_colors + if color_used[ci] + offsets[ci] = num_colors_useless + else + num_colors_useless += 1 + end + end + + # assign the neutral color to every vertex with a useless color and remap the colors + for i in eachindex(color) + ci = color[i] + if !color_used[ci] + # assign the neutral color + color[i] = 0 + else + # remap the color to not have any gap + color[i] -= offsets[ci] + end + end + end + return color +end From b7c02072da957f72f9197ea1562e7fc41c278cdd Mon Sep 17 00:00:00 2001 From: Alexis Montoison <35051714+amontoison@users.noreply.github.com> Date: Sat, 1 Nov 2025 02:48:34 -0500 Subject: [PATCH 12/39] Rename has_diagonal into augmented_graph (#273) * Rename has_diagonal into augmented_graph * Fix a typo in postprocessing.jl * Update src/graph.jl Co-authored-by: Guillaume Dalle <22795598+gdalle@users.noreply.github.com> --------- Co-authored-by: Guillaume Dalle <22795598+gdalle@users.noreply.github.com> --- src/coloring.jl | 22 +++++++++++----------- src/decompression.jl | 4 ++-- src/graph.jl | 32 ++++++++++++++++---------------- src/interface.jl | 8 ++++---- src/order.jl | 2 +- src/postprocessing.jl | 2 +- src/result.jl | 2 +- test/graph.jl | 2 +- test/order.jl | 18 ++++++++---------- 9 files changed, 45 insertions(+), 47 deletions(-) diff --git a/src/coloring.jl b/src/coloring.jl index 2746b30b..cf3d8f79 100644 --- a/src/coloring.jl +++ b/src/coloring.jl @@ -124,7 +124,7 @@ function star_coloring( for v in vertices_in_order for (w, index_vw) in neighbors_with_edge_indices(g, v) - !has_diagonal(g) || (v == w && continue) + augmented_graph(g) || (v == w && continue) color_w = color[w] iszero(color_w) && continue forbidden_colors[color_w] = v @@ -139,7 +139,7 @@ function star_coloring( else first_neighbor[color[w]] = (v, w, index_vw) for (x, index_wx) in neighbors_with_edge_indices(g, w) - !has_diagonal(g) || (w == x && continue) + augmented_graph(g) || (w == x && continue) color_x = color[x] (x == v || iszero(color_x)) && continue if x == hub[star[index_wx]] # potential Case 2 (which is always false for trivial stars with two vertices, since the associated hub is negative) @@ -184,7 +184,7 @@ function _treat!( color::AbstractVector{<:Integer}, ) for x in neighbors(g, w) - !has_diagonal(g) || (w == x && continue) + augmented_graph(g) || (w == x && continue) color_x = color[x] iszero(color_x) && continue forbidden_colors[color_x] = v @@ -204,12 +204,12 @@ function _update_stars!( first_neighbor::AbstractVector{<:Tuple}, ) for (w, index_vw) in neighbors_with_edge_indices(g, v) - !has_diagonal(g) || (v == w && continue) + augmented_graph(g) || (v == w && continue) color_w = color[w] iszero(color_w) && continue x_exists = false for (x, index_wx) in neighbors_with_edge_indices(g, w) - !has_diagonal(g) || (w == x && continue) + augmented_graph(g) || (w == x && continue) if x != v && color[x] == color[v] # vw, wx ∈ E star_wx = star[index_wx] hub[star_wx] = w # this may already be true @@ -286,16 +286,16 @@ function acyclic_coloring( for v in vertices_in_order for w in neighbors(g, v) - !has_diagonal(g) || (v == w && continue) + augmented_graph(g) || (v == w && continue) color_w = color[w] iszero(color_w) && continue forbidden_colors[color_w] = v end for w in neighbors(g, v) - !has_diagonal(g) || (v == w && continue) + augmented_graph(g) || (v == w && continue) iszero(color[w]) && continue for (x, index_wx) in neighbors_with_edge_indices(g, w) - !has_diagonal(g) || (w == x && continue) + augmented_graph(g) || (w == x && continue) color_x = color[x] iszero(color_x) && continue if forbidden_colors[color_x] != v @@ -320,16 +320,16 @@ function acyclic_coloring( end end for (w, index_vw) in neighbors_with_edge_indices(g, v) # grow two-colored stars around the vertex v - !has_diagonal(g) || (v == w && continue) + augmented_graph(g) || (v == w && continue) color_w = color[w] iszero(color_w) && continue _grow_star!(v, w, index_vw, color_w, first_neighbor, forest) end for (w, index_vw) in neighbors_with_edge_indices(g, v) - !has_diagonal(g) || (v == w && continue) + augmented_graph(g) || (v == w && continue) iszero(color[w]) && continue for (x, index_wx) in neighbors_with_edge_indices(g, w) - !has_diagonal(g) || (w == x && continue) + augmented_graph(g) || (w == x && continue) color_x = color[x] (x == v || iszero(color_x)) && continue if color_x == color[v] diff --git a/src/decompression.jl b/src/decompression.jl index a48c7e34..42bf51eb 100644 --- a/src/decompression.jl +++ b/src/decompression.jl @@ -547,7 +547,7 @@ function decompress!( end # Recover the diagonal coefficients of A - if has_diagonal(ag) + if !augmented_graph(ag) for i in axes(S, 1) if !iszero(S[i, i]) A[i, i] = B[i, color[i]] @@ -616,7 +616,7 @@ function decompress!( end # Recover the diagonal coefficients of A - if has_diagonal(ag) + if !augmented_graph(ag) if uplo == :L for i in diagonal_indices # A[i, i] is the first element in column i diff --git a/src/graph.jl b/src/graph.jl index 9b53a36d..84c3353f 100644 --- a/src/graph.jl +++ b/src/graph.jl @@ -202,7 +202,7 @@ end ## Adjacency graph """ - AdjacencyGraph{T,has_diagonal} + AdjacencyGraph{T,augmented_graph} Undirected graph without self-loops representing the nonzeros of a symmetric matrix (typically a Hessian matrix). @@ -213,18 +213,18 @@ The adjacency graph of a symmetric matrix `A ∈ ℝ^{n × n}` is `G(A) = (V, E) # Constructors - AdjacencyGraph(A::SparseMatrixCSC; has_diagonal::Bool=true) + AdjacencyGraph(A::SparseMatrixCSC; augmented_graph::Bool=false) # Fields -- `S::SparsityPatternCSC{T}`: Underlying sparsity pattern, whose diagonal is empty whenever `has_diagonal` is `false` +- `S::SparsityPatternCSC{T}`: Underlying sparsity pattern, which represents an augmented graph whenever `augmented_graph` is `true`. Here, "augmented graph" means the sparsity pattern of the augmented matrix `H = [0 Jᵀ; J 0]`. - `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) # References > [_What Color Is Your Jacobian? SparsityPatternCSC Coloring for Computing Derivatives_](https://epubs.siam.org/doi/10.1137/S0036144504444711), Gebremedhin et al. (2005) """ -struct AdjacencyGraph{T<:Integer,has_diagonal} +struct AdjacencyGraph{T<:Integer,augmented_graph} S::SparsityPatternCSC{T} edge_to_index::Vector{T} end @@ -234,24 +234,24 @@ Base.eltype(::AdjacencyGraph{T}) where {T} = T function AdjacencyGraph( S::SparsityPatternCSC{T}, edge_to_index::Vector{T}=build_edge_to_index(S); - has_diagonal::Bool=true, + augmented_graph::Bool=false, ) where {T} - return AdjacencyGraph{T,has_diagonal}(S, edge_to_index) + return AdjacencyGraph{T,augmented_graph}(S, edge_to_index) end -function AdjacencyGraph(A::SparseMatrixCSC; has_diagonal::Bool=true) - return AdjacencyGraph(SparsityPatternCSC(A); has_diagonal) +function AdjacencyGraph(A::SparseMatrixCSC; augmented_graph::Bool=false) + return AdjacencyGraph(SparsityPatternCSC(A); augmented_graph) end -function AdjacencyGraph(A::AbstractMatrix; has_diagonal::Bool=true) - return AdjacencyGraph(SparseMatrixCSC(A); has_diagonal) +function AdjacencyGraph(A::AbstractMatrix; augmented_graph::Bool=false) + return AdjacencyGraph(SparseMatrixCSC(A); augmented_graph) end pattern(g::AdjacencyGraph) = g.S edge_indices(g::AdjacencyGraph) = g.edge_to_index nb_vertices(g::AdjacencyGraph) = pattern(g).n vertices(g::AdjacencyGraph) = 1:nb_vertices(g) -has_diagonal(::AdjacencyGraph{T,hd}) where {T,hd} = hd +augmented_graph(::AdjacencyGraph{T,ag}) where {T,ag} = ag function neighbors(g::AdjacencyGraph, v::Integer) S = pattern(g) @@ -267,17 +267,17 @@ function neighbors_with_edge_indices(g::AdjacencyGraph, v::Integer) return zip(neighbors_v, edges_indices_v) end -degree(g::AdjacencyGraph{T,false}, v::Integer) where {T} = g.S.colptr[v + 1] - g.S.colptr[v] +degree(g::AdjacencyGraph{T,true}, v::Integer) where {T} = g.S.colptr[v + 1] - g.S.colptr[v] -function degree(g::AdjacencyGraph{T,true}, v::Integer) where {T} +function degree(g::AdjacencyGraph{T,false}, v::Integer) where {T} neigh = neighbors(g, v) has_selfloop = insorted(v, neigh) return g.S.colptr[v + 1] - g.S.colptr[v] - has_selfloop end -nb_edges(g::AdjacencyGraph{T,false}) where {T} = nnz(g.S) ÷ 2 +nb_edges(g::AdjacencyGraph{T,true}) where {T} = nnz(g.S) ÷ 2 -function nb_edges(g::AdjacencyGraph{T,true}) where {T} +function nb_edges(g::AdjacencyGraph{T,false}) where {T} ne = 0 for v in vertices(g) ne += degree(g, v) @@ -290,7 +290,7 @@ minimum_degree(g::AdjacencyGraph) = minimum(Base.Fix1(degree, g), vertices(g)) function has_neighbor(g::AdjacencyGraph, v::Integer, u::Integer) for w in neighbors(g, v) - !has_diagonal(g) || (v == w && continue) + augmented_graph(g) || (v == w && continue) if w == u return true end diff --git a/src/interface.jl b/src/interface.jl index b8d4ce54..cc57bcc7 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -279,7 +279,7 @@ function _coloring( symmetric_pattern::Bool; forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing, ) - ag = AdjacencyGraph(A; has_diagonal=true) + ag = AdjacencyGraph(A; augmented_graph=false) color_and_star_set_by_order = map(algo.orders) do order vertices_in_order = vertices(ag, order) return star_coloring(ag, vertices_in_order, algo.postprocessing; forced_colors) @@ -300,7 +300,7 @@ function _coloring( decompression_eltype::Type{R}, symmetric_pattern::Bool, ) where {R} - ag = AdjacencyGraph(A; has_diagonal=true) + ag = AdjacencyGraph(A; augmented_graph=false) color_and_tree_set_by_order = map(algo.orders) do order vertices_in_order = vertices(ag, order) return acyclic_coloring(ag, vertices_in_order, algo.postprocessing) @@ -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; has_diagonal=false) + ag = AdjacencyGraph(A_and_Aᵀ, edge_to_index; augmented_graph=true) outputs_by_order = map(algo.orders) do order vertices_in_order = vertices(ag, order) _color, _star_set = star_coloring( @@ -370,7 +370,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; has_diagonal=false) + ag = AdjacencyGraph(A_and_Aᵀ, edge_to_index; 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) diff --git a/src/order.jl b/src/order.jl index aab5bb23..20884c6a 100644 --- a/src/order.jl +++ b/src/order.jl @@ -334,7 +334,7 @@ function vertices( π[index] = u for v in neighbors(g, u) - !has_diagonal(g) || (u == v && continue) + augmented_graph(g) || (u == v && continue) dv = degrees[v] dv == -1 && continue update_bucket!(db, v, dv; degtype, direction) diff --git a/src/postprocessing.jl b/src/postprocessing.jl index 217b663b..8174d863 100644 --- a/src/postprocessing.jl +++ b/src/postprocessing.jl @@ -13,7 +13,7 @@ function postprocess!( color_used = zeros(Bool, nb_colors) # nonzero diagonal coefficients force the use of their respective color (there can be no neutral colors if the diagonal is fully nonzero) - if has_diagonal(g) + if !augmented_graph(g) for i in axes(S, 1) if !iszero(S[i, i]) color_used[color[i]] = true diff --git a/src/result.jl b/src/result.jl index 7a3e3b68..89da22d6 100644 --- a/src/result.jl +++ b/src/result.jl @@ -411,7 +411,7 @@ function TreeSetColoringResult( diagonal_nzind = T[] ndiag = 0 - if has_diagonal(ag) + if !augmented_graph(ag) for j in axes(S, 2) for k in nzrange(S, j) i = rv[k] diff --git a/test/graph.jl b/test/graph.jl index 9c8f66c1..e9784e56 100644 --- a/test/graph.jl +++ b/test/graph.jl @@ -166,7 +166,7 @@ end; @test degree(g, 7) == 6 @test degree(g, 8) == 6 - g = AdjacencyGraph(transpose(A) * A; has_diagonal=false) + g = AdjacencyGraph(transpose(A) * A; augmented_graph=true) # wrong degree @test degree(g, 1) == 4 @test degree(g, 2) == 4 diff --git a/test/order.jl b/test/order.jl index 838c32a9..9cae16f1 100644 --- a/test/order.jl +++ b/test/order.jl @@ -61,16 +61,14 @@ end; end; @testset "LargestFirst" begin - for has_diagonal in (false, true) - A = sparse([ - 0 1 0 0 - 1 0 1 1 - 0 1 0 1 - 0 1 1 0 - ]) - ag = AdjacencyGraph(A; has_diagonal) - @test vertices(ag, LargestFirst()) == [2, 3, 4, 1] - end + A = sparse([ + 0 1 0 0 + 1 0 1 1 + 0 1 0 1 + 0 1 1 0 + ]) + ag = AdjacencyGraph(A) + @test vertices(ag, LargestFirst()) == [2, 3, 4, 1] A = sparse([ 1 1 0 0 From d4b2e00277a63b3db198132322dbafd382498867 Mon Sep 17 00:00:00 2001 From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com> Date: Sat, 1 Nov 2025 10:34:56 +0100 Subject: [PATCH 13/39] Try tests on Julia 1.12 (#276) --- .github/workflows/Test-GPU.yml | 2 +- .github/workflows/Test.yml | 2 +- Project.toml | 6 +++--- test/runtests.jl | 17 ++++++++++------- 4 files changed, 15 insertions(+), 12 deletions(-) diff --git a/.github/workflows/Test-GPU.yml b/.github/workflows/Test-GPU.yml index f5a4fbc8..1795d222 100644 --- a/.github/workflows/Test-GPU.yml +++ b/.github/workflows/Test-GPU.yml @@ -21,7 +21,7 @@ jobs: env: CUDA_VISIBLE_DEVICES: 1 JULIA_DEPOT_PATH: /scratch/github-actions/julia_depot_alexis - JULIA_SMC_TEST_GROUP: "GPU" + JULIA_SMC_TEST_GROUP: 'GPU' strategy: matrix: julia-version: ['1.10', '1.11'] diff --git a/.github/workflows/Test.yml b/.github/workflows/Test.yml index 3e17ae97..e582efa3 100644 --- a/.github/workflows/Test.yml +++ b/.github/workflows/Test.yml @@ -20,7 +20,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - julia-version: ['1.10', '1.11'] + julia-version: ['1.10', '1.11', '1.12'] steps: - uses: actions/checkout@v5 diff --git a/Project.toml b/Project.toml index 864fae83..e34ce604 100644 --- a/Project.toml +++ b/Project.toml @@ -31,9 +31,9 @@ CliqueTrees = "1" Colors = "0.12.11, 0.13" DocStringExtensions = "0.8,0.9" JuMP = "1.29.1" -LinearAlgebra = "<0.0.1, 1" +LinearAlgebra = "1" MathOptInterface = "1.45.0" PrecompileTools = "1.2.1" -Random = "<0.0.1, 1" -SparseArrays = "<0.0.1, 1" +Random = "1" +SparseArrays = "1" julia = "1.10" diff --git a/test/runtests.jl b/test/runtests.jl index fa5cc7ff..7a070756 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -19,16 +19,16 @@ include("utils.jl") else @testset verbose = true "Code quality" begin @testset "Aqua" begin - Aqua.test_all(SparseMatrixColorings; stale_deps=(; ignore=[:Requires],)) + Aqua.test_all(SparseMatrixColorings; stale_deps = (; ignore = [:Requires])) end @testset "JET" begin - JET.test_package(SparseMatrixColorings; target_defined_modules=true) + JET.test_package(SparseMatrixColorings; target_defined_modules = true) end # @testset "JuliaFormatter" begin - # TODO: switch to Runic (temporarily deactivated) - # @test JuliaFormatter.format( - # SparseMatrixColorings; verbose=false, overwrite=false - # ) + # TODO: switch to Runic (temporarily deactivated) + # @test JuliaFormatter.format( + # SparseMatrixColorings; verbose=false, overwrite=false + # ) # end @testset "Doctests" begin Documenter.doctest(SparseMatrixColorings) @@ -88,7 +88,10 @@ include("utils.jl") end @testset verbose = true "Performance" begin @testset "Type stability" begin - include("type_stability.jl") + if VERSION < v"1.12" + # TODO: fix JET misbehaving + include("type_stability.jl") + end end @testset "Allocations" begin include("allocations.jl") From 7a831815d4812a00a8537d791507fb6a3b316bfc Mon Sep 17 00:00:00 2001 From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com> Date: Sat, 1 Nov 2025 13:08:50 +0100 Subject: [PATCH 14/39] Bump SMC to v0.4.23 (#277) --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index e34ce604..0e01b3d4 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SparseMatrixColorings" uuid = "0a514795-09f3-496d-8182-132a7b665d35" authors = ["Guillaume Dalle", "Alexis Montoison"] -version = "0.4.22" +version = "0.4.23" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From 85e2d057f71c18592f395b45446cc916cbc1e4f9 Mon Sep 17 00:00:00 2001 From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com> Date: Sat, 1 Nov 2025 14:15:43 +0100 Subject: [PATCH 15/39] Patch failing alloc test (#278) --- test/allocations.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/allocations.jl b/test/allocations.jl index caed10f8..31305e63 100644 --- a/test/allocations.jl +++ b/test/allocations.jl @@ -154,7 +154,7 @@ end b64 = @b fast_coloring(A64, problem, algo) b32 = @b fast_coloring(A32, problem, algo) # check that we allocate no more than 50% + epsilon with Int32 - @test b32.bytes < 0.6 * b64.bytes + @test b32.bytes < 0.7 * b64.bytes end end end; From 2cf5f8bc8a937df357b46d2830352496c1824d5e Mon Sep 17 00:00:00 2001 From: Alexis Montoison <35051714+amontoison@users.noreply.github.com> Date: Mon, 3 Nov 2025 01:46:48 -0600 Subject: [PATCH 16/39] Add postprocess_with_star_set! and postprocess_with_tree_set! (#279) --- src/postprocessing.jl | 208 +++++++++++++++++++++++------------------- 1 file changed, 116 insertions(+), 92 deletions(-) diff --git a/src/postprocessing.jl b/src/postprocessing.jl index 8174d863..9b0443b2 100644 --- a/src/postprocessing.jl +++ b/src/postprocessing.jl @@ -22,99 +22,11 @@ function postprocess!( end if star_or_tree_set isa StarSet - # only the colors of the hubs are used - (; star, hub) = star_or_tree_set - nb_trivial_stars = 0 - - # Iterate through all non-trivial stars - for s in eachindex(hub) - h = hub[s] - if h > 0 - color_used[color[h]] = true - else - nb_trivial_stars += 1 - end - end - - # Process the trivial stars (if any) - if nb_trivial_stars > 0 - rvS = rowvals(S) - for j in axes(S, 2) - for k in nzrange(S, j) - i = rvS[k] - if i > j - index_ij = edge_to_index[k] - s = star[index_ij] - h = hub[s] - if h < 0 - h = abs(h) - spoke = h == j ? i : j - if color_used[color[spoke]] - # Switch the hub and the spoke to possibly avoid adding one more used color - hub[s] = spoke - else - # Keep the current hub - color_used[color[h]] = true - end - end - end - end - end - end + # star_or_tree_set is a StarSet + postprocess_with_star_set!(g, color_used, color, star_or_tree_set) else - # only the colors of non-leaf vertices are used - (; reverse_bfs_orders, is_star, tree_edge_indices, nt) = star_or_tree_set - nb_trivial_trees = 0 - - # Iterate through all non-trivial trees - for k in 1:nt - # Position of the first edge in the tree - first = tree_edge_indices[k] - - # Total number of edges in the tree - ne_tree = tree_edge_indices[k + 1] - first - - # Check if we have more than one edge in the tree (non-trivial tree) - if ne_tree > 1 - # Determine if the tree is a star - if is_star[k] - # It is a non-trivial star and only the color of the hub is needed - (_, hub) = reverse_bfs_orders[first] - color_used[color[hub]] = true - else - # It is not a star and both colors are needed during the decompression - (i, j) = reverse_bfs_orders[first] - color_used[color[i]] = true - color_used[color[j]] = true - end - else - nb_trivial_trees += 1 - end - end - - # Process the trivial trees (if any) - if nb_trivial_trees > 0 - for k in 1:nt - # Position of the first edge in the tree - first = tree_edge_indices[k] - - # Total number of edges in the tree - ne_tree = tree_edge_indices[k + 1] - first - - # Check if we have exactly one edge in the tree - if ne_tree == 1 - (i, j) = reverse_bfs_orders[first] - if color_used[color[i]] - # Make i the root to avoid possibly adding one more used color - # Switch it with the (only) leaf - reverse_bfs_orders[first] = (j, i) - else - # Keep j as the root - color_used[color[j]] = true - end - end - end - end + # star_or_tree_set is a TreeSet + postprocess_with_tree_set!(color_used, color, star_or_tree_set) end # if at least one of the colors is useless, modify the color assignments of vertices @@ -144,3 +56,115 @@ function postprocess!( end return color end + +function postprocess_with_star_set!( + g::AdjacencyGraph, + color_used::Vector{Bool}, + color::AbstractVector{<:Integer}, + star_set::StarSet, +) + S = pattern(g) + edge_to_index = edge_indices(g) + + # only the colors of the hubs are used + (; star, hub) = star_set + nb_trivial_stars = 0 + + # Iterate through all non-trivial stars + for s in eachindex(hub) + h = hub[s] + if h > 0 + color_used[color[h]] = true + else + nb_trivial_stars += 1 + end + end + + # Process the trivial stars (if any) + if nb_trivial_stars > 0 + rvS = rowvals(S) + for j in axes(S, 2) + for k in nzrange(S, j) + i = rvS[k] + if i > j + index_ij = edge_to_index[k] + s = star[index_ij] + h = hub[s] + if h < 0 + h = abs(h) + spoke = h == j ? i : j + if color_used[color[spoke]] + # Switch the hub and the spoke to possibly avoid adding one more used color + hub[s] = spoke + else + # Keep the current hub + color_used[color[h]] = true + end + end + end + end + end + end + return color_used +end + +function postprocess_with_tree_set!( + color_used::Vector{Bool}, + color::AbstractVector{<:Integer}, + tree_set::TreeSet, +) + # only the colors of non-leaf vertices are used + (; reverse_bfs_orders, is_star, tree_edge_indices, nt) = tree_set + nb_trivial_trees = 0 + + # Iterate through all non-trivial trees + for k in 1:nt + # Position of the first edge in the tree + first = tree_edge_indices[k] + + # Total number of edges in the tree + ne_tree = tree_edge_indices[k + 1] - first + + # Check if we have more than one edge in the tree (non-trivial tree) + if ne_tree > 1 + # Determine if the tree is a star + if is_star[k] + # It is a non-trivial star and only the color of the hub is needed + (_, hub) = reverse_bfs_orders[first] + color_used[color[hub]] = true + else + # It is not a star and both colors are needed during the decompression + (i, j) = reverse_bfs_orders[first] + color_used[color[i]] = true + color_used[color[j]] = true + end + else + nb_trivial_trees += 1 + end + end + + # Process the trivial trees (if any) + if nb_trivial_trees > 0 + for k in 1:nt + # Position of the first edge in the tree + first = tree_edge_indices[k] + + # Total number of edges in the tree + ne_tree = tree_edge_indices[k + 1] - first + + # Check if we have exactly one edge in the tree + if ne_tree == 1 + (i, j) = reverse_bfs_orders[first] + if color_used[color[i]] + # Make i the root to avoid possibly adding one more used color + # Switch it with the (only) leaf + reverse_bfs_orders[first] = (j, i) + else + # Keep j as the root + color_used[color[j]] = true + end + end + end + end + return color_used +end From 9b52faccdaae41d3ce27158434cc5597d1a61a36 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 24 Nov 2025 15:44:30 +1300 Subject: [PATCH 17/39] Bump actions/checkout from 5 to 6 (#282) Bumps [actions/checkout](https://github.com/actions/checkout) from 5 to 6. - [Release notes](https://github.com/actions/checkout/releases) - [Changelog](https://github.com/actions/checkout/blob/main/CHANGELOG.md) - [Commits](https://github.com/actions/checkout/compare/v5...v6) --- updated-dependencies: - dependency-name: actions/checkout dependency-version: '6' dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/Benchmark.yml | 2 +- .github/workflows/Documentation.yml | 2 +- .github/workflows/Test-GPU.yml | 2 +- .github/workflows/Test.yml | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/.github/workflows/Benchmark.yml b/.github/workflows/Benchmark.yml index 23dd4973..30e2e1e5 100644 --- a/.github/workflows/Benchmark.yml +++ b/.github/workflows/Benchmark.yml @@ -12,7 +12,7 @@ jobs: runs-on: ubuntu-latest if: contains(github.event.pull_request.labels.*.name, 'benchmark') steps: - - uses: actions/checkout@v5 + - uses: actions/checkout@v6 - uses: julia-actions/setup-julia@v2 with: version: "1" diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml index c5527ae9..b73255ea 100644 --- a/.github/workflows/Documentation.yml +++ b/.github/workflows/Documentation.yml @@ -19,7 +19,7 @@ jobs: statuses: write runs-on: ubuntu-latest steps: - - uses: actions/checkout@v5 + - uses: actions/checkout@v6 - uses: julia-actions/setup-julia@v2 with: version: '1' diff --git a/.github/workflows/Test-GPU.yml b/.github/workflows/Test-GPU.yml index 1795d222..508da480 100644 --- a/.github/workflows/Test-GPU.yml +++ b/.github/workflows/Test-GPU.yml @@ -26,7 +26,7 @@ jobs: matrix: julia-version: ['1.10', '1.11'] steps: - - uses: actions/checkout@v5 + - uses: actions/checkout@v6 - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.julia-version }} diff --git a/.github/workflows/Test.yml b/.github/workflows/Test.yml index e582efa3..d979b583 100644 --- a/.github/workflows/Test.yml +++ b/.github/workflows/Test.yml @@ -23,7 +23,7 @@ jobs: julia-version: ['1.10', '1.11', '1.12'] steps: - - uses: actions/checkout@v5 + - uses: actions/checkout@v6 - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.julia-version }} From e2f039da1d511dac1785d8247e8735a40d367e2d Mon Sep 17 00:00:00 2001 From: Alexis Montoison <35051714+amontoison@users.noreply.github.com> Date: Mon, 15 Dec 2025 02:37:48 -0600 Subject: [PATCH 18/39] Replace is_solved_and_feasible with assert_is_solved_and_feasible (#284) --- ext/SparseMatrixColoringsJuMPExt.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ext/SparseMatrixColoringsJuMPExt.jl b/ext/SparseMatrixColoringsJuMPExt.jl index ca4f5e7b..77b74dbd 100644 --- a/ext/SparseMatrixColoringsJuMPExt.jl +++ b/ext/SparseMatrixColoringsJuMPExt.jl @@ -3,7 +3,7 @@ module SparseMatrixColoringsJuMPExt using ADTypes: ADTypes using JuMP: Model, - is_solved_and_feasible, + assert_is_solved_and_feasible, optimize!, primal_status, set_silent, @@ -44,7 +44,7 @@ function optimal_distance2_coloring( optimize!(model) if assert_solved # assert feasibility and optimality - @assert is_solved_and_feasible(model) + assert_is_solved_and_feasible(model) else # only assert feasibility @assert primal_status(model) == MOI.FEASIBLE_POINT From 0702578b9ba69bb98dfe9203b7ebcbc90b4a8225 Mon Sep 17 00:00:00 2001 From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com> Date: Mon, 15 Dec 2025 10:41:16 +0100 Subject: [PATCH 19/39] Fix JET tests (#285) * Fix JET tests * Move JuliaFormatter to action * Fix formatting --- .github/workflows/Format.yml | 12 ++++++++++++ ext/SparseMatrixColoringsColorsExt.jl | 5 ++--- src/postprocessing.jl | 4 +--- test/Project.toml | 1 - test/runtests.jl | 15 +++++---------- 5 files changed, 20 insertions(+), 17 deletions(-) create mode 100644 .github/workflows/Format.yml diff --git a/.github/workflows/Format.yml b/.github/workflows/Format.yml new file mode 100644 index 00000000..dd2b2c50 --- /dev/null +++ b/.github/workflows/Format.yml @@ -0,0 +1,12 @@ +name: Format suggestions +on: + pull_request: + # this argument is not required if you don't use the `suggestion-label` input + types: [opened, reopened, synchronize, labeled, unlabeled] +jobs: + code-style: + runs-on: ubuntu-latest + steps: + - uses: julia-actions/julia-format@v4 + with: + version: '1' # Set `version` to '1.0.54' if you need to use JuliaFormatter.jl v1.0.54 (default: '1') diff --git a/ext/SparseMatrixColoringsColorsExt.jl b/ext/SparseMatrixColoringsColorsExt.jl index 7f844e89..c934b5e9 100644 --- a/ext/SparseMatrixColoringsColorsExt.jl +++ b/ext/SparseMatrixColoringsColorsExt.jl @@ -271,9 +271,8 @@ function show_colors!( A_ccolor_indices = mod1.(column_colors(res), length(colorscheme)) A_rcolor_indices = mod1.(row_shift .+ row_colors(res), length(colorscheme)) B_ccolor_indices = mod1.(1:maximum(column_colors(res)), length(colorscheme)) - B_rcolor_indices = mod1.( - (row_shift + 1):(row_shift + maximum(row_colors(res))), length(colorscheme) - ) + B_rcolor_indices = + mod1.((row_shift + 1):(row_shift + maximum(row_colors(res))), length(colorscheme)) A_ccolors = colorscheme[A_ccolor_indices] A_rcolors = colorscheme[A_rcolor_indices] B_ccolors = colorscheme[B_ccolor_indices] diff --git a/src/postprocessing.jl b/src/postprocessing.jl index 9b0443b2..7b2b58f5 100644 --- a/src/postprocessing.jl +++ b/src/postprocessing.jl @@ -109,9 +109,7 @@ function postprocess_with_star_set!( end function postprocess_with_tree_set!( - color_used::Vector{Bool}, - color::AbstractVector{<:Integer}, - tree_set::TreeSet, + color_used::Vector{Bool}, color::AbstractVector{<:Integer}, tree_set::TreeSet ) # only the colors of non-leaf vertices are used (; reverse_bfs_orders, is_star, tree_edge_indices, nt) = tree_set diff --git a/test/Project.toml b/test/Project.toml index 217d62c9..284f410b 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -15,7 +15,6 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" -JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MatrixDepot = "b51810bb-c9f3-55da-ae3c-350fc1fbce05" MiniZinc = "a7f392d2-6c35-496e-b8cc-0974fbfcbf91" diff --git a/test/runtests.jl b/test/runtests.jl index 7a070756..9494e5e2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,7 +1,6 @@ using Aqua using Documenter using JET -using JuliaFormatter using SparseMatrixColorings using Test @@ -19,17 +18,13 @@ include("utils.jl") else @testset verbose = true "Code quality" begin @testset "Aqua" begin - Aqua.test_all(SparseMatrixColorings; stale_deps = (; ignore = [:Requires])) + Aqua.test_all(SparseMatrixColorings; undocumented_names=true) end @testset "JET" begin - JET.test_package(SparseMatrixColorings; target_defined_modules = true) - end - # @testset "JuliaFormatter" begin - # TODO: switch to Runic (temporarily deactivated) - # @test JuliaFormatter.format( - # SparseMatrixColorings; verbose=false, overwrite=false - # ) - # end + JET.test_package( + SparseMatrixColorings; target_modules=(SparseMatrixColorings,) + ) + end @testset "Doctests" begin Documenter.doctest(SparseMatrixColorings) end From e04bd4e815967bab2364222b7546e0f8329298db Mon Sep 17 00:00:00 2001 From: Alexis Montoison <35051714+amontoison@users.noreply.github.com> Date: Wed, 31 Dec 2025 11:17:41 +0100 Subject: [PATCH 20/39] Use diagonal_indices in the general decompress! for acyclic coloring (#287) --- src/decompression.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/decompression.jl b/src/decompression.jl index 42bf51eb..f8fd610a 100644 --- a/src/decompression.jl +++ b/src/decompression.jl @@ -534,7 +534,8 @@ end function decompress!( A::AbstractMatrix, B::AbstractMatrix, result::TreeSetColoringResult, uplo::Symbol=:F ) - (; ag, color, reverse_bfs_orders, tree_edge_indices, nt, buffer) = result + (; ag, color, reverse_bfs_orders, tree_edge_indices, nt, diagonal_indices, buffer) = + result (; S) = ag uplo == :F && check_same_pattern(A, S) R = eltype(A) @@ -548,10 +549,8 @@ function decompress!( # Recover the diagonal coefficients of A if !augmented_graph(ag) - for i in axes(S, 1) - if !iszero(S[i, i]) - A[i, i] = B[i, color[i]] - end + for i in diagonal_indices + A[i, i] = B[i, color[i]] end end From 662df63ee283a6ddd774cba7a7613ccc971a1d1c Mon Sep 17 00:00:00 2001 From: Alexis Montoison <35051714+amontoison@users.noreply.github.com> Date: Wed, 31 Dec 2025 11:18:34 +0100 Subject: [PATCH 21/39] Enhance decompress! for bicoloring (#288) --- src/decompression.jl | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/src/decompression.jl b/src/decompression.jl index f8fd610a..1f8e88db 100644 --- a/src/decompression.jl +++ b/src/decompression.jl @@ -769,10 +769,22 @@ end function decompress!( A::AbstractMatrix, Br::AbstractMatrix, Bc::AbstractMatrix, result::BicoloringResult ) + (; large_colptr, large_rowval, symmetric_result) = result m, n = size(A) + R = eltype(A) + fill!(A, zero(R)) + nzval = Vector{R}(undef, length(large_rowval)) + A_and_noAᵀ = SparseMatrixCSC(m + n, m + n, large_colptr, large_rowval, nzval) Br_and_Bc = _join_compressed!(result, Br, Bc) - A_and_Aᵀ = decompress(Br_and_Bc, result.symmetric_result) - copyto!(A, A_and_Aᵀ[(n + 1):(n + m), 1:n]) # original matrix in bottom left corner + decompress!(A_and_noAᵀ, Br_and_Bc, symmetric_result, :L) + rvA = rowvals(A_and_noAᵀ) + nzA = nonzeros(A_and_noAᵀ) + for j in 1:n + for k in nzrange(A_and_noAᵀ, j) + i = rvA[k] + A[i - n, j] = nzA[k] + end + end return A end @@ -781,10 +793,10 @@ function decompress!( ) (; large_colptr, large_rowval, symmetric_result) = result m, n = size(A) - Br_and_Bc = _join_compressed!(result, Br, Bc) # pretend A is larger A_and_noAᵀ = SparseMatrixCSC(m + n, m + n, large_colptr, large_rowval, A.nzval) # decompress lower triangle only + Br_and_Bc = _join_compressed!(result, Br, Bc) decompress!(A_and_noAᵀ, Br_and_Bc, symmetric_result, :L) return A end From a5f32f1708f4d65ce52a6d1bca38a60dcb56dfe8 Mon Sep 17 00:00:00 2001 From: Alexis Montoison <35051714+amontoison@users.noreply.github.com> Date: Sun, 4 Jan 2026 22:20:14 +0100 Subject: [PATCH 22/39] Add a compat entry for CUDA.jl in test/Project.toml (#293) --- .github/workflows/Test-GPU.yml | 4 ++-- test/Project.toml | 3 +++ 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/.github/workflows/Test-GPU.yml b/.github/workflows/Test-GPU.yml index 508da480..aa30a5fc 100644 --- a/.github/workflows/Test-GPU.yml +++ b/.github/workflows/Test-GPU.yml @@ -24,7 +24,7 @@ jobs: JULIA_SMC_TEST_GROUP: 'GPU' strategy: matrix: - julia-version: ['1.10', '1.11'] + julia-version: ['lts', '1'] steps: - uses: actions/checkout@v6 - uses: julia-actions/setup-julia@v2 @@ -32,7 +32,7 @@ jobs: version: ${{ matrix.julia-version }} arch: x64 - uses: julia-actions/julia-downgrade-compat@v2 - if: ${{ matrix.version == '1.10' }} + if: ${{ matrix.version == 'lts' }} with: skip: LinearAlgebra, Random, SparseArrays - uses: julia-actions/cache@v2 diff --git a/test/Project.toml b/test/Project.toml index 284f410b..4c505efa 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -23,3 +23,6 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[compat] +CUDA = "=5.9.5" From 31999f1409f5398ccd6224fb76f8dccdfe2d50c6 Mon Sep 17 00:00:00 2001 From: Alexis Montoison <35051714+amontoison@users.noreply.github.com> Date: Sun, 4 Jan 2026 22:46:40 +0100 Subject: [PATCH 23/39] Add nb_self_loops in AdjacencyGraph (#290) --- docs/src/dev.md | 2 +- src/decompression.jl | 59 ++++++++++++++++++++++---------------------- src/graph.jl | 25 ++++++++++--------- src/interface.jl | 4 +-- src/matrices.jl | 59 ++++++++++++++++++++++++++++++++++---------- src/result.jl | 16 ++++++------ test/matrices.jl | 41 +++++++++++++++++++++++++----- 7 files changed, 135 insertions(+), 71 deletions(-) diff --git a/docs/src/dev.md b/docs/src/dev.md index 2f253dbe..82579dd6 100644 --- a/docs/src/dev.md +++ b/docs/src/dev.md @@ -58,7 +58,7 @@ SparseMatrixColorings.valid_dynamic_order ```@docs SparseMatrixColorings.respectful_similar SparseMatrixColorings.matrix_versions -SparseMatrixColorings.same_pattern +SparseMatrixColorings.compatible_pattern ``` ## Visualization diff --git a/src/decompression.jl b/src/decompression.jl index 1f8e88db..1528c885 100644 --- a/src/decompression.jl +++ b/src/decompression.jl @@ -339,9 +339,9 @@ end ## ColumnColoringResult function decompress!(A::AbstractMatrix, B::AbstractMatrix, result::ColumnColoringResult) - (; color) = result - S = result.bg.S2 - check_same_pattern(A, S) + (; bg, color) = result + S = bg.S2 + check_compatible_pattern(A, bg) fill!(A, zero(eltype(A))) rvS = rowvals(S) for j in axes(S, 2) @@ -357,9 +357,9 @@ end function decompress_single_color!( A::AbstractMatrix, b::AbstractVector, c::Integer, result::ColumnColoringResult ) - (; group) = result - S = result.bg.S2 - check_same_pattern(A, S) + (; bg, group) = result + S = bg.S2 + check_compatible_pattern(A, bg) rvS = rowvals(S) for j in group[c] for k in nzrange(S, j) @@ -371,9 +371,9 @@ function decompress_single_color!( end function decompress!(A::SparseMatrixCSC, B::AbstractMatrix, result::ColumnColoringResult) - (; compressed_indices) = result - S = result.bg.S2 - check_same_pattern(A, S) + (; bg, compressed_indices) = result + S = bg.S2 + check_compatible_pattern(A, bg) nzA = nonzeros(A) for k in eachindex(nzA, compressed_indices) nzA[k] = B[compressed_indices[k]] @@ -384,9 +384,9 @@ end function decompress_single_color!( A::SparseMatrixCSC, b::AbstractVector, c::Integer, result::ColumnColoringResult ) - (; group) = result - S = result.bg.S2 - check_same_pattern(A, S) + (; bg, group) = result + S = bg.S2 + check_compatible_pattern(A, bg) rvS = rowvals(S) nzA = nonzeros(A) for j in group[c] @@ -401,9 +401,9 @@ end ## RowColoringResult function decompress!(A::AbstractMatrix, B::AbstractMatrix, result::RowColoringResult) - (; color) = result - S = result.bg.S2 - check_same_pattern(A, S) + (; bg, color) = result + S = bg.S2 + check_compatible_pattern(A, bg) fill!(A, zero(eltype(A))) rvS = rowvals(S) for j in axes(S, 2) @@ -419,9 +419,9 @@ end function decompress_single_color!( A::AbstractMatrix, b::AbstractVector, c::Integer, result::RowColoringResult ) - (; group) = result - S, Sᵀ = result.bg.S2, result.bg.S1 - check_same_pattern(A, S) + (; bg, group) = result + S, Sᵀ = bg.S2, bg.S1 + check_compatible_pattern(A, bg) rvSᵀ = rowvals(Sᵀ) for i in group[c] for k in nzrange(Sᵀ, i) @@ -433,9 +433,9 @@ function decompress_single_color!( end function decompress!(A::SparseMatrixCSC, B::AbstractMatrix, result::RowColoringResult) - (; compressed_indices) = result - S = result.bg.S2 - check_same_pattern(A, S) + (; bg, compressed_indices) = result + S = bg.S2 + check_compatible_pattern(A, bg) nzA = nonzeros(A) for k in eachindex(nzA, compressed_indices) nzA[k] = B[compressed_indices[k]] @@ -450,7 +450,7 @@ function decompress!( ) (; ag, compressed_indices) = result (; S) = ag - uplo == :F && check_same_pattern(A, S) + check_compatible_pattern(A, ag, uplo) fill!(A, zero(eltype(A))) rvS = rowvals(S) @@ -474,7 +474,7 @@ function decompress_single_color!( ) (; ag, compressed_indices, group) = result (; S) = ag - uplo == :F && check_same_pattern(A, S) + check_compatible_pattern(A, ag, uplo) lower_index = (c - 1) * S.n + 1 upper_index = c * S.n @@ -508,8 +508,8 @@ function decompress!( (; ag, compressed_indices) = result (; S) = ag nzA = nonzeros(A) + check_compatible_pattern(A, ag, uplo) if uplo == :F - check_same_pattern(A, S) for k in eachindex(nzA, compressed_indices) nzA[k] = B[compressed_indices[k]] end @@ -537,7 +537,7 @@ function decompress!( (; ag, color, reverse_bfs_orders, tree_edge_indices, nt, diagonal_indices, buffer) = result (; S) = ag - uplo == :F && check_same_pattern(A, S) + check_compatible_pattern(A, ag, uplo) R = eltype(A) fill!(A, zero(R)) @@ -606,7 +606,7 @@ function decompress!( (; S) = ag A_colptr = A.colptr nzA = nonzeros(A) - uplo == :F && check_same_pattern(A, S) + check_compatible_pattern(A, ag, uplo) if eltype(buffer) == R buffer_right_type = buffer @@ -707,9 +707,10 @@ function decompress!( result::LinearSystemColoringResult, uplo::Symbol=:F, ) - (; color, strict_upper_nonzero_inds, M_factorization, strict_upper_nonzeros_A) = result - S = result.ag.S - uplo == :F && check_same_pattern(A, S) + (; ag, color, strict_upper_nonzero_inds, M_factorization, strict_upper_nonzeros_A) = + result + S = ag.S + check_compatible_pattern(A, ag, uplo) # TODO: for some reason I cannot use ldiv! with a sparse QR strict_upper_nonzeros_A = M_factorization \ vec(B) diff --git a/src/graph.jl b/src/graph.jl index 84c3353f..b93a864b 100644 --- a/src/graph.jl +++ b/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) @@ -275,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)) diff --git a/src/interface.jl b/src/interface.jl index cc57bcc7..936b974c 100644 --- a/src/interface.jl +++ b/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( @@ -370,7 +370,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) diff --git a/src/matrices.jl b/src/matrices.jl index eaf79882..453e3062 100644 --- a/src/matrices.jl +++ b/src/matrices.jl @@ -62,24 +62,57 @@ function respectful_similar(A::Union{Symmetric,Hermitian}, ::Type{T}) where {T} end """ - same_pattern(A, B) + compatible_pattern(A::AbstractMatrix, bg::BipartiteGraph) + compatible_pattern(A::AbstractMatrix, ag::AdjacencyGraph, uplo::Symbol) -Perform a partial equality check on the sparsity patterns of `A` and `B`: +Perform a coarse compatibility check between the sparsity pattern of `A` +and the reference sparsity pattern encoded in a graph structure. -- if the return is `true`, they might have the same sparsity pattern but we're not sure -- if the return is `false`, they definitely don't have the same sparsity pattern +This function only checks necessary conditions for the two sparsity patterns to match. +In particular, it may return `true` even if the patterns are not identical. + +When A is a `SparseMatrixCSC`, additional checks on the sparsity structure are performed. + +# Return value +- `true` : the sparsity patterns are potentially compatible +- `false` : the sparsity patterns are definitely incompatible """ -same_pattern(A, B) = size(A) == size(B) +compatible_pattern(A::AbstractMatrix, bg::BipartiteGraph) = size(A) == size(bg.S2) +function compatible_pattern(A::SparseMatrixCSC, bg::BipartiteGraph) + return size(A) == size(bg.S2) && nnz(A) == nnz(bg.S2) +end + +function compatible_pattern(A::AbstractMatrix, ag::AdjacencyGraph, uplo::Symbol) + return size(A) == size(ag.S) +end + +function compatible_pattern(A::SparseMatrixCSC, ag::AdjacencyGraph, uplo::Symbol) + nnzS = (uplo == :L || uplo == :U) ? (nb_edges(ag) + ag.nb_self_loops) : nnz(ag.S) + return size(A) == size(ag.S) && nnz(A) == nnzS +end -function same_pattern( - A::Union{SparseMatrixCSC,SparsityPatternCSC}, - B::Union{SparseMatrixCSC,SparsityPatternCSC}, -) - return size(A) == size(B) && nnz(A) == nnz(B) +function check_compatible_pattern(A::AbstractMatrix, bg::BipartiteGraph) + if !compatible_pattern(A, bg) + throw(DimensionMismatch("`A` and `bg.S2` must have the same sparsity pattern.")) + end end -function check_same_pattern(A, S) - if !same_pattern(A, S) - throw(DimensionMismatch("`A` and `S` must have the same sparsity pattern.")) +function check_compatible_pattern(A::AbstractMatrix, ag::AdjacencyGraph, uplo::Symbol) + if !compatible_pattern(A, ag, uplo) + if uplo == :L + throw( + DimensionMismatch( + "`A` and `tril(ag.S)` must have the same sparsity pattern." + ), + ) + elseif uplo == :U + throw( + DimensionMismatch( + "`A` and `triu(ag.S)` must have the same sparsity pattern." + ), + ) + else # uplo == :F + throw(DimensionMismatch("`A` and `ag.S` must have the same sparsity pattern.")) + end end end diff --git a/src/result.jl b/src/result.jl index 89da22d6..c74ceaa1 100644 --- a/src/result.jl +++ b/src/result.jl @@ -401,31 +401,31 @@ function TreeSetColoringResult( decompression_eltype::Type{R}, ) 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) # Vector for the decompression of the diagonal coefficients - diagonal_indices = T[] - diagonal_nzind = T[] - ndiag = 0 + diagonal_indices = Vector{T}(undef, nb_self_loops) + diagonal_nzind = Vector{T}(undef, nb_self_loops) if !augmented_graph(ag) + l = 0 for j in axes(S, 2) for k in nzrange(S, j) i = rv[k] if i == j - push!(diagonal_indices, i) - push!(diagonal_nzind, k) - ndiag += 1 + l += 1 + diagonal_indices[l] = i + diagonal_nzind[l] = k end end end end # Vectors for the decompression of the off-diagonal coefficients - nedges = (nnz(S) - ndiag) ÷ 2 + nedges = nb_edges(ag) lower_triangle_offsets = Vector{T}(undef, nedges) upper_triangle_offsets = Vector{T}(undef, nedges) diff --git a/test/matrices.jl b/test/matrices.jl index 1571497b..fa76946c 100644 --- a/test/matrices.jl +++ b/test/matrices.jl @@ -1,7 +1,12 @@ using LinearAlgebra using SparseArrays using SparseMatrixColorings: - check_same_pattern, matrix_versions, respectful_similar, same_pattern + BipartiteGraph, + AdjacencyGraph, + matrix_versions, + respectful_similar, + compatible_pattern, + check_compatible_pattern using StableRNGs using Test @@ -33,7 +38,7 @@ same_view(::Adjoint, ::Adjoint) = true end end -@testset "Sparsity pattern" begin +@testset "Compatible sparsity pattern -- BipartiteGraph" begin S = sparse([ 0 1 1 0 1 0 @@ -44,9 +49,33 @@ end A2 = copy(S) A2[1, 1] = 1 - @test same_pattern(A1, S) - @test !same_pattern(A2, S) - @test same_pattern(Matrix(A2), S) + bg1 = BipartiteGraph(A1) + bg2 = BipartiteGraph(A2) + @test compatible_pattern(S, bg1) + @test !compatible_pattern(S, bg2) + @test compatible_pattern(Matrix(S), bg1) - @test_throws DimensionMismatch check_same_pattern(A2, S) + @test_throws DimensionMismatch check_compatible_pattern(S, bg2) +end + +@testset "Compatible sparsity pattern -- AdjacencyGraph" begin + S = sparse([ + 1 0 1 + 0 1 1 + 1 1 0 + ]) + + A1 = copy(S) + A2 = copy(S) + A2[3, 3] = 1 + + ag1 = AdjacencyGraph(A1) + ag2 = AdjacencyGraph(A2) + for (op, uplo) in ((tril, :L), (triu, :U), (identity, :F)) + @test compatible_pattern(op(S), ag1, uplo) + @test !compatible_pattern(op(S), ag2, uplo) + @test compatible_pattern(Matrix(op(S)), ag1, uplo) + + @test_throws DimensionMismatch check_compatible_pattern(op(S), ag2, uplo) + end end From a314d12cc02e37962e914045ac808d8e33f6519a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 9 Jan 2026 18:23:06 +0100 Subject: [PATCH 24/39] Remove unused arguments from internal functions (#295) --- src/coloring.jl | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/coloring.jl b/src/coloring.jl index cf3d8f79..d4b80d99 100644 --- a/src/coloring.jl +++ b/src/coloring.jl @@ -302,7 +302,6 @@ function acyclic_coloring( _prevent_cycle!( v, w, - x, index_wx, color_x, first_visit_to_tree, @@ -333,7 +332,7 @@ function acyclic_coloring( color_x = color[x] (x == v || iszero(color_x)) && continue if color_x == color[v] - _merge_trees!(v, w, x, index_vw, index_wx, forest) # merge trees T₁ ∋ vw and T₂ ∋ wx if T₁ != T₂ + _merge_trees!(index_vw, index_wx, forest) # merge trees T₁ ∋ vw and T₂ ∋ wx if T₁ != T₂ end end end @@ -354,7 +353,6 @@ function _prevent_cycle!( # not modified v::Integer, w::Integer, - x::Integer, index_wx::Integer, color_x::Integer, # modified @@ -396,9 +394,6 @@ end function _merge_trees!( # not modified - v::Integer, - w::Integer, - x::Integer, index_vw::Integer, index_wx::Integer, # modified From 3b9fe4bbffc7e0c17a6042445eaaa03f466e0455 Mon Sep 17 00:00:00 2001 From: Alexis Montoison <35051714+amontoison@users.noreply.github.com> Date: Sun, 18 Jan 2026 21:03:03 -0600 Subject: [PATCH 25/39] Add test functions substitutable_columns and substitutable_bidirectional (#297) --- docs/src/dev.md | 3 + src/check.jl | 250 +++++++++++++++++++++++++++++++++++++++++++++++- test/check.jl | 124 ++++++++++++++++++++++++ test/utils.jl | 24 ++++- 4 files changed, 397 insertions(+), 4 deletions(-) diff --git a/docs/src/dev.md b/docs/src/dev.md index 82579dd6..a9e0d303 100644 --- a/docs/src/dev.md +++ b/docs/src/dev.md @@ -50,6 +50,9 @@ SparseMatrixColorings.directly_recoverable_columns SparseMatrixColorings.symmetrically_orthogonal_columns SparseMatrixColorings.structurally_orthogonal_columns SparseMatrixColorings.structurally_biorthogonal +SparseMatrixColorings.substitutable_columns +SparseMatrixColorings.substitutable_bidirectional +SparseMatrixColorings.rank_nonzeros_from_trees SparseMatrixColorings.valid_dynamic_order ``` diff --git a/src/check.jl b/src/check.jl index f47f0958..db5c0a34 100644 --- a/src/check.jl +++ b/src/check.jl @@ -86,11 +86,15 @@ A partition of the columns of a symmetric matrix `A` is _symmetrically orthogona 1. the group containing the column `A[:, j]` has no other column with a nonzero in row `i` 2. the group containing the column `A[:, i]` has no other column with a nonzero in row `j` +It is equivalent to a __star coloring__. + !!! warning This function is not coded with efficiency in mind, it is designed for small-scale tests. # References +> [_On the Estimation of Sparse Hessian Matrices_](https://doi.org/10.1137/0716078), Powell and Toint (1979) +> [_Estimation of sparse hessian matrices and graph coloring problems_](https://doi.org/10.1007/BF02612334), Coleman and Moré (1984) > [_What Color Is Your Jacobian? Graph Coloring for Computing Derivatives_](https://epubs.siam.org/doi/10.1137/S0036144504444711), Gebremedhin et al. (2005) """ function symmetrically_orthogonal_columns( @@ -102,7 +106,7 @@ function symmetrically_orthogonal_columns( end issymmetric(A) || return false group = group_by_color(color) - for i in axes(A, 2), j in axes(A, 2) + for i in axes(A, 1), j in axes(A, 2) iszero(A[i, j]) && continue ci, cj = color[i], color[j] check = _bilateral_check( @@ -126,6 +130,8 @@ A bipartition of the rows and columns of a matrix `A` is _structurally biorthogo 1. the group containing the column `A[:, j]` has no other column with a nonzero in row `i` 2. the group containing the row `A[i, :]` has no other row with a nonzero in column `j` +It is equivalent to a __star bicoloring__. + !!! warning This function is not coded with efficiency in mind, it is designed for small-scale tests. """ @@ -138,8 +144,8 @@ function structurally_biorthogonal( if !proper_length_bicoloring(A, row_color, column_color; verbose) return false end - column_group = group_by_color(column_color) row_group = group_by_color(row_color) + column_group = group_by_color(column_color) for i in axes(A, 1), j in axes(A, 2) iszero(A[i, j]) && continue ci, cj = row_color[i], column_color[j] @@ -261,6 +267,174 @@ function directly_recoverable_columns( return true end +""" + substitutable_columns( + A::AbstractMatrix, rank_nonzeros::AbstractMatrix, color::AbstractVector{<:Integer}; + verbose=false + ) + +Return `true` if coloring the columns of the symmetric matrix `A` with the vector `color` results in a partition that is substitutable, and `false` otherwise. +For all nonzeros `A[i, j]`, `rank_nonzeros[i, j]` provides its rank of recovery. + +A partition of the columns of a symmetric matrix `A` is _substitutable_ if, for every nonzero element `A[i, j]`, either of the following statements holds: + +1. the group containing the column `A[:, j]` has all nonzeros in row `i` ordered before `A[i, j]` +2. the group containing the column `A[:, i]` has all nonzeros in row `j` ordered before `A[i, j]` + +It is equivalent to an __acyclic coloring__. + +!!! warning + This function is not coded with efficiency in mind, it is designed for small-scale tests. + +# References + +> [_On the Estimation of Sparse Hessian Matrices_](https://doi.org/10.1137/0716078), Powell and Toint (1979) +> [_The Cyclic Coloring Problem and Estimation of Sparse Hessian Matrices_](https://doi.org/10.1137/0607026), Coleman and Cai (1986) +> [_What Color Is Your Jacobian? Graph Coloring for Computing Derivatives_](https://epubs.siam.org/doi/10.1137/S0036144504444711), Gebremedhin et al. (2005) +""" +function substitutable_columns( + A::AbstractMatrix, + rank_nonzeros::AbstractMatrix, + color::AbstractVector{<:Integer}; + verbose::Bool=false, +) + checksquare(A) + if !proper_length_coloring(A, color; verbose) + return false + end + issymmetric(A) || return false + group = group_by_color(color) + for i in axes(A, 1), j in axes(A, 2) + iszero(A[i, j]) && continue + ci, cj = color[i], color[j] + check = _substitutable_check( + A, rank_nonzeros; i, j, ci, cj, row_group=group, column_group=group, verbose + ) + !check && return false + end + return true +end + +""" + substitutable_bidirectional( + A::AbstractMatrix, rank_nonzeros::AbstractMatrix, row_color::AbstractVector{<:Integer}, column_color::AbstractVector{<:Integer}; + verbose=false + ) + +Return `true` if bicoloring of the matrix `A` with the vectors `row_color` and `column_color` results in a bipartition that is substitutable, and `false` otherwise. +For all nonzeros `A[i, j]`, `rank_nonzeros[i, j]` provides its rank of recovery. + +A bipartition of the rows and columns of a matrix `A` is _substitutable_ if, for every nonzero element `A[i, j]`, either of the following statements holds: + +1. the group containing the column `A[:, j]` has all nonzeros in row `i` ordered before `A[i, j]` +2. the group containing the row `A[i, :]` has all nonzeros in column `j` ordered before `A[i, j]` + +It is equivalent to an __acyclic bicoloring__. + +!!! warning + This function is not coded with efficiency in mind, it is designed for small-scale tests. +""" +function substitutable_bidirectional( + A::AbstractMatrix, + rank_nonzeros::AbstractMatrix, + row_color::AbstractVector{<:Integer}, + column_color::AbstractVector{<:Integer}; + verbose::Bool=false, +) + if !proper_length_bicoloring(A, row_color, column_color; verbose) + return false + end + row_group = group_by_color(row_color) + column_group = group_by_color(column_color) + for i in axes(A, 1), j in axes(A, 2) + iszero(A[i, j]) && continue + ci, cj = row_color[i], column_color[j] + check = _substitutable_check( + A, rank_nonzeros; i, j, ci, cj, row_group, column_group, verbose + ) + !check && return false + end + return true +end + +function _substitutable_check( + A::AbstractMatrix, + rank_nonzeros::AbstractMatrix; + i::Integer, + j::Integer, + ci::Integer, + cj::Integer, + row_group::AbstractVector, + column_group::AbstractVector, + verbose::Bool, +) + order_ij = rank_nonzeros[i, j] + k_row = 0 + k_column = 0 + if ci != 0 + for k in row_group[ci] + (k == i) && continue + if !iszero(A[k, j]) + order_kj = rank_nonzeros[k, j] + @assert !iszero(order_kj) + if order_kj > order_ij + k_row = k + end + end + end + end + if cj != 0 + for k in column_group[cj] + (k == j) && continue + if !iszero(A[i, k]) + order_ik = rank_nonzeros[i, k] + @assert !iszero(order_ik) + if order_ik > order_ij + k_column = k + end + end + end + end + if ci == 0 && cj == 0 + if verbose + @warn """ + For coefficient (i=$i, j=$j) with colors (ci=$ci, cj=$cj): + - Row color ci=$ci is neutral. + - Column color cj=$cj is neutral. + """ + end + return false + elseif ci == 0 && !iszero(k_column) + if verbose + @warn """ + For coefficient (i=$i, j=$j) with colors (ci=$ci, cj=$cj): + - Row color ci=$ci is neutral. + - For the column $k_column in column color cj=$cj, A[$i, $k_column] is ordered after A[$i, $j]. + """ + end + return false + elseif cj == 0 && !iszero(k_row) + if verbose + @warn """ + For coefficient (i=$i, j=$j) with colors (ci=$ci, cj=$cj): + - For the row $k_row in row color ci=$ci, A[$k_row, $j] is ordered after A[$i, $j]. + - Column color cj=$cj is neutral. + """ + end + return false + elseif !iszero(k_row) && !iszero(k_column) + if verbose + @warn """ + For coefficient (i=$i, j=$j) with colors (ci=$ci, cj=$cj): + - For the row $k_row in row color ci=$ci, A[$k_row, $j] is ordered after A[$i, $j]. + - For the column $k_column in column color cj=$cj, A[$i, $k_column] is ordered after A[$i, $j]. + """ + end + return false + end + return true +end + """ valid_dynamic_order(g::AdjacencyGraph, π::AbstractVector{<:Integer}, order::DynamicDegreeBasedOrder) valid_dynamic_order(bg::AdjacencyGraph, ::Val{side}, π::AbstractVector{<:Integer}, order::DynamicDegreeBasedOrder) @@ -326,3 +500,75 @@ function valid_dynamic_order( end return true end + +""" + rank_nonzeros_from_trees(result::TreeSetColoringResult) + rank_nonzeros_from_trees(result::BicoloringResult) + +Construct a sparse matrix `rank_nonzeros` that assigns a unique recovery rank +to each nonzero coefficient associated with an acyclic coloring or bicoloring. + +For every nonzero entry `result.A[i, j]`, `rank_nonzeros[i, j]` stores a strictly positive +integer representing the order in which this coefficient is recovered during the decompression. +A larger value means the coefficient is recovered later. + +This ranking is used to test substitutability (acyclicity) of colorings: +for a given nonzero `result.A[i, j]`, the ranks allow one to check whether all competing +nonzeros in the same row or column (within a color group) are recovered before it. +""" +function rank_nonzeros_from_trees end + +function rank_nonzeros_from_trees(result::TreeSetColoringResult) + (; A, ag, reverse_bfs_orders, diagonal_indices, tree_edge_indices, nt) = result + (; S) = ag + m, n = size(A) + nnzS = nnz(S) + nzval = zeros(Int, nnzS) + rank_nonzeros = SparseMatrixCSC(n, n, S.colptr, S.rowval, nzval) + counter = 0 + for i in diagonal_indices + counter += 1 + rank_nonzeros[i, i] = counter + end + for k in 1:nt + first = tree_edge_indices[k] + last = tree_edge_indices[k + 1] - 1 + for pos in first:last + (i, j) = reverse_bfs_orders[pos] + counter += 1 + rank_nonzeros[i, j] = counter + rank_nonzeros[j, i] = counter + end + end + return rank_nonzeros +end + +function rank_nonzeros_from_trees(result::BicoloringResult) + (; A, abg, row_color, column_color, symmetric_result, large_colptr, large_rowval) = + result + @assert symmetric_result isa TreeSetColoringResult + (; ag, reverse_bfs_orders, tree_edge_indices, nt) = symmetric_result + (; S) = ag + m, n = size(A) + nnzA = nnz(S) ÷ 2 + nzval = zeros(Int, nnzA) + colptr = large_colptr[1:(n + 1)] + rowval = large_rowval[1:nnzA] + rowval .-= n + rank_nonzeros = SparseMatrixCSC(m, n, colptr, rowval, nzval) + counter = 0 + for k in 1:nt + first = tree_edge_indices[k] + last = tree_edge_indices[k + 1] - 1 + for pos in first:last + (i, j) = reverse_bfs_orders[pos] + counter += 1 + if i > j + rank_nonzeros[i - n, j] = counter + else + rank_nonzeros[j - n, i] = counter + end + end + end + return rank_nonzeros +end diff --git a/test/check.jl b/test/check.jl index 4e9e5a61..b2701095 100644 --- a/test/check.jl +++ b/test/check.jl @@ -4,6 +4,8 @@ using SparseMatrixColorings: symmetrically_orthogonal_columns, structurally_biorthogonal, directly_recoverable_columns, + substitutable_columns, + substitutable_bidirectional, what_fig_41, efficient_fig_1 using Test @@ -184,3 +186,125 @@ For coefficient (i=1, j=1) with colors (ci=0, cj=0): """, ) structurally_biorthogonal(A, [0, 2, 2, 3], [0, 2, 2, 2, 3], verbose=true) end + +@testset "Substitutable columns" begin + A1 = [ + 1 1 1 1 1 + 1 1 0 0 0 + 1 0 1 0 0 + 1 0 0 1 0 + 1 0 0 0 1 + ] + B1 = [ + 1 6 7 8 9 + 6 2 0 0 0 + 7 0 3 0 0 + 8 0 0 4 0 + 9 0 0 0 5 + ] + A2 = [ + 1 1 0 0 0 + 1 1 1 0 0 + 0 1 1 1 0 + 0 0 1 1 1 + 0 0 0 1 1 + ] + B2 = [ + 5 1 0 0 0 + 1 6 2 0 0 + 0 2 7 3 0 + 0 0 3 8 4 + 0 0 0 4 9 + ] + A3 = [ + 0 1 1 1 1 + 1 0 1 1 1 + 1 1 0 1 1 + 1 1 1 0 1 + 1 1 1 1 0 + ] + B3 = [ + 0 1 2 3 4 + 1 0 5 6 7 + 2 5 0 8 9 + 3 6 8 0 10 + 4 7 9 10 0 + ] + + # success + + @test substitutable_columns(A1, B1, [1, 2, 2, 2, 2]) + @test substitutable_columns(A2, B2, [1, 2, 3, 1, 2]) + @test substitutable_columns(A3, B3, [1, 2, 3, 4, 0]) + + # failure + + @test !substitutable_columns(A1, B1, [1, 1, 1, 1]) + log = (:warn, "4 colors provided for 5 columns.") + @test_logs log substitutable_columns(A1, B1, [1, 1, 1, 1]; verbose=true) + + @test !substitutable_columns(A1, B1, [1, 1, 1, 1, 1]) + @test_logs ( + :warn, + """ +For coefficient (i=1, j=1) with colors (ci=1, cj=1): +- For the row 5 in row color ci=1, A[5, 1] is ordered after A[1, 1]. +- For the column 5 in column color cj=1, A[1, 5] is ordered after A[1, 1]. +""", + ) substitutable_columns(A1, B1, [1, 1, 1, 1, 1]; verbose=true) + + @test !substitutable_columns(A2, B2, [1, 2, 0, 1, 2]) + @test_logs ( + :warn, + """ +For coefficient (i=3, j=3) with colors (ci=0, cj=0): +- Row color ci=0 is neutral. +- Column color cj=0 is neutral. +""", + ) substitutable_columns(A2, B2, [1, 2, 0, 1, 2]; verbose=true) + + @test !substitutable_columns(A3, B3, [0, 1, 2, 3, 3]) + @test_logs ( + :warn, + """ +For coefficient (i=1, j=4) with colors (ci=0, cj=3): +- Row color ci=0 is neutral. +- For the column 5 in column color cj=3, A[1, 5] is ordered after A[1, 4]. +""", + ) substitutable_columns(A3, B3, [0, 1, 2, 3, 3]; verbose=true) + + @test !substitutable_columns(A3, B3, [1, 2, 3, 3, 0]) + @test_logs ( + :warn, + """ +For coefficient (i=3, j=5) with colors (ci=3, cj=0): +- For the row 4 in row color ci=3, A[4, 5] is ordered after A[3, 5]. +- Column color cj=0 is neutral. +""", + ) substitutable_columns(A3, B3, [1, 2, 3, 3, 0]; verbose=true) +end + +@testset "Substitutable bidirectional" begin + A = [ + 1 0 0 + 0 1 0 + 0 0 1 + ] + B = [ + 1 0 0 + 0 2 0 + 0 0 3 + ] + + # success + + @test substitutable_bidirectional(A, B, [1, 0, 0], [0, 1, 1]) + + # failure + + log = (:warn, "2 colors provided for 3 columns.") + @test_logs log !substitutable_bidirectional(A, B, [1, 0, 0], [0, 1]; verbose=true) + + log = (:warn, "4 colors provided for 3 rows.") + @test_logs log !substitutable_bidirectional(A, B, [1, 0, 0, 1], [0, 1, 1]; verbose=true) +end diff --git a/test/utils.jl b/test/utils.jl index 2462f8c1..77c88610 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -12,7 +12,10 @@ using SparseMatrixColorings: respectful_similar, structurally_orthogonal_columns, symmetrically_orthogonal_columns, - structurally_biorthogonal + structurally_biorthogonal, + substitutable_columns, + substitutable_bidirectional, + rank_nonzeros_from_trees using Test const _ALL_ORDERS = ( @@ -77,7 +80,6 @@ function test_coloring_decompression( end @testset "Recoverability" begin - # TODO: find tests for recoverability for substitution decompression if decompression == :direct if structure == :nonsymmetric if partition == :column @@ -97,6 +99,15 @@ function test_coloring_decompression( end end + @testset "Substitutable" begin + if decompression == :substitution + if structure == :symmetric + rank_nonzeros = rank_nonzeros_from_trees(result) + @test substitutable_columns(A0, rank_nonzeros, color) + end + end + end + @testset "Single-color decompression" begin if decompression == :direct # TODO: implement for :substitution too A2 = respectful_similar(A, eltype(B)) @@ -235,6 +246,15 @@ function test_bicoloring_decompression( @test structurally_biorthogonal(A0, row_color, column_color) end end + + if decompression == :substitution + @testset "Substitutable" begin + rank_nonzeros = rank_nonzeros_from_trees(result) + @test substitutable_bidirectional( + A0, rank_nonzeros, row_color, column_color + ) + end + end end @testset "More orders is better" begin From 921c7a06f467d005187276314d7d4adb0e1ef1bd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 29 Jan 2026 09:49:03 +0100 Subject: [PATCH 26/39] Document fields of TreeSet (#296) --- src/coloring.jl | 25 ++++++++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/src/coloring.jl b/src/coloring.jl index d4b80d99..699dbc64 100644 --- a/src/coloring.jl +++ b/src/coloring.jl @@ -417,9 +417,23 @@ Encode a set of 2-colored trees resulting from the [`acyclic_coloring`](@ref) al $TYPEDFIELDS """ struct TreeSet{T} + """ + For a tree index `1 <= k <= nt`, the list + `list = reverse_bfs_order[tree_edge_indices[k]:(tree_edge_indices[k+1]-1)]` is a list of edges + `list[i] = (leaf, inner)` of the `k`th tree such that `leaf` is a leaf of the tree containing + the edges `list[i:end]`. + From an other point of view, `reverse(list)` contains the edges in the order of a breadth first + (BFS) traversal order of the `k`th tree starting from a depth-minimizing root. + """ reverse_bfs_orders::Vector{Tuple{T,T}} + "For a tree index `1 <= k <= nt`, `is_star[k]` indicates whether the `k`th three is a star." is_star::Vector{Bool} + """ + `tree_edge_indices[1]` is one and `tree_edge_indices[k+1] - tree_edge_indices[k]` is the number of edges in the `k`th tree. + One can think of it as a kind of fused vector of offsets (similar to the `colptr` field of `SparseMatrixCSC`) of all trees together. + """ tree_edge_indices::Vector{T} + "numbers of 2-colored trees for which trees sharing the same 2 colors have disjoint edges" nt::T end @@ -427,6 +441,7 @@ function TreeSet( g::AdjacencyGraph{T}, forest::Forest{T}, buffer::AbstractVector{T}, + # The value of `reverse_bfs_orders` is ignored, we just provide the storage for it (or reuse memory allocated during acyclic coloring) reverse_bfs_orders::Vector{Tuple{T,T}}, ne::Integer, ) where {T} @@ -561,7 +576,15 @@ function TreeSet( # Number of edges treated num_edges_treated = zero(T) - # reverse_bfs_orders contains the reverse breadth first (BFS) traversal order for each tree in the forest + # The `rank` of the `k`th tree encoded in `forest` does not correspond + # to the depth of the tree rooted as the root encoded in `forest` because + # `forest.parents[u] = v` only needs a path to exists from `u` to `v` but + # there may not be an edge `(u, v)`. + # We also want a root `r` that minimizes the depth of the tree rooted at + # `r`. To achieve this, we start from each leaf and remove the corresponding + # edges. We then look at all leaves of the corresponding graphs and repeat + # the process until there is only one vertex left. This vertex will then be + # a depth-minimizing root. for k in 1:nt # Initialize the queue to store the leaves queue_start = 1 From a7cb1547e240a68afcfead3d8e5c88441c587128 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 29 Jan 2026 18:36:23 +0100 Subject: [PATCH 27/39] Fix respectful_similar with SparsityPatternCSC (#299) * Fix respectful_similar with SparsityPatternCSC * Add test * Fix format * Add test --------- Co-authored-by: Alexis Montoison <35051714+amontoison@users.noreply.github.com> --- src/graph.jl | 5 +++++ test/structured.jl | 17 +++++++++++++++++ 2 files changed, 22 insertions(+) diff --git a/src/graph.jl b/src/graph.jl index b93a864b..39df611c 100644 --- a/src/graph.jl +++ b/src/graph.jl @@ -32,6 +32,11 @@ SparseArrays.nnz(S::SparsityPatternCSC) = length(S.rowval) SparseArrays.rowvals(S::SparsityPatternCSC) = S.rowval SparseArrays.nzrange(S::SparsityPatternCSC, j::Integer) = S.colptr[j]:(S.colptr[j + 1] - 1) +# Needed if using `coloring(::SparsityPatternCSC, ...)` +function Base.similar(A::SparsityPatternCSC, ::Type{T}) where {T} + return SparseArrays.SparseMatrixCSC(A.m, A.n, A.colptr, A.rowval, similar(A.rowval, T)) +end + """ transpose(S::SparsityPatternCSC) diff --git a/test/structured.jl b/test/structured.jl index 23ba967f..5c7b5874 100644 --- a/test/structured.jl +++ b/test/structured.jl @@ -2,6 +2,7 @@ using ArrayInterface: ArrayInterface using BandedMatrices: BandedMatrix, brand using BlockBandedMatrices: BandedBlockBandedMatrix, BlockBandedMatrix using LinearAlgebra +using SparseArrays using SparseMatrixColorings using Test @@ -56,3 +57,19 @@ end; test_structured_coloring_decompression(A) end end; + +# See https://github.com/gdalle/SparseMatrixColorings.jl/pull/299 +@testset "SparsityPatternCSC $T" for T in [Int, Float32] + S = sparse(T[ + 0 0 1 1 0 1 + 1 0 0 0 1 0 + 0 1 0 0 1 0 + 0 1 1 0 0 0 + ]) + P = SparseMatrixColorings.SparsityPatternCSC(S) + problem = ColoringProblem() + algo = GreedyColoringAlgorithm() + result = coloring(P, problem, algo) + B = compress(S, result) + @test decompress(B, result) isa SparseMatrixCSC{T,Int} +end; From 47899442d0fcb3a0b4b0b26d21b0f7fd15b44950 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 18 Feb 2026 12:00:33 +0100 Subject: [PATCH 28/39] Fix coloring with empty matrix as input (#300) * Fix coloring with empty matrix as input * Update src/result.jl Co-authored-by: Guillaume Dalle <22795598+gdalle@users.noreply.github.com> * Add check * Fix format * Change check * Apply suggestion from @gdalle --------- Co-authored-by: Guillaume Dalle <22795598+gdalle@users.noreply.github.com> --- src/interface.jl | 8 +++++++- src/result.jl | 3 +++ test/check.jl | 11 +++++++++++ 3 files changed, 21 insertions(+), 1 deletion(-) diff --git a/src/interface.jl b/src/interface.jl index 936b974c..0d183c9c 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -305,7 +305,13 @@ function _coloring( vertices_in_order = vertices(ag, order) return acyclic_coloring(ag, vertices_in_order, algo.postprocessing) end - color, tree_set = argmin(maximum ∘ first, color_and_tree_set_by_order) + # if `color` is empty, `maximum` will fail but `color_and_tree_set_by_order` + # is also one so we can just add a special case for this + if length(color_and_tree_set_by_order) == 1 + color, tree_set = only(color_and_tree_set_by_order) + else + color, tree_set = argmin(maximum ∘ first, color_and_tree_set_by_order) + end if speed_setting isa WithResult return TreeSetColoringResult(A, ag, color, tree_set, R) else diff --git a/src/result.jl b/src/result.jl index c74ceaa1..0d47cf13 100644 --- a/src/result.jl +++ b/src/result.jl @@ -82,6 +82,9 @@ Create a color-indexed vector `group` such that `i ∈ group[c]` iff `color[i] = Assumes the colors are contiguously numbered from `0` to some `cmax`. """ function group_by_color(::Type{T}, color::AbstractVector) where {T<:Integer} + if isempty(color) + return typeof(view(T[], 1:0))[] + end cmin, cmax = extrema(color) @assert cmin >= 0 # Compute group sizes and offsets for a joint storage diff --git a/test/check.jl b/test/check.jl index b2701095..246261d0 100644 --- a/test/check.jl +++ b/test/check.jl @@ -1,4 +1,5 @@ using LinearAlgebra +using SparseArrays using SparseMatrixColorings: structurally_orthogonal_columns, symmetrically_orthogonal_columns, @@ -308,3 +309,13 @@ end log = (:warn, "4 colors provided for 3 rows.") @test_logs log !substitutable_bidirectional(A, B, [1, 0, 0, 1], [0, 1, 1]; verbose=true) end + +# See https://github.com/gdalle/SparseMatrixColorings.jl/pull/300 +@testset "Empty matrix" begin + problem = ColoringProblem(; structure=:symmetric, partition=:column) + algo = GreedyColoringAlgorithm(; decompression=:substitution) + S = spzeros(Int, 0, 0) + result = coloring(S, problem, algo) + @test isempty(result.color) + @test isempty(result.group) +end From 1c3c9da10304047708695e46d03072fde765ae0e Mon Sep 17 00:00:00 2001 From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com> Date: Fri, 27 Feb 2026 09:17:27 +0100 Subject: [PATCH 29/39] test: skip failing MiniZinc tests (#305) --- test/optimal.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/optimal.jl b/test/optimal.jl index 10fc32ef..36bf91c6 100644 --- a/test/optimal.jl +++ b/test/optimal.jl @@ -15,13 +15,14 @@ asymmetric_params = vcat( algo = GreedyColoringAlgorithm() optalgo = OptimalColoringAlgorithm(() -> MiniZinc.Optimizer{Float64}("highs"); silent=false) +# TODO: reactivate tests once https://github.com/jump-dev/MiniZinc.jl/issues/103 is fixed + @testset "Column coloring" begin problem = ColoringProblem(; structure=:nonsymmetric, partition=:column) for (m, n, p) in asymmetric_params A = sprand(rng, m, n, p) result = coloring(A, problem, algo) - optresult = coloring(A, problem, optalgo) - @test ncolors(result) >= ncolors(optresult) + @test_skip ncolors(result) >= ncolors(coloring(A, problem, optalgo)) end end @@ -30,8 +31,7 @@ end for (m, n, p) in asymmetric_params A = sprand(rng, m, n, p) result = coloring(A, problem, algo) - optresult = coloring(A, problem, optalgo) - @test ncolors(result) >= ncolors(optresult) + @test_skip ncolors(result) >= ncolors(coloring(A, problem, optalgo)) end end From 408d8f34ead90522dd277daf47818e165ea5d9ea Mon Sep 17 00:00:00 2001 From: Chris Rackauckas - Beep Boop Edition Date: Fri, 27 Feb 2026 00:55:27 -0800 Subject: [PATCH 30/39] Fix eltype invalidation for SparsityPatternCSC (#304) SparsityPatternCSC{Ti} <: AbstractMatrix{Bool} but the custom Base.eltype method was returning Ti (the index type) instead of Bool. This caused 24,906 method invalidations when loading the package, as it invalidated the backedge from Base.eltype(::AbstractArray). Change the custom method to SparseArrays.indtype instead, since that's what the type parameter Ti actually represents. The eltype is now correctly inherited from AbstractMatrix{Bool}. Co-authored-by: ChrisRackauckas-Claude Co-authored-by: Guillaume Dalle <22795598+gdalle@users.noreply.github.com> --- src/graph.jl | 2 +- test/graph.jl | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/graph.jl b/src/graph.jl index 39df611c..94954602 100644 --- a/src/graph.jl +++ b/src/graph.jl @@ -23,7 +23,7 @@ end SparsityPatternCSC(A::SparseMatrixCSC) = SparsityPatternCSC(A.m, A.n, A.colptr, A.rowval) -Base.eltype(::SparsityPatternCSC{T}) where {T} = T +SparseArrays.indtype(::SparsityPatternCSC{T}) where {T} = T Base.size(S::SparsityPatternCSC) = (S.m, S.n) Base.size(S::SparsityPatternCSC, d::Integer) = d::Integer <= 2 ? size(S)[d] : 1 Base.axes(S::SparsityPatternCSC, d::Integer) = Base.OneTo(size(S, d)) diff --git a/test/graph.jl b/test/graph.jl index e9784e56..773a0206 100644 --- a/test/graph.jl +++ b/test/graph.jl @@ -15,7 +15,8 @@ using Test ## SparsityPatternCSC @testset "SparsityPatternCSC" begin - @test eltype(SparsityPatternCSC(sprand(10, 10, 0.1))) == Int + @test eltype(SparsityPatternCSC(sprand(10, 10, 0.1))) == Bool + @test SparseArrays.indtype(SparsityPatternCSC(sprand(10, 10, 0.1))) == Int @testset "Transpose" begin for _ in 1:1000 m, n = rand(100:1000), rand(100:1000) From 099603049103f40918eb607351d85d3a721f4af0 Mon Sep 17 00:00:00 2001 From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com> Date: Fri, 27 Feb 2026 09:58:42 +0100 Subject: [PATCH 31/39] chore: bump version (#306) --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 0e01b3d4..3e20fc96 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SparseMatrixColorings" uuid = "0a514795-09f3-496d-8182-132a7b665d35" authors = ["Guillaume Dalle", "Alexis Montoison"] -version = "0.4.23" +version = "0.4.24" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From 0d3501aff3bc3e0d27bd18f4b070d77d73baaf7e Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 9 Mar 2026 15:11:09 +0100 Subject: [PATCH 32/39] Bump julia-actions/cache from 2 to 3 (#307) --- .github/workflows/Benchmark.yml | 2 +- .github/workflows/Documentation.yml | 2 +- .github/workflows/Test-GPU.yml | 2 +- .github/workflows/Test.yml | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/.github/workflows/Benchmark.yml b/.github/workflows/Benchmark.yml index 30e2e1e5..fec67016 100644 --- a/.github/workflows/Benchmark.yml +++ b/.github/workflows/Benchmark.yml @@ -16,7 +16,7 @@ jobs: - uses: julia-actions/setup-julia@v2 with: version: "1" - - uses: julia-actions/cache@v2 + - uses: julia-actions/cache@v3 - name: Extract Package Name from Project.toml id: extract-package-name run: | diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml index b73255ea..4d85c43d 100644 --- a/.github/workflows/Documentation.yml +++ b/.github/workflows/Documentation.yml @@ -23,7 +23,7 @@ jobs: - uses: julia-actions/setup-julia@v2 with: version: '1' - - uses: julia-actions/cache@v2 + - uses: julia-actions/cache@v3 - name: Install dependencies run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' - name: Build and deploy diff --git a/.github/workflows/Test-GPU.yml b/.github/workflows/Test-GPU.yml index aa30a5fc..cc390241 100644 --- a/.github/workflows/Test-GPU.yml +++ b/.github/workflows/Test-GPU.yml @@ -35,7 +35,7 @@ jobs: if: ${{ matrix.version == 'lts' }} with: skip: LinearAlgebra, Random, SparseArrays - - uses: julia-actions/cache@v2 + - uses: julia-actions/cache@v3 - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 diff --git a/.github/workflows/Test.yml b/.github/workflows/Test.yml index d979b583..46601d04 100644 --- a/.github/workflows/Test.yml +++ b/.github/workflows/Test.yml @@ -32,7 +32,7 @@ jobs: if: ${{ matrix.version == '1.10' }} with: skip: LinearAlgebra, Random, SparseArrays - - uses: julia-actions/cache@v2 + - uses: julia-actions/cache@v3 - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 From 6948b0afdcba7f57443e73a47ecf9a6519f6db71 Mon Sep 17 00:00:00 2001 From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com> Date: Mon, 16 Mar 2026 13:05:51 +0100 Subject: [PATCH 33/39] Update repo owner to JuliaDiff org (#309) * Update repo owner * Bump version * Switch GPU CI to buildkite * Add badge --- .buildkite/pipeline.yml | 18 +++++++++++++ .github/workflows/Test-GPU.yml | 46 ---------------------------------- CITATION.cff | 2 +- Project.toml | 2 +- README.md | 10 +++++--- docs/make.jl | 4 ++- test/check.jl | 2 +- test/structured.jl | 2 +- 8 files changed, 31 insertions(+), 55 deletions(-) create mode 100644 .buildkite/pipeline.yml delete mode 100644 .github/workflows/Test-GPU.yml diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml new file mode 100644 index 00000000..a1a04cb0 --- /dev/null +++ b/.buildkite/pipeline.yml @@ -0,0 +1,18 @@ +env: + SECRET_CODECOV_TOKEN: "l7PxpMgxHEiz2l3UDDj5122Xirv5du1QFQiTf09YISC8VOz/nXeQC2yYY6iTMbccQG3xrumfWke9vtLbAc9LPH75kw3zQODuWMbz0eGSf9wtw7foqO6KAo5neMTeZJ0ZFgFMa2Q89T5SNIQsi4b5Zs7BVto6qB9Z/3QEs/BpsR25cYkY4Y6JBU6XuqZ6GRAc6BtlB4OmBEp3BBsXatx6y64zF0qbp4rocmcBQEoeQkcXtxf0dfA0KNAHpweWPrzIAMZ0aUYp6iEikxwLY5TjVFhOIcpVXUlxIdhl2qaaDF6b6WlgGXiGLpBjsLQfvlOgIlXH59Ddg0IVxboF2a37OA==;U2FsdGVkX19MveWMe1aasYZoJwaeiKO5XqwMZ/utOFj7h1CE9UX4nduBMXTyJ77tpNmCBYsQNJ/PadbX2224hQ==" + +steps: + - label: "Julia v1" + plugins: + - JuliaCI/julia#v1: + version: "1" + - JuliaCI/julia-test#v1: ~ + - JuliaCI/julia-coverage#v1: + codecov: true + agents: + queue: "juliagpu" + cuda: "*" + if: build.message !~ /\[skip tests\]/ + timeout_in_minutes: 60 + env: + JULIA_SMC_TEST_GROUP: 'GPU' \ No newline at end of file diff --git a/.github/workflows/Test-GPU.yml b/.github/workflows/Test-GPU.yml deleted file mode 100644 index cc390241..00000000 --- a/.github/workflows/Test-GPU.yml +++ /dev/null @@ -1,46 +0,0 @@ -name: Test-GPU - -on: - push: - branches: - - main - pull_request: - -concurrency: - group: ${{ github.workflow }}-${{ github.ref }} - cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} - -# needed to allow julia-actions/cache to delete old caches that it has created -permissions: - actions: write - contents: read - -jobs: - test: - runs-on: self-hosted - env: - CUDA_VISIBLE_DEVICES: 1 - JULIA_DEPOT_PATH: /scratch/github-actions/julia_depot_alexis - JULIA_SMC_TEST_GROUP: 'GPU' - strategy: - matrix: - julia-version: ['lts', '1'] - steps: - - uses: actions/checkout@v6 - - uses: julia-actions/setup-julia@v2 - with: - version: ${{ matrix.julia-version }} - arch: x64 - - uses: julia-actions/julia-downgrade-compat@v2 - if: ${{ matrix.version == 'lts' }} - with: - skip: LinearAlgebra, Random, SparseArrays - - uses: julia-actions/cache@v3 - - uses: julia-actions/julia-buildpkg@v1 - - uses: julia-actions/julia-runtest@v1 - - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v5 - with: - files: lcov.info - token: ${{ secrets.CODECOV_TOKEN }} - fail_ci_if_error: false diff --git a/CITATION.cff b/CITATION.cff index c2c0a8fb..a36bc54e 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -18,7 +18,7 @@ identifiers: - type: doi value: 10.5281/zenodo.11314275 description: Zenodo -repository-code: 'https://github.com/gdalle/SparseMatrixColorings.jl' +repository-code: 'https://github.com/JuliaDiff/SparseMatrixColorings.jl' abstract: >- Coloring algorithms for sparse Jacobian and Hessian matrices diff --git a/Project.toml b/Project.toml index 3e20fc96..4fb07535 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SparseMatrixColorings" uuid = "0a514795-09f3-496d-8182-132a7b665d35" authors = ["Guillaume Dalle", "Alexis Montoison"] -version = "0.4.24" +version = "0.4.25" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" diff --git a/README.md b/README.md index 66b10dd3..1e376a70 100644 --- a/README.md +++ b/README.md @@ -1,9 +1,11 @@ # SparseMatrixColorings.jl -[![Build Status](https://github.com/gdalle/SparseMatrixColorings.jl/actions/workflows/Test.yml/badge.svg?branch=main)](https://github.com/gdalle/SparseMatrixColorings.jl/actions/workflows/Test.yml?query=branch%3Amain) -[![Stable Documentation](https://img.shields.io/badge/docs-stable-blue.svg)](https://gdalle.github.io/SparseMatrixColorings.jl/stable/) -[![Dev Documentation](https://img.shields.io/badge/docs-dev-blue.svg)](https://gdalle.github.io/SparseMatrixColorings.jl/dev/) -[![Coverage](https://codecov.io/gh/gdalle/SparseMatrixColorings.jl/branch/main/graph/badge.svg)](https://app.codecov.io/gh/gdalle/SparseMatrixColorings.jl) +[![Build Status](https://github.com/juliadiff/SparseMatrixColorings.jl/actions/workflows/Test.yml/badge.svg?branch=main)](https://github.com/juliadiff/SparseMatrixColorings.jl/actions/workflows/Test.yml?query=branch%3Amain) +[![GPU build status](https://badge.buildkite.com/7d8ed289d7bdb5a25ae48b2c778a202ce4990b7ee558cdfef8.svg)](https://buildkite.com/julialang/sparsematrixcolorings-dot-jl) +[![Coverage](https://codecov.io/gh/juliadiff/SparseMatrixColorings.jl/branch/main/graph/badge.svg)](https://app.codecov.io/gh/juliadiff/SparseMatrixColorings.jl) + +[![Stable Documentation](https://img.shields.io/badge/docs-stable-blue.svg)](https://juliadiff.org/SparseMatrixColorings.jl/stable/) +[![Dev Documentation](https://img.shields.io/badge/docs-dev-blue.svg)](https://juliadiff.org/SparseMatrixColorings.jl/dev/) [![Code Style: Blue](https://img.shields.io/badge/code%20style-blue-4495d1.svg)](https://github.com/JuliaDiff/BlueStyle) [![arXiv](https://img.shields.io/badge/arXiv-2505.07308-b31b1b.svg)](https://arxiv.org/abs/2505.07308) [![DOI](https://zenodo.org/badge/801999408.svg)](https://zenodo.org/doi/10.5281/zenodo.11314275) diff --git a/docs/make.jl b/docs/make.jl index 78fba268..80d2a23d 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -21,5 +21,7 @@ makedocs(; ) deploydocs(; - repo="github.com/gdalle/SparseMatrixColorings.jl", push_preview=true, devbranch="main" + repo="github.com/JuliaDiff/SparseMatrixColorings.jl", + push_preview=true, + devbranch="main", ) diff --git a/test/check.jl b/test/check.jl index 246261d0..8cf2c495 100644 --- a/test/check.jl +++ b/test/check.jl @@ -310,7 +310,7 @@ end @test_logs log !substitutable_bidirectional(A, B, [1, 0, 0, 1], [0, 1, 1]; verbose=true) end -# See https://github.com/gdalle/SparseMatrixColorings.jl/pull/300 +# See https://github.com/JuliaDiff/SparseMatrixColorings.jl/pull/300 @testset "Empty matrix" begin problem = ColoringProblem(; structure=:symmetric, partition=:column) algo = GreedyColoringAlgorithm(; decompression=:substitution) diff --git a/test/structured.jl b/test/structured.jl index 5c7b5874..0d2a7564 100644 --- a/test/structured.jl +++ b/test/structured.jl @@ -58,7 +58,7 @@ end; end end; -# See https://github.com/gdalle/SparseMatrixColorings.jl/pull/299 +# See https://github.com/JuliaDiff/SparseMatrixColorings.jl/pull/299 @testset "SparsityPatternCSC $T" for T in [Int, Float32] S = sparse(T[ 0 0 1 1 0 1 From f56da2909ddd2319b16953b395a7b049c33f5c33 Mon Sep 17 00:00:00 2001 From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com> Date: Mon, 16 Mar 2026 21:07:51 +0100 Subject: [PATCH 34/39] Fix buildkite badge [skip tests] (#310) --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 1e376a70..b37b0ca8 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,7 @@ # SparseMatrixColorings.jl [![Build Status](https://github.com/juliadiff/SparseMatrixColorings.jl/actions/workflows/Test.yml/badge.svg?branch=main)](https://github.com/juliadiff/SparseMatrixColorings.jl/actions/workflows/Test.yml?query=branch%3Amain) -[![GPU build status](https://badge.buildkite.com/7d8ed289d7bdb5a25ae48b2c778a202ce4990b7ee558cdfef8.svg)](https://buildkite.com/julialang/sparsematrixcolorings-dot-jl) +[![GPU build status](https://badge.buildkite.com/7d8ed289d7bdb5a25ae48b2c778a202ce4990b7ee558cdfef8.svg?branch=main)](https://buildkite.com/julialang/sparsematrixcolorings-dot-jl) [![Coverage](https://codecov.io/gh/juliadiff/SparseMatrixColorings.jl/branch/main/graph/badge.svg)](https://app.codecov.io/gh/juliadiff/SparseMatrixColorings.jl) [![Stable Documentation](https://img.shields.io/badge/docs-stable-blue.svg)](https://juliadiff.org/SparseMatrixColorings.jl/stable/) From 9aebbe6d5b2299c8af2b14bb9f45a9cf65616ea7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 20 Mar 2026 10:00:24 +0100 Subject: [PATCH 35/39] Accept AbstractSparseMatrixCSC for decompression (#298) * Accept AbstractSparseMatrixCSC for decompression * Define decompress_csc * Remove unused import of AbstractSparseMatrixCSC * Apply suggestion from @blegat * Change return statement to return nothing * Apply suggestion from @blegat * Add to docs * Apply suggestion from @gdalle Co-authored-by: Guillaume Dalle <22795598+gdalle@users.noreply.github.com> * Move to internals --------- Co-authored-by: Guillaume Dalle <22795598+gdalle@users.noreply.github.com> --- docs/src/dev.md | 1 + src/decompression.jl | 34 ++++++++++++++++++++++++++++------ 2 files changed, 29 insertions(+), 6 deletions(-) diff --git a/docs/src/dev.md b/docs/src/dev.md index a9e0d303..0f7b86fa 100644 --- a/docs/src/dev.md +++ b/docs/src/dev.md @@ -41,6 +41,7 @@ SparseMatrixColorings.TreeSetColoringResult SparseMatrixColorings.LinearSystemColoringResult SparseMatrixColorings.BicoloringResult SparseMatrixColorings.remap_colors +SparseMatrixColorings.decompress_csc! ``` ## Testing diff --git a/src/decompression.jl b/src/decompression.jl index 1528c885..849f0f6b 100644 --- a/src/decompression.jl +++ b/src/decompression.jl @@ -585,8 +585,23 @@ function decompress!( return A end -function decompress!( - A::SparseMatrixCSC{R}, +""" + decompress_csc!( + nzA::AbstractVector{R}, + A_colptr::AbstractVector, + B::AbstractMatrix{R}, + result::TreeSetColoringResult, + uplo::Symbol=:F, + ) where {R<:Real} + +Decompress the values of `B` into the vector of nonzero entries `nzA` of a +sparse matrix with column pointers `colptr`. This function assumes that the +row indices are sorted in increasing order and are the same as those of the +sparse matrix given to the `coloring` function that returned `result`. +""" +function decompress_csc!( + nzA::AbstractVector{R}, + A_colptr::AbstractVector{<:Integer}, B::AbstractMatrix{R}, result::TreeSetColoringResult, uplo::Symbol=:F, @@ -603,10 +618,6 @@ function decompress!( upper_triangle_offsets, buffer, ) = result - (; S) = ag - A_colptr = A.colptr - nzA = nonzeros(A) - check_compatible_pattern(A, ag, uplo) if eltype(buffer) == R buffer_right_type = buffer @@ -696,6 +707,17 @@ function decompress!( #! format: on end end + return nothing +end + +function decompress!( + A::SparseMatrixCSC{R}, + B::AbstractMatrix{R}, + result::TreeSetColoringResult, + uplo::Symbol=:F, +) where {R<:Real} + check_compatible_pattern(A, result.ag, uplo) + decompress_csc!(nonzeros(A), A.colptr, B, result, uplo) return A end From dd8aa54e6f33f0388885aa93176ba5ae9bb1d90c Mon Sep 17 00:00:00 2001 From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com> Date: Fri, 20 Mar 2026 10:16:13 +0100 Subject: [PATCH 36/39] Bump version (#311) --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 4fb07535..17a66304 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SparseMatrixColorings" uuid = "0a514795-09f3-496d-8182-132a7b665d35" authors = ["Guillaume Dalle", "Alexis Montoison"] -version = "0.4.25" +version = "0.4.26" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From 09fc114a3eb222faeb31d79bc35005acc182ea3e Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 30 Mar 2026 16:51:45 +0200 Subject: [PATCH 37/39] Bump codecov/codecov-action from 5 to 6 (#312) Bumps [codecov/codecov-action](https://github.com/codecov/codecov-action) from 5 to 6. - [Release notes](https://github.com/codecov/codecov-action/releases) - [Changelog](https://github.com/codecov/codecov-action/blob/main/CHANGELOG.md) - [Commits](https://github.com/codecov/codecov-action/compare/v5...v6) --- updated-dependencies: - dependency-name: codecov/codecov-action dependency-version: '6' dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/Test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/Test.yml b/.github/workflows/Test.yml index 46601d04..6b27e932 100644 --- a/.github/workflows/Test.yml +++ b/.github/workflows/Test.yml @@ -36,7 +36,7 @@ jobs: - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v5 + - uses: codecov/codecov-action@v6 with: files: lcov.info token: ${{ secrets.CODECOV_TOKEN }} From 5d1ae0abe0a56d331909d89ceae1c9b83522c005 Mon Sep 17 00:00:00 2001 From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com> Date: Thu, 16 Apr 2026 14:17:07 +0200 Subject: [PATCH 38/39] chore: bump CUDA compat to v6 (cuSPARSE now a separate package) (#314) * chore: bump CUDA compat to v6 (cuSPARSE now a separate package) * More renaming * Format * Fix compression? * Generic compression and result * Import * No generic result * Fix ambiguity --- Project.toml | 11 ++++++--- ext/SparseMatrixColoringsCUDAExt.jl | 18 +-------------- ext/SparseMatrixColoringsGPUArraysExt.jl | 29 ++++++++++++++++++++++++ test/Project.toml | 4 +--- test/cuda.jl | 4 +--- test/runtests.jl | 2 +- 6 files changed, 41 insertions(+), 27 deletions(-) create mode 100644 ext/SparseMatrixColoringsGPUArraysExt.jl diff --git a/Project.toml b/Project.toml index 17a66304..9ef7826d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SparseMatrixColorings" uuid = "0a514795-09f3-496d-8182-132a7b665d35" +version = "0.4.27" authors = ["Guillaume Dalle", "Alexis Montoison"] -version = "0.4.26" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" @@ -15,25 +15,30 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8" Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" +GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" +cuSPARSE = "b26da814-b3bc-49ef-b0ee-c816305aa060" [extensions] -SparseMatrixColoringsCUDAExt = "CUDA" +SparseMatrixColoringsCUDAExt = ["CUDA", "cuSPARSE"] SparseMatrixColoringsCliqueTreesExt = "CliqueTrees" SparseMatrixColoringsColorsExt = "Colors" +SparseMatrixColoringsGPUArraysExt = "GPUArrays" SparseMatrixColoringsJuMPExt = ["JuMP", "MathOptInterface"] [compat] ADTypes = "1.2.1" -CUDA = "5.8.2" +CUDA = "6.0.0" CliqueTrees = "1" Colors = "0.12.11, 0.13" DocStringExtensions = "0.8,0.9" +GPUArrays = "11.5.0" JuMP = "1.29.1" LinearAlgebra = "1" MathOptInterface = "1.45.0" PrecompileTools = "1.2.1" Random = "1" SparseArrays = "1" +cuSPARSE = "6.0.0" julia = "1.10" diff --git a/ext/SparseMatrixColoringsCUDAExt.jl b/ext/SparseMatrixColoringsCUDAExt.jl index 2750964d..ed33eece 100644 --- a/ext/SparseMatrixColoringsCUDAExt.jl +++ b/ext/SparseMatrixColoringsCUDAExt.jl @@ -3,23 +3,7 @@ module SparseMatrixColoringsCUDAExt import SparseMatrixColorings as SMC using SparseArrays: SparseMatrixCSC, rowvals, nnz, nzrange using CUDA: CuVector, CuMatrix -using CUDA.CUSPARSE: AbstractCuSparseMatrix, CuSparseMatrixCSC, CuSparseMatrixCSR - -SMC.matrix_versions(A::AbstractCuSparseMatrix) = (A,) - -## Compression (slow, through CPU) - -function SMC.compress( - A::AbstractCuSparseMatrix, result::SMC.AbstractColoringResult{structure,:column} -) where {structure} - return CuMatrix(SMC.compress(SparseMatrixCSC(A), result)) -end - -function SMC.compress( - A::AbstractCuSparseMatrix, result::SMC.AbstractColoringResult{structure,:row} -) where {structure} - return CuMatrix(SMC.compress(SparseMatrixCSC(A), result)) -end +using cuSPARSE: AbstractCuSparseMatrix, CuSparseMatrixCSC, CuSparseMatrixCSR ## CSC Result diff --git a/ext/SparseMatrixColoringsGPUArraysExt.jl b/ext/SparseMatrixColoringsGPUArraysExt.jl new file mode 100644 index 00000000..37574e73 --- /dev/null +++ b/ext/SparseMatrixColoringsGPUArraysExt.jl @@ -0,0 +1,29 @@ +module SparseMatrixColoringsGPUArraysExt + +using GPUArrays: AbstractGPUSparseMatrix, dense_array_type +using SparseArrays: SparseMatrixCSC +import SparseMatrixColorings as SMC + +SMC.matrix_versions(A::AbstractGPUSparseMatrix) = (A,) + +## Compression (slow, through CPU) + +function SMC.compress( + A::AbstractGPUSparseMatrix, result::SMC.AbstractColoringResult{structure,:column} +) where {structure} + A_cpu = SparseMatrixCSC(A) + B_cpu = SMC.compress(A_cpu, result) + B = dense_array_type(A)(B_cpu) + return B +end + +function SMC.compress( + A::AbstractGPUSparseMatrix, result::SMC.AbstractColoringResult{structure,:row} +) where {structure} + A_cpu = SparseMatrixCSC(A) + B_cpu = SMC.compress(A_cpu, result) + B = dense_array_type(A)(B_cpu) + return B +end + +end diff --git a/test/Project.toml b/test/Project.toml index 4c505efa..5e41ceb2 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -23,6 +23,4 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" - -[compat] -CUDA = "=5.9.5" +cuSPARSE = "b26da814-b3bc-49ef-b0ee-c816305aa060" diff --git a/test/cuda.jl b/test/cuda.jl index 6549b168..8b27e07b 100644 --- a/test/cuda.jl +++ b/test/cuda.jl @@ -1,4 +1,4 @@ -using CUDA.CUSPARSE: CuSparseMatrixCSC, CuSparseMatrixCSR +using cuSPARSE: CuSparseMatrixCSC, CuSparseMatrixCSR using LinearAlgebra using SparseArrays using SparseMatrixColorings @@ -6,8 +6,6 @@ import SparseMatrixColorings as SMC using StableRNGs using Test -include("utils.jl") - rng = StableRNG(63) asymmetric_params = vcat( diff --git a/test/runtests.jl b/test/runtests.jl index 9494e5e2..ec679d4e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -12,7 +12,7 @@ include("utils.jl") @testset verbose = true "SparseMatrixColorings" begin if get(ENV, "JULIA_SMC_TEST_GROUP", nothing) == "GPU" @testset "CUDA" begin - using CUDA + using CUDA, cuSPARSE include("cuda.jl") end else From d318be536f49025720932ca300d0648ddaa4044c Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Fri, 17 Apr 2026 09:49:23 +0200 Subject: [PATCH 39/39] Bump julia-actions/setup-julia from 2 to 3 (#315) Bumps [julia-actions/setup-julia](https://github.com/julia-actions/setup-julia) from 2 to 3. - [Release notes](https://github.com/julia-actions/setup-julia/releases) - [Commits](https://github.com/julia-actions/setup-julia/compare/v2...v3) --- updated-dependencies: - dependency-name: julia-actions/setup-julia dependency-version: '3' dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/Benchmark.yml | 2 +- .github/workflows/CompatHelper.yml | 2 +- .github/workflows/Documentation.yml | 2 +- .github/workflows/Test.yml | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/.github/workflows/Benchmark.yml b/.github/workflows/Benchmark.yml index fec67016..4a4a444b 100644 --- a/.github/workflows/Benchmark.yml +++ b/.github/workflows/Benchmark.yml @@ -13,7 +13,7 @@ jobs: if: contains(github.event.pull_request.labels.*.name, 'benchmark') steps: - uses: actions/checkout@v6 - - uses: julia-actions/setup-julia@v2 + - uses: julia-actions/setup-julia@v3 with: version: "1" - uses: julia-actions/cache@v3 diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml index 1f74cd19..492a6c0b 100644 --- a/.github/workflows/CompatHelper.yml +++ b/.github/workflows/CompatHelper.yml @@ -15,7 +15,7 @@ jobs: run: which julia continue-on-error: true - name: Install Julia, but only if it is not already available in the PATH - uses: julia-actions/setup-julia@v2 + uses: julia-actions/setup-julia@v3 with: version: '1' arch: ${{ runner.arch }} diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml index 4d85c43d..2f2e57df 100644 --- a/.github/workflows/Documentation.yml +++ b/.github/workflows/Documentation.yml @@ -20,7 +20,7 @@ jobs: runs-on: ubuntu-latest steps: - uses: actions/checkout@v6 - - uses: julia-actions/setup-julia@v2 + - uses: julia-actions/setup-julia@v3 with: version: '1' - uses: julia-actions/cache@v3 diff --git a/.github/workflows/Test.yml b/.github/workflows/Test.yml index 6b27e932..6e8c6af5 100644 --- a/.github/workflows/Test.yml +++ b/.github/workflows/Test.yml @@ -24,7 +24,7 @@ jobs: steps: - uses: actions/checkout@v6 - - uses: julia-actions/setup-julia@v2 + - uses: julia-actions/setup-julia@v3 with: version: ${{ matrix.julia-version }} arch: x64