From 8b748ad0ea8a149f42e5858da4dec75f83389231 Mon Sep 17 00:00:00 2001
From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com>
Date: Tue, 21 May 2024 15:38:26 +0200
Subject: [PATCH 1/8] Test against reference papers

---
 Project.toml                           |   2 +
 README.md                              |   2 +-
 src/SparseMatrixColorings.jl           |  22 ++---
 src/adtypes.jl                         |   8 +-
 src/coloring.jl                        |  18 ++--
 src/graph.jl                           | 118 +++++++++++++++++--------
 src/order.jl                           |  18 ++--
 test/Project.toml                      |   3 +
 test/graph.jl                          |  28 +++---
 test/order.jl                          |  36 ++++----
 test/reference/colpack_table_3.csv     |  11 +++
 test/reference/colpack_table_6_7.csv   |  11 +++
 test/reference/what_table_31_32_33.csv |  23 +++++
 test/reference/what_table_41_42.csv    |   7 ++
 test/runtests.jl                       |   3 +
 test/suitesparse.jl                    | 106 ++++++++++++++++++++++
 16 files changed, 318 insertions(+), 98 deletions(-)
 create mode 100644 test/reference/colpack_table_3.csv
 create mode 100644 test/reference/colpack_table_6_7.csv
 create mode 100644 test/reference/what_table_31_32_33.csv
 create mode 100644 test/reference/what_table_41_42.csv
 create mode 100644 test/suitesparse.jl

diff --git a/Project.toml b/Project.toml
index fa9b8b72..a826b2d6 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,12 +5,14 @@ version = "0.1.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
+DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
 ADTypes = "1.2.1"
+DocStringExtensions = "0.9"
 LinearAlgebra = "1"
 Random = "1"
 SparseArrays = "1"
