diff --git a/Project.toml b/Project.toml
index ac81ff2b40..abbf1bbc07 100644
--- a/Project.toml
+++ b/Project.toml
@@ -14,6 +14,7 @@ OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
@@ -42,6 +43,7 @@ ParallelTestRunner = "2.4.1"
 PrecompileTools = "1.0.1"
 Printf = "1"
 SparseArrays = "1"
+SparseMatrixColorings = "0.4.26"
 SpecialFunctions = "1, 2"
 Test = "1"
 julia = "1.10"
diff --git a/src/Nonlinear/Nonlinear.jl b/src/Nonlinear/Nonlinear.jl
index 9a5e4bd8d9..bd0f7a4443 100644
--- a/src/Nonlinear/Nonlinear.jl
+++ b/src/Nonlinear/Nonlinear.jl
@@ -10,6 +10,7 @@ import ForwardDiff
 import ..MathOptInterface as MOI
 import OrderedCollections: OrderedDict
 import SparseArrays
+import SparseMatrixColorings as SMC
 
 using SpecialFunctions
 
diff --git a/src/Nonlinear/ReverseAD/Coloring/Coloring.jl b/src/Nonlinear/ReverseAD/Coloring/Coloring.jl
index 52e5d6f431..ae563784a2 100644
--- a/src/Nonlinear/ReverseAD/Coloring/Coloring.jl
+++ b/src/Nonlinear/ReverseAD/Coloring/Coloring.jl
@@ -6,9 +6,6 @@
 
 module Coloring
 
