diff --git a/README.md b/README.md
index 5965509a..9056936b 100644
--- a/README.md
+++ b/README.md
@@ -18,11 +18,12 @@ pkg> add https://github.com/gdalle/SparseMatrixColorings.jl
 
 ## Background
 
-The algorithms implemented in this package are mainly taken from the following articles:
+The algorithms implemented in this package are 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)
+- [_What Color Is Your Jacobian? Graph Coloring for Computing Derivatives_](https://epubs.siam.org/doi/10.1137/S0036144504444711), Gebremedhin et al. (2005)
+- [_New Acyclic and Star Coloring Algorithms with Application to Computing Hessians_](https://epubs.siam.org/doi/abs/10.1137/050639879), Gebremedhin et al. (2007)
+- [_Efficient Computation of Sparse Hessians Using Coloring and Automatic Differentiation_](https://pubsonline.informs.org/doi/abs/10.1287/ijoc.1080.0286), Gebremedhin et al. (2009)
+- [_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/docs/src/api.md b/docs/src/api.md
index 51be2bc5..c9734413 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -10,6 +10,9 @@ CurrentModule = SparseMatrixColorings
 ```@docs
 SparseMatrixColorings
 GreedyColoringAlgorithm
+column_coloring
+row_coloring
+symmetric_coloring
 ```
 
 ## Public, not exported
@@ -60,13 +63,12 @@ vertices
 
 ```@docs
 partial_distance2_coloring
-star_coloring1
+star_coloring
 ```
 
 ### Testing
 
 ```@docs
 check_structurally_orthogonal_columns
-check_structurally_orthogonal_rows
-check_symmetrically_orthogonal
+check_symmetrically_orthogonal_columns
 ```
diff --git a/src/SparseMatrixColorings.jl b/src/SparseMatrixColorings.jl
index 3d0bb4d2..70c22c79 100644
--- a/src/SparseMatrixColorings.jl
+++ b/src/SparseMatrixColorings.jl
@@ -5,7 +5,8 @@ $README
 """
 module SparseMatrixColorings
 
-using ADTypes: ADTypes, AbstractColoringAlgorithm
+using ADTypes:
+    ADTypes, AbstractColoringAlgorithm, column_coloring, row_coloring, symmetric_coloring
 using Compat: @compat
 using DocStringExtensions: README
 using LinearAlgebra:
@@ -46,7 +47,9 @@ include("check.jl")
 @compat public decompress_columns, decompress_columns!
 @compat public decompress_rows, decompress_rows!
 @compat public decompress_symmetric, decompress_symmetric!
+@compat public column_coloring, row_coloring, symmetric_coloring
 
 export GreedyColoringAlgorithm
+export column_coloring, row_coloring, symmetric_coloring
 
 end
diff --git a/src/adtypes.jl b/src/adtypes.jl
index 49fd55ff..a1f949f9 100644
--- a/src/adtypes.jl
+++ b/src/adtypes.jl
@@ -9,24 +9,42 @@ Compatible with the [ADTypes.jl coloring framework](https://sciml.github.io/ADTy
 
     GreedyColoringAlgorithm(order::AbstractOrder=NaturalOrder())
 
-# Implements
+# See also
+
+- [`AbstractOrder`](@ref)
+"""
+struct GreedyColoringAlgorithm{O<:AbstractOrder} <: ADTypes.AbstractColoringAlgorithm
+    order::O
+end
+
+GreedyColoringAlgorithm() = GreedyColoringAlgorithm(NaturalOrder())
+
+function Base.show(io::IO, algo::GreedyColoringAlgorithm)
+    return print(io, "GreedyColoringAlgorithm($(algo.order))")
+end
+
+"""
+    column_coloring(A::AbstractMatrix, algo::GreedyColoringAlgorithm)
+
+Compute a partial distance-2 coloring of the columns in the bipartite graph of the matrix `A`.
 
-- [`ADTypes.column_coloring`](@extref ADTypes) and [`ADTypes.row_coloring`](@extref ADTypes) with a partial distance-2 coloring of the bipartite graph
-- [`ADTypes.symmetric_coloring`](@extref ADTypes) with a star coloring of the adjacency graph
+Function defined by ADTypes, re-exported by SparseMatrixColorings.
 
-# Example use
+# Example
 
 ```jldoctest
-using ADTypes, SparseMatrixColorings, SparseArrays
+using SparseMatrixColorings, SparseArrays
 
 algo = GreedyColoringAlgorithm(SparseMatrixColorings.LargestFirst())
+
 A = sparse([
     0 0 1 1 0
     1 0 0 0 1
     0 1 1 0 0
     0 1 1 0 1
 ])
-ADTypes.column_coloring(A, algo)
+
+column_coloring(A, algo)
 
 # output
 
@@ -37,32 +55,62 @@ ADTypes.column_coloring(A, algo)
  2
  3
 ```
-
-# See also
-
-- [`AbstractOrder`](@ref)
 """
-struct GreedyColoringAlgorithm{O<:AbstractOrder} <: ADTypes.AbstractColoringAlgorithm
-    order::O
-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)
 end
 
+"""
+    row_coloring(A::AbstractMatrix, algo::GreedyColoringAlgorithm)
+
+Compute a partial distance-2 coloring of the rows in the bipartite graph of the matrix `A`.
+
+Function defined by ADTypes, re-exported by SparseMatrixColorings.
+
+# Example
+
+```jldoctest
+using SparseMatrixColorings, SparseArrays
+
+algo = GreedyColoringAlgorithm(SparseMatrixColorings.LargestFirst())
+
+A = sparse([
+    0 0 1 1 0
+    1 0 0 0 1
+    0 1 1 0 0
+    0 1 1 0 1
+])
+
+row_coloring(A, algo)
+
+# output
+
+4-element Vector{Int64}:
+ 2
+ 2
+ 3
+ 1
+```
+"""
 function ADTypes.row_coloring(A::AbstractMatrix, algo::GreedyColoringAlgorithm)
     bg = bipartite_graph(A)
     return partial_distance2_coloring(bg, Val(1), algo.order)
 end
 
+"""
+    symmetric_coloring(A::AbstractMatrix, algo::GreedyColoringAlgorithm)
+
+Compute a star coloring of the columns in the adjacency graph of the symmetric matrix `A`.
+
+Function defined by ADTypes, re-exported by SparseMatrixColorings.
+
+# Example
+
+!!! warning
+    Work in progress.
+"""
 function ADTypes.symmetric_coloring(A::AbstractMatrix, algo::GreedyColoringAlgorithm)
     ag = adjacency_graph(A)
-    return star_coloring1(ag, algo.order)
+    return star_coloring(ag, algo.order)
 end
diff --git a/src/check.jl b/src/check.jl
index 0cb6d34a..8e2516f1 100644
--- a/src/check.jl
+++ b/src/check.jl
@@ -1,55 +1,39 @@
 """
     check_structurally_orthogonal_columns(
-        A::AbstractMatrix, colors::AbstractVector{<:Integer}
+        A::AbstractMatrix, color::AbstractVector{<:Integer}
         verbose=false
     )
 
-Return `true` if coloring the columns of the matrix `A` with the vector `colors` results in a partition that is structurally orthogonal, and `false` otherwise.
+Return `true` if coloring the columns of the matrix `A` with the vector `color` results in a partition that is structurally orthogonal, and `false` otherwise.
     
 A partition of the columns of a matrix `A` is _structurally orthogonal_ if, for every nonzero element `A[i, j]`, the group containing column `A[:, j]` has no other column with a nonzero in row `i`.
 
 !!! warning
     This function is not coded with efficiency in mind, it is designed for small-scale tests.
+
+# References
+
+> [_What Color Is Your Jacobian? Graph Coloring for Computing Derivatives_](https://epubs.siam.org/doi/10.1137/S0036144504444711), Gebremedhin et al. (2005)
 """
 function check_structurally_orthogonal_columns(
-    A::AbstractMatrix, colors::AbstractVector{<:Integer}; verbose::Bool=true
+    A::AbstractMatrix, color::AbstractVector{<:Integer}; verbose::Bool=false
 )
-    groups = color_groups(colors)
-    for (c, g) in enumerate(groups)
-        Ag = @view A[:, g]
+    if length(color) != size(A, 2)
+        if verbose
+            @warn "$(length(color)) colors provided for $(size(A, 2)) columns"
+        end
+        return false
+    end
+    group = color_groups(color)
+    for (c, g) in enumerate(group)
+        Ag = view(A, :, g)
         nonzeros_per_row = dropdims(count(!iszero, Ag; dims=2); dims=2)
         max_nonzeros_per_row, i = findmax(nonzeros_per_row)
         if max_nonzeros_per_row > 1
-            verbose && @warn "Columns $g (with color $c) share nonzeros in row $i"
-            return false
-        end
-    end
-    return true
-end
-
-"""
-    check_structurally_orthogonal_rows(
-        A::AbstractMatrix, colors::AbstractVector{<:Integer};
-        verbose=false
-    )
-
-Return `true` if coloring the rows of the matrix `A` with the vector `colors` results in a partition that is structurally orthogonal, and `false` otherwise.
-    
-A partition of the rows of a matrix `A` is _structurally orthogonal_ if, for every nonzero element `A[i, j]`, the group containing row `A[i, :]` has no other row with a nonzero in column `j`.
-
-!!! warning
-    This function is not coded with efficiency in mind, it is designed for small-scale tests.
-"""
-function check_structurally_orthogonal_rows(
-    A::AbstractMatrix, colors::AbstractVector{<:Integer}; verbose::Bool=true
-)
-    groups = color_groups(colors)
-    for (c, g) in enumerate(groups)
-        Ag = @view A[g, :]
-        nonzeros_per_col = dropdims(count(!iszero, Ag; dims=1); dims=1)
-        max_nonzeros_per_col, j = findmax(nonzeros_per_col)
-        if max_nonzeros_per_col > 1
-            verbose && @warn "Rows $g (with color $c) share nonzeros in column $j"
+            if verbose
+                incompatible_columns = g[findall(!iszero, view(Ag, i, :))]
+                @warn "In color $c, columns $incompatible_columns all have nonzeros in row $i"
+            end
             return false
         end
     end
@@ -57,41 +41,55 @@ function check_structurally_orthogonal_rows(
 end
 
 """
-    check_symmetrically_orthogonal(
-        A::AbstractMatrix, colors::AbstractVector{<:Integer};
+    check_symmetrically_orthogonal_columns(
+        A::AbstractMatrix, color::AbstractVector{<:Integer};
         verbose=false
     )
 
-Return `true` if coloring the columns of the symmetric matrix `A` with the vector `colors` results in a partition that is symmetrically orthogonal, and `false` otherwise.
+Return `true` if coloring the columns of the symmetric matrix `A` with the vector `color` results in a partition that is symmetrically orthogonal, and `false` otherwise.
     
-A partition of the columns of a symmetrix matrix `A` is _symmetrically orthogonal_ if, for every nonzero element `A[i, j]`, either
+A partition of the columns of a symmetrix matrix `A` is _symmetrically orthogonal_ if, for every nonzero element `A[i, j]`, either of the following statements holds:
 
 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`
 
 !!! warning
     This function is not coded with efficiency in mind, it is designed for small-scale tests.
+
+# References
+
+> [_What Color Is Your Jacobian? Graph Coloring for Computing Derivatives_](https://epubs.siam.org/doi/10.1137/S0036144504444711), Gebremedhin et al. (2005)
 """
-function check_symmetrically_orthogonal(
-    A::AbstractMatrix, colors::AbstractVector{<:Integer}; verbose::Bool=true
+function check_symmetrically_orthogonal_columns(
+    A::AbstractMatrix, color::AbstractVector{<:Integer}; verbose::Bool=false
 )
     checksquare(A)
+    if length(color) != size(A, 2)
+        if verbose
+            @warn "$(length(color)) colors provided for $(size(A, 2)) columns"
+        end
+        return false
+    end
     issymmetric(A) || return false
-    groups = color_groups(colors)
+    group = color_groups(color)
     for i in axes(A, 2), j in axes(A, 2)
         iszero(A[i, j]) && continue
-        ki, kj = colors[i], colors[j]
-        gi, gj = groups[ki], groups[kj]
+        ci, cj = color[i], color[j]
+        gi, gj = group[ci], group[cj]
         A_gj_rowi = view(A, i, gj)
         A_gi_rowj = view(A, j, gi)
         nonzeros_gj_rowi = count(!iszero, A_gj_rowi)
         nonzeros_gi_rowj = count(!iszero, A_gi_rowj)
         if nonzeros_gj_rowi > 1 && nonzeros_gi_rowj > 1
-            verbose && @warn """
-            For coefficient $((i, j)):
-            - columns $gj (with color $kj) share nonzeros in row $i
-            - columns $gi (with color $ki) share nonzeros in row $j
-            """
+            if verbose
+                gj_incompatible_columns = gj[findall(!iszero, A_gj_rowi)]
+                gi_incompatible_columns = gi[findall(!iszero, A_gi_rowj)]
+                @warn """
+                For coefficient (i=$i, j=$j) with column colors (ci=$ci, cj=$cj):
+                - in color ci=$ci, columns $gi_incompatible_columns all have nonzeros in row j=$j
+                - in color cj=$cj, columns $gj_incompatible_columns all have nonzeros in row i=$i
+                """
+            end
             return false
         end
     end
diff --git a/src/coloring.jl b/src/coloring.jl
index 1c68b518..fad778f8 100644
--- a/src/coloring.jl
+++ b/src/coloring.jl
@@ -11,38 +11,42 @@ The vertices are colored in a greedy fashion, following the `order` supplied.
 
 - [`BipartiteGraph`](@ref)
 - [`AbstractOrder`](@ref)
+
+# References
+
+> [_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, ::Val{side}, order::AbstractOrder
 ) where {side}
-    colors = Vector{Int}(undef, length(bg, Val(side)))
+    color = 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
+    partial_distance2_coloring!(color, forbidden_colors, bg, Val(side), vertices_in_order)
+    return color
 end
 
 function partial_distance2_coloring!(
-    colors::Vector{Int},
+    color::Vector{Int},
     forbidden_colors::Vector{Int},
     bg::BipartiteGraph,
     ::Val{side},
     vertices_in_order::AbstractVector{<:Integer},
 ) where {side}
-    colors .= 0
+    color .= 0
     forbidden_colors .= 0
     other_side = 3 - side
     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])
-                    forbidden_colors[colors[x]] = v
+                if !iszero(color[x])
+                    forbidden_colors[color[x]] = v
                 end
             end
         end
-        for c in eachindex(forbidden_colors)
-            if forbidden_colors[c] != v
-                colors[v] = c
+        for i in eachindex(forbidden_colors)
+            if forbidden_colors[i] != v
+                color[v] = i
                 break
             end
         end
@@ -50,7 +54,7 @@ function partial_distance2_coloring!(
 end
 
 """
-    star_coloring1(g::Graph, order::AbstractOrder)
+    star_coloring(g::Graph, order::AbstractOrder)
 
 Compute a star coloring of all vertices in the adjacency graph `g` and return a vector of integer colors.
 
@@ -62,48 +66,109 @@ The vertices are colored in a greedy fashion, following the `order` supplied.
 
 - [`Graph`](@ref)
 - [`AbstractOrder`](@ref)
+
+# References
+
+> [_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_coloring1(g::Graph, order::AbstractOrder)
-    colors = Vector{Int}(undef, length(g))
-    forbidden_colors = Vector{Int}(undef, length(g))
+function star_coloring(g::Graph, order::AbstractOrder)
+    # Initialize data structures
+    color = zeros(Int, length(g))
+    forbidden_colors = zeros(Int, length(g))
+    first_neighbor = Vector{Tuple{Int,Int}}(undef, length(g))
+    treated = zeros(Int, length(g))
+    star = Dict{Tuple{Int,Int},Int}()
+    hub = Int[]
     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
-            end
-            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(g, x)
-                        if !iszero(colors[y]) && y != w
-                            if colors[y] == colors[w]
-                                forbidden_colors[colors[x]] = v
-                                break
-                            end
-                        end
+            iszero(color[w]) && continue
+            forbidden_colors[color[w]] = v
+            (p, q) = first_neighbor[color[w]]
+            if p == v  # Case 1
+                if treated[q] != v
+                    # forbid colors of neighbors of q
+                    _treat!(treated, forbidden_colors, g, v, q, color)
+                end
+                # forbid colors of neighbors of w
+                _treat!(treated, forbidden_colors, g, v, w, color)
+            else
+                first_neighbor[color[w]] = (v, w)
+                for x in neighbors(g, w)
+                    (x == v || iszero(color[x])) && continue
+                    wx = _sort(w, x)
+                    if x == hub[star[wx]]  # potential Case 2
+                        forbidden_colors[color[x]] = v
                     end
                 end
             end
         end
-        for c in eachindex(forbidden_colors)
-            if forbidden_colors[c] != v
-                colors[v] = c
+        for i in eachindex(forbidden_colors)
+            if forbidden_colors[i] != v
+                color[v] = i
+                break
+            end
+        end
+        _update_stars!(star, hub, g, v, color, first_neighbor)
+    end
+    return color
+end
+
+_sort(u, v) = (min(u, v), max(u, v))
+
+function _treat!(
+    # modified
+    treated::AbstractVector{<:Integer},
+    forbidden_colors::AbstractVector{<:Integer},
+    # not modified
+    g::Graph,
+    v::Integer,
+    w::Integer,
+    color::AbstractVector{<:Integer},
+)
+    for x in neighbors(g, w)
+        iszero(color[x]) && continue
+        forbidden_colors[color[x]] = v
+    end
+    treated[w] = v
+    return nothing
+end
+
+function _update_stars!(
+    # modified
+    star,
+    hub::AbstractVector{<:Integer},
+    # not modified
+    g::Graph,
+    v::Integer,
+    color::AbstractVector{<:Integer},
+    first_neighbor::AbstractVector{<:Tuple{<:Integer,<:Integer}},
+)
+    for w in neighbors(g, v)
+        iszero(color[w]) && continue
+        vw = _sort(v, w)
+        x_exists = false
+        for x in neighbors(g, w)
+            if x != v && color[x] == color[v]  # vw, wx ∈ E
+                wx = _sort(w, x)
+                hub[star[wx]] = w  # this may already be true
+                star[vw] = star[wx]
+                x_exists = true
                 break
             end
         end
+        if !x_exists
+            (p, q) = first_neighbor[color[w]]
+            if p == v && q != w  # vw, vq ∈ E and color[w] = color[q]
+                vq = _sort(v, q)
+                hub[star[vq]] = v  # this may already be true
+                star[vw] = star[vq]
+            else  # vw forms a new star
+                push!(hub, 0)  # hub is yet undefined
+                star[vw] = length(hub)
+            end
+        end
     end
+    return nothing
 end
diff --git a/src/decompression.jl b/src/decompression.jl
index 4a4716f1..d8db9275 100644
--- a/src/decompression.jl
+++ b/src/decompression.jl
@@ -4,32 +4,32 @@
     decompress_columns!(
         A::AbstractMatrix{R},
         S::AbstractMatrix{Bool},
-        C::AbstractMatrix{R},
-        colors::AbstractVector{<:Integer}
+        B::AbstractMatrix{R},
+        color::AbstractVector{<:Integer}
     ) where {R<:Real}
 
-Decompress the thin matrix `C` into the fat matrix `A` which must have the same sparsity pattern as `S`.
+Decompress the thin matrix `B` into the fat matrix `A` which must have the same sparsity pattern as `S`.
 
-Here, `colors` is a column coloring of `S`, while `C` is a compressed representation of matrix `A` obtained by summing the columns that share the same color.
+Here, `color` is a column coloring of `S`, while `B` is a compressed representation of matrix `A` obtained by summing the columns that share the same color.
 """
 function decompress_columns! end
 
 function decompress_columns!(
     A::AbstractMatrix{R},
     S::AbstractMatrix{Bool},
-    C::AbstractMatrix{R},
-    colors::AbstractVector{<:Integer},
+    B::AbstractMatrix{R},
+    color::AbstractVector{<:Integer},
 ) where {R<:Real}
     if !same_sparsity_pattern(A, S)
         throw(DimensionMismatch("`A` and `S` must have the same sparsity pattern."))
     end
     A .= zero(R)
     for j in axes(A, 2)
-        k = colors[j]
+        cj = color[j]
         rows_j = (!iszero).(view(S, :, j))
         Aj = view(A, rows_j, j)
-        Cj = view(C, rows_j, k)
-        copyto!(Aj, Cj)
+        Bj = view(B, rows_j, cj)
+        copyto!(Aj, Bj)
     end
     return A
 end
@@ -37,8 +37,8 @@ end
 function decompress_columns!(
     A::SparseMatrixCSC{R},
     S::SparseMatrixCSC{Bool},
-    C::AbstractMatrix{R},
-    colors::AbstractVector{<:Integer},
+    B::AbstractMatrix{R},
+    color::AbstractVector{<:Integer},
 ) where {R<:Real}
     if !same_sparsity_pattern(A, S)
         throw(DimensionMismatch("`A` and `S` must have the same sparsity pattern."))
@@ -46,12 +46,12 @@ function decompress_columns!(
     Anz, Arv = nonzeros(A), rowvals(A)
     Anz .= zero(R)
     for j in axes(A, 2)
-        k = colors[j]
+        cj = color[j]
         nzrange_j = nzrange(A, j)
         rows_j = view(Arv, nzrange_j)
         Aj = view(Anz, nzrange_j)
-        Cj = view(C, rows_j, k)
-        copyto!(Aj, Cj)
+        Bj = view(B, rows_j, cj)
+        copyto!(Aj, Bj)
     end
     return A
 end
@@ -59,19 +59,19 @@ end
 """
     decompress_columns(
         S::AbstractMatrix{Bool},
-        C::AbstractMatrix{R},
-        colors::AbstractVector{<:Integer}
+        B::AbstractMatrix{R},
+        color::AbstractVector{<:Integer}
     ) where {R<:Real}
 
-Decompress the thin matrix `C` into a new fat matrix `A` with the same sparsity pattern as `S`.
+Decompress the thin matrix `B` into a new fat matrix `A` with the same sparsity pattern as `S`.
 
-Here, `colors` is a column coloring of `S`, while `C` is a compressed representation of matrix `A` obtained by summing the columns that share the same color.
+Here, `color` is a column coloring of `S`, while `B` is a compressed representation of matrix `A` obtained by summing the columns that share the same color.
 """
 function decompress_columns(
-    S::AbstractMatrix{Bool}, C::AbstractMatrix{R}, colors::AbstractVector{<:Integer}
+    S::AbstractMatrix{Bool}, B::AbstractMatrix{R}, color::AbstractVector{<:Integer}
 ) where {R<:Real}
     A = respectful_similar(S, R)
-    return decompress_columns!(A, S, C, colors)
+    return decompress_columns!(A, S, B, color)
 end
 
 ## Row decompression
@@ -80,32 +80,32 @@ end
     decompress_rows!(
         A::AbstractMatrix{R},
         S::AbstractMatrix{Bool},
-        C::AbstractMatrix{R},
-        colors::AbstractVector{<:Integer}
+        B::AbstractMatrix{R},
+        color::AbstractVector{<:Integer}
     ) where {R<:Real}
 
-Decompress the small matrix `C` into the tall matrix `A` which must have the same sparsity pattern as `S`.
+Decompress the small matrix `B` into the tall matrix `A` which must have the same sparsity pattern as `S`.
 
-Here, `colors` is a row coloring of `S`, while `C` is a compressed representation of matrix `A` obtained by summing the columns that share the same color.
+Here, `color` is a row coloring of `S`, while `B` is a compressed representation of matrix `A` obtained by summing the columns that share the same color.
 """
 function decompress_rows! end
 
 function decompress_rows!(
     A::AbstractMatrix{R},
     S::AbstractMatrix{Bool},
-    C::AbstractMatrix{R},
-    colors::AbstractVector{<:Integer},
+    B::AbstractMatrix{R},
+    color::AbstractVector{<:Integer},
 ) where {R<:Real}
     if !same_sparsity_pattern(A, S)
         throw(DimensionMismatch("`A` and `S` must have the same sparsity pattern."))
     end
     A .= zero(R)
     for i in axes(A, 1)
-        k = colors[i]
+        ci = color[i]
         cols_i = (!iszero).(view(S, i, :))
         Ai = view(A, i, cols_i)
-        Ci = view(C, k, cols_i)
-        copyto!(Ai, Ci)
+        Bi = view(B, ci, cols_i)
+        copyto!(Ai, Bi)
     end
     return A
 end
@@ -113,8 +113,8 @@ end
 function decompress_rows!(
     A::TransposeOrAdjoint{R,<:SparseMatrixCSC{R}},
     S::TransposeOrAdjoint{Bool,<:SparseMatrixCSC{Bool}},
-    C::AbstractMatrix{R},
-    colors::AbstractVector{<:Integer},
+    B::AbstractMatrix{R},
+    color::AbstractVector{<:Integer},
 ) where {R<:Real}
     if !same_sparsity_pattern(A, S)
         throw(DimensionMismatch("`A` and `S` must have the same sparsity pattern."))
@@ -123,12 +123,12 @@ function decompress_rows!(
     PAnz, PArv = nonzeros(PA), rowvals(PA)
     PAnz .= zero(R)
     for i in axes(A, 1)
-        k = colors[i]
+        ci = color[i]
         nzrange_i = nzrange(PA, i)
         cols_i = view(PArv, nzrange_i)
         Ai = view(PAnz, nzrange_i)
-        Ci = view(C, k, cols_i)
-        copyto!(Ai, Ci)
+        Bi = view(B, ci, cols_i)
+        copyto!(Ai, Bi)
     end
     return A
 end
@@ -136,19 +136,19 @@ end
 """
     decompress_rows(
         S::AbstractMatrix{Bool},
-        C::AbstractMatrix{R},
-        colors::AbstractVector{<:Integer}
+        B::AbstractMatrix{R},
+        color::AbstractVector{<:Integer}
     ) where {R<:Real}
 
-Decompress the small matrix `C` into a new tall matrix `A` with the same sparsity pattern as `S`.
+Decompress the small matrix `B` into a new tall matrix `A` with the same sparsity pattern as `S`.
 
-Here, `colors` is a row coloring of `S`, while `C` is a compressed representation of matrix `A` obtained by summing the columns that share the same color.
+Here, `color` is a row coloring of `S`, while `B` is a compressed representation of matrix `A` obtained by summing the columns that share the same color.
 """
 function decompress_rows(
-    S::AbstractMatrix{Bool}, C::AbstractMatrix{R}, colors::AbstractVector{<:Integer}
+    S::AbstractMatrix{Bool}, B::AbstractMatrix{R}, color::AbstractVector{<:Integer}
 ) where {R<:Real}
     A = respectful_similar(S, R)
-    return decompress_rows!(A, S, C, colors)
+    return decompress_rows!(A, S, B, color)
 end
 
 ## Symmetric decompression
@@ -157,35 +157,45 @@ end
     decompress_symmetric!(
         A::AbstractMatrix{R},
         S::AbstractMatrix{Bool},
-        C::AbstractMatrix{R},
-        colors::AbstractVector{<:Integer}
+        B::AbstractMatrix{R},
+        color::AbstractVector{<:Integer}
     ) where {R<:Real}
 
-Decompress the thin matrix `C` into the symmetric matrix `A` which must have the same sparsity pattern as `S`.
+Decompress the thin matrix `B` into the symmetric matrix `A` which must have the same sparsity pattern as `S`.
+
+Here, `color` is a symmetric coloring of `S`, while `B` is a compressed representation of matrix `A` obtained by summing the columns that share the same color.
 
-Here, `colors` is a symmetric coloring of `S`, while `C` is a compressed representation of matrix `A` obtained by summing the columns that share the same color.
+# References
+
+> [_Efficient Computation of Sparse Hessians Using Coloring and Automatic Differentiation_](https://pubsonline.informs.org/doi/abs/10.1287/ijoc.1080.0286), Gebremedhin et al. (2009), Figure 2
 """
 function decompress_symmetric! end
 
 function decompress_symmetric!(
     A::AbstractMatrix{R},
     S::AbstractMatrix{Bool},
-    C::AbstractMatrix{R},
-    colors::AbstractVector{<:Integer},
+    B::AbstractMatrix{R},
+    color::AbstractVector{<:Integer},
 ) where {R<:Real}
-    A .= zero(R)
-    groups = color_groups(colors)
     checksquare(A)
-    for i in axes(A, 1), j in axes(A, 2)
-        iszero(S[i, j]) && continue
-        ki, kj = colors[i], colors[j]
-        gi, gj = groups[ki], groups[kj]
-        if sum(!iszero, view(S, i, gj)) == 1
-            A[i, j] = C[i, kj]
-        elseif sum(!iszero, view(S, j, gi)) == 1
-            A[i, j] = C[j, ki]
-        else
-            error("Symmetric coloring is not valid")
+    if !same_sparsity_pattern(A, S)
+        throw(DimensionMismatch("`A` and `S` must have the same sparsity pattern."))
+    end
+    A .= zero(R)
+    group = color_groups(color)
+    for ij in findall(!iszero, S)
+        i, j = Tuple(ij)
+        j2_exists = false
+        for j2 in group[color[j]]
+            j2 == j && continue
+            if !iszero(S[i, j2])
+                A[i, j] = A[j, i] = B[j, color[i]]
+                j2_exists = true
+                break
+            end
+        end
+        if !j2_exists
+            A[i, j] = A[j, i] = B[i, color[j]]
         end
     end
     return A
@@ -194,28 +204,32 @@ end
 function decompress_symmetric!(
     A::Symmetric{R},
     S::AbstractMatrix{Bool},
-    C::AbstractMatrix{R},
+    B::AbstractMatrix{R},
     colors::AbstractVector{<:Integer},
 ) where {R<:Real}
     # requires parent decompression to handle both upper and lower triangles
-    decompress_symmetric!(parent(A), S, C, colors)
+    decompress_symmetric!(parent(A), S, B, colors)
     return A
 end
 
 """
     decompress_symmetric(
         S::AbstractMatrix{Bool},
-        C::AbstractMatrix{R},
+        B::AbstractMatrix{R},
         colors::AbstractVector{<:Integer}
     ) where {R<:Real}
 
-Decompress the thin matrix `C` into a new symmetric matrix `A` with the same sparsity pattern as `S`.
+Decompress the thin matrix `B` into a new symmetric matrix `A` with the same sparsity pattern as `S`.
+
+Here, `colors` is a symmetric coloring of `S`, while `B` is a compressed representation of matrix `A` obtained by summing the columns that share the same color.
+
+# References
 
-Here, `colors` is a symmetric coloring of `S`, while `C` is a compressed representation of matrix `A` obtained by summing the columns that share the same color.
+> [_Efficient Computation of Sparse Hessians Using Coloring and Automatic Differentiation_](https://pubsonline.informs.org/doi/abs/10.1287/ijoc.1080.0286), Gebremedhin et al. (2009)
 """
 function decompress_symmetric(
-    S::AbstractMatrix{Bool}, C::AbstractMatrix{R}, colors::AbstractVector{<:Integer}
+    S::AbstractMatrix{Bool}, B::AbstractMatrix{R}, colors::AbstractVector{<:Integer}
 ) where {R<:Real}
     A = respectful_similar(S, R)
-    return decompress_symmetric!(A, S, C, colors)
+    return decompress_symmetric!(A, S, B, colors)
 end
diff --git a/src/graph.jl b/src/graph.jl
index 6502c391..a46d654e 100644
--- a/src/graph.jl
+++ b/src/graph.jl
@@ -88,6 +88,10 @@ The adjacency graph of a symmetrix matric `A ∈ ℝ^{n × n}` is `G(A) = (V, E)
 
 - `V = 1:n` is the set of rows or columns `i`/`j`
 - `(i, j) ∈ E` whenever `A[i, j] ≠ 0` and `i ≠ j`
+
+# References
+
+> [_What Color Is Your Jacobian? Graph Coloring for Computing Derivatives_](https://epubs.siam.org/doi/10.1137/S0036144504444711), Gebremedhin et al. (2005)
 """
 adjacency_graph(H::SparseMatrixCSC) = Graph(H - Diagonal(H))
 adjacency_graph(H::AbstractMatrix) = adjacency_graph(sparse(H))
@@ -102,6 +106,10 @@ The bipartite graph of a matrix `A ∈ ℝ^{m × n}` is `Gb(A) = (V₁, V₂, E)
 - `V₁ = 1:m` is the set of rows `i`
 - `V₂ = 1:n` is the set of columns `j`
 - `(i, j) ∈ E` whenever `A[i, j] ≠ 0`
+
+# References
+
+> [_What Color Is Your Jacobian? Graph Coloring for Computing Derivatives_](https://epubs.siam.org/doi/10.1137/S0036144504444711), Gebremedhin et al. (2005)
 """
 function bipartite_graph(J::SparseMatrixCSC)
     g1 = Graph(transpose(J))  # rows to columns
diff --git a/src/groups.jl b/src/groups.jl
index 729b8ab5..578f3ded 100644
--- a/src/groups.jl
+++ b/src/groups.jl
@@ -1,16 +1,16 @@
 """
-    color_groups(colors)
+    color_groups(color)
 
-Return `groups::Vector{Vector{Int}}` such that `i ∈ groups[c]` iff `colors[i] == c`.
+Return `group::Vector{Vector{Int}}` such that `i ∈ group[c]` iff `color[i] == c`.
 
 Assumes the colors are contiguously numbered from `1` to some `cmax`.
 """
-function color_groups(colors::AbstractVector{<:Integer})
-    cmin, cmax = extrema(colors)
+function color_groups(color::AbstractVector{<:Integer})
+    cmin, cmax = extrema(color)
     @assert cmin == 1
-    groups = [Int[] for c in 1:cmax]
-    for (k, c) in enumerate(colors)
-        push!(groups[c], k)
+    group = [Int[] for c in 1:cmax]
+    for (k, c) in enumerate(color)
+        push!(group[c], k)
     end
-    return groups
+    return group
 end
diff --git a/test/Project.toml b/test/Project.toml
index 820476d6..c685ee33 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,5 +1,4 @@
 [deps]
-ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 Chairmarks = "0ca39b1e-fe0b-4e98-acfc-b1656634c4de"
diff --git a/test/check.jl b/test/check.jl
index df451ac6..6959136a 100644
--- a/test/check.jl
+++ b/test/check.jl
@@ -1,7 +1,5 @@
 using SparseMatrixColorings:
-    check_structurally_orthogonal_columns,
-    check_structurally_orthogonal_rows,
-    check_symmetrically_orthogonal
+    check_structurally_orthogonal_columns, check_symmetrically_orthogonal_columns
 using Test
 
 @testset "Structurally orthogonal columns" begin
@@ -10,10 +8,24 @@ using Test
         0 1 0
         0 1 1
     ]
+
+    # success
+
     @test check_structurally_orthogonal_columns(A, [1, 2, 3])
     @test check_structurally_orthogonal_columns(A, [1, 2, 1])
     @test check_structurally_orthogonal_columns(A, [1, 1, 2])
+
+    # failure
+
+    @test !check_structurally_orthogonal_columns(A, [1, 2])
+    @test_logs (:warn, "2 colors provided for 3 columns") check_structurally_orthogonal_columns(
+        A, [1, 2]; verbose=true
+    )
+
     @test !check_structurally_orthogonal_columns(A, [1, 2, 2])
+    @test_logs (:warn, "In color 2, columns [2, 3] all have nonzeros in row 3") check_structurally_orthogonal_columns(
+        A, [1, 2, 2]; verbose=true
+    )
 end
 
 @testset "Structurally orthogonal rows" begin
@@ -22,14 +34,29 @@ end
         0 1 0
         0 1 1
     ]
-    @test check_structurally_orthogonal_rows(A, [1, 2, 3])
-    @test check_structurally_orthogonal_rows(A, [1, 2, 1])
-    @test check_structurally_orthogonal_rows(A, [1, 1, 2])
-    @test !check_structurally_orthogonal_rows(A, [1, 2, 2])
+
+    # success
+
+    @test check_structurally_orthogonal_columns(transpose(A), [1, 2, 3])
+    @test check_structurally_orthogonal_columns(transpose(A), [1, 2, 1])
+    @test check_structurally_orthogonal_columns(transpose(A), [1, 1, 2])
+
+    # failure
+
+    @test !check_structurally_orthogonal_columns(transpose(A), [1, 2, 2, 3])
+    @test_logs (:warn, "4 colors provided for 3 columns") check_structurally_orthogonal_columns(
+        transpose(A), [1, 2, 2, 3]; verbose=true
+    )
+
+    @test !check_structurally_orthogonal_columns(transpose(A), [1, 2, 2])
+    @test_logs (:warn, "In color 2, columns [2, 3] all have nonzeros in row 2") !check_structurally_orthogonal_columns(
+        transpose(A), [1, 2, 2]; verbose=true
+    )
 end
 
 @testset "Symmetrically orthogonal" begin
     # Fig 4.1 of "What color is your Jacobian?"
+
     A = [
         1 1 0 0 0 0
         1 1 1 0 1 1
@@ -38,6 +65,52 @@ end
         0 1 0 0 1 0
         0 1 0 1 0 1
     ]
-    @test check_symmetrically_orthogonal(A, [1, 2, 1, 3, 1, 1])
-    @test !check_symmetrically_orthogonal(A, [1, 3, 1, 3, 1, 1])
+    @test issymmetric(A)
+
+    # success
+
+    @test check_symmetrically_orthogonal_columns(A, [1, 2, 1, 3, 1, 1])
+
+    # failure
+
+    @test !check_symmetrically_orthogonal_columns(A, [1, 2, 1, 3, 1])
+    @test_logs (:warn, "5 colors provided for 6 columns") check_symmetrically_orthogonal_columns(
+        A, [1, 2, 1, 3, 1]; verbose=true
+    )
+
+    @test !check_symmetrically_orthogonal_columns(A, [1, 3, 1, 3, 1, 1])
+    @test_logs (
+        :warn,
+        """
+For coefficient (i=2, j=3) with column colors (ci=3, cj=1):
+- in color ci=3, columns [2, 4] all have nonzeros in row j=3
+- in color cj=1, columns [1, 3, 5, 6] all have nonzeros in row i=2
+""",
+    ) check_symmetrically_orthogonal_columns(A, [1, 3, 1, 3, 1, 1]; verbose=true)
+
+    # Fig 1 of "Efficient computation of sparse hessians using coloring and AD"
+
+    A = [
+        1 1 0 0 0 0 1 0 0 0
+        1 1 1 0 1 0 0 0 0 0
+        0 1 1 1 0 1 0 0 0 0
+        0 0 1 1 0 0 0 0 0 1
+        0 1 0 0 1 1 0 1 0 0
+        0 0 1 0 1 1 0 0 1 0
+        1 0 0 0 0 0 1 1 0 0
+        0 0 0 0 1 0 1 1 1 0
+        0 0 0 0 0 1 0 1 1 1
+        0 0 0 1 0 0 0 0 1 1
+    ]
+    @test issymmetric(A)
+
+    # success
+
+    @test check_symmetrically_orthogonal_columns(A, [1, 2, 1, 3, 1, 4, 3, 5, 1, 2])
+
+    # failure
+
+    @test !check_symmetrically_orthogonal_columns(A, [1, 2, 1, 3, 1, 4, 3, 4, 1, 2])
+    @test !check_symmetrically_orthogonal_columns(A, [1, 2, 1, 3, 1, 4, 2, 5, 1, 2])
+    @test !check_symmetrically_orthogonal_columns(A, [1, 2, 1, 4, 1, 4, 3, 5, 1, 2])
 end
diff --git a/test/coloring_correctness.jl b/test/coloring_correctness.jl
index 8c98b29e..a1acb21d 100644
--- a/test/coloring_correctness.jl
+++ b/test/coloring_correctness.jl
@@ -1,11 +1,10 @@
-using ADTypes: column_coloring, row_coloring, symmetric_coloring
 using LinearAlgebra: I, Symmetric
 using SparseArrays: sprand
+using SparseMatrixColorings
 using SparseMatrixColorings:
     GreedyColoringAlgorithm,
     check_structurally_orthogonal_columns,
-    check_structurally_orthogonal_rows,
-    check_symmetrically_orthogonal,
+    check_symmetrically_orthogonal_columns,
     matrix_versions
 using StableRNGs
 using Test
@@ -17,36 +16,39 @@ algo = GreedyColoringAlgorithm()
 @test startswith(string(algo), "GreedyColoringAlgorithm(")
 
 @testset "Column coloring" begin
-    @testset "A::$(typeof(A)) of size $(size(A))" for A in vcat(
-        matrix_versions(sprand(rng, Bool, 100, 200, 0.05)),
-        matrix_versions(sprand(rng, Bool, 200, 100, 0.05)),
-    )
-        column_colors = column_coloring(A, algo)
-        @test check_structurally_orthogonal_columns(A, column_colors)
-        @test minimum(column_colors) == 1
-        @test maximum(column_colors) < size(A, 2) ÷ 2
+    @testset "Size ($m, $n) - sparsity $p" for (m, n, p) in [
+        (10, 20, 0.05), (20, 10, 0.05), (100, 200, 0.05), (200, 100, 0.05)
+    ]
+        @testset "$(typeof(A))" for A in matrix_versions(sprand(rng, Bool, m, n, p))
+            column_color = column_coloring(A, algo)
+            @test check_structurally_orthogonal_columns(A, column_color)
+            @test minimum(column_color) == 1
+            @test maximum(column_color) < size(A, 2) ÷ 2
+        end
     end
 end;
 
 @testset "Row coloring" begin
-    @testset "A::$(typeof(A)) of size $(size(A))" for A in vcat(
-        matrix_versions(sprand(rng, Bool, 100, 200, 0.05)),
-        matrix_versions(sprand(rng, Bool, 200, 100, 0.05)),
-    )
-        row_colors = row_coloring(A, algo)
-        @test check_structurally_orthogonal_rows(A, row_colors)
-        @test minimum(row_colors) == 1
-        @test maximum(row_colors) < size(A, 1) ÷ 2
+    @testset "Size ($m, $n) - sparsity $p" for (m, n, p) in [
+        (10, 20, 0.05), (20, 10, 0.05), (100, 200, 0.05), (200, 100, 0.05)
+    ]
+        @testset "$(typeof(A))" for A in matrix_versions(sprand(rng, Bool, m, n, 0.05))
+            row_color = row_coloring(A, algo)
+            @test check_structurally_orthogonal_columns(transpose(A), row_color)
+            @test minimum(row_color) == 1
+            @test maximum(row_color) < size(A, 1) ÷ 2
+        end
     end
 end;
 
 @testset "Symmetric coloring" begin
-    @testset "A::$(typeof(A)) of size $(size(A))" for A in matrix_versions(
-        Symmetric(sprand(rng, Bool, 100, 100, 0.05))
-    )
-        symmetric_colors = symmetric_coloring(A, algo)
-        @test check_symmetrically_orthogonal(A, symmetric_colors)
-        @test minimum(symmetric_colors) == 1
-        @test maximum(symmetric_colors) < size(A, 2) ÷ 2
+    @testset "Size ($n, $n) - sparsity $p" for (n, p) in [(10, 0.05), (100, 0.05)]
+        @testset "$(typeof(A))" for A in
+                                    matrix_versions(Symmetric(sprand(rng, Bool, n, n, p)))
+            symmetric_color = symmetric_coloring(A, algo)
+            @test check_symmetrically_orthogonal_columns(A, symmetric_color)
+            @test minimum(symmetric_color) == 1
+            @test maximum(symmetric_color) < size(A, 2) ÷ 2
+        end
     end
 end;
diff --git a/test/coloring_performance.jl b/test/coloring_performance.jl
index 132c62aa..5f55d0cc 100644
--- a/test/coloring_performance.jl
+++ b/test/coloring_performance.jl
@@ -1,10 +1,9 @@
-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 SparseMatrixColorings
+using SparseMatrixColorings: adjacency_graph, bipartite_graph, partial_distance2_coloring!
 using StableRNGs
 using Test
 
@@ -17,27 +16,18 @@ rng = StableRNG(63)
     @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
+        Symmetric(A), algo
     )
 end
 
 function benchmark_distance2_coloring(n)
     @be (;
         bg=bipartite_graph(sprand(Bool, n, n, 3 / n)),
-        colors=Vector{Int}(undef, n),
+        color=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
+    ) partial_distance2_coloring!(_.color, _.forbidden_colors, _.bg, Val(1), 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/decompression_correctness.jl b/test/decompression_correctness.jl
index 26cc90da..286f9bae 100644
--- a/test/decompression_correctness.jl
+++ b/test/decompression_correctness.jl
@@ -1,5 +1,4 @@
 using Base.Iterators: product
-using ADTypes: column_coloring, row_coloring, symmetric_coloring
 using Compat
 using LinearAlgebra
 using SparseArrays
@@ -33,16 +32,16 @@ algo = GreedyColoringAlgorithm()
             0 1 1
             1 0 0
         ]
-        C = [
+        B = [
             1 2
             3 4
             5 0
         ]
-        colors = [1, 1, 2]
+        color = [1, 1, 2]
         @testset "A::$(typeof(A)) - S::$(typeof(S))" for (A, S) in product(
             matrix_versions(A0), matrix_versions(S0)
         )
-            @test decompress_columns(S, C, colors) == A
+            @test decompress_columns(S, B, color) == A
         end
     end
 
@@ -50,16 +49,16 @@ algo = GreedyColoringAlgorithm()
         m, n = 18, 20
         A0 = sprand(rng, Bool, m, n, 0.2)
         S0 = map(!iszero, A0)
-        colors = column_coloring(A0, algo)
-        groups = color_groups(colors)
-        C = stack(groups; dims=2) do group
-            dropdims(sum(A0[:, group]; dims=2); dims=2)
+        color = column_coloring(A0, algo)
+        group = color_groups(color)
+        B = stack(group; dims=2) do g
+            dropdims(sum(A0[:, g]; dims=2); dims=2)
         end
-        @test size(C) == (size(A0, 1), length(groups))
+        @test size(B) == (size(A0, 1), length(group))
         @testset "A::$(typeof(A)) - S::$(typeof(S))" for (A, S) in product(
             matrix_versions(A0), matrix_versions(S0)
         )
-            @test decompress_columns(S, C, colors) == A
+            @test decompress_columns(S, B, color) == A
         end
     end
 end;
@@ -76,15 +75,15 @@ end;
             0 1 0
             1 1 0
         ]
-        C = [
+        B = [
             1 2 3
             4 5 0
         ]
-        colors = [1, 1, 2]
+        color = [1, 1, 2]
         @testset "A::$(typeof(A)) - S::$(typeof(S))" for (A, S) in product(
             matrix_versions(A0), matrix_versions(S0)
         )
-            @test decompress_rows(S, C, colors) == A
+            @test decompress_rows(S, B, color) == A
         end
     end
 
@@ -92,16 +91,16 @@ end;
         m, n = 18, 20
         A0 = sprand(rng, Bool, m, n, 0.2)
         S0 = map(!iszero, A0)
-        colors = row_coloring(A0, algo)
-        groups = color_groups(colors)
-        C = stack(groups; dims=1) do group
-            dropdims(sum(A0[group, :]; dims=1); dims=1)
+        color = row_coloring(A0, algo)
+        group = color_groups(color)
+        B = stack(group; dims=1) do g
+            dropdims(sum(A0[g, :]; dims=1); dims=1)
         end
-        @test size(C) == (length(groups), size(A0, 2))
+        @test size(B) == (length(group), size(A0, 2))
         @testset "A::$(typeof(A)) - S::$(typeof(S))" for (A, S) in product(
             matrix_versions(A0), matrix_versions(S0)
         )
-            @test decompress_rows(S, C, colors) == A
+            @test decompress_rows(S, B, color) == A
         end
     end
 end;
@@ -118,7 +117,7 @@ end;
             0 6 0 10 0 12
         ]
         S0 = (!iszero).(A0)
-        colors = [
+        color = [
             1,  # green
             2,  # red
             1,  # green
@@ -126,14 +125,14 @@ end;
             1,  # green
             1,  # green
         ]
-        groups = color_groups(colors)
-        C = stack(groups; dims=2) do group
-            dropdims(sum(A0[:, group]; dims=2); dims=2)
+        group = color_groups(color)
+        B = stack(group; dims=2) do g
+            dropdims(sum(A0[:, g]; dims=2); dims=2)
         end
         @testset "A::$(typeof(A)) - S::$(typeof(S))" for (A, S) in product(
             matrix_versions(A0), matrix_versions(S0)
         )
-            @test decompress_symmetric(S, C, colors) == A
+            @test decompress_symmetric(S, B, color) == A
         end
     end
 end;
diff --git a/test/reference/what_table_41_42.csv b/test/reference/what_table_41_42.csv
deleted file mode 100644
index ed67137f..00000000
--- a/test/reference/what_table_41_42.csv
+++ /dev/null
@@ -1,7 +0,0 @@
-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,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/suitesparse.jl b/test/suitesparse.jl
index 8dab06a5..07b441df 100644
--- a/test/suitesparse.jl
+++ b/test/suitesparse.jl
@@ -14,7 +14,7 @@ using SparseMatrixColorings:
     maximum_degree,
     neighbors,
     partial_distance2_coloring,
-    star_coloring1,
+    star_coloring,
     vertices
 using Test
 
@@ -39,10 +39,10 @@ colpack_table_6_7 = CSV.read(
         @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]
+        color_N1 = partial_distance2_coloring(bg, Val(1), NaturalOrder())
+        color_N2 = partial_distance2_coloring(bg, Val(2), NaturalOrder())
+        @test length(unique(color_N1)) == row[:N1]
+        @test length(unique(color_N2)) == row[:N2]
     end
 end;
 
@@ -68,40 +68,13 @@ what_table_31_32 = CSV.read(
         @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]
+        color_Nb = partial_distance2_coloring(bg, Val(2), NaturalOrder())
+        if length(unique(color_Nb)) == row[:K]
+            @test length(unique(color_Nb)) == row[:K]
         else
-            @test_broken length(unique(colors_Nb)) == row[:K]
+            @test_broken length(unique(color_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]) against survey paper"
-        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]
-        @test maximum_degree(g) == row[:Δ]
-        @test minimum_degree(g) == row[:δ]
-        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;