diff --git a/README.md b/README.md
index 0e58b7d3..527ac0ff 100644
--- a/README.md
+++ b/README.md
@@ -21,7 +21,7 @@ The algorithms implemented in this package are mainly taken from the following a
 
 > [_What Color Is Your Jacobian? Graph Coloring for Computing Derivatives_](https://epubs.siam.org/doi/10.1137/S0036144504444711), Gebremedhin et al. (2005)
 
-> [ColPack: Software for graph coloring and related problems in scientific computing](https://dl.acm.org/doi/10.1145/2513109.2513110), Gebremedhin et al. (2013)
+> [_ColPack: Software for graph coloring and related problems in scientific computing_](https://dl.acm.org/doi/10.1145/2513109.2513110), Gebremedhin et al. (2013)
 
 Some parts of the articles (like definitions) are thus copied verbatim in the documentation.
 
diff --git a/src/SparseMatrixColorings.jl b/src/SparseMatrixColorings.jl
index 899ddbc5..05203bf1 100644
--- a/src/SparseMatrixColorings.jl
+++ b/src/SparseMatrixColorings.jl
@@ -1,22 +1,24 @@
 """
     SparseMatrixColorings
 
-Coloring algorithms for sparse Jacobian and Hessian matrices.
-
-The algorithms implemented in this package are mainly taken from the following articles:
-
-> [_What Color Is Your Jacobian? Graph Coloring for Computing Derivatives_](https://epubs.siam.org/doi/10.1137/S0036144504444711), Gebremedhin et al. (2005)
-
-> [ColPack: Software for graph coloring and related problems in scientific computing](https://dl.acm.org/doi/10.1145/2513109.2513110), Gebremedhin et al. (2013)
-
-Some parts of the articles (like definitions) are thus copied verbatim in the documentation.
+$README
 """
 module SparseMatrixColorings
 
 using ADTypes: ADTypes, AbstractColoringAlgorithm
+using DocStringExtensions
 using LinearAlgebra: Diagonal, Transpose, checksquare, parent, transpose
 using Random: AbstractRNG, default_rng, randperm
-using SparseArrays: SparseArrays, SparseMatrixCSC, nzrange, rowvals, spzeros
+using SparseArrays:
+    SparseArrays,
+    SparseMatrixCSC,
+    dropzeros,
+    dropzeros!,
+    nnz,
+    nzrange,
+    rowvals,
+    sparse,
+    spzeros
 
 include("graph.jl")
 include("order.jl")
diff --git a/src/adtypes.jl b/src/adtypes.jl
index e9bece08..e41f6a79 100644
--- a/src/adtypes.jl
+++ b/src/adtypes.jl
@@ -21,16 +21,16 @@ end
 GreedyColoringAlgorithm() = GreedyColoringAlgorithm(NaturalOrder())
 
 function ADTypes.column_coloring(A::AbstractMatrix, algo::GreedyColoringAlgorithm)
-    bg = BipartiteGraph(A)
+    bg = bipartite_graph(A)
     return partial_distance2_coloring(bg, Val(2), algo.order)
 end
 
 function ADTypes.row_coloring(A::AbstractMatrix, algo::GreedyColoringAlgorithm)
-    bg = BipartiteGraph(A)
+    bg = bipartite_graph(A)
     return partial_distance2_coloring(bg, Val(1), algo.order)
 end
 
 function ADTypes.symmetric_coloring(A::AbstractMatrix, algo::GreedyColoringAlgorithm)
-    ag = AdjacencyGraph(A)
-    return star_coloring(ag, algo.order)
+    ag = adjacency_graph(A)
+    return star_coloring1(ag, algo.order)
 end
diff --git a/src/coloring.jl b/src/coloring.jl
index e194e3f4..dfd4de70 100644
--- a/src/coloring.jl
+++ b/src/coloring.jl
@@ -37,9 +37,9 @@ function partial_distance2_coloring(
 end
 
 """
-    star_coloring(ag::AdjacencyGraph, order::AbstractOrder)
+    star_coloring1(g::Graph, order::AbstractOrder)
 
-Compute a star coloring of all vertices in the adjacency graph `ag` and return a vector of integer colors.
+Compute a star coloring of all vertices in the adjacency graph `g` and return a vector of integer colors.
 
 A _star coloring_ is a distance-1 coloring such that every path on 4 vertices uses at least 3 colors.
 
@@ -47,23 +47,23 @@ The vertices are colored in a greedy fashion, following the `order` supplied.
 
 # See also
 
-- [`AdjacencyGraph`](@ref)
+- [`Graph`](@ref)
 - [`AbstractOrder`](@ref)
 """
-function star_coloring(ag::AdjacencyGraph, order::AbstractOrder)
-    n = length(ag)
+function star_coloring1(g::Graph, order::AbstractOrder)
+    n = length(g)
     colors = zeros(Int, n)
     forbidden_colors = zeros(Int, n)
-    for v in vertices(ag, order)
-        for w in neighbors(ag, v)
+    for v in vertices(g, order)
+        for w in neighbors(g, v)
             if !iszero(colors[w])  # w is colored
                 forbidden_colors[colors[w]] = v
             end
-            for x in neighbors(ag, w)
+            for x in neighbors(g, w)
                 if !iszero(colors[x]) && iszero(colors[w])  # w is not colored
                     forbidden_colors[colors[x]] = v
                 else
-                    for y in neighbors(ag, x)
+                    for y in neighbors(g, x)
                         if !iszero(colors[y]) && y != w
                             if colors[y] == colors[w]
                                 forbidden_colors[colors[x]] = v
diff --git a/src/graph.jl b/src/graph.jl
index 68c0f025..e62fde1e 100644
--- a/src/graph.jl
+++ b/src/graph.jl
@@ -1,52 +1,101 @@
+## Standard graph
+
+"""
+    Graph{T}
+
+Undirected graph structure stored in Compressed Sparse Column (CSC) format.
+
+# Fields
+
+- `colptr::Vector{T}`: same as for `SparseMatrixCSC`
+- `rowval::Vector{T}`: same as for `SparseMatrixCSC`
+"""
 struct Graph{T<:Integer}
     colptr::Vector{T}
     rowval::Vector{T}
 end
 
 Graph(A::SparseMatrixCSC) = Graph(A.colptr, A.rowval)
+Graph(A::AbstractMatrix) = Graph(sparse(A))
 
 Base.length(g::Graph) = length(g.colptr) - 1
+SparseArrays.nnz(g::Graph) = length(g.rowval)
 
+vertices(g::Graph) = 1:length(g)
 neighbors(g::Graph, v::Integer) = view(g.rowval, g.colptr[v]:(g.colptr[v + 1] - 1))
+degree(g::Graph, v::Integer) = length(g.colptr[v]:(g.colptr[v + 1] - 1))
 
-## Adjacency graph
+maximum_degree(g::Graph) = maximum(Base.Fix1(degree, g), vertices(g))
+minimum_degree(g::Graph) = minimum(Base.Fix1(degree, g), vertices(g))
+
+## Bipartite graph
 
 """
-    AdjacencyGraph
+    BipartiteGraph{T}
 
-Undirected graph representing the nonzeros of a symmetrix matrix (typically a Hessian matrix).
+Undirected bipartite graph structure stored in bidirectional Compressed Sparse Column format (redundancy allows for faster access).
 
-The adjacency graph of a symmetrix matrix `A ∈ ℝ^{n × n}` is `G(A) = (V, E)` where
+A bipartite graph has two "sides", which we number `1` and `2`.
 
-- `V = 1:n` is the set of rows or columns
-- `(i, j) ∈ E` whenever `A[i, j] ≠ 0` and `i ≠ j`
+# Fields
+
+- `g1::Graph{T}`: contains the neighbors for vertices on side `1`
+- `g2::Graph{T}`: contains the neighbors for vertices on side `2`
 """
-struct AdjacencyGraph{T}
-    g::Graph{T}
+struct BipartiteGraph{T<:Integer}
+    g1::Graph{T}
+    g2::Graph{T}
 end
 
-function AdjacencyGraph(H::SparseMatrixCSC)
-    g = Graph(H - Diagonal(H))
-    return AdjacencyGraph(g)
-end
+Base.length(bg::BipartiteGraph, ::Val{1}) = length(bg.g1)
+Base.length(bg::BipartiteGraph, ::Val{2}) = length(bg.g2)
+SparseArrays.nnz(bg::BipartiteGraph) = nnz(bg.g1)
 
-Base.length(ag::AdjacencyGraph) = length(ag.g)
+"""
+    vertices(bg::BipartiteGraph, Val(side))
 
+Return the list of vertices of `bg` from the specified `side` as a range `1:n`.
 """
-    neighbors(ag::AdjacencyGraph, v::Integer)
+vertices(bg::BipartiteGraph, ::Val{side}) where {side} = 1:length(bg, Val(side))
 
-Return the neighbors of `v` in the graph `ag`.
 """
-neighbors(ag::AdjacencyGraph, v::Integer) = neighbors(ag.g, v)
+    neighbors(bg::BipartiteGraph, Val(side), v::Integer)
 
-degree(ag::AdjacencyGraph, v::Integer) = length(neighbors(ag, v))
+Return the neighbors of `v` (a vertex from the specified `side`, `1` or `2`), in the graph `bg`.
+"""
+neighbors(bg::BipartiteGraph, ::Val{1}, v::Integer) = neighbors(bg.g1, v)
+neighbors(bg::BipartiteGraph, ::Val{2}, v::Integer) = neighbors(bg.g2, v)
 
-## Bipartite graph
+degree(bg::BipartiteGraph, ::Val{1}, v::Integer) = degree(bg.g1, v)
+degree(bg::BipartiteGraph, ::Val{2}, v::Integer) = degree(bg.g2, v)
+
+function maximum_degree(bg::BipartiteGraph, ::Val{side}) where {side}
+    return maximum(v -> degree(bg, Val(side), v), vertices(bg, Val(side)))
+end
+
+function minimum_degree(bg::BipartiteGraph, ::Val{side}) where {side}
+    return minimum(v -> degree(bg, Val(side), v), vertices(bg, Val(side)))
+end
+
+## Construct from matrices
+
+"""
+    adjacency_graph(H::AbstractMatrix)
+
+Return a [`Graph`](@ref) representing the nonzeros of a symmetric matrix (typically a Hessian matrix).
+
+The adjacency graph of a symmetrix matric `A ∈ ℝ^{n × n}` is `G(A) = (V, E)` where
+
+- `V = 1:n` is the set of rows or columns `i`/`j`
+- `(i, j) ∈ E` whenever `A[i, j] ≠ 0` and `i ≠ j`
+"""
+adjacency_graph(H::SparseMatrixCSC) = Graph(H - Diagonal(H))
+adjacency_graph(H::AbstractMatrix) = adjacency_graph(sparse(H))
 
 """
-    BipartiteGraph
+    bipartite_graph(J::AbstractMatrix)
 
-Undirected bipartite graph representing the nonzeros of a non-symmetric matrix (typically a Jacobian matrix).
+Return a [`BipartiteGraph`](@ref) representing the nonzeros of a non-symmetric matrix (typically a Jacobian matrix).
 
 The bipartite graph of a matrix `A ∈ ℝ^{m × n}` is `Gb(A) = (V₁, V₂, E)` where
 
@@ -54,28 +103,27 @@ The bipartite graph of a matrix `A ∈ ℝ^{m × n}` is `Gb(A) = (V₁, V₂, E)
 - `V₂ = 1:n` is the set of columns `j`
 - `(i, j) ∈ E` whenever `A[i, j] ≠ 0`
 """
-struct BipartiteGraph{T}
-    g1::Graph{T}
-    g2::Graph{T}
-end
-
-function BipartiteGraph(J::SparseMatrixCSC)
+function bipartite_graph(J::SparseMatrixCSC)
     g1 = Graph(SparseMatrixCSC(transpose(J)))  # rows to columns
     g2 = Graph(J)  # columns to rows
     return BipartiteGraph(g1, g2)
 end
 
-Base.length(bg::BipartiteGraph, ::Val{1}) = length(bg.g1)
-Base.length(bg::BipartiteGraph, ::Val{2}) = length(bg.g2)
+bipartite_graph(J::AbstractMatrix) = bipartite_graph(sparse(J))
 
 """
-    neighbors(bg::BipartiteGraph, Val(side), v::Integer)
+    column_intersection_graph(J::AbstractMatrix)
 
-Return the neighbors of `v`, which is a vertex from the specified `side` (`1` or `2`), in the graph `bg`.
-"""
-neighbors(bg::BipartiteGraph, ::Val{1}, v::Integer) = neighbors(bg.g1, v)
-neighbors(bg::BipartiteGraph, ::Val{2}, v::Integer) = neighbors(bg.g2, v)
+Return a [`Graph`](@ref) representing the column intersections of a non-symmetric matrix (typically a Jacobian matrix).
+
+The column intersection graph of a matrix `A ∈ ℝ^{m × n}` is `Gc(A) = (V, E)` where
 
-function degree(bg::BipartiteGraph, ::Val{side}, v::Integer) where {side}
-    return length(neighbors(bg, Val(side), v))
+- `V = 1:n` is the set of columns `j`
+- `(j1, j2) ∈ E` whenever `A[:, j1] ∩ A[:, j2] ≠ ∅`
+"""
+function column_intersection_graph(J::SparseMatrixCSC)
+    A = transpose(J) * J
+    return adjacency_graph(A - Diagonal(A))
 end
+
+column_intersection_graph(J::AbstractMatrix) = column_intersection_graph(sparse(J))
diff --git a/src/order.jl b/src/order.jl
index b7a9c1e1..b49973eb 100644
--- a/src/order.jl
+++ b/src/order.jl
@@ -18,12 +18,12 @@ Order vertices as they come in the graph.
 """
 struct NaturalOrder <: AbstractOrder end
 
-function vertices(ag::AdjacencyGraph, ::NaturalOrder)
-    return 1:length(ag)
+function vertices(g::Graph, ::NaturalOrder)
+    return vertices(g)
 end
 
 function vertices(bg::BipartiteGraph, ::Val{side}, ::NaturalOrder) where {side}
-    return 1:length(bg, Val(side))
+    return vertices(bg, Val(side))
 end
 
 """
@@ -37,8 +37,8 @@ end
 
 RandomOrder() = RandomOrder(default_rng())
 
-function vertices(ag::AdjacencyGraph, order::RandomOrder)
-    return randperm(order.rng, length(ag))
+function vertices(g::Graph, order::RandomOrder)
+    return randperm(order.rng, length(g))
 end
 
 function vertices(bg::BipartiteGraph, ::Val{side}, order::RandomOrder) where {side}
@@ -52,12 +52,12 @@ Order vertices by decreasing degree.
 """
 struct LargestFirst <: AbstractOrder end
 
-function vertices(ag::AdjacencyGraph, ::LargestFirst)
-    criterion(v) = degree(ag, v)
-    return sort(1:length(ag); by=criterion, rev=true)
+function vertices(g::Graph, ::LargestFirst)
+    criterion(v) = degree(g, v)
+    return sort(vertices(g); by=criterion, rev=true)
 end
 
 function vertices(bg::BipartiteGraph, ::Val{side}, ::LargestFirst) where {side}
     criterion(v) = degree(bg, Val(side), v)
-    return sort(1:length(bg, Val(side)); by=criterion, rev=true)
+    return sort(vertices(bg, Val(side)); by=criterion, rev=true)
 end
diff --git a/test/Project.toml b/test/Project.toml
index 617bcb56..f71c8c9e 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,10 +1,13 @@
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MatrixDepot = "b51810bb-c9f3-55da-ae3c-350fc1fbce05"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/graph.jl b/test/graph.jl
index e4f6a703..7c67ca42 100644
--- a/test/graph.jl
+++ b/test/graph.jl
@@ -1,15 +1,15 @@
 using SparseArrays
-using SparseMatrixColorings: Graph, BipartiteGraph, AdjacencyGraph, neighbors
+using SparseMatrixColorings: Graph, adjacency_graph, bipartite_graph, neighbors
 using Test
 
 ## Standard graph
 
 @testset "Graph" begin
-    g = Graph(sparse([
+    g = Graph([
         1 0 1
         1 1 0
         0 0 0
-    ]))
+    ])
 
     @test length(g) == 3
     @test neighbors(g, 1) == [1, 2]
@@ -27,7 +27,7 @@ end;
         0 0 0 1 1 1 1 0
     ])
 
-    bg = BipartiteGraph(A)
+    bg = bipartite_graph(A)
     @test length(bg, Val(1)) == 4
     @test length(bg, Val(2)) == 8
     # neighbors of rows
@@ -56,14 +56,14 @@ end;
         0 0 0 1 1 1 1 0
     ])
 