-include("IntDisjointSet.jl")
-include("topological_sort.jl")
-
 """
     mutable struct IndexedSet
         nzidx::Vector{Int}
@@ -66,500 +63,4 @@ function Base.union!(v::IndexedSet, s)
     return
 end
 
-"""
-    struct UndirectedGraph
-        adjlist::Vector{Int}
-        edgeindex::Vector{Int}
-        offsets::Vector{Int}
-        edges::Vector{Tuple{Int,Int}}
-    end
-
-Compact storage for an undirected graph. The number of nodes is given by
-`length(offsets) - 1`. The edges of node `u` are given by `edges[e]` for
-`e in edgeindex[offsets[u]:(offsets[u+1] - 1)]`. The neighbors are also
-stored at `adjlist[offsets[u]:(offsets[u+1] - 1)]`.
-"""
-struct UndirectedGraph
-    adjlist::Vector{Int}
-    edgeindex::Vector{Int} # corresponding edge number, indexed as adjlist
-    offsets::Vector{Int}
-    edges::Vector{Tuple{Int,Int}}
-end
-
-_num_vertices(g::UndirectedGraph) = length(g.offsets) - 1
-
-_num_edges(g::UndirectedGraph) = length(g.edges)
-
-_num_neighbors(i::Int, g::UndirectedGraph) = g.offsets[i+1] - g.offsets[i]
-
-_start_neighbors(i::Int, g::UndirectedGraph) = g.offsets[i]
-
-function UndirectedGraph(I, J, nel)
-    adjcount = zeros(Int, nel)
-    n_edges = 0
-    for k in eachindex(I)
-        i = I[k]
-        j = J[k]
-        if i == j
-            continue
-        end
-        n_edges += 1
-        adjcount[i] += 1
-        adjcount[j] += 1
-    end
-    offsets = Vector{Int}(undef, nel + 1)
-    offsets[1] = 1
-    for k in 1:nel
-        offsets[k+1] = offsets[k] + adjcount[k]
-    end
-    fill!(adjcount, 0)
-    edges = Vector{Tuple{Int,Int}}(undef, n_edges)
-    adjlist = Vector{Int}(undef, offsets[nel+1] - 1)
-    edgeindex = Vector{Int}(undef, length(adjlist))
-    edge_count = 0
-    for k in eachindex(I)
-        i = I[k]
-        j = J[k]
-        if i == j
-            continue
-        end
-        edge_count += 1
-        adjlist[offsets[i]+adjcount[i]] = j
-        edgeindex[offsets[i]+adjcount[i]] = edge_count
-        adjcount[i] += 1
-        adjlist[offsets[j]+adjcount[j]] = i
-        edgeindex[offsets[j]+adjcount[j]] = edge_count
-        adjcount[j] += 1
-        edges[edge_count] = (i, j)
-    end
-    @assert edge_count == n_edges
-    return UndirectedGraph(adjlist, edgeindex, offsets, edges)
-end
-
-struct _Edge
-    index::Int
-    source::Int
-    target::Int
-end
-
-function _prevent_cycle(
-    v,
-    w,
-    x,
-    e_idx1,
-    e_idx2,
-    S,
-    firstVisitToTree,
-    forbiddenColors,
-    color,
-)
-    er = _find_root!(S, e_idx2)
-    @inbounds first = firstVisitToTree[er]
-    p = first.source # but this depends on the order?
-    q = first.target
-    @inbounds if p != v
-        firstVisitToTree[er] = _Edge(e_idx1, v, w)
-    elseif q != w
-        forbiddenColors[color[x]] = v
-    end
-    return
-end
-
-function _grow_star(v, w, e_idx, firstNeighbor, color, S)
-    @inbounds e = firstNeighbor[color[w]]
-    p = e.source
-    @inbounds if p != v
-        firstNeighbor[color[w]] = _Edge(e_idx, v, w)
-    else
-        root1 = _find_root!(S, e_idx)
-        root2 = _find_root!(S, e.index)
-        _root_union!(S, root1, root2)
-    end
-    return
-end
-
-"""
-    acyclic_coloring(g::UndirectedGraph)
-
-Implement the acyclic coloring algorithm of Gebremdehin, Tarafdar, Manne, and
-Pothen, "New Acyclic and Star Coloring Algorithms with Application to Computing
-Hessians." SIAM J. Sci. Comput. 2007.
-
-Returns `Tuple{Vector{Int},Int}` giving the color index of each node in `g`, as
-well as the total number of colors used.
-"""
-function acyclic_coloring(g::UndirectedGraph)
-    if _num_edges(g) == 0
-        return fill(1, _num_vertices(g)), 1
-    end
-    num_colors = 0
-    forbiddenColors = Int[]
-    firstNeighbor = _Edge[]
-    firstVisitToTree = fill(_Edge(0, 0, 0), _num_edges(g))
-    color = fill(0, _num_vertices(g))
-    S = _IntDisjointSet(_num_edges(g))
-    @inbounds for v in 1:_num_vertices(g)
-        n_neighbor = _num_neighbors(v, g)
-        start_neighbor = _start_neighbors(v, g)
-        for k in 0:(n_neighbor-1)
-            w = g.adjlist[start_neighbor+k]
-            if color[w] == 0
-                continue
-            end
-            forbiddenColors[color[w]] = v
-        end
-        for k in 0:(n_neighbor-1)
-            w = g.adjlist[start_neighbor+k]
-            e_idx = g.edgeindex[start_neighbor+k]
-            if color[w] == 0
-                continue
-            end
-            n_neighbor_w = _num_neighbors(w, g)
-            start_neighbor_w = _start_neighbors(w, g)
-            for k2 in 0:(n_neighbor_w-1)
-                x = g.adjlist[start_neighbor_w+k2]
-                e2_idx = g.edgeindex[start_neighbor_w+k2]
-                if color[x] == 0
-                    continue
-                end
-                if forbiddenColors[color[x]] != v
-                    _prevent_cycle(
-                        v,
-                        w,
-                        x,
-                        e_idx,
-                        e2_idx,
-                        S,
-                        firstVisitToTree,
-                        forbiddenColors,
-                        color,
-                    )
-                end
-            end
-        end
-        # find feasible color
-        found = false
-        for k in 1:num_colors
-            if forbiddenColors[k] != v
-                color[v] = k
-                found = true
-                break
-            end
-        end
-        if !found
-            num_colors += 1
-            push!(forbiddenColors, 0)
-            push!(firstNeighbor, _Edge(0, 0, 0))
-            color[v] = num_colors
-        end
-        for k in 0:(n_neighbor-1)
-            w = g.adjlist[start_neighbor+k]
-            e_idx = g.edgeindex[start_neighbor+k]
-            if color[w] == 0
-                continue
-            end
-            _grow_star(v, w, e_idx, firstNeighbor, color, S)
-        end
-        for k in 0:(n_neighbor-1)
-            w = g.adjlist[start_neighbor+k]
-            e_idx = g.edgeindex[start_neighbor+k]
-            if color[w] == 0
-                continue
-            end
-            n_neighbor_w = _num_neighbors(w, g)
-            start_neighbor_w = _start_neighbors(w, g)
-            for k2 in 0:(n_neighbor_w-1)
-                x = g.adjlist[start_neighbor_w+k2]
-                e2_idx = g.edgeindex[start_neighbor_w+k2]
-                if color[x] == 0 || x == v
-                    continue
-                end
-                if color[x] == color[v]
-                    _union!(S, e_idx, e2_idx)
-                end
-            end
-        end
-    end
-    return color, num_colors
-end
-
-struct RecoveryInfo
-    vertexmap::Vector{Vector{Int}}
-    postorder::Vector{Vector{Int}}
-    parents::Vector{Vector{Int}}
-    color::Vector{Int}
-    num_colors::Int
-    nnz::Int # number of off-diagonal nonzeros
-    local_indices::Vector{Int} # map back to global indices
-end
-
-function RecoveryInfo()
-    return RecoveryInfo(
-        Vector{Vector{Int}}(undef, 0),
-        Vector{Vector{Int}}(undef, 0),
-        Vector{Vector{Int}}(undef, 0),
-        Vector{Int}(undef, 0),
-        0,
-        0,
-        Vector{Int}(undef, 0),
-    )
-end
-
-function recovery_preprocess(
-    g::UndirectedGraph,
-    color,
-    num_colors,
-    local_indices,
-)
-    # represent two-color subgraph as:
-    # list of vertices (with map to global indices)
-    # adjacency list in a single vector (with list of offsets)
-    # linear index of pair of colors
-    twocolorindex = zeros(Int32, num_colors, num_colors)
-    seen_twocolors = 0
-    # count of edges in each subgraph
-    edge_count = Int[]
-    for k in eachindex(g.edges)
-        u, v = g.edges[k]
-        i = min(color[u], color[v])
-        j = max(color[u], color[v])
-        if twocolorindex[i, j] == 0
-            seen_twocolors += 1
-            twocolorindex[i, j] = seen_twocolors
-            push!(edge_count, 0)
-        end
-        idx = twocolorindex[i, j]
-        edge_count[idx] += 1
-    end
-    # edges sorted by twocolor subgraph
-    sorted_edges = Vector{Vector{Tuple{Int,Int}}}(undef, seen_twocolors)
-    for idx in 1:seen_twocolors
-        sorted_edges[idx] = Tuple{Int,Int}[]
-        sizehint!(sorted_edges[idx], edge_count[idx])
-    end
-    for k in eachindex(g.edges)
-        u, v = g.edges[k]
-        i = min(color[u], color[v])
-        j = max(color[u], color[v])
-        idx = twocolorindex[i, j]
-        push!(sorted_edges[idx], (u, v))
-    end
-    # list of unique vertices in each twocolor subgraph
-    vertexmap = Vector{Vector{Int}}(undef, seen_twocolors)
-    postorder = Vector{Vector{Int}}(undef, seen_twocolors)
-    parents = Vector{Vector{Int}}(undef, seen_twocolors)
-    # temporary lookup map from global index to subgraph index
-    revmap = zeros(Int, _num_vertices(g))
-    adjcount = zeros(Int, _num_vertices(g))
-    cmap = zeros(Int, 0) # shared storage for DFS
-    for idx in 1:seen_twocolors
-        my_edges = sorted_edges[idx]
-        vlist = Int[]
-        # build up the vertex list and adjacency count
-        for k in eachindex(my_edges)
-            u, v = my_edges[k]
-            # seen these vertices yet?
-            if revmap[u] == 0
-                push!(vlist, u)
-                revmap[u] = length(vlist)
-                adjcount[u] = 0
-            end
-            if revmap[v] == 0
-                push!(vlist, v)
-                revmap[v] = length(vlist)
-                adjcount[v] = 0
-            end
-            adjcount[u] += 1
-            adjcount[v] += 1
-        end
-        # set up offsets for adjlist
-        offset = Vector{Int}(undef, length(vlist) + 1)
-        offset[1] = 1
-        for k in eachindex(vlist)
-            offset[k+1] = offset[k] + adjcount[vlist[k]]
-            adjcount[vlist[k]] = 0
-        end
-        # adjlist for node u in twocolor idx starts at
-        # vec[offset[u]]
-        # u has global index vlist[u]
-        vec = Vector{Int}(undef, offset[length(vlist)+1] - 1)
-        # now fill in
-        for k in eachindex(my_edges)
-            u, v = my_edges[k]
-            u_rev = revmap[u] # indices in the subgraph
-            v_rev = revmap[v]
-            vec[offset[u_rev]+adjcount[u]] = v_rev
-            vec[offset[v_rev]+adjcount[v]] = u_rev
-            adjcount[u] += 1
-            adjcount[v] += 1
-        end
-        resize!(cmap, length(vlist))
-        order, parent =
-            reverse_topological_sort_by_dfs(vec, offset, length(vlist), cmap)
-        for k in eachindex(vlist)
-            revmap[vlist[k]] = 0  # clear for reuse
-        end
-        postorder[idx] = order
-        parents[idx] = parent
-        vertexmap[idx] = vlist
-    end
-    return RecoveryInfo(
-        vertexmap,
-        postorder,
-        parents,
-        color,
-        num_colors,
-        _num_edges(g),
-        local_indices,
-    )
-end
-
-_normalize(i, j) = (j > i) ? (j, i) : (i, j)
-
-function _indirect_recover_structure(rinfo::RecoveryInfo)
-    N = length(rinfo.color)
-    I = zeros(Int, rinfo.nnz + N)
-    J = zeros(Int, rinfo.nnz + N)
-    k = 0
-    for i in 1:N
-        k += 1
-        I[k] = J[k] = i
-    end
-    for t in eachindex(rinfo.postorder)
-        vmap = rinfo.vertexmap[t]
-        order = rinfo.postorder[t]
-        parent = rinfo.parents[t]
-        for z in eachindex(order)
-            v = order[z]
-            p = parent[v]
-            if p == 0
-                continue
-            end
-            k += 1
-            I[k], J[k] = _normalize(vmap[v], vmap[p])
-        end
-    end
-    @assert k == rinfo.nnz + N
-    return I, J
-end
-
-"""
-    hessian_color_preprocess(
-        edgelist,
-        num_total_var,
-        seen_idx = IndexedSet(0),
-    )
-
-`edgelist` contains the nonzeros in the Hessian, *including* nonzeros on the
-diagonal.
-"""
-function hessian_color_preprocess(
-    edgelist,
-    num_total_var,
-    seen_idx = IndexedSet(0),
-)
-    resize!(seen_idx, num_total_var)
-    I, J = Int[], Int[]
-    for (i, j) in edgelist
-        push!(seen_idx, i)
-        push!(seen_idx, j)
-        push!(I, i)
-        push!(J, j)
-    end
-    local_indices = sort!(collect(seen_idx))
-    empty!(seen_idx)
-    global_to_local_idx = seen_idx.nzidx # steal for storage
-    for k in eachindex(local_indices)
-        global_to_local_idx[local_indices[k]] = k
-    end
-    # only do the coloring on the local indices
-    for k in eachindex(I)
-        I[k] = global_to_local_idx[I[k]]
-        J[k] = global_to_local_idx[J[k]]
-    end
-    g = UndirectedGraph(I, J, length(local_indices))
-    color, num_colors = acyclic_coloring(g)
-    @assert length(color) == _num_vertices(g)
-    rinfo = recovery_preprocess(g, color, num_colors, local_indices)
-    I, J = _indirect_recover_structure(rinfo)
-    # convert back to global indices
-    for k in eachindex(I)
-        I[k] = local_indices[I[k]]
-        J[k] = local_indices[J[k]]
-    end
-    return I, J, rinfo
-end
-
-"""
-    seed_matrix(rinfo::RecoveryInfo)
-
-Allocate a seed matrix for the Coloring.
-"""
-function seed_matrix(rinfo::RecoveryInfo)
-    return Matrix{Float64}(undef, length(rinfo.local_indices), rinfo.num_colors)
-end
-
-"""
-    prepare_seed_matrix!(R, rinfo::RecoveryInfo)
-"""
-function prepare_seed_matrix!(R, rinfo::RecoveryInfo)
-    N = length(rinfo.color)
-    @assert N == size(R, 1) == length(rinfo.local_indices)
-    @assert size(R, 2) == rinfo.num_colors
-    fill!(R, 0.0)
-    for i in 1:N
-        R[i, rinfo.color[i]] = 1
-    end
-    return
-end
-
-"""
-    recover_from_matmat!(
-        V::AbstractVector{T},
-        R::AbstractMatrix{T},
-        rinfo::RecoveryInfo,
-        stored_values::AbstractVector{T},
-    ) where {T}
-
-recover the hessian values from the hessian-matrix solution
-stored_values is a temporary vector of length >= length(rinfo.local_indices)
-"""
-function recover_from_matmat!(
-    V::AbstractVector{T},
-    R::AbstractMatrix{T},
-    rinfo::RecoveryInfo,
-    stored_values::AbstractVector{T},
-) where {T}
-    N = length(rinfo.color) # number of local variables
-    @assert length(V) == rinfo.nnz + N
-    @assert length(stored_values) >= length(rinfo.local_indices)
-    k = 0
-    for i in 1:N
-        k += 1
-        V[k] = R[i, rinfo.color[i]]
-    end
-    for t in eachindex(rinfo.vertexmap)
-        vmap = rinfo.vertexmap[t]
-        order = rinfo.postorder[t]
-        parent = rinfo.parents[t]
-        for z in eachindex(order)
-            @inbounds stored_values[z] = 0.0
-        end
-        @inbounds for z in eachindex(order)
-            v = order[z]
-            p = parent[v]
-            if p == 0
-                continue
-            end
-            i, j = vmap[v], vmap[p]
-            value = R[i, rinfo.color[j]] - stored_values[v]
-            stored_values[p] += value
-            k += 1
-            V[k] = value
-        end
-    end
-    @assert k == rinfo.nnz + N
-    return
-end
-
-end  # module
+end # module Coloring
diff --git a/src/Nonlinear/ReverseAD/Coloring/IntDisjointSet.jl b/src/Nonlinear/ReverseAD/Coloring/IntDisjointSet.jl
deleted file mode 100644
index 4fb6ea261d..0000000000
--- a/src/Nonlinear/ReverseAD/Coloring/IntDisjointSet.jl
+++ /dev/null
@@ -1,56 +0,0 @@
-# Copyright (c) 2017: Miles Lubin and contributors
-# Copyright (c) 2017: Google Inc.
-# Copyright (c) 2024: Guillaume Dalle and Alexis Montoison
-#
-# Use of this source code is governed by an MIT-style license that can be found
-# in the LICENSE.md file or at https://opensource.org/licenses/MIT.
-
-# The code in this file was taken from
-# https://github.com/gdalle/SparseMatrixColorings.jl/blob/main/src/Forest.jl
-#
-# It was copied at the suggestion of Alexis in his JuMP-dev 2025 talk.
-#
-# @odow made minor changes to match MOI coding styles.
-#
-# x-ref https://github.com/gdalle/SparseMatrixColorings.jl/pull/190
-
-mutable struct _IntDisjointSet
-    # current number of distinct trees in the S
-    number_of_trees::Int
-    # vector storing the index of a parent in the tree for each edge, used in
-    # union-find operations
-    parents::Vector{Int}
-    # vector approximating the depth of each tree to optimize path compression
-    ranks::Vector{Int}
-
-    _IntDisjointSet(n::Integer) = new(n, collect(1:n), zeros(Int, n))
-end
-
-function _find_root!(S::_IntDisjointSet, x::Integer)
-    p = S.parents[x]
-    if S.parents[p] != p
-        S.parents[x] = p = _find_root!(S, p)
-    end
-    return p
-end
-
-function _root_union!(S::_IntDisjointSet, x::Int, y::Int)
-    rank1, rank2 = S.ranks[x], S.ranks[y]
-    if rank1 < rank2
-        x, y = y, x
-    elseif rank1 == rank2
-        S.ranks[x] += 1
-    end
-    S.parents[y] = x
-    S.number_of_trees -= 1
-    return
-end
-
-function _union!(S, x::Int, y::Int)
-    root_x = _find_root!(S, x)
-    root_y = _find_root!(S, y)
-    if root_x != root_y
-        _root_union!(S, root_x, root_y)
-    end
-    return
-end
diff --git a/src/Nonlinear/ReverseAD/Coloring/topological_sort.jl b/src/Nonlinear/ReverseAD/Coloring/topological_sort.jl
deleted file mode 100644
index 0b8d0ca2eb..0000000000
--- a/src/Nonlinear/ReverseAD/Coloring/topological_sort.jl
+++ /dev/null
@@ -1,106 +0,0 @@
-# Modified from Graphs.jl.
-# Workaround for edge_index not being O(1) on SimpleGraph.
-# edge_index was only needed to test for cycles, so
-# this implementation skips that check.
-
-# Graphs.jl is licensed under the MIT License:
-#
-# Copyright (c) 2012: John Myles White and other contributors.
-#
-# Permission is hereby granted, free of charge, to any person obtaining
-# a copy of this software and associated documentation files (the
-# "Software"), to deal in the Software without restriction, including
-# without limitation the rights to use, copy, modify, merge, publish,
-# distribute, sublicense, and/or sell copies of the Software, and to
-# permit persons to whom the Software is furnished to do so, subject to
-# the following conditions:
-#
-# The above copyright notice and this permission notice shall be
-# included in all copies or substantial portions of the Software.
-#
-# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
-# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
-# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
-# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
-# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
-# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
-# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
-
-mutable struct _TopologicalSortVisitor
-    vertices::Vector{Int}
-    parents::Vector{Int}
-
-    function _TopologicalSortVisitor(n::Int)
-        vs = Int[]
-        sizehint!(vs, n)
-        return new(vs, zeros(Int, n))
-    end
-end
-
-function _traverse_graph(
-    adjlist,
-    offsets,
-    s,
-    visitor::_TopologicalSortVisitor,
-    vertexcolormap,
-    vertex_stack,
-    index_stack,
-)
-    vertexcolormap[s] = 1
-    resize!(vertex_stack, 1)
-    vertex_stack[1] = s
-    resize!(index_stack, 1)
-    index_stack[1] = 1
-    while !isempty(vertex_stack)
-        u = pop!(vertex_stack)
-        out_idx = pop!(index_stack)
-        len_uegs = offsets[u+1] - offsets[u]
-        found_new_vertex = false
-        while out_idx <= len_uegs && !found_new_vertex
-            v = adjlist[offsets[u]+out_idx-1]
-            out_idx += 1
-            v_color = vertexcolormap[v]
-            if v_color == 0
-                found_new_vertex = true
-                vertexcolormap[v] = 1
-                push!(vertex_stack, u)
-                push!(index_stack, out_idx)
-                visitor.parents[v] = u
-                push!(vertex_stack, v)
-                push!(index_stack, 1)
-            end
-        end
-        if !found_new_vertex
-            push!(visitor.vertices, u)
-            vertexcolormap[u] = 2
-        end
-    end
-    return
-end
-
-function reverse_topological_sort_by_dfs(
-    adjlist,
-    offsets,
-    num_vertices,
-    cmap::Vector{Int},
-    vertex_stack::Vector{Int} = Int[],
-    index_stack::Vector{Int} = Int[],
-)
-    @assert length(cmap) == num_vertices
-    fill!(cmap, 0)
-    visitor = _TopologicalSortVisitor(num_vertices)
-    for s in 1:num_vertices
-        if cmap[s] == 0
-            _traverse_graph(
-                adjlist,
-                offsets,
-                s,
-                visitor,
-                cmap,
-                vertex_stack,
-                index_stack,
-            )
-        end
-    end
-    return visitor.vertices, visitor.parents
-end
diff --git a/src/Nonlinear/ReverseAD/ReverseAD.jl b/src/Nonlinear/ReverseAD/ReverseAD.jl
index c874f3e7f7..94c4c50df5 100644
--- a/src/Nonlinear/ReverseAD/ReverseAD.jl
+++ b/src/Nonlinear/ReverseAD/ReverseAD.jl
@@ -7,9 +7,11 @@
 module ReverseAD
 
 import ForwardDiff
+import LinearAlgebra
 import MathOptInterface as MOI
 import ..Nonlinear
 import SparseArrays
+import SparseMatrixColorings as SMC
 
 # Override basic math functions to return NaN instead of throwing errors.
 # This is what NLP solvers expect, and sometimes the results aren't needed
@@ -31,6 +33,7 @@ import NaNMath:
     sqrt
 
 include("Coloring/Coloring.jl")
+include("coloring.jl")
 include("graph_tools.jl")
 include("types.jl")
 include("utils.jl")
diff --git a/src/Nonlinear/ReverseAD/coloring.jl b/src/Nonlinear/ReverseAD/coloring.jl
new file mode 100644
index 0000000000..037448dd22
--- /dev/null
+++ b/src/Nonlinear/ReverseAD/coloring.jl
@@ -0,0 +1,161 @@
+# Copyright (c) 2017: Miles Lubin and contributors
+# Copyright (c) 2017: Google Inc.
+#
+# Use of this source code is governed by an MIT-style license that can be found
+# in the LICENSE.md file or at https://opensource.org/licenses/MIT.
+
+"""
+    struct ColoringResult{R<:SMC.AbstractColoringResult}
+        result::R
+        local_indices::Vector{Int}  # map from local to global indices
+    end
+
+Wrapper around AbstractColoringResult that maps local Hessian indices back
+to global variable indices.
+"""
+struct ColoringResult{R<:SMC.AbstractColoringResult}
+    result::R
+    local_indices::Vector{Int}   # map from local to global indices
+end
+
+"""
+    _hessian_color_preprocess(
+        edgelist,
+        num_total_var,
+        algo::SMC.GreedyColoringAlgorithm,
+        seen_idx = MOI.Nonlinear.ReverseAD.Coloring.IndexedSet(0),
+    )
+
+`edgelist` contains the nonzeros in the Hessian, *including* nonzeros on the
+diagonal. Diagonal entries absent from `edgelist` are not added to the pattern.
+
+Returns `(colptr, I, J, result)` where `colptr` contains the local lower-
+triangular CSC column pointers and `I`, `J` give the corresponding entries
+in global variable indices. `result` wraps an SMC coloring result, or is
+`nothing` if the pattern is empty.
+"""
+function _hessian_color_preprocess(
+    edgelist,
+    num_total_var,
+    algo::SMC.GreedyColoringAlgorithm,
+    seen_idx = MOI.Nonlinear.ReverseAD.Coloring.IndexedSet(0),
+)
+    resize!(seen_idx, num_total_var)
+    I, J = Int[], Int[]
+    for (i, j) in edgelist
+        push!(seen_idx, i)
+        push!(seen_idx, j)
+        push!(I, max(i, j))
+        push!(J, min(i, j))
+    end
+    local_indices = sort!(collect(seen_idx))
+    empty!(seen_idx)
+    # SMC's star coloring does not support an empty pattern.
+    # TODO: remove when https://github.com/JuliaDiff/SparseMatrixColorings.jl/pull/352 lands
+    if isempty(local_indices)
+        return [1], Int[], Int[], nothing
+    end
+
+    global_to_local_idx = seen_idx.nzidx # steal for storage
+    for k in eachindex(local_indices)
+        global_to_local_idx[local_indices[k]] = k
+    end
+    for k in eachindex(I)
+        I[k] = global_to_local_idx[I[k]]
+        J[k] = global_to_local_idx[J[k]]
+    end
+
+    n = length(local_indices)
+    mat = SparseArrays.sparse(I, J, trues(length(I)), n, n, |)
+    # Color the full symmetric pattern, but recover only its lower triangle.
+    S = SMC.SparsityPatternCSC(
+        SparseArrays.sparse(LinearAlgebra.Symmetric(mat, :L)),
+    )
+    problem = SMC.ColoringProblem(; structure = :symmetric, partition = :column)
+    tree_result = SMC.coloring(S, problem, algo)
+
+    I, J, _ = SparseArrays.findnz(mat)
+    I_global = local_indices[I]
+    J_global = local_indices[J]
+    result = ColoringResult(tree_result, local_indices)
+    return mat.colptr, I_global, J_global, result
+end
+
+"""
+    _seed_matrix(result::ColoringResult)
+
+Allocate a seed matrix for the coloring result.
+"""
+function _seed_matrix(result::ColoringResult)
+    n = length(result.local_indices)
+    ncolors = SMC.ncolors(result.result)
+    return Matrix{Float64}(undef, n, ncolors)
+end
+
+function _seed_matrix(::Nothing)
+    return Matrix{Float64}(undef, 0, 0)
+end
+
+"""
+    _prepare_seed_matrix!(R, result::ColoringResult)
+
+Prepare the seed matrix R for Hessian computation.
+"""
+function _prepare_seed_matrix!(R, result::ColoringResult)
+    color = SMC.column_colors(result.result)
+    N = length(result.local_indices)
+    @assert N == size(R, 1)
+    @assert size(R, 2) == SMC.ncolors(result.result)
+    fill!(R, 0.0)
+    for i in 1:N
+        if color[i] > 0
+            R[i, color[i]] = 1
+        end
+    end
+    return
+end
+
+"""
+    _recover_from_matmat!(
+        colptr::AbstractVector,
+        V::AbstractVector{T},
+        R::AbstractMatrix{T},
+        result::ColoringResult,
+    ) where {T}
+
+Recover the Hessian values from the Hessian-matrix product H*R_seed.
+R is the result of H*R_seed where R_seed is the seed matrix.
+"""
+function _recover_from_matmat!(
+    colptr::AbstractVector,
+    V::AbstractVector{T},
+    R::AbstractMatrix{T},
+    result::ColoringResult,
+) where {T}
+    SMC.decompress_csc!(V, colptr, R, result.result, :L)
+    return
+end
+
+function _recover_from_matmat!(
+    colptr::AbstractVector,
+    V::AbstractVector{T},
+    R::AbstractMatrix{T},
+    result::ColoringResult{<:SMC.StarSetColoringResult},
+) where {T}
+    # `SMC.decompress_csc!` only supports substitution decompression. Mirror
+    # `SMC.decompress!` for a lower-triangular CSC matrix without constructing
+    # it.
+    S = SMC.sparsity_pattern(result.result)
+    rows = SparseArrays.rowvals(S)
+    compressed_indices = result.result.compressed_indices
+    l = 0
+    for j in axes(S, 2)
+        for k in SparseArrays.nzrange(S, j)
+            if rows[k] >= j
+                l += 1
+                V[l] = R[compressed_indices[k]]
+            end
+        end
+    end
+    return
+end
diff --git a/src/Nonlinear/ReverseAD/forward_over_reverse.jl b/src/Nonlinear/ReverseAD/forward_over_reverse.jl
index 8c4b6feb13..db765707b1 100644
--- a/src/Nonlinear/ReverseAD/forward_over_reverse.jl
+++ b/src/Nonlinear/ReverseAD/forward_over_reverse.jl
@@ -38,12 +38,11 @@ function _eval_hessian(
     scale::Float64,
     nzcount::Int,
 )::Int
-    if ex.linearity == LINEAR
-        @assert length(ex.hess_I) == 0
+    if isempty(ex.hess_I)
         return 0
     end
     chunk = min(size(ex.seed_matrix, 2), d.max_chunk)
-    Coloring.prepare_seed_matrix!(ex.seed_matrix, ex.rinfo)
+    _prepare_seed_matrix!(ex.seed_matrix, ex.rinfo)
     # Compute hessian-vector products
     num_products = size(ex.seed_matrix, 2) # number of hessian-vector products
     num_chunks = div(num_products, chunk)
@@ -66,11 +65,11 @@ function _eval_hessian(
     end
     # TODO(odow): consider reverting to a view.
     output_slice = _UnsafeVectorView(nzcount, length(ex.hess_I), pointer(H))
-    Coloring.recover_from_matmat!(
+    _recover_from_matmat!(
+        ex.hess_colptr,
         output_slice,
         ex.seed_matrix,
         ex.rinfo,
-        d.output_ϵ,
     )
     for i in 1:length(output_slice)
         output_slice[i] *= scale
diff --git a/src/Nonlinear/ReverseAD/mathoptinterface_api.jl b/src/Nonlinear/ReverseAD/mathoptinterface_api.jl
index f67d9caf8e..87af71c3fe 100644
--- a/src/Nonlinear/ReverseAD/mathoptinterface_api.jl
+++ b/src/Nonlinear/ReverseAD/mathoptinterface_api.jl
@@ -19,7 +19,10 @@ function MOI.features_available(d::NLPEvaluator)
     return [:Grad, :Jac, :JacVec, :Hess, :HessVec]
 end
 
-function MOI.initialize(d::NLPEvaluator, requested_features::Vector{Symbol})
+function MOI.initialize(
+    d::NLPEvaluator{R},
+    requested_features::Vector{Symbol},
+) where {R}
     # Check that we support the features requested by the user.
     available_features = MOI.features_available(d)
     for feature in requested_features
@@ -38,7 +41,7 @@ function MOI.initialize(d::NLPEvaluator, requested_features::Vector{Symbol})
     d.objective = nothing
     d.user_output_buffer = zeros(largest_user_input_dimension)
     d.jac_storage = zeros(max(N, largest_user_input_dimension))
-    d.constraints = _FunctionStorage[]
+    d.constraints = _FunctionStorage{R}[]
     d.last_x = fill(NaN, N)
     d.want_hess = :Hess in requested_features
     want_hess_storage = (:HessVec in requested_features) || d.want_hess
@@ -101,12 +104,12 @@ function MOI.initialize(d::NLPEvaluator, requested_features::Vector{Symbol})
     end
     max_chunk = 1
     if d.data.objective !== nothing
-        objective = _FunctionStorage(
+        objective = _FunctionStorage{R}(
             main_expressions[1],
             something(d.data.objective).values,
             N,
             coloring_storage,
-            d.want_hess,
+            d.want_hess ? d.coloring_algorithm : nothing,
             d.subexpressions,
             individual_order[1],
             d.subexpression_linearity,
@@ -122,12 +125,12 @@ function MOI.initialize(d::NLPEvaluator, requested_features::Vector{Symbol})
         idx = d.data.objective !== nothing ? k + 1 : k
         push!(
             d.constraints,
-            _FunctionStorage(
+            _FunctionStorage{R}(
                 main_expressions[idx],
                 constraint.expression.values,
                 N,
                 coloring_storage,
-                d.want_hess,
+                d.want_hess ? d.coloring_algorithm : nothing,
                 d.subexpressions,
                 individual_order[idx],
                 d.subexpression_linearity,
diff --git a/src/Nonlinear/ReverseAD/types.jl b/src/Nonlinear/ReverseAD/types.jl
index ce3e807b94..bd652f91bd 100644
--- a/src/Nonlinear/ReverseAD/types.jl
+++ b/src/Nonlinear/ReverseAD/types.jl
@@ -42,7 +42,7 @@ struct _SubexpressionStorage
     end
 end
 
-struct _FunctionStorage
+struct _FunctionStorage{R<:SMC.AbstractColoringResult}
     nodes::Vector{Nonlinear.Node}
     adj::SparseArrays.SparseMatrixCSC{Bool,Int}
     const_values::Vector{Float64}
@@ -51,27 +51,28 @@ struct _FunctionStorage
     reverse_storage::Vector{Float64}
     grad_sparsity::Vector{Int}
     # Nonzero pattern of Hessian matrix
+    hess_colptr::Vector{Int}
     hess_I::Vector{Int}
     hess_J::Vector{Int}
-    rinfo::Coloring.RecoveryInfo # coloring info for hessians
+    rinfo::Union{Nothing,ColoringResult{R}}
     seed_matrix::Matrix{Float64}
     linearity::Linearity
     # subexpressions which this function depends on, ordered for forward pass.
     dependent_subexpressions::Vector{Int}
 
-    function _FunctionStorage(
+    function _FunctionStorage{R}(
         nodes::Vector{Nonlinear.Node},
         const_values,
         num_variables,
         coloring_storage::Coloring.IndexedSet,
-        want_hess::Bool,
+        coloring_algorithm::Union{Nothing,SMC.GreedyColoringAlgorithm},
         subexpressions::Vector{_SubexpressionStorage},
         dependent_subexpressions,
         subexpression_linearity,
         subexpression_edgelist,
         subexpression_variables,
         moi_index_to_consecutive_index,
-    )
+    ) where {R<:SMC.AbstractColoringResult}
         nodes = _replace_moi_variables(nodes, moi_index_to_consecutive_index)
         adj = Nonlinear.adjacency_matrix(nodes)
         N = length(nodes)
@@ -85,7 +86,7 @@ struct _FunctionStorage
         end
         grad_sparsity = sort!(collect(coloring_storage))
         empty!(coloring_storage)
-        if want_hess
+        if !isnothing(coloring_algorithm)
             linearity = _classify_linearity(nodes, adj, subexpression_linearity)
             edgelist = _compute_hessian_sparsity(
                 nodes,
@@ -94,13 +95,14 @@ struct _FunctionStorage
                 subexpression_edgelist,
                 subexpression_variables,
             )
-            hess_I, hess_J, rinfo = Coloring.hessian_color_preprocess(
+            hess_colptr, hess_I, hess_J, rinfo = _hessian_color_preprocess(
                 edgelist,
                 num_variables,
+                coloring_algorithm,
                 coloring_storage,
             )
-            seed_matrix = Coloring.seed_matrix(rinfo)
-            return new(
+            seed_matrix = _seed_matrix(rinfo)
+            return new{R}(
                 nodes,
                 adj,
                 const_values,
@@ -108,6 +110,7 @@ struct _FunctionStorage
                 zeros(N),  # partials_storage,
                 zeros(N),  # reverse_storage,
                 grad_sparsity,
+                hess_colptr,
                 hess_I,
                 hess_J,
                 rinfo,
@@ -116,7 +119,7 @@ struct _FunctionStorage
                 dependent_subexpressions,
             )
         else
-            return new(
+            return new{R}(
                 nodes,
                 adj,
                 const_values,
@@ -126,7 +129,8 @@ struct _FunctionStorage
                 grad_sparsity,
                 Int[],
                 Int[],
-                Coloring.RecoveryInfo(),
+                Int[],
+                nothing,
                 Array{Float64}(undef, 0, 0),
                 NONLINEAR,
                 dependent_subexpressions,
@@ -147,12 +151,14 @@ interface.
 !!! warning
     Before using, you must initialize the evaluator using `MOI.initialize`.
 """
