diff --git a/README.md b/README.md
index 407683ea..0e58b7d3 100644
--- a/README.md
+++ b/README.md
@@ -17,11 +17,13 @@ pkg> add https://github.com/gdalle/SparseMatrixColorings.jl
 
 ## Background
 
-The algorithms implemented in this package are all taken from the following survey:
+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)
 
-Some parts of the survey (like definitions and theorems) are also copied verbatim or referred to by their number in the documentation.
+> [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.
 
 ## Alternatives
 
diff --git a/docs/Project.toml b/docs/Project.toml
index 555dcc97..1dca586a 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,3 +1,4 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656"
 SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35"
diff --git a/docs/make.jl b/docs/make.jl
index 534e3953..4bdd90ee 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,6 +1,9 @@
 using Documenter
+using DocumenterInterLinks
 using SparseMatrixColorings
 
+links = InterLinks("ADTypes" => "https://sciml.github.io/ADTypes.jl/stable/")
+
 cp(joinpath(@__DIR__, "..", "README.md"), joinpath(@__DIR__, "src", "index.md"); force=true)
 
 makedocs(;
@@ -9,6 +12,7 @@ makedocs(;
     sitename="SparseMatrixColorings.jl",
     format=Documenter.HTML(),
     pages=["Home" => "index.md", "API reference" => "api.md"],
+    plugins=[links],
 )
 
 deploydocs(; repo="github.com/gdalle/SparseMatrixColorings.jl", devbranch="main")
diff --git a/src/SparseMatrixColorings.jl b/src/SparseMatrixColorings.jl
index f2d1e04e..899ddbc5 100644
--- a/src/SparseMatrixColorings.jl
+++ b/src/SparseMatrixColorings.jl
@@ -3,27 +3,24 @@
 
 Coloring algorithms for sparse Jacobian and Hessian matrices.
 
-The algorithms implemented in this package are all taken from the following survey:
+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)
 
-Some parts of the survey (like definitions and theorems) are also copied verbatim or referred to by their number in the documentation.
+> [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.
 """
 module SparseMatrixColorings
 
 using ADTypes: ADTypes, AbstractColoringAlgorithm
-using LinearAlgebra: Transpose, parent, transpose
-using Random: AbstractRNG
-using SparseArrays: SparseMatrixCSC, nzrange, rowvals
-
-include("utils.jl")
-
-include("bipartite_graph.jl")
-include("adjacency_graph.jl")
-
-include("distance2_coloring.jl")
-include("star_coloring.jl")
+using LinearAlgebra: Diagonal, Transpose, checksquare, parent, transpose
+using Random: AbstractRNG, default_rng, randperm
+using SparseArrays: SparseArrays, SparseMatrixCSC, nzrange, rowvals, spzeros
 
+include("graph.jl")
+include("order.jl")
+include("coloring.jl")
 include("adtypes.jl")
 include("check.jl")
 
diff --git a/src/adjacency_graph.jl b/src/adjacency_graph.jl
deleted file mode 100644
index 1f4142da..00000000
--- a/src/adjacency_graph.jl
+++ /dev/null
@@ -1,26 +0,0 @@
-"""
-    AdjacencyGraph
-
-Represent a graph between the columns of a symmetric `n × n` matrix `A` with nonzero diagonal elements.
-
-This graph is defined as `G = (C, E)` where `C = 1:n` is the set of columns and `(i, j) ∈ E` whenever `A[i, j]` is nonzero for some `j ∈ 1:m`, `j ≠ i`.
-
-# Fields
-
-- `A_colmajor::AbstractMatrix`: output of [`col_major`](@ref) applied to `A`
-"""
-struct AdjacencyGraph{M<:AbstractMatrix}
-    A_colmajor::M
-
-    function AdjacencyGraph(A::AbstractMatrix)
-        A_colmajor = col_major(A)
-        return new{typeof(A_colmajor)}(A_colmajor)
-    end
-end
-
-rows(g::AdjacencyGraph) = axes(g.A_colmajor, 1)
-columns(g::AdjacencyGraph) = axes(g.A_colmajor, 2)
-
-function neighbors(g::AdjacencyGraph, j::Integer)
-    return filter(!isequal(j), nz_in_col(g.A_colmajor, j))
-end
diff --git a/src/adtypes.jl b/src/adtypes.jl
index d9c87775..e9bece08 100644
--- a/src/adtypes.jl
+++ b/src/adtypes.jl
@@ -5,25 +5,32 @@ Matrix coloring algorithm for sparse Jacobians and Hessians.
 
 Compatible with the [ADTypes.jl coloring framework](https://sciml.github.io/ADTypes.jl/stable/#Coloring-algorithm).
 
+# Constructor
+
+    GreedyColoringAlgorithm(order::AbstractOrder=NaturalOrder())
+
 # Implements
 
-- `ADTypes.column_coloring` with a partial distance-2 coloring of the bipartite graph
-- `ADTypes.row_coloring` with a partial distance-2 coloring of the bipartite graph
-- `ADTypes.symmetric_coloring` with a star coloring of the adjacency graph
+- [`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
 """
-struct GreedyColoringAlgorithm <: ADTypes.AbstractColoringAlgorithm end
+struct GreedyColoringAlgorithm{O<:AbstractOrder} <: ADTypes.AbstractColoringAlgorithm
+    order::O
+end
+
+GreedyColoringAlgorithm() = GreedyColoringAlgorithm(NaturalOrder())
 
-function ADTypes.column_coloring(A::AbstractMatrix, ::GreedyColoringAlgorithm)
-    g = BipartiteGraph(A)
-    return distance2_column_coloring(g)
+function ADTypes.column_coloring(A::AbstractMatrix, algo::GreedyColoringAlgorithm)
+    bg = BipartiteGraph(A)
+    return partial_distance2_coloring(bg, Val(2), algo.order)
 end
 
-function ADTypes.row_coloring(A::AbstractMatrix, ::GreedyColoringAlgorithm)
-    g = BipartiteGraph(A)
-    return distance2_row_coloring(g)
+function ADTypes.row_coloring(A::AbstractMatrix, algo::GreedyColoringAlgorithm)
+    bg = BipartiteGraph(A)
+    return partial_distance2_coloring(bg, Val(1), algo.order)
 end
 
-function ADTypes.symmetric_coloring(A::AbstractMatrix, ::GreedyColoringAlgorithm)
-    g = AdjacencyGraph(A)
-    return star_coloring(g)
+function ADTypes.symmetric_coloring(A::AbstractMatrix, algo::GreedyColoringAlgorithm)
+    ag = AdjacencyGraph(A)
+    return star_coloring(ag, algo.order)
 end
diff --git a/src/bipartite_graph.jl b/src/bipartite_graph.jl
deleted file mode 100644
index 45d7a6f1..00000000
--- a/src/bipartite_graph.jl
+++ /dev/null
@@ -1,28 +0,0 @@
-"""
-    BipartiteGraph
-
-Represent a bipartite graph between the rows and the columns of a non-symmetric `m × n` matrix `A`.
-
-This graph is defined as `G = (R, C, E)` where `R = 1:m` is the set of row indices, `C = 1:n` is the set of column indices and `(i, j) ∈ E` whenever `A[i, j]` is nonzero.
-
-# Fields
-
-- `A_colmajor::AbstractMatrix`: output of [`col_major`](@ref) applied to `A` (useful to get neighbors of a column)
-- `A_rowmajor::AbstractMatrix`: output of [`row_major`](@ref) applied to `A` (useful to get neighbors of a row)
-"""
-struct BipartiteGraph{M1<:AbstractMatrix,M2<:AbstractMatrix}
-    A_colmajor::M1
-    A_rowmajor::M2
-
-    function BipartiteGraph(A::AbstractMatrix)
-        A_colmajor = col_major(A)
-        A_rowmajor = row_major(A)
-        return new{typeof(A_colmajor),typeof(A_rowmajor)}(A_colmajor, A_rowmajor)
-    end
-end
-
-rows(g::BipartiteGraph) = axes(g.A_colmajor, 1)
-columns(g::BipartiteGraph) = axes(g.A_colmajor, 2)
-
-neighbors_of_column(g::BipartiteGraph, j::Integer) = nz_in_col(g.A_colmajor, j)
-neighbors_of_row(g::BipartiteGraph, i::Integer) = nz_in_row(g.A_rowmajor, i)
diff --git a/src/check.jl b/src/check.jl
index af0d640f..3f8deba5 100644
--- a/src/check.jl
+++ b/src/check.jl
@@ -1,14 +1,18 @@
 """
-    check_structurally_orthogonal_columns(A, colors; verbose=false)
+    check_structurally_orthogonal_columns(
+        A::AbstractMatrix, colors::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.
     
-Def 3.2: 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`.
+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`.
 
-Thm 3.5: The function [`distance2_column_coloring`](@ref) applied to the [`BipartiteGraph`](@ref) of `A` should return a suitable coloring.
+!!! warning
+    This function is not coded with efficiency in mind, it is designed for small-scale tests.
 """
 function check_structurally_orthogonal_columns(
-    A::AbstractMatrix, colors::AbstractVector{<:Integer}; verbose=false
+    A::AbstractMatrix, colors::AbstractVector{<:Integer}; verbose::Bool=false
 )
     for c in unique(colors)
         js = filter(j -> colors[j] == c, axes(A, 2))
@@ -23,16 +27,20 @@ function check_structurally_orthogonal_columns(
 end
 
 """
-    check_structurally_orthogonal_rows(A, colors; verbose=false)
+    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.
     
-Def 3.2: 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`.
+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`.
 
-Thm 3.5: The function [`distance2_row_coloring`](@ref) applied to the [`BipartiteGraph`](@ref) of `A` should return a suitable coloring.
+!!! 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=false
+    A::AbstractMatrix, colors::AbstractVector{<:Integer}; verbose::Bool=false
 )
     for c in unique(colors)
         is = filter(i -> colors[i] == c, axes(A, 1))
@@ -47,17 +55,23 @@ function check_structurally_orthogonal_rows(
 end
 
 """
-    check_symmetrically_orthogonal(A, colors; verbose=false)
+    check_symmetrically_orthogonal(
+        A::AbstractMatrix, colors::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.
     
-Def 4.2: 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
 
 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.
 """
 function check_symmetrically_orthogonal(
-    A::AbstractMatrix, colors::AbstractVector{<:Integer}; verbose=false
+    A::AbstractMatrix, colors::AbstractVector{<:Integer}; verbose::Bool=false
 )
     for i in axes(A, 2), j in axes(A, 2)
         if !iszero(A[i, j])
diff --git a/src/coloring.jl b/src/coloring.jl
new file mode 100644
index 00000000..e194e3f4
--- /dev/null
+++ b/src/coloring.jl
@@ -0,0 +1,85 @@
+"""
+    partial_distance2_coloring(bg::BipartiteGraph, ::Val{side}, order::AbstractOrder)
+
+Compute a distance-2 coloring of the given `side` (`1` or `2`) in the bipartite graph `bg` and return a vector of integer colors.
+
+A _distance-2 coloring_ is such that two vertices have different colors if they are at distance at most 2.
+
+The vertices are colored in a greedy fashion, following the `order` supplied.
+
+# See also
+
+- [`BipartiteGraph`](@ref)
+- [`AbstractOrder`](@ref)
+"""
+function partial_distance2_coloring(
+    bg::BipartiteGraph, ::Val{side}, order::AbstractOrder
+) where {side}
+    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 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
+                end
+            end
+        end
+        for c in eachindex(forbidden_colors)
+            if forbidden_colors[c] != v
+                colors[v] = c
+                break
+            end
+        end
+    end
+    return colors
+end
+
+"""
+    star_coloring(ag::AdjacencyGraph, order::AbstractOrder)
+
+Compute a star coloring of all vertices in the adjacency graph `ag` 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.
+
+The vertices are colored in a greedy fashion, following the `order` supplied.
+
+# See also
+
+- [`AdjacencyGraph`](@ref)
+- [`AbstractOrder`](@ref)
+"""
+function star_coloring(ag::AdjacencyGraph, order::AbstractOrder)
+    n = length(ag)
+    colors = zeros(Int, n)
+    forbidden_colors = zeros(Int, n)
+    for v in vertices(ag, order)
+        for w in neighbors(ag, v)
+            if !iszero(colors[w])  # w is colored
+                forbidden_colors[colors[w]] = v
+            end
+            for x in neighbors(ag, w)
+                if !iszero(colors[x]) && iszero(colors[w])  # w is not colored
+                    forbidden_colors[colors[x]] = v
+                else
+                    for y in neighbors(ag, x)
+                        if !iszero(colors[y]) && y != w
+                            if colors[y] == colors[w]
+                                forbidden_colors[colors[x]] = v
+                                break
+                            end
+                        end
+                    end
+                end
+            end
+        end
+        for c in eachindex(forbidden_colors)
+            if forbidden_colors[c] != v
+                colors[v] = c
+                break
+            end
+        end
+    end
+    return colors
+end
diff --git a/src/distance2_coloring.jl b/src/distance2_coloring.jl
deleted file mode 100644
index 88e0d8a5..00000000
--- a/src/distance2_coloring.jl
+++ /dev/null
@@ -1,61 +0,0 @@
-"""
-    distance2_column_coloring(g::BipartiteGraph)
-
-Compute a distance-2 coloring of the column vertices in the [`BipartiteGraph`](@ref) `g` and return a vector of integer colors.
-
-A _distance-2 coloring_ is such that two vertices have different colors if they are at distance at most 2.
-
-The algorithm used is the greedy Algorithm 3.2.
-"""
-function distance2_column_coloring(g::BipartiteGraph)
-    n = length(columns(g))
-    colors = zeros(Int, n)
-    forbidden_colors = zeros(Int, n)
-    for v in sort(columns(g); by=j -> length(neighbors_of_column(g, j)), rev=true)
-        for w in neighbors_of_column(g, v)
-            for x in neighbors_of_row(g, w)
-                if !iszero(colors[x])
-                    forbidden_colors[colors[x]] = v
-                end
-            end
-        end
-        for c in eachindex(forbidden_colors)
-            if forbidden_colors[c] != v
-                colors[v] = c
-                break
-            end
-        end
-    end
-    return colors
-end
-
-"""
-    distance2_row_coloring(g::BipartiteGraph)
-
-Compute a distance-2 coloring of the row vertices in the [`BipartiteGraph`](@ref) `g` and return a vector of integer colors.
-
-A _distance-2 coloring_ is such that two vertices have different colors if they are at distance at most 2.
-
-The algorithm used is the greedy Algorithm 3.2.
-"""
-function distance2_row_coloring(g::BipartiteGraph)
-    m = length(rows(g))
-    colors = zeros(Int, m)
-    forbidden_colors = zeros(Int, m)
-    for v in sort(rows(g); by=i -> length(neighbors_of_row(g, i)), rev=true)
-        for w in neighbors_of_row(g, v)
-            for x in neighbors_of_column(g, w)
-                if !iszero(colors[x])
-                    forbidden_colors[colors[x]] = v
-                end
-            end
-        end
-        for c in eachindex(forbidden_colors)
-            if forbidden_colors[c] != v
-                colors[v] = c
-                break
-            end
-        end
-    end
-    return colors
-end
diff --git a/src/graph.jl b/src/graph.jl
new file mode 100644
index 00000000..68c0f025
--- /dev/null
+++ b/src/graph.jl
@@ -0,0 +1,81 @@
+struct Graph{T<:Integer}
+    colptr::Vector{T}
+    rowval::Vector{T}
+end
+
+Graph(A::SparseMatrixCSC) = Graph(A.colptr, A.rowval)
+
+Base.length(g::Graph) = length(g.colptr) - 1
+
+neighbors(g::Graph, v::Integer) = view(g.rowval, g.colptr[v]:(g.colptr[v + 1] - 1))
+
+## Adjacency graph
+
+"""
+    AdjacencyGraph
+
+Undirected graph representing the nonzeros of a symmetrix matrix (typically a Hessian matrix).
+
+The adjacency graph of a symmetrix matrix `A ∈ ℝ^{n × n}` is `G(A) = (V, E)` where
+
+- `V = 1:n` is the set of rows or columns
+- `(i, j) ∈ E` whenever `A[i, j] ≠ 0` and `i ≠ j`
+"""
+struct AdjacencyGraph{T}
+    g::Graph{T}
+end
+
+function AdjacencyGraph(H::SparseMatrixCSC)
+    g = Graph(H - Diagonal(H))
+    return AdjacencyGraph(g)
+end
+
+Base.length(ag::AdjacencyGraph) = length(ag.g)
+
+"""
+    neighbors(ag::AdjacencyGraph, v::Integer)
+
+Return the neighbors of `v` in the graph `ag`.
+"""
+neighbors(ag::AdjacencyGraph, v::Integer) = neighbors(ag.g, v)
+
+degree(ag::AdjacencyGraph, v::Integer) = length(neighbors(ag, v))
+
+## Bipartite graph
+
+"""
+    BipartiteGraph
+
+Undirected bipartite graph 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
+
+- `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`
+"""
+struct BipartiteGraph{T}
+    g1::Graph{T}
+    g2::Graph{T}
+end
+
+function BipartiteGraph(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)
+
+"""
+    neighbors(bg::BipartiteGraph, Val(side), v::Integer)
+
+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)
+
+function degree(bg::BipartiteGraph, ::Val{side}, v::Integer) where {side}
+    return length(neighbors(bg, Val(side), v))
+end
diff --git a/src/order.jl b/src/order.jl
new file mode 100644
index 00000000..b7a9c1e1
--- /dev/null
+++ b/src/order.jl
@@ -0,0 +1,63 @@
+"""
+    AbstractOrder
+
+Abstract supertype for vertex ordering schemes.
+
+# Subtypes
+
+- [`NaturalOrder`](@ref)
+- [`RandomOrder`](@ref)
+- [`LargestFirst`](@ref)
+"""
+abstract type AbstractOrder end
+
+"""
+    NaturalOrder()
+
+Order vertices as they come in the graph.
+"""
+struct NaturalOrder <: AbstractOrder end
+
+function vertices(ag::AdjacencyGraph, ::NaturalOrder)
+    return 1:length(ag)
+end
+
+function vertices(bg::BipartiteGraph, ::Val{side}, ::NaturalOrder) where {side}
+    return 1:length(bg, Val(side))
+end
+
+"""
+    RandomOrder(rng=default_rng())
+
+Order vertices with a random permutation.
+"""
+struct RandomOrder{R<:AbstractRNG} <: AbstractOrder
+    rng::R
+end
+
+RandomOrder() = RandomOrder(default_rng())
+
+function vertices(ag::AdjacencyGraph, order::RandomOrder)
+    return randperm(order.rng, length(ag))
+end
+
+function vertices(bg::BipartiteGraph, ::Val{side}, order::RandomOrder) where {side}
+    return randperm(order.rng, length(bg, Val(side)))
+end
+
+"""
+    LargestFirst()
+
+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)
+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)
+end
diff --git a/src/star_coloring.jl b/src/star_coloring.jl
deleted file mode 100644
index 116e951e..00000000
--- a/src/star_coloring.jl
+++ /dev/null
@@ -1,42 +0,0 @@
-"""
-    star_coloring(g::BipartiteGraph)
-
-Compute a star coloring of the column vertices in the [`AdjacencyGraph`](@ref) `g` and return a vector of integer colors.
-
-Def 4.5: A _star coloring_ is a distance-1 coloring such that every path on four vertices uses at least three colors.
-
-The algorithm used is the greedy Algorithm 4.1.
-"""
-function star_coloring(g::AdjacencyGraph)
-    n = length(columns(g))
-    colors = zeros(Int, n)
-    forbidden_colors = zeros(Int, n)
-    for v in sort(columns(g); by=j -> length(neighbors(g, j)), rev=true)
-        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
-                    end
-                end
-            end
-        end
-        for c in eachindex(forbidden_colors)
-            if forbidden_colors[c] != v
-                colors[v] = c
-                break
-            end
-        end
-    end
-    return colors
-end
diff --git a/src/utils.jl b/src/utils.jl
deleted file mode 100644
index e7c010c7..00000000
--- a/src/utils.jl
+++ /dev/null
@@ -1,42 +0,0 @@
-## Conversion between row- and col-major
-
-"""
-    col_major(A::AbstractMatrix)
-
-Construct a column-major representation of the matrix `A`.
-"""
-col_major(A::M) where {M<:AbstractMatrix} = A
-col_major(A::Transpose{<:Any,M}) where {M<:AbstractMatrix} = M(A)
-
-"""
-    row_major(A::AbstractMatrix)
-
-Construct a row-major representation of the matrix `A`.
-"""
-row_major(A::M) where {M<:AbstractMatrix} = transpose(M(transpose(A)))
-row_major(A::Transpose{<:Any,M}) where {M<:AbstractMatrix} = A
-
-## Generic nz
-
-function nz_in_col(A_colmajor::AbstractMatrix, j::Integer)
-    return filter(i -> !iszero(A_colmajor[i, j]), axes(A_colmajor, 1))
-end
-
-function nz_in_row(A_rowmajor::AbstractMatrix, i::Integer)
-    return filter(j -> !iszero(A_rowmajor[i, j]), axes(A_rowmajor, 2))
-end
-
-## Sparse nz
-
-function nz_in_col(A_colmajor::SparseMatrixCSC{T}, j::Integer) where {T}
-    rv = rowvals(A_colmajor)
-    ind = nzrange(A_colmajor, j)
-    return view(rv, ind)
-end
-
-function nz_in_row(A_rowmajor::Transpose{T,<:SparseMatrixCSC{T}}, i::Integer) where {T}
-    A_transpose_colmajor = parent(A_rowmajor)
-    rv = rowvals(A_transpose_colmajor)
-    ind = nzrange(A_transpose_colmajor, i)
-    return view(rv, ind)
-end
diff --git a/test/Project.toml b/test/Project.toml
index 201e7a86..617bcb56 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -6,4 +6,5 @@ JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/adtypes.jl b/test/adtypes.jl
index b3f0e01f..36041bed 100644
--- a/test/adtypes.jl
+++ b/test/adtypes.jl
@@ -29,7 +29,7 @@ end
 end
 
 @testset "Symmetric coloring" begin
-    S = Symmetric(sprand(Bool, 100, 100, 0.05)) + I
+    S = sparse(Symmetric(sprand(Bool, 100, 100, 0.05)))
     symmetric_colors = symmetric_coloring(S, alg)
     @test check_symmetrically_orthogonal(S, symmetric_colors)
     @test minimum(symmetric_colors) == 1
diff --git a/test/check.jl b/test/check.jl
index 11b20ea1..0bea965e 100644
--- a/test/check.jl
+++ b/test/check.jl
@@ -4,7 +4,7 @@ using SparseMatrixColorings:
     check_symmetrically_orthogonal
 using Test
 
-@testset "Column coloring" begin
+@testset "Structurally orthogonal columns" begin
     A = [
         1 0 0
         0 1 0
@@ -16,7 +16,7 @@ using Test
     @test !check_structurally_orthogonal_columns(A, [1, 2, 2])
 end
 
-@testset "Row coloring" begin
+@testset "Structurally orthogonal rows" begin
     A = [
         1 0 0
         0 1 0
@@ -28,8 +28,8 @@ end
     @test !check_structurally_orthogonal_rows(A, [1, 2, 2])
 end
 
-@testset "Symmetric coloring" begin
-    # example from "What color is your Jacobian", fig 4.1
+@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
diff --git a/test/graph.jl b/test/graph.jl
new file mode 100644
index 00000000..e4f6a703
--- /dev/null
+++ b/test/graph.jl
@@ -0,0 +1,69 @@
+using SparseArrays
+using SparseMatrixColorings: Graph, BipartiteGraph, AdjacencyGraph, neighbors
+using Test
+
+## Standard graph
+
+@testset "Graph" begin
+    g = Graph(sparse([
+        1 0 1
+        1 1 0
+        0 0 0
+    ]))
+
+    @test length(g) == 3
+    @test neighbors(g, 1) == [1, 2]
+    @test neighbors(g, 2) == [2]
+    @test neighbors(g, 3) == [1]
+end;
+
+## Bipartite graph (fig 3.1 of "What color is your Jacobian?")
+
+@testset "BipartiteGraph" begin
+    A = sparse([
+        1 0 0 0 0 1 1 1
+        0 1 0 0 1 0 1 1
+        0 0 1 0 1 1 0 1
+        0 0 0 1 1 1 1 0
+    ])
+
+    bg = BipartiteGraph(A)
+    @test length(bg, Val(1)) == 4
+    @test length(bg, Val(2)) == 8
+    # neighbors of rows
+    @test neighbors(bg, Val(1), 1) == [1, 6, 7, 8]
+    @test neighbors(bg, Val(1), 2) == [2, 5, 7, 8]
+    @test neighbors(bg, Val(1), 3) == [3, 5, 6, 8]
+    @test neighbors(bg, Val(1), 4) == [4, 5, 6, 7]
+    # neighbors of columns
+    @test neighbors(bg, Val(2), 1) == [1]
+    @test neighbors(bg, Val(2), 2) == [2]
+    @test neighbors(bg, Val(2), 3) == [3]
+    @test neighbors(bg, Val(2), 4) == [4]
+    @test neighbors(bg, Val(2), 5) == [2, 3, 4]
+    @test neighbors(bg, Val(2), 6) == [1, 3, 4]
+    @test neighbors(bg, Val(2), 7) == [1, 2, 4]
+    @test neighbors(bg, Val(2), 8) == [1, 2, 3]
+end;
+
+## Adjacency graph (fig 3.1 of "What color is your Jacobian?")
+
+@testset "AdjacencyGraph" begin
+    A = sparse([
+        1 0 0 0 0 1 1 1
+        0 1 0 0 1 0 1 1
+        0 0 1 0 1 1 0 1
+        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]
+end
diff --git a/test/order.jl b/test/order.jl
new file mode 100644
index 00000000..e89ad448
--- /dev/null
+++ b/test/order.jl
@@ -0,0 +1,64 @@
+using SparseArrays
+using SparseMatrixColorings: Graph, LargestFirst, NaturalOrder, RandomOrder, vertices
+using StableRNGs
+using Test
+
+rng = StableRNG(63)
+
+@testset "NaturalOrder" begin
+    A = sprand(rng, Bool, 5, 5, 0.5)
+    ag = AdjacencyGraph(A)
+    @test vertices(ag, NaturalOrder()) == 1:5
+
+    A = sprand(rng, Bool, 5, 4, 0.5)
+    bg = BipartiteGraph(A)
+    @test vertices(bg, Val(1), NaturalOrder()) == 1:5
+
+    A = sprand(rng, Bool, 5, 4, 0.5)
+    bg = BipartiteGraph(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
+
+    A = sprand(rng, Bool, 5, 4, 0.5)
+    bg = BipartiteGraph(A)
+    @test sort(vertices(bg, Val(1), RandomOrder(rng))) == 1:5
+
+    A = sprand(rng, Bool, 5, 4, 0.5)
+    bg = BipartiteGraph(A)
+    @test sort(vertices(bg, Val(2), RandomOrder(rng))) == 1:4
+end;
+
+@testset "LargestFirst" begin
+    A = sparse([
+        0 1 0
+        1 0 0
+        0 1 0
+    ])
+    ag = AdjacencyGraph(A)
+
+    @test vertices(ag, LargestFirst()) == [2, 1, 3]
+
+    A = sparse([
+        1 1 1
+        1 0 0
+        0 1 1
+        0 0 0
+    ])
+    bg = BipartiteGraph(A)
+
+    @test vertices(bg, Val(1), LargestFirst()) == [1, 3, 2, 4]
+
+    A = sparse([
+        1 1 0 0
+        1 0 1 0
+        1 0 1 0
+    ])
+    bg = BipartiteGraph(A)
+
+    @test vertices(bg, Val(2), LargestFirst()) == [1, 3, 2, 4]
+end;
diff --git a/test/runtests.jl b/test/runtests.jl
index 703e14e6..1486dce2 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -24,8 +24,11 @@ using Test
     @testset "Doctests" begin
         Documenter.doctest(SparseMatrixColorings)
     end
-    @testset "Utils" begin
-        include("utils.jl")
+    @testset "Graph" begin
+        include("graph.jl")
+    end
+    @testset "Order" begin
+        include("order.jl")
     end
     @testset "Check" begin
         include("check.jl")