-    ag = AdjacencyGraph(transpose(A) * A)
-    @test length(ag) == 8
-    @test neighbors(ag, 1) == [6, 7, 8]
-    @test neighbors(ag, 2) == [5, 7, 8]
-    @test neighbors(ag, 3) == [5, 6, 8]
-    @test neighbors(ag, 4) == [5, 6, 7]
-    @test neighbors(ag, 5) == [2, 3, 4, 6, 7, 8]
-    @test neighbors(ag, 6) == [1, 3, 4, 5, 7, 8]
-    @test neighbors(ag, 7) == [1, 2, 4, 5, 6, 8]
-    @test neighbors(ag, 8) == [1, 2, 3, 5, 6, 7]
+    g = adjacency_graph(transpose(A) * A)
+    @test length(g) == 8
+    @test neighbors(g, 1) == [6, 7, 8]
+    @test neighbors(g, 2) == [5, 7, 8]
+    @test neighbors(g, 3) == [5, 6, 8]
+    @test neighbors(g, 4) == [5, 6, 7]
+    @test neighbors(g, 5) == [2, 3, 4, 6, 7, 8]
+    @test neighbors(g, 6) == [1, 3, 4, 5, 7, 8]
+    @test neighbors(g, 7) == [1, 2, 4, 5, 6, 8]
+    @test neighbors(g, 8) == [1, 2, 3, 5, 6, 7]
 end
diff --git a/test/order.jl b/test/order.jl
index e89ad448..77c5cb25 100644
--- a/test/order.jl
+++ b/test/order.jl
@@ -1,36 +1,40 @@
 using SparseArrays
-using SparseMatrixColorings: Graph, LargestFirst, NaturalOrder, RandomOrder, vertices
-using StableRNGs
+using SparseMatrixColorings:
+    Graph,
+    adjacency_graph,
+    bipartite_graph,
+    LargestFirst,
+    NaturalOrder,
+    RandomOrder,
+    vertices
 using Test
 
-rng = StableRNG(63)
-
 @testset "NaturalOrder" begin
     A = sprand(rng, Bool, 5, 5, 0.5)
-    ag = AdjacencyGraph(A)
+    ag = adjacency_graph(A)
     @test vertices(ag, NaturalOrder()) == 1:5
 
     A = sprand(rng, Bool, 5, 4, 0.5)
-    bg = BipartiteGraph(A)
+    bg = bipartite_graph(A)
     @test vertices(bg, Val(1), NaturalOrder()) == 1:5
 
     A = sprand(rng, Bool, 5, 4, 0.5)