-mutable struct NLPEvaluator <: MOI.AbstractNLPEvaluator
+mutable struct NLPEvaluator{R,C<:SMC.GreedyColoringAlgorithm} <:
+               MOI.AbstractNLPEvaluator
     data::Nonlinear.Model
     ordered_variables::Vector{MOI.VariableIndex}
+    coloring_algorithm::C
 
-    objective::Union{Nothing,_FunctionStorage}
-    constraints::Vector{_FunctionStorage}
+    objective::Union{Nothing,_FunctionStorage{R}}
+    constraints::Vector{_FunctionStorage{R}}
     subexpressions::Vector{_SubexpressionStorage}
     subexpression_order::Vector{Int}
     # Storage for the subexpressions in reverse-mode automatic differentiation.
@@ -187,7 +193,17 @@ mutable struct NLPEvaluator <: MOI.AbstractNLPEvaluator
     function NLPEvaluator(
         data::Nonlinear.Model,
         ordered_variables::Vector{MOI.VariableIndex},
+        coloring_algorithm::SMC.GreedyColoringAlgorithm = SMC.GreedyColoringAlgorithm(;
+            decompression = :substitution,
+        ),
     )
-        return new(data, ordered_variables)
+        problem =
+            SMC.ColoringProblem(; structure = :symmetric, partition = :column)
+        C = typeof(coloring_algorithm)
+        _S = SMC.SparsityPatternCSC(
+            SparseArrays.sparse([1], [1], Bool[true], 1, 1),
+        )
+        R = typeof(SMC.coloring(_S, problem, coloring_algorithm))
+        return new{R,C}(data, ordered_variables, coloring_algorithm)
     end
 end
diff --git a/src/Nonlinear/types.jl b/src/Nonlinear/types.jl
index 36ac83a245..42a35418b3 100644
--- a/src/Nonlinear/types.jl
+++ b/src/Nonlinear/types.jl
@@ -276,20 +276,39 @@ function Evaluator(model::Model, ::ExprGraphOnly, ::Vector{MOI.VariableIndex})
 end
 
 """
-    SparseReverseMode() <: AbstractAutomaticDifferentiation
+    SparseReverseMode([coloring_algorithm]) <: AbstractAutomaticDifferentiation
 
 An implementation of `AbstractAutomaticDifferentiation` that uses sparse
 reverse-mode automatic differentiation to compute derivatives. Supports all
 features in the MOI nonlinear interface.
+
+The optional `coloring_algorithm` defaults to SparseMatrixColorings' greedy
+coloring with substitution decompression.
 """
-struct SparseReverseMode <: AbstractAutomaticDifferentiation end
+struct SparseReverseMode{C<:SMC.GreedyColoringAlgorithm} <:
+       AbstractAutomaticDifferentiation
+    coloring_algorithm::C
+end
+
+function SparseReverseMode()
+    return SparseReverseMode(
+        SMC.GreedyColoringAlgorithm(; decompression = :substitution),
+    )
+end
 
 function Evaluator(
     model::Model,
-    ::SparseReverseMode,
+    mode::SparseReverseMode,
     ordered_variables::Vector{MOI.VariableIndex},
 )