-    bg = BipartiteGraph(A)
+    bg = bipartite_graph(A)
     @test vertices(bg, Val(2), NaturalOrder()) == 1:4
 end;
 
 @testset "RandomOrder" begin
     A = sprand(rng, Bool, 5, 5, 0.5)
-    ag = AdjacencyGraph(A)
-    @test sort(vertices(ag, RandomOrder(rng))) == 1:5
+    ag = adjacency_graph(A)
+    @test sort(vertices(ag, RandomOrder())) == 1:5
 
     A = sprand(rng, Bool, 5, 4, 0.5)
-    bg = BipartiteGraph(A)
-    @test sort(vertices(bg, Val(1), RandomOrder(rng))) == 1:5
+    bg = bipartite_graph(A)
+    @test sort(vertices(bg, Val(1), RandomOrder())) == 1:5
 
     A = sprand(rng, Bool, 5, 4, 0.5)
-    bg = BipartiteGraph(A)
-    @test sort(vertices(bg, Val(2), RandomOrder(rng))) == 1:4
+    bg = bipartite_graph(A)
+    @test sort(vertices(bg, Val(2), RandomOrder())) == 1:4
 end;
 
 @testset "LargestFirst" begin
@@ -39,7 +43,7 @@ end;
         1 0 0
         0 1 0
     ])
-    ag = AdjacencyGraph(A)
+    ag = adjacency_graph(A)
 
     @test vertices(ag, LargestFirst()) == [2, 1, 3]
 
@@ -49,7 +53,7 @@ end;
         0 1 1
         0 0 0
     ])
-    bg = BipartiteGraph(A)
+    bg = bipartite_graph(A)
 
     @test vertices(bg, Val(1), LargestFirst()) == [1, 3, 2, 4]
 
@@ -58,7 +62,7 @@ end;
         1 0 1 0
         1 0 1 0
     ])
-    bg = BipartiteGraph(A)
+    bg = bipartite_graph(A)
 
     @test vertices(bg, Val(2), LargestFirst()) == [1, 3, 2, 4]
 end;
diff --git a/test/reference/colpack_table_3.csv b/test/reference/colpack_table_3.csv
new file mode 100644
index 00000000..b75e6b1a
--- /dev/null
+++ b/test/reference/colpack_table_3.csv
@@ -0,0 +1,11 @@
+row,modified,group,name,V,E,Δ,LF,N
+1,false,INPRO,msdoor,415863,9912536,76,42,42
+2,false,GHS_psdef,ldoor,952203,22785136,76,42,42
+3,false,DNVS,shipsec1,140874,3836265,101,54,48
+4,false,DNVS,shipsec5,179860,4966618,125,48,50
+5,false,Chen,pkustk11,87804,2565054,131,54,66
+6,false,Boeing,ct20stif,52329,1323067,206,47,49
+7,false,Boeing,pwtk,217918,5708253,179,42,48
+8,false,Chen,pkustk13,94893,3260967,299,45,57
+9,false,Nasa,nasasrb,54870,1311227,275,40,41
+10,false,GHS_indef,bmw3_2,227362,5530634,335,48,48
diff --git a/test/reference/colpack_table_6_7.csv b/test/reference/colpack_table_6_7.csv
new file mode 100644
index 00000000..7348ef07
--- /dev/null
+++ b/test/reference/colpack_table_6_7.csv
@@ -0,0 +1,11 @@
+row,modified,group,name,V1,V2,E,Δ1,Δ2,LF2,N2,LF1,N1
+1,false,LPnetlib,lp_cre_a,3516,7248,18168,360,14,360,360,14,16
+2,false,LPnetlib,lp_ken_11,14694,21349,49058,122,3,128,130,4,5
+3,false,LPnetlib,lp_ken_13,28632,42659,97246,170,3,174,176,5,4
+4,false,LPnetlib,lp_maros_r7,3136,9408,144848,48,46,70,74,100,72
+5,false,LPnetlib,lp_cre_d,8926,73948,246614,808,13,808,813,13,15
+6,false,LPnetlib,lp_ken_18,105127,154699,358171,325,3,325,330,5,5
+7,false,Bai,af23560,23560,23560,460598,21,21,59,44,43,43
+8,false,Shen,e40r0100,17281,17281,553562,62,62,85,95,85,87
+9,false,vanHeukelum,cage11,39082,39082,559722,31,31,70,81,67,81
+10,false,vanHeukelum,cage12,130228,130228,2032536,33,33,79,96,73,96
diff --git a/test/reference/what_table_31_32_33.csv b/test/reference/what_table_31_32_33.csv
new file mode 100644
index 00000000..379122e3
--- /dev/null
+++ b/test/reference/what_table_31_32_33.csv
@@ -0,0 +1,23 @@
+row,modified,group,name,m,n,nnz,κmax,κmin,ρmax,ρmin,K,Ec
+1,false,LPnetlib,lp_cre_a,3516,7248,18168,14,1,360,0,360,253411
+2,false,LPnetlib,lp_dfl001,6071,12230,35632,14,1,228,2,228,250976
+3,false,LPnetlib,lp_ken_11,14694,21349,49058,3,1,122,1,130,459921
+4,false,LPnetlib,lp_stocfor3,16675,23541,76473,18,1,15,1,16,125969
+5,false,LPnetlib,lp_ken_13,28632,42659,97246,3,1,170,1,176,1158664
+6,false,LPnetlib,lp_pds_10,16558,49932,107605,3,1,96,1,96,594681
+7,false,LPnetlib,lp_maros_r7,3136,9408,144848,46,1,48,5,74,610760
+8,false,Mallya,lhr10,10672,10672,232633,36,1,63,1,65,431411
+9,false,LPnetlib,lp_pds_20,33874,108175,232647,3,1,96,0,96,1325891
+10,false,LPnetlib,lp_cre_d,8926,73948,246614,13,1,808,0,813,21347885
+11,false,LPnetlib,lp_cre_b,9648,77137,260785,14,1,844,0,845,20852569
+12,false,,e30r2000,9661,9661,306356,62,8,62,8,65,688848
+13,false,Mallya,lhr14,14270,14270,307858,36,1,63,1,65,572463
+14,false,LPnetlib,lp_ken_18,105127,154699,358171,3,1,325,1,330,8412174
+15,false,Bai,af23560,23560,23560,484256,21,10,21,11,32,1210004
+16,["nnz"],Shen,e40r0100,17281,17281,553562,62,8,62,8,66,1254328
+17,false,vanHeukelum,cage11,39082,39082,559722,31,3,31,3,81,1887384
+18,false,Mallya,lhr34,35152,35152,764014,36,1,63,1,65,1417888
+19,["m"],Mallya,lhr71c,70304,70304,1528092,36,1,63,1,65,2835968
+20,false,vanHeukelum,cage12,130228,130228,2032536,33,5,33,5,96,7550823
+21,false,LPnetlib,lp_osa_07,1118,25067,144812,6,1,17613,18,17613,
+22,false,LPnetlib,lp_fit2d,25,10524,129042,17,1,10500,1427,10501,
\ No newline at end of file
diff --git a/test/reference/what_table_41_42.csv b/test/reference/what_table_41_42.csv
new file mode 100644
index 00000000..34050f35
--- /dev/null
+++ b/test/reference/what_table_41_42.csv
@@ -0,0 +1,7 @@
+row,modified,group,name,V,E,∆,δ,K_d2,K_star1,K_star2
+1,false,,mrng1,257000,505048,4,2,12,8,10
+2,false,,mrng2,1017253,2015714,4,2,12,9,10
+3,false,DIMACS10,598a,110971,741934,26,5 13,38,27,32
+4,false,DIMACS10,144,144649,1074393,26,4,41,28,35
+5,false,DIMACS10,m14b,214765,1679018,40,4,42,29,34
+6,false,DIMACS10,auto,448695,3314611,37,4,42,29,36
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 1486dce2..d3041c89 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -36,4 +36,7 @@ using Test
     @testset "ADTypes" begin
         include("adtypes.jl")
     end