-    return Evaluator(model, ReverseAD.NLPEvaluator(model, ordered_variables))
+    return Evaluator(
+        model,
+        ReverseAD.NLPEvaluator(
+            model,
+            ordered_variables,
+            mode.coloring_algorithm,
+        ),
+    )
 end
 
 """
diff --git a/test/Nonlinear/test_ReverseAD.jl b/test/Nonlinear/test_ReverseAD.jl
index 688dede800..79cf165b0f 100644
--- a/test/Nonlinear/test_ReverseAD.jl
+++ b/test/Nonlinear/test_ReverseAD.jl
@@ -104,15 +104,15 @@ function test_objective_quadratic_multivariate()
     g = [NaN, NaN]
     MOI.eval_objective_gradient(evaluator, g, [1.2, 2.3])
     @test g == [2 * 1.2 + 2.3, 1.2 + 2 * 2.3]
-    @test MOI.hessian_objective_structure(evaluator) == [(1, 1), (2, 2), (2, 1)]
+    @test MOI.hessian_objective_structure(evaluator) == [(1, 1), (2, 1), (2, 2)]
     H = [NaN, NaN, NaN]
     MOI.eval_hessian_objective(evaluator, H, [1.2, 2.3])
-    @test H == [2.0, 2.0, 1.0]
+    @test H == [2.0, 1.0, 2.0]
     @test MOI.hessian_lagrangian_structure(evaluator) ==
-          [(1, 1), (2, 2), (2, 1)]
+          [(1, 1), (2, 1), (2, 2)]
     H = [NaN, NaN, NaN]
     MOI.eval_hessian_lagrangian(evaluator, H, [1.2, 2.3], 1.5, Float64[])
-    @test H == 1.5 .* [2.0, 2.0, 1.0]
+    @test H == 1.5 .* [2.0, 1.0, 2.0]
     v = [0.3, 0.4]
     hv = [NaN, NaN]
     MOI.eval_hessian_lagrangian_product(
@@ -145,19 +145,19 @@ function test_objective_quadratic_multivariate_subexpressions()
     MOI.eval_objective_gradient(evaluator, g, val)
     @test g == [2 * 1.2 + 2.3, 1.2 + 2 * 2.3]
     @test 0 == @allocated MOI.eval_objective_gradient(evaluator, g, val)
-    @test MOI.hessian_objective_structure(evaluator) == [(1, 1), (2, 2), (2, 1)]
+    @test MOI.hessian_objective_structure(evaluator) == [(1, 1), (2, 1), (2, 2)]
     H = [NaN, NaN, NaN]
     MOI.eval_hessian_objective(evaluator, H, val)
-    @test H == [2.0, 2.0, 1.0]
+    @test H == [2.0, 1.0, 2.0]
     @test evaluator.backend.max_chunk == 2
     @test 0 == @allocated MOI.eval_hessian_objective(evaluator, H, val)
     @test MOI.hessian_lagrangian_structure(evaluator) ==
-          [(1, 1), (2, 2), (2, 1)]
+          [(1, 1), (2, 1), (2, 2)]
     H = [NaN, NaN, NaN]
     μ = Float64[]
     MOI.eval_hessian_lagrangian(evaluator, H, val, 1.5, μ)
     @test 0 == @allocated MOI.eval_hessian_lagrangian(evaluator, H, val, 1.5, μ)
-    @test H == 1.5 .* [2.0, 2.0, 1.0]
+    @test H == 1.5 .* [2.0, 1.0, 2.0]
     v = [0.3, 0.4]
     hv = [NaN, NaN]
     MOI.eval_hessian_lagrangian_product(evaluator, hv, val, v, 1.5, μ)
@@ -281,10 +281,10 @@ function test_constraint_quadratic_multivariate()
     MOI.eval_constraint_jacobian(evaluator, J, x_val)
     @test J == [2 * x_val[1] + x_val[2], x_val[1] + 2 * x_val[2]]
     @test MOI.hessian_lagrangian_structure(evaluator) ==
-          [(1, 1), (2, 2), (2, 1)]
+          [(1, 1), (2, 1), (2, 2)]
     H = [NaN, NaN, NaN]
     MOI.eval_hessian_lagrangian(evaluator, H, x_val, 0.0, [1.5])
-    @test H == 1.5 .* [2.0, 2.0, 1.0]
+    @test H == 1.5 .* [2.0, 1.0, 2.0]
     return
 end
 
@@ -319,10 +319,10 @@ function test_constraint_quadratic_multivariate_subexpressions()
     @test y ≈ wJ[:]
     # Hessian-Lagrangian
     @test MOI.hessian_lagrangian_structure(evaluator) ==
-          [(1, 1), (2, 2), (2, 1)]
+          [(1, 1), (2, 1), (2, 2)]
     H = [NaN, NaN, NaN]
     MOI.eval_hessian_lagrangian(evaluator, H, x_val, 0.0, [1.5])
-    @test H ≈ 1.5 .* [2.0, 2.0, 1.0]
+    @test H ≈ 1.5 .* [2.0, 1.0, 2.0]
     return
 end
 
@@ -349,10 +349,10 @@ function test_hessian_sparsity_registered_function()
     @test :Hess in MOI.features_available(evaluator)
     MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
     @test MOI.hessian_lagrangian_structure(evaluator) ==
-          [(1, 1), (2, 2), (3, 3), (3, 1)]
+          [(1, 1), (3, 1), (2, 2), (3, 3)]
     H = fill(NaN, 4)
     MOI.eval_hessian_lagrangian(evaluator, H, rand(3), 1.5, Float64[])
-    @test H == 1.5 .* [2.0, 2.0, 2.0, 0.0]
+    @test H == 1.5 .* [2.0, 0.0, 2.0, 2.0]
     return
 end
 
@@ -380,10 +380,10 @@ function test_hessian_sparsity_registered_rosenbrock()
     @test :Hess in MOI.features_available(evaluator)
     MOI.initialize(evaluator, [:Grad, :Jac, :Hess])
     @test MOI.hessian_lagrangian_structure(evaluator) ==
-          [(1, 1), (2, 2), (2, 1)]
+          [(1, 1), (2, 1), (2, 2)]
     H = fill(NaN, 3)
     MOI.eval_hessian_lagrangian(evaluator, H, [1.0, 1.0], 1.5, Float64[])
-    @test H == 1.5 .* [802.0, 200.0, -400.0]
+    @test H == 1.5 .* [802.0, -400.0, 200.0]
     return
 end
 
@@ -417,85 +417,187 @@ function test_hessian_registered_error()
     return
 end
 
-struct _ColoringGraph
-    num_vertices::Int
-    edges::Vector{Tuple{Int,Int}}
-end
-
-function to_adjlist(graph::_ColoringGraph)
-    I = [i for (i, _) in graph.edges]
-    J = [j for (_, j) in graph.edges]
-    return Coloring.UndirectedGraph(I, J, graph.num_vertices)
-end
-
-function test_coloring_edge_free_graph()
-    graph = _ColoringGraph(10, [])
-    _, numcolors = Coloring.acyclic_coloring(to_adjlist(graph))
-    @test numcolors == 1
-    return
-end
-
-function test_coloring_one_edge_graph()
-    graph = _ColoringGraph(10, [(2, 4)])
-    color, numcolors = Coloring.acyclic_coloring(to_adjlist(graph))
-    @test numcolors == 2
-    @test color[2] != color[4]
+function test_mode_coloring()
+    algorithm = ReverseAD.SMC.GreedyColoringAlgorithm(
+        ReverseAD.SMC.LargestFirst();
+        decompression = :substitution,
+    )
+    for mode in
+        (Nonlinear.SparseReverseMode(), Nonlinear.SparseReverseMode(algorithm))
+        model = Nonlinear.Model()
+        x, y = MOI.VariableIndex.(1:2)
+        Nonlinear.set_objective(model, :($x * $y + $y^3))
+        evaluator = Nonlinear.Evaluator(model, mode, [x, y])
+        MOI.initialize(evaluator, [:Grad])
+        @test evaluator.backend.coloring_algorithm === mode.coloring_algorithm
+        @test evaluator.backend.objective.rinfo === nothing
+        @test MOI.eval_objective(evaluator, [3.0, 2.0]) == 14.0
+        g = zeros(2)
+        MOI.eval_objective_gradient(evaluator, g, [3.0, 2.0])
+        @test g == [2.0, 15.0]
+        MOI.initialize(evaluator, [:Hess])
+        @test MOI.hessian_objective_structure(evaluator) == [(2, 1), (2, 2)]
+        H = zeros(2)
+        MOI.eval_hessian_objective(evaluator, H, [3.0, 2.0])
+        @test H == [1.0, 12.0]
+    end
     return
 end
 
-function test_coloring_two_edge_graph()
-    graph = _ColoringGraph(10, [(2, 4), (2, 3)])
-    color, numcolors = Coloring.acyclic_coloring(to_adjlist(graph))
-    @test numcolors == 2
-    @test color[3] == color[4]
+function test_coloring_end_to_end_hessian_coloring_and_recovery()
+    # Test the new coloring API through the compatibility layer
+    coloring_algorithm =
+        ReverseAD.SMC.GreedyColoringAlgorithm(; decompression = :substitution)
+    for (edges, hess, expected_I, expected_J, expected_V) in (
+        (
+            Set([(1, 1), (1, 2), (2, 2)]),
+            [3.4 2.1; 2.1 1.3],
+            [1, 2, 2],
+            [1, 1, 2],
+            [3.4, 2.1, 1.3],
+        ),
+        (Set([(1, 2)]), [0.0 2.1; 2.1 0.0], [2], [1], [2.1]),
+        (Set([(1, 2), (2, 2)]), [0.0 2.1; 2.1 1.3], [2, 2], [1, 2], [2.1, 1.3]),
+    )
+        colptr, I, J, rinfo = ReverseAD._hessian_color_preprocess(
+            edges,
+            2,
+            coloring_algorithm,
+            MOI.Nonlinear.ReverseAD.Coloring.IndexedSet(0),
+        )
+        R = ReverseAD._seed_matrix(rinfo)
+        ReverseAD._prepare_seed_matrix!(R, rinfo)
+        @test I == expected_I
+        @test J == expected_J
+        @test R == [1.0 0.0; 0.0 1.0]
+        matmat = hess * R
+        V = zeros(length(I))
+        ReverseAD._recover_from_matmat!(colptr, V, matmat, rinfo)
+        @test V == expected_V
+    end
     return
 end
 
-function test_coloring_three_edge_graph()
-    graph = _ColoringGraph(10, [(2, 4), (2, 3), (3, 4)])
-    color, numcolors = Coloring.acyclic_coloring(to_adjlist(graph))
-    @test numcolors == 3
-    # TODO: What is this testing?
-    Coloring.recovery_preprocess(to_adjlist(graph), color, numcolors, Int[])
-    return
-end
+function test_hessian_structurally_zero_diagonal()
+    variables = MOI.VariableIndex.(1:4)
+    x, y = variables[2], variables[4]
+    for decompression in (:direct, :substitution),
+        (expr, structure) in (
+            (:($x * $y), [(4, 2)]),
+            (:($x * $y + $y^3), [(4, 2), (4, 4)]),
+            (:($x + $y), Tuple{Int,Int}[]),
+            (:(3.0), Tuple{Int,Int}[]),
+            (:(ifelse($x < 0, $x, $y)), Tuple{Int,Int}[]),
+        )
 
-function test_coloring_two_edge_three_vertex_graph()
-    graph = _ColoringGraph(3, [(1, 3), (2, 3)])
-    _, numcolors = Coloring.acyclic_coloring(to_adjlist(graph))
-    @test numcolors == 2
+        model = Nonlinear.Model()
+        Nonlinear.set_objective(model, expr)
+        algorithm = ReverseAD.SMC.GreedyColoringAlgorithm(; decompression)
+        evaluator = Nonlinear.Evaluator(
+            model,
+            Nonlinear.SparseReverseMode(algorithm),
+            variables,
+        )
+        MOI.initialize(evaluator, [:Hess])
+        @test MOI.hessian_objective_structure(evaluator) == structure
+        @test MOI.hessian_lagrangian_structure(evaluator) == structure
+        for y_value in (0.0, 2.0)
+            values = [0.0, 3.0, 0.0, y_value]
+            expected = if isempty(structure)
+                Float64[]
+            elseif length(structure) == 1
+                [1.0]
+            else
+                [1.0, 6 * y_value]
+            end
+            H = fill(NaN, length(structure))
+            MOI.eval_hessian_objective(evaluator, H, values)
+            @test H == expected
+            MOI.eval_hessian_lagrangian(evaluator, H, values, 1.5, Float64[])
+            @test H == 1.5 .* expected
+        end
+    end
     return
 end
 
-function test_coloring_four_edge_four_vertex_graph()
-    graph = _ColoringGraph(4, [(1, 2), (2, 3), (3, 4), (4, 1)])
-    _, numcolors = Coloring.acyclic_coloring(to_adjlist(graph))
-    @test numcolors == 3
-    return
-end
+function test_coloring_lower_triangular_recovery()
+    for (T, decompression) in
+        ((Float32, :direct), (Float64, :direct), (Float64, :substitution)),
+        order in (ReverseAD.SMC.NaturalOrder(), ReverseAD.SMC.LargestFirst())
 
-function test_coloring_topological_sort()
-    # graph = _ColoringGraph(6, [(1, 2), (1, 3), (1, 6), (2, 4), (2, 5)])
-    vec = [3, 6, 2, 1, 4, 5, 1, 2, 2, 1]
-    offset = [1, 4, 7, 8, 9, 10, 11]
-    v = Coloring.reverse_topological_sort_by_dfs(vec, offset, 6, zeros(Int, 6))
-    @test v[1] == [3, 6, 4, 5, 2, 1]
-    @test v[2] == [0, 1, 1, 2, 2, 1]
+        algorithm = ReverseAD.SMC.GreedyColoringAlgorithm(order; decompression)
+        for edges in (
+            Set([(2, 4), (4, 7), (7, 9)]),
+            Set([(4, 2), (2, 7), (9, 2), (7, 7)]),
+            Set([(4, 2), (4, 7), (7, 9), (9, 2), (4, 4), (9, 9)]),
+            Set([(2, 4), (4, 2), (2, 2)]),
+            Set([(2, 2), (9, 9)]),
+            Set{Tuple{Int,Int}}(),
+        )
+            H = zeros(T, 9, 9)
+            for (i, j) in edges
+                H[i, j] = H[j, i] = 10 * min(i, j) + max(i, j)
+            end
+            colptr, I, J, result =
+                ReverseAD._hessian_color_preprocess(edges, 9, algorithm)
+            expected_structure = sort!(
+                unique([(max(i, j), min(i, j)) for (i, j) in edges]);
+                by = reverse,
+            )
+            @test collect(zip(I, J)) == expected_structure
+            R = Matrix{T}(ReverseAD._seed_matrix(result))
+            if isempty(expected_structure)
+                @test result === nothing
+                @test colptr == [1]
+                @test size(R) == (0, 0)
+                continue
+            end
+            ReverseAD._prepare_seed_matrix!(R, result)
+            local_indices = result.local_indices
+            compressed = H[local_indices, local_indices] * R
+            V = fill(T(NaN), length(I))
+            ReverseAD._recover_from_matmat!(colptr, V, compressed, result)
+            @test V == [H[i, j] for (i, j) in expected_structure]
+            @test 0 == @allocated ReverseAD._recover_from_matmat!(
+                colptr,
+                V,
+                compressed,
+                result,
+            )
+        end
+    end
     return
 end
 
-function test_coloring_end_to_end_hessian_coloring_and_recovery()
-    I, J, rinfo = Coloring.hessian_color_preprocess(Set([(1, 2)]), 2)
-    R = Coloring.seed_matrix(rinfo)
-    Coloring.prepare_seed_matrix!(R, rinfo)
-    @test I == [1, 2, 2]
-    @test J == [1, 2, 1]
-    @test R == [1.0 0.0; 0.0 1.0]
-    hess = [3.4 2.1; 2.1 1.3]
-    matmat = hess * R
-    V = zeros(3)
-    Coloring.recover_from_matmat!(V, matmat, rinfo, zeros(3))
-    @test V == [3.4, 1.3, 2.1]
+function test_hessian_star_coloring()
+    x, y = MOI.VariableIndex.(1:2)
+    model = Nonlinear.Model()
+    Nonlinear.set_objective(model, :((1 - $x)^2 + 100 * ($y - $x^2)^2))
+    Nonlinear.add_constraint(model, :($x * $y), MOI.LessThan(2.0))
+    Nonlinear.add_constraint(model, :($y^2), MOI.LessThan(2.0))
+    algorithm = ReverseAD.SMC.GreedyColoringAlgorithm(; decompression = :direct)
+    evaluator = Nonlinear.Evaluator(
+        model,
+        Nonlinear.SparseReverseMode(algorithm),
+        [x, y],
+    )
+    MOI.initialize(evaluator, [:Hess])
+    @test MOI.hessian_lagrangian_structure(evaluator) ==
+          [(1, 1), (2, 1), (2, 2), (2, 1), (2, 2)]
+    H = zeros(5)
+    μ = [2.0, 3.0]
+    for xv in ([1.0, 1.0], [-1.2, 1.0])
+        MOI.eval_hessian_lagrangian(evaluator, H, xv, 1.5, μ)
+        @test H ≈ [
+            1.5 * (2 + 1200 * xv[1]^2 - 400 * xv[2]),
+            1.5 * (-400 * xv[1]),
+            300.0,
+            2.0,
+            6.0,
+        ]
+        @test 0 ==
+              @allocated MOI.eval_hessian_lagrangian(evaluator, H, xv, 1.5, μ)
+    end
     return
 end
 
@@ -1415,46 +1517,27 @@ function test_hessian_reinterpret_unsafe()
     H = zeros(length(H_s))
     x_v = ones(5)
     MOI.eval_hessian_lagrangian(evaluator, H, x_v, 0.0, [1.0, 1.0])
-    @test count(isapprox.(H, 1.0; atol = 1e-8)) == 3
-    @test count(isapprox.(H, 0.0; atol = 1e-8)) == 5
-    @test sort(H_s[round.(Bool, H)]) == [(3, 1), (3, 2), (5, 4)]
-    return
-end
-
-function test_IntDisjointSet()
-    for case in [
-        [(1, 2) => [1, 1, 3], (1, 3) => [1, 1, 1]],
-        [(1, 2) => [1, 1, 3], (3, 1) => [1, 1, 1]],
-        [(2, 1) => [2, 2, 3], (1, 3) => [2, 2, 2]],
-        [(2, 1) => [2, 2, 3], (3, 1) => [3, 2, 3]],
-        [(1, 3) => [1, 2, 1], (2, 3) => [1, 2, 2]],
-        [(1, 3) => [1, 2, 1], (3, 2) => [1, 1, 1]],
-        [(3, 1) => [3, 2, 3], (2, 3) => [3, 3, 3]],
-        [(3, 1) => [3, 2, 3], (3, 2) => [3, 3, 3]],
-        [(2, 3) => [1, 2, 2], (1, 3) => [1, 2, 1]],
-        [(2, 3) => [1, 2, 2], (3, 1) => [2, 2, 2]],
-        [(3, 2) => [1, 3, 3], (1, 3) => [3, 3, 3]],
-        [(3, 2) => [1, 3, 3], (3, 1) => [3, 3, 3]],
-    ]
-        S = Coloring._IntDisjointSet(3)
-        @test Coloring._find_root!.((S,), [1, 2, 3]) == [1, 2, 3]
-        @test S.number_of_trees == 3
-        for (i, (union, result)) in enumerate(case)
-            Coloring._root_union!(S, union[1], union[2])
-            @test Coloring._find_root!.((S,), [1, 2, 3]) == result
-            @test S.number_of_trees == 3 - i
-        end
-    end
+    @test H ≈ ones(3)
+    @test sort(H_s) == [(3, 1), (3, 2), (5, 4)]
     return
 end
 
 function test_issue_2897()
     I = [4, 5, 4, 6, 5, 6]
     J = [2, 1, 1, 2, 3, 3]
-    g = Coloring.UndirectedGraph(I, J, length(I))
-    color, num_colors = Coloring.acyclic_coloring(g)
-    @test color == [1, 1, 1, 2, 2, 3]
-    @test num_colors == 3
+    algorithm =
+        ReverseAD.SMC.GreedyColoringAlgorithm(; decompression = :substitution)
+    colptr, I, J, result =
+        ReverseAD._hessian_color_preprocess(Set(zip(I, J)), 6, algorithm)
+    H = zeros(6, 6)
+    for (i, j) in zip(I, J)
+        H[i, j] = H[j, i] = 10 * j + i
+    end
+    R = ReverseAD._seed_matrix(result)
+    ReverseAD._prepare_seed_matrix!(R, result)
+    V = zeros(length(I))
+    ReverseAD._recover_from_matmat!(colptr, V, H * R, result)
+    @test V == [H[i, j] for (i, j) in zip(I, J)]
     return
 end
 
@@ -1574,11 +1657,11 @@ function test_hessian_pow_constant_exponent_nonpositive_base()
     )
     MOI.initialize(evaluator, [:Grad, :Hess])
     @test MOI.hessian_lagrangian_structure(evaluator) ==
-          [(1, 1), (2, 2), (2, 1)]
+          [(1, 1), (2, 1), (2, 2)]
     H = zeros(3)
     for x_val in (-1.5, 0.0)
         MOI.eval_hessian_lagrangian(evaluator, H, [x_val, 2.0], 1.0, Float64[])
-        @test isapprox(H, [4.0, 2 * x_val, 2 * x_val + 4.0])
+        @test isapprox(H, [4.0, 2 * x_val + 4.0, 2 * x_val])
     end
     return
 end