+    @testset "SuiteSparse" begin
+        include("suitesparse.jl")
+    end
 end
diff --git a/test/suitesparse.jl b/test/suitesparse.jl
new file mode 100644
index 00000000..ebc0e23b
--- /dev/null
+++ b/test/suitesparse.jl
@@ -0,0 +1,106 @@
+using CSV
+using DataFrames
+using LinearAlgebra
+using MatrixDepot
+using SparseArrays
+using SparseMatrixColorings:
+    Graph,
+    adjacency_graph,
+    bipartite_graph,
+    column_intersection_graph,
+    LargestFirst,
+    NaturalOrder,
+    degree,
+    minimum_degree,
+    maximum_degree,
+    neighbors,
+    partial_distance2_coloring,
+    star_coloring1,
+    vertices
+using Test
+
+## Distance-2 coloring
+
+#=
+Comparison with Tables VI and VII of the ColPack paper
+=#
+
+colpack_table_6_7 = CSV.read(
+    joinpath(@__DIR__, "reference", "colpack_table_6_7.csv"), DataFrame
+)
+
+@testset "Distance-2 coloring (ColPack paper)" begin
+    @testset "$(row[:name])" for row in eachrow(colpack_table_6_7)
+        @info "Testing distance-2 coloring for $(row[:name])"
+        original_mat = matrixdepot("$(row[:group])/$(row[:name])")
+        mat = dropzeros(original_mat)
+        bg = bipartite_graph(mat)
+        @test length(bg, Val(1)) == row[:V1]
+        @test length(bg, Val(2)) == row[:V2]
+        @test nnz(bg) == row[:E]
+        @test maximum_degree(bg, Val(1)) == row[:Δ1]
+        @test maximum_degree(bg, Val(2)) == row[:Δ2]
+        colors_N1 = partial_distance2_coloring(bg, Val(1), NaturalOrder())
+        colors_N2 = partial_distance2_coloring(bg, Val(2), NaturalOrder())
+        @test length(unique(colors_N1)) == row[:N1]
+        @test length(unique(colors_N2)) == row[:N2]
+    end
+end;
+
+#=
+Comparison with Tables 3.1 and 3.2 of "What color is your Jacobian?"
+=#
+
+what_table_31_32_33 = CSV.read(
+    joinpath(@__DIR__, "reference", "what_table_31_32_33.csv"), DataFrame
+)
+
+@testset "Distance-2 coloring (survey paper)" begin
+    @testset "$(row[:name])" for row in eachrow(what_table_31_32_33)
+        ismissing(row[:group]) && continue
+        @info "Testing distance-2 coloring for $(row[:name])"
+        original_mat = matrixdepot("$(row[:group])/$(row[:name])")
+        mat = original_mat  # no dropzeros
+        bg = bipartite_graph(mat)
+        @test length(bg, Val(1)) == row[:m]
+        @test length(bg, Val(2)) == row[:n]
+        @test nnz(bg) == row[:nnz]
+        @test minimum_degree(bg, Val(1)) == row[:ρmin]
+        @test maximum_degree(bg, Val(1)) == row[:ρmax]
+        @test minimum_degree(bg, Val(2)) == row[:κmin]
+        @test maximum_degree(bg, Val(2)) == row[:κmax]
+        colors_Nb = partial_distance2_coloring(bg, Val(2), NaturalOrder())
+        if length(unique(colors_Nb)) == row[:K]
+            @test length(unique(colors_Nb)) == row[:K]
+        else
+            @test_broken length(unique(colors_Nb)) == row[:K]
+        end
+    end
+end;
+
+## Star coloring
+
+#=
+Comparison with Tables 4.1 and 4.2 of "What color is your Jacobian?"
+=#
+
+what_table_41_42 = CSV.read(
+    joinpath(@__DIR__, "reference", "what_table_41_42.csv"), DataFrame
+)
+
+@testset "Star coloring (survey paper)" begin
+    @testset "$(row[:name])" for row in eachrow(what_table_41_42)
+        ismissing(row[:group]) && continue
+        @info "Testing star coloring for $(row[:name])"
+        original_mat = matrixdepot("$(row[:group])/$(row[:name])")
+        mat = original_mat  # no dropzeros
+        bg = bipartite_graph(mat)
+        g = adjacency_graph(mat)
+        @test length(g) == row[:V]
+        @test nnz(g) ÷ 2 == row[:E]
+        colors_Nb = partial_distance2_coloring(bg, Val(2), NaturalOrder())
+        colors_Ns1 = star_coloring1(g, NaturalOrder())
+        @test length(unique(colors_Nb)) == row[:K_d2]
+        @test length(unique(colors_Ns1)) == row[:K_star1]
+    end
+end;

From 4458272dab0ffe913a7ce586802f249f831a2ba1 Mon Sep 17 00:00:00 2001
From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com>
Date: Tue, 21 May 2024 15:38:57 +0200
Subject: [PATCH 2/8] Rng

---
 test/order.jl | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/test/order.jl b/test/order.jl
index 77c5cb25..7054e18e 100644
--- a/test/order.jl
+++ b/test/order.jl
@@ -9,6 +9,8 @@ using SparseMatrixColorings:
     vertices
 using Test
 
+rng = StableRNG(63)
+
 @testset "NaturalOrder" begin
     A = sprand(rng, Bool, 5, 5, 0.5)
     ag = adjacency_graph(A)
@@ -26,14 +28,17 @@ end;
 @testset "RandomOrder" begin
     A = sprand(rng, Bool, 5, 5, 0.5)
     ag = adjacency_graph(A)
+    @test sort(vertices(ag, RandomOrder(rng))) == 1:5
     @test sort(vertices(ag, RandomOrder())) == 1:5
 
     A = sprand(rng, Bool, 5, 4, 0.5)
     bg = bipartite_graph(A)
+    @test sort(vertices(bg, Val(1), RandomOrder(rng))) == 1:5
     @test sort(vertices(bg, Val(1), RandomOrder())) == 1:5
 
     A = sprand(rng, Bool, 5, 4, 0.5)
     bg = bipartite_graph(A)
+    @test sort(vertices(bg, Val(2), RandomOrder(rng))) == 1:4
     @test sort(vertices(bg, Val(2), RandomOrder())) == 1:4
 end;
 

From 193b203d7483d2aa4c25073df935a48408d8887f Mon Sep 17 00:00:00 2001
From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com>
Date: Tue, 21 May 2024 15:40:35 +0200
Subject: [PATCH 3/8] Using

---
 test/order.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/order.jl b/test/order.jl
index 7054e18e..5b872bfd 100644
--- a/test/order.jl
+++ b/test/order.jl
@@ -7,6 +7,7 @@ using SparseMatrixColorings:
     NaturalOrder,
     RandomOrder,
     vertices
+using StableRNGs
 using Test
 
 rng = StableRNG(63)

From 3a38cd1067df35f39d5777fa49f1aac9e860e625 Mon Sep 17 00:00:00 2001
From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com>
Date: Tue, 21 May 2024 15:43:11 +0200
Subject: [PATCH 4/8] Rename

---
 test/reference/colpack_table_3.csv     | 11 -----------
 test/reference/what_table_31_32.csv    | 23 +++++++++++++++++++++++
 test/reference/what_table_31_32_33.csv | 23 -----------------------
 test/suitesparse.jl                    | 12 ++++++------
 4 files changed, 29 insertions(+), 40 deletions(-)
 delete mode 100644 test/reference/colpack_table_3.csv
 create mode 100644 test/reference/what_table_31_32.csv
 delete mode 100644 test/reference/what_table_31_32_33.csv

diff --git a/test/reference/colpack_table_3.csv b/test/reference/colpack_table_3.csv
deleted file mode 100644
index b75e6b1a..00000000
--- a/test/reference/colpack_table_3.csv
+++ /dev/null
@@ -1,11 +0,0 @@
-row,modified,group,name,V,E,Δ,LF,N
-1,false,INPRO,msdoor,415863,9912536,76,42,42
-2,false,GHS_psdef,ldoor,952203,22785136,76,42,42
-3,false,DNVS,shipsec1,140874,3836265,101,54,48
-4,false,DNVS,shipsec5,179860,4966618,125,48,50
-5,false,Chen,pkustk11,87804,2565054,131,54,66
-6,false,Boeing,ct20stif,52329,1323067,206,47,49
-7,false,Boeing,pwtk,217918,5708253,179,42,48
-8,false,Chen,pkustk13,94893,3260967,299,45,57
-9,false,Nasa,nasasrb,54870,1311227,275,40,41
-10,false,GHS_indef,bmw3_2,227362,5530634,335,48,48
diff --git a/test/reference/what_table_31_32.csv b/test/reference/what_table_31_32.csv
new file mode 100644
index 00000000..e518d9d9
--- /dev/null
+++ b/test/reference/what_table_31_32.csv
@@ -0,0 +1,23 @@
+row,modified,group,name,m,n,nnz,κmax,κmin,ρmax,ρmin,K
+1,false,LPnetlib,lp_cre_a,3516,7248,18168,14,1,360,0,360
+2,false,LPnetlib,lp_dfl001,6071,12230,35632,14,1,228,2,228
+3,false,LPnetlib,lp_ken_11,14694,21349,49058,3,1,122,1,130
+4,false,LPnetlib,lp_stocfor3,16675,23541,76473,18,1,15,1,16
+5,false,LPnetlib,lp_ken_13,28632,42659,97246,3,1,170,1,176
+6,false,LPnetlib,lp_pds_10,16558,49932,107605,3,1,96,1,96
+7,false,LPnetlib,lp_maros_r7,3136,9408,144848,46,1,48,5,74
+8,false,Mallya,lhr10,10672,10672,232633,36,1,63,1,65
+9,false,LPnetlib,lp_pds_20,33874,108175,232647,3,1,96,0,96
+10,false,LPnetlib,lp_cre_d,8926,73948,246614,13,1,808,0,813
+11,false,LPnetlib,lp_cre_b,9648,77137,260785,14,1,844,0,845
+12,false,,e30r2000,9661,9661,306356,62,8,62,8,65
+13,false,Mallya,lhr14,14270,14270,307858,36,1,63,1,65
+14,false,LPnetlib,lp_ken_18,105127,154699,358171,3,1,325,1,330
+15,false,Bai,af23560,23560,23560,484256,21,10,21,11,32
+16,["nnz"],Shen,e40r0100,17281,17281,553562,62,8,62,8,66
+17,false,vanHeukelum,cage11,39082,39082,559722,31,3,31,3,81
+18,false,Mallya,lhr34,35152,35152,764014,36,1,63,1,65
+19,["m"],Mallya,lhr71c,70304,70304,1528092,36,1,63,1,65
+20,false,vanHeukelum,cage12,130228,130228,2032536,33,5,33,5,96
+21,false,LPnetlib,lp_osa_07,1118,25067,144812,6,1,17613,18,17613,
+22,false,LPnetlib,lp_fit2d,25,10524,129042,17,1,10500,1427,10501,
\ No newline at end of file
diff --git a/test/reference/what_table_31_32_33.csv b/test/reference/what_table_31_32_33.csv
deleted file mode 100644
index 379122e3..00000000
--- a/test/reference/what_table_31_32_33.csv
+++ /dev/null
@@ -1,23 +0,0 @@
-row,modified,group,name,m,n,nnz,κmax,κmin,ρmax,ρmin,K,Ec
-1,false,LPnetlib,lp_cre_a,3516,7248,18168,14,1,360,0,360,253411
-2,false,LPnetlib,lp_dfl001,6071,12230,35632,14,1,228,2,228,250976
-3,false,LPnetlib,lp_ken_11,14694,21349,49058,3,1,122,1,130,459921
-4,false,LPnetlib,lp_stocfor3,16675,23541,76473,18,1,15,1,16,125969
-5,false,LPnetlib,lp_ken_13,28632,42659,97246,3,1,170,1,176,1158664
-6,false,LPnetlib,lp_pds_10,16558,49932,107605,3,1,96,1,96,594681
-7,false,LPnetlib,lp_maros_r7,3136,9408,144848,46,1,48,5,74,610760
-8,false,Mallya,lhr10,10672,10672,232633,36,1,63,1,65,431411
-9,false,LPnetlib,lp_pds_20,33874,108175,232647,3,1,96,0,96,1325891
-10,false,LPnetlib,lp_cre_d,8926,73948,246614,13,1,808,0,813,21347885
-11,false,LPnetlib,lp_cre_b,9648,77137,260785,14,1,844,0,845,20852569
-12,false,,e30r2000,9661,9661,306356,62,8,62,8,65,688848
-13,false,Mallya,lhr14,14270,14270,307858,36,1,63,1,65,572463
-14,false,LPnetlib,lp_ken_18,105127,154699,358171,3,1,325,1,330,8412174
-15,false,Bai,af23560,23560,23560,484256,21,10,21,11,32,1210004
-16,["nnz"],Shen,e40r0100,17281,17281,553562,62,8,62,8,66,1254328
-17,false,vanHeukelum,cage11,39082,39082,559722,31,3,31,3,81,1887384
-18,false,Mallya,lhr34,35152,35152,764014,36,1,63,1,65,1417888
-19,["m"],Mallya,lhr71c,70304,70304,1528092,36,1,63,1,65,2835968
-20,false,vanHeukelum,cage12,130228,130228,2032536,33,5,33,5,96,7550823
-21,false,LPnetlib,lp_osa_07,1118,25067,144812,6,1,17613,18,17613,
-22,false,LPnetlib,lp_fit2d,25,10524,129042,17,1,10500,1427,10501,
\ No newline at end of file
diff --git a/test/suitesparse.jl b/test/suitesparse.jl
index ebc0e23b..aff3da82 100644
--- a/test/suitesparse.jl
+++ b/test/suitesparse.jl
@@ -31,7 +31,7 @@ colpack_table_6_7 = CSV.read(
 
 @testset "Distance-2 coloring (ColPack paper)" begin
     @testset "$(row[:name])" for row in eachrow(colpack_table_6_7)
-        @info "Testing distance-2 coloring for $(row[:name])"
+        @info "Testing distance-2 coloring for $(row[:name]) against ColPack paper"
         original_mat = matrixdepot("$(row[:group])/$(row[:name])")
         mat = dropzeros(original_mat)
         bg = bipartite_graph(mat)
@@ -51,14 +51,14 @@ end;
 Comparison with Tables 3.1 and 3.2 of "What color is your Jacobian?"
 =#
 
-what_table_31_32_33 = CSV.read(
-    joinpath(@__DIR__, "reference", "what_table_31_32_33.csv"), DataFrame
+what_table_31_32 = CSV.read(
+    joinpath(@__DIR__, "reference", "what_table_31_32.csv"), DataFrame
 )
 
 @testset "Distance-2 coloring (survey paper)" begin
-    @testset "$(row[:name])" for row in eachrow(what_table_31_32_33)
+    @testset "$(row[:name])" for row in eachrow(what_table_31_32)
         ismissing(row[:group]) && continue
-        @info "Testing distance-2 coloring for $(row[:name])"
+        @info "Testing distance-2 coloring for $(row[:name]) against survey paper"
         original_mat = matrixdepot("$(row[:group])/$(row[:name])")
         mat = original_mat  # no dropzeros
         bg = bipartite_graph(mat)
@@ -91,7 +91,7 @@ what_table_41_42 = CSV.read(
 @testset "Star coloring (survey paper)" begin
     @testset "$(row[:name])" for row in eachrow(what_table_41_42)
         ismissing(row[:group]) && continue
-        @info "Testing star coloring for $(row[:name])"
+        @info "Testing star coloring for $(row[:name]) against survey paper"
         original_mat = matrixdepot("$(row[:group])/$(row[:name])")
         mat = original_mat  # no dropzeros
         bg = bipartite_graph(mat)

From c029befb789e98c031b7f830e7df904ffa300314 Mon Sep 17 00:00:00 2001
From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com>
Date: Tue, 21 May 2024 15:44:03 +0200
Subject: [PATCH 5/8] CSV

---
 test/reference/what_table_31_32.csv | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/reference/what_table_31_32.csv b/test/reference/what_table_31_32.csv
index e518d9d9..733a3343 100644
--- a/test/reference/what_table_31_32.csv
+++ b/test/reference/what_table_31_32.csv
@@ -19,5 +19,5 @@ row,modified,group,name,m,n,nnz,κmax,κmin,ρmax,ρmin,K
 18,false,Mallya,lhr34,35152,35152,764014,36,1,63,1,65
 19,["m"],Mallya,lhr71c,70304,70304,1528092,36,1,63,1,65
 20,false,vanHeukelum,cage12,130228,130228,2032536,33,5,33,5,96
-21,false,LPnetlib,lp_osa_07,1118,25067,144812,6,1,17613,18,17613,
-22,false,LPnetlib,lp_fit2d,25,10524,129042,17,1,10500,1427,10501,
\ No newline at end of file
+21,false,LPnetlib,lp_osa_07,1118,25067,144812,6,1,17613,18,17613
+22,false,LPnetlib,lp_fit2d,25,10524,129042,17,1,10500,1427,10501
\ No newline at end of file

From 80d24b8e5b619449088d7c3a029f1ccdd2aa679a Mon Sep 17 00:00:00 2001
From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com>
Date: Tue, 21 May 2024 16:23:03 +0200
Subject: [PATCH 6/8] Add performance tests

---
 src/adtypes.jl      |  4 ++++
 src/coloring.jl     | 42 +++++++++++++++++++++++++++++++++---------
 src/graph.jl        |  2 +-
 test/Project.toml   |  1 +
 test/performance.jl | 43 +++++++++++++++++++++++++++++++++++++++++++
 test/runtests.jl    | 43 ++++++++++++++++++++++++++-----------------
 6 files changed, 108 insertions(+), 27 deletions(-)
 create mode 100644 test/performance.jl

diff --git a/src/adtypes.jl b/src/adtypes.jl
index e41f6a79..4469dfba 100644
--- a/src/adtypes.jl
+++ b/src/adtypes.jl
@@ -20,6 +20,10 @@ end
 
 GreedyColoringAlgorithm() = GreedyColoringAlgorithm(NaturalOrder())
 
+function Base.show(io::IO, algo::GreedyColoringAlgorithm)
+    return print(io, "GreedyColoringAlgorithm($(algo.order))")
+end
+
 function ADTypes.column_coloring(A::AbstractMatrix, algo::GreedyColoringAlgorithm)
     bg = bipartite_graph(A)
     return partial_distance2_coloring(bg, Val(2), algo.order)
diff --git a/src/coloring.jl b/src/coloring.jl
index dfd4de70..1c68b518 100644
--- a/src/coloring.jl
+++ b/src/coloring.jl
@@ -15,10 +15,24 @@ The vertices are colored in a greedy fashion, following the `order` supplied.
 function partial_distance2_coloring(
     bg::BipartiteGraph, ::Val{side}, order::AbstractOrder
 ) where {side}
+    colors = Vector{Int}(undef, length(bg, Val(side)))
+    forbidden_colors = Vector{Int}(undef, length(bg, Val(side)))
+    vertices_in_order = vertices(bg, Val(side), order)
+    partial_distance2_coloring!(colors, forbidden_colors, bg, Val(side), vertices_in_order)
+    return colors
+end
+
+function partial_distance2_coloring!(
+    colors::Vector{Int},
+    forbidden_colors::Vector{Int},
+    bg::BipartiteGraph,
+    ::Val{side},
+    vertices_in_order::AbstractVector{<:Integer},
+) where {side}
+    colors .= 0
+    forbidden_colors .= 0
     other_side = 3 - side
-    colors = zeros(Int, length(bg, Val(side)))
-    forbidden_colors = zeros(Int, length(bg, Val(side)))
-    for v in vertices(bg, Val(side), order)
+    for v in vertices_in_order
         for w in neighbors(bg, Val(side), v)
             for x in neighbors(bg, Val(other_side), w)
                 if !iszero(colors[x])
@@ -33,7 +47,6 @@ function partial_distance2_coloring(
             end
         end
     end
-    return colors
 end
 
 """
@@ -51,10 +64,22 @@ The vertices are colored in a greedy fashion, following the `order` supplied.
 - [`AbstractOrder`](@ref)
 """
 function star_coloring1(g::Graph, order::AbstractOrder)
-    n = length(g)
-    colors = zeros(Int, n)
-    forbidden_colors = zeros(Int, n)
-    for v in vertices(g, order)
+    colors = Vector{Int}(undef, length(g))
+    forbidden_colors = Vector{Int}(undef, length(g))
+    vertices_in_order = vertices(g, order)
+    star_coloring1!(colors, forbidden_colors, g, vertices_in_order)
+    return colors
+end
+
+function star_coloring1!(
+    colors::Vector{Int},
+    forbidden_colors::Vector{Int},
+    g::Graph,
+    vertices_in_order::AbstractVector{<:Integer},
+)
+    colors .= 0
+    forbidden_colors .= 0
+    for v in vertices_in_order
         for w in neighbors(g, v)
             if !iszero(colors[w])  # w is colored
                 forbidden_colors[colors[w]] = v
@@ -81,5 +106,4 @@ function star_coloring1(g::Graph, order::AbstractOrder)
             end
         end
     end
-    return colors
 end
diff --git a/src/graph.jl b/src/graph.jl
index e62fde1e..de3fb725 100644
--- a/src/graph.jl
+++ b/src/graph.jl
@@ -104,7 +104,7 @@ The bipartite graph of a matrix `A ∈ ℝ^{m × n}` is `Gb(A) = (V₁, V₂, E)
 - `(i, j) ∈ E` whenever `A[i, j] ≠ 0`
 """
 function bipartite_graph(J::SparseMatrixCSC)
-    g1 = Graph(SparseMatrixCSC(transpose(J)))  # rows to columns
+    g1 = Graph(transpose(J))  # rows to columns
     g2 = Graph(J)  # columns to rows
     return BipartiteGraph(g1, g2)
 end
diff --git a/test/Project.toml b/test/Project.toml
index f71c8c9e..261cbc8a 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -2,6 +2,7 @@
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+Chairmarks = "0ca39b1e-fe0b-4e98-acfc-b1656634c4de"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
diff --git a/test/performance.jl b/test/performance.jl
new file mode 100644
index 00000000..132c62aa
--- /dev/null
+++ b/test/performance.jl
@@ -0,0 +1,43 @@
+using ADTypes: column_coloring, row_coloring, symmetric_coloring
+using Chairmarks
+using JET
+using LinearAlgebra
+using SparseArrays
+using SparseMatrixColorings:
+    adjacency_graph, bipartite_graph, partial_distance2_coloring!, star_coloring1!
+using StableRNGs
+using Test
+
+rng = StableRNG(63)
+
+@testset "Type stability" begin
+    n = 10
+    A = sprand(rng, Bool, n, n, 3 / n)
+    algo = GreedyColoringAlgorithm()
+    @test_opt target_modules = (SparseMatrixColorings,) column_coloring(A, algo)
+    @test_opt target_modules = (SparseMatrixColorings,) row_coloring(A, algo)
+    @test_opt target_modules = (SparseMatrixColorings,) symmetric_coloring(
+        A + transpose(A), algo
+    )
+end
+
+function benchmark_distance2_coloring(n)
+    @be (;
+        bg=bipartite_graph(sprand(Bool, n, n, 3 / n)),
+        colors=Vector{Int}(undef, n),
+        forbidden_colors=Vector{Int}(undef, n),
+    ) partial_distance2_coloring!(_.colors, _.forbidden_colors, _.bg, Val(1), 1:n) evals = 1
+end
+
+function benchmark_star_coloring(n)
+    @be (;
+        g=adjacency_graph(Symmetric(sprand(Bool, n, n, 3 / n))),
+        colors=Vector{Int}(undef, n),
+        forbidden_colors=Vector{Int}(undef, n),
+    ) star_coloring1!(_.colors, _.forbidden_colors, _.g, 1:n) evals = 1
+end
+
+@testset "Allocations" begin
+    @test minimum(benchmark_distance2_coloring(10)).allocs == 0
+    @test minimum(benchmark_star_coloring(10)).allocs == 0
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index d3041c89..b236f025 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -6,8 +6,8 @@ using SparseMatrixColorings
 using Test
 
 @testset verbose = true "SparseMatrixColorings" begin
-    if VERSION >= v"1.10"
-        @testset verbose = true "Code quality" begin
+    @testset verbose = true "Code quality" begin
+        if VERSION >= v"1.10"
             @testset "Aqua" begin
                 Aqua.test_all(SparseMatrixColorings)
             end
@@ -20,23 +20,32 @@ using Test
                 )
             end
         end
+        @testset "Doctests" begin
+            Documenter.doctest(SparseMatrixColorings)
+        end
     end
-    @testset "Doctests" begin
-        Documenter.doctest(SparseMatrixColorings)
-    end
-    @testset "Graph" begin
-        include("graph.jl")
-    end
-    @testset "Order" begin
-        include("order.jl")
-    end
-    @testset "Check" begin
-        include("check.jl")
+    @testset verbose = true "Internals" begin
+        @testset "Graph" begin
+            include("graph.jl")
+        end
+        @testset "Order" begin
+            include("order.jl")
+        end
+        @testset "Check" begin
+            include("check.jl")
+        end
     end
-    @testset "ADTypes" begin
-        include("adtypes.jl")
+    @testset verbose = true "Correctness" begin
+        @testset "ADTypes" begin
+            include("adtypes.jl")
+        end
+        @testset "SuiteSparse" begin
+            include("suitesparse.jl")
+        end
     end
-    @testset "SuiteSparse" begin
-        include("suitesparse.jl")
+    @testset "Performance" begin
+        if VERSION >= v"1.10"
+            include("performance.jl")
+        end
     end
 end

From 8800d8a8294fb6d2cb6fc47169ffbfd362984176 Mon Sep 17 00:00:00 2001
From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com>
Date: Tue, 21 May 2024 16:24:55 +0200
Subject: [PATCH 7/8] Correct column intersection graph

---
 test/graph.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/test/graph.jl b/test/graph.jl
index 7c67ca42..29f59ab8 100644
--- a/test/graph.jl
+++ b/test/graph.jl
@@ -56,7 +56,8 @@ end;
         0 0 0 1 1 1 1 0
     ])
 
-    g = adjacency_graph(transpose(A) * A)
+    B = transpose(A) * A
+    g = adjacency_graph(B - Diagonal(B))
     @test length(g) == 8
     @test neighbors(g, 1) == [6, 7, 8]
     @test neighbors(g, 2) == [5, 7, 8]

From 04137178dcd977b40105ec894549497778fde8af Mon Sep 17 00:00:00 2001
From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com>
Date: Tue, 21 May 2024 16:26:09 +0200
Subject: [PATCH 8/8] Missing using

---
 test/graph.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/graph.jl b/test/graph.jl
index 29f59ab8..6a46984b 100644
--- a/test/graph.jl
+++ b/test/graph.jl
@@ -1,3 +1,4 @@
+using LinearAlgebra
 using SparseArrays
 using SparseMatrixColorings: Graph, adjacency_graph, bipartite_graph, neighbors
 using Test