diff --git a/Project.toml b/Project.toml
index 178e656..71cab14 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "SparseMatrixColorings"
 uuid = "0a514795-09f3-496d-8182-132a7b665d35"
 authors = ["Guillaume Dalle", "Alexis Montoison"]
-version = "0.4.8"
+version = "0.4.9"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
diff --git a/docs/src/api.md b/docs/src/api.md
index 2912efd..f171ce6 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -26,6 +26,7 @@ ConstantColoringAlgorithm
 AbstractColoringResult
 column_colors
 row_colors
+ncolors
 column_groups
 row_groups
 sparsity_pattern
diff --git a/docs/src/dev.md b/docs/src/dev.md
index 5e1907a..84234f5 100644
--- a/docs/src/dev.md
+++ b/docs/src/dev.md
@@ -38,6 +38,8 @@ SparseMatrixColorings.RowColoringResult
 SparseMatrixColorings.StarSetColoringResult
 SparseMatrixColorings.TreeSetColoringResult
 SparseMatrixColorings.LinearSystemColoringResult
+SparseMatrixColorings.BicoloringResult
+SparseMatrixColorings.remap_colors
 ```
 
 ## Testing
diff --git a/ext/SparseMatrixColoringsColorsExt.jl b/ext/SparseMatrixColoringsColorsExt.jl
index 0984c19..a2827f7 100644
--- a/ext/SparseMatrixColoringsColorsExt.jl
+++ b/ext/SparseMatrixColoringsColorsExt.jl
@@ -21,7 +21,8 @@ using SparseMatrixColorings:
     AbstractColoringResult,
     sparsity_pattern,
     column_colors,
-    row_colors
+    row_colors,
+    ncolors
 using Colors: Colorant, RGB, RGBA, distinguishable_colors
 
 const DEFAULT_BACKGROUND = RGBA(0, 0, 0, 0)
@@ -31,9 +32,6 @@ const DEFAULT_PAD = 0 # update docstring in src/images.jl when changing this def
 # Sample n distinguishable colors, excluding the background color
 default_colorscheme(n, background) = distinguishable_colors(n, background; dropseed=true)
 
-ncolors(res::AbstractColoringResult{s,:column}) where {s} = maximum(column_colors(res))
-ncolors(res::AbstractColoringResult{s,:row}) where {s} = maximum(row_colors(res))
-
 ## Top-level function that handles argument errors, eagerly promotes types and allocates output buffer
 
 function SparseMatrixColorings.show_colors(
@@ -119,4 +117,30 @@ function show_colors!(
     return out
 end
 
+function show_colors!(
+    out, res::AbstractColoringResult{s,:bidirectional}, colorscheme, scale, pad
+) where {s}
+    scale < 3 && throw(ArgumentError("`scale` has to be ≥ 3 to visualize bicoloring"))
+    ccolor_indices = mod1.(column_colors(res), length(colorscheme)) # cycle color indices if necessary
+    row_shift = maximum(column_colors(res))
+    rcolor_indices = mod1.(row_shift .+ row_colors(res), length(colorscheme)) # cycle color indices if necessary
+    ccolors = colorscheme[ccolor_indices]
+    rcolors = colorscheme[rcolor_indices]
+    pattern = sparsity_pattern(res)
+    for I in CartesianIndices(pattern)
+        if !iszero(pattern[I])
+            r, c = Tuple(I)
+            area = matrix_entry_area(I, scale, pad)
+            for i in axes(area, 1), j in axes(area, 2)
+                if j > i
+                    out[area[i, j]] = ccolors[c]
+                elseif i > j
+                    out[area[i, j]] = rcolors[r]
+                end
+            end
+        end
+    end
+    return out
+end
+
 end # module
diff --git a/src/SparseMatrixColorings.jl b/src/SparseMatrixColorings.jl
index d4f5f3e..ea2e3b4 100644
--- a/src/SparseMatrixColorings.jl
+++ b/src/SparseMatrixColorings.jl
@@ -60,7 +60,7 @@ export DynamicDegreeBasedOrder, SmallestLast, IncidenceDegree, DynamicLargestFir
 export ColoringProblem, GreedyColoringAlgorithm, AbstractColoringResult
 export ConstantColoringAlgorithm
 export coloring
-export column_colors, row_colors
+export column_colors, row_colors, ncolors
 export column_groups, row_groups
 export sparsity_pattern
 export compress, decompress, decompress!, decompress_single_color!
diff --git a/src/decompression.jl b/src/decompression.jl
index b1452d4..6cba735 100644
--- a/src/decompression.jl
+++ b/src/decompression.jl
@@ -9,9 +9,6 @@ Compress `A` given a coloring `result` of the sparsity pattern of `A`.
 Compression means summing either the columns or the rows of `A` which share the same color.
 It is undone by calling [`decompress`](@ref) or [`decompress!`](@ref).
 
-!!! warning
-    At the moment, `:bidirectional` partitions are not implemented.
-
 # Example
 
 ```jldoctest
@@ -63,10 +60,25 @@ function compress(A, result::AbstractColoringResult{structure,:row}) where {stru
     return B
 end
 
+function compress(
+    A, result::AbstractColoringResult{structure,:bidirectional}
+) where {structure}
+    row_group = row_groups(result)
+    column_group = column_groups(result)
+    Br = stack(row_group; dims=1) do g
+        dropdims(sum(A[g, :]; dims=1); dims=1)
+    end
+    Bc = stack(column_group; dims=2) do g
+        dropdims(sum(A[:, g]; dims=2); dims=2)
+    end
+    return Br, Bc
+end
+
 """
-    decompress(B::AbstractMatrix, result::AbstractColoringResult)
+    decompress(B::AbstractMatrix, result::AbstractColoringResult{_,:column/:row})
+    decompress(Br::AbstractMatrix, Bc::AbstractMatrix, result::AbstractColoringResult{_,:bidirectional})
 
-Decompress `B` into a new matrix `A`, given a coloring `result` of the sparsity pattern of `A`.
+Decompress `B` (or the tuple `(Br,Bc)`) into a new matrix `A`, given a coloring `result` of the sparsity pattern of `A`.
 The in-place alternative is [`decompress!`](@ref).
 
 Compression means summing either the columns or the rows of `A` which share the same color.
@@ -120,13 +132,27 @@ function decompress(B::AbstractMatrix, result::AbstractColoringResult)
     return decompress!(A, B, result)
 end
 
+function decompress(
+    Br::AbstractMatrix,
+    Bc::AbstractMatrix,
+    result::AbstractColoringResult{structure,:bidirectional},
+) where {structure}
+    A = respectful_similar(result.A, Base.promote_eltype(Br, Bc))
+    return decompress!(A, Br, Bc, result)
+end
+
 """
     decompress!(
         A::AbstractMatrix, B::AbstractMatrix,
-        result::AbstractColoringResult, [uplo=:F]
+        result::AbstractColoringResult{_,:column/:row}, [uplo=:F]
+    )
+
+    decompress!(
+        A::AbstractMatrix, Br::AbstractMatrix, Bc::AbstractMatrix
+        result::AbstractColoringResult{_,:bidirectional}
     )
 
-Decompress `B` in-place into `A`, given a coloring `result` of the sparsity pattern of `A`.
+Decompress `B` (or the tuple `(Br,Bc)`) in-place into `A`, given a coloring `result` of the sparsity pattern of `A`.
 The out-of-place alternative is [`decompress`](@ref).
 
 !!! note
@@ -632,3 +658,60 @@ function decompress!(
     end
     return A
 end
+
+## BicoloringResult
+
+function _join_compressed!(result::BicoloringResult, Br::AbstractMatrix, Bc::AbstractMatrix)
+    #=
+    Say we have an original matrix `A` of size `(n, m)` and we build an augmented matrix `A_and_Aᵀ = [zeros(n, n) Aᵀ; A zeros(m, m)]`.
+    Its first `1:n` columns have the form `[zeros(n); A[:, j]]` and its following `n+1:n+m` columns have the form `[A[i, :]; zeros(m)]`.
+    The symmetric column coloring is performed on `A_and_Aᵀ` and the column-wise compression of `A_and_Aᵀ` should return a matrix `Br_and_Bc`.
+    But in reality, `Br_and_Bc` is computed as two partial compressions: the row-wise compression `Br` (corresponding to `Aᵀ`) and the columnwise compression `Bc` (corresponding to `A`).
+    Before symmetric decompression, we must reconstruct `Br_and_Bc` from `Br` and `Bc`, knowing that the symmetric colors (those making up `Br_and_Bc`) are present in either a row of `Br`, a column of `Bc`, or both.
+    Therefore, the column indices in `Br_and_Bc` don't necessarily match with the row indices in `Br` or the column indices in `Bc` since some colors may be missing in the partial compressions.
+    The columns of the top part of `Br_and_Bc` (rows `1:n`) are the rows of `Br`, interlaced with zero columns whenever the current color hasn't been used to color any row.
+    The columns of the bottom part of `Br_and_Bc` (rows `n+1:n+m`) are the columns of `Bc`, interlaced with zero columns whenever the current color hasn't been used to color any column.
+    We use the dictionaries `col_color_ind` and `row_color_ind` to map from symmetric colors to row/column colors.
+    =#
+    (; A, col_color_ind, row_color_ind) = result
+    m, n = size(A)
+    R = Base.promote_eltype(Br, Bc)
+    if eltype(result.Br_and_Bc) == R
+        Br_and_Bc = result.Br_and_Bc
+    else
+        Br_and_Bc = similar(result.Br_and_Bc, R)
+    end
+    fill!(Br_and_Bc, zero(R))
+    for c in axes(Br_and_Bc, 2)
+        if haskey(row_color_ind, c)  # some rows were colored with symmetric color c
+            copyto!(view(Br_and_Bc, 1:n, c), view(Br, row_color_ind[c], :))
+        end
+        if haskey(col_color_ind, c)  # some columns were colored with symmetric c
+            copyto!(view(Br_and_Bc, (n + 1):(n + m), c), view(Bc, :, col_color_ind[c]))
+        end
+    end
+    return Br_and_Bc
+end
+
+function decompress!(
+    A::AbstractMatrix, Br::AbstractMatrix, Bc::AbstractMatrix, result::BicoloringResult
+)
+    m, n = size(A)
+    Br_and_Bc = _join_compressed!(result, Br, Bc)
+    A_and_Aᵀ = decompress(Br_and_Bc, result.symmetric_result)
+    copyto!(A, A_and_Aᵀ[(n + 1):(n + m), 1:n])  # original matrix in bottom left corner
+    return A
+end
+
+function decompress!(
+    A::SparseMatrixCSC, Br::AbstractMatrix, Bc::AbstractMatrix, result::BicoloringResult
+)
+    (; large_colptr, large_rowval, symmetric_result) = result
+    m, n = size(A)
+    Br_and_Bc = _join_compressed!(result, Br, Bc)
+    # pretend A is larger
+    A_and_noAᵀ = SparseMatrixCSC(m + n, m + n, large_colptr, large_rowval, A.nzval)
+    # decompress lower triangle only
+    decompress!(A_and_noAᵀ, Br_and_Bc, symmetric_result, :L)
+    return A
+end
diff --git a/src/graph.jl b/src/graph.jl
index 65854f3..f4af4ae 100644
--- a/src/graph.jl
+++ b/src/graph.jl
@@ -145,17 +145,13 @@ function degree(g::AdjacencyGraph, v::Integer)
 end
 
 function nb_edges(g::AdjacencyGraph)
-    S = pattern(g)
     ne = 0
-    for j in vertices(g)
-        for k in nzrange(S, j)
-            i = rowvals(S)[k]
-            if i > j
-                ne += 1
-            end
+    for v in vertices(g)
+        for u in neighbors(g, v)
+            ne += 1
         end
     end
-    return ne
+    return ne ÷ 2
 end
 
 maximum_degree(g::AdjacencyGraph) = maximum(Base.Fix1(degree, g), vertices(g))
diff --git a/src/interface.jl b/src/interface.jl
index 9bc9363..f2b601f 100644
--- a/src/interface.jl
+++ b/src/interface.jl
@@ -1,6 +1,6 @@
 function check_valid_problem(structure::Symbol, partition::Symbol)
     valid = (
-        (structure == :nonsymmetric && partition in (:column, :row)) ||
+        (structure == :nonsymmetric && partition in (:column, :row, :bidirectional)) ||
         (structure == :symmetric && partition == :column)
     )
     if !valid
@@ -49,7 +49,7 @@ Matrix coloring is often used in automatic differentiation, and here is the tran
 | -------- | ------- | --------------- | ---------------- | ----------- |
 | Jacobian | forward | `:nonsymmetric` | `:column`        | yes         |
 | Jacobian | reverse | `:nonsymmetric` | `:row`           | yes         |
-| Jacobian | mixed   | `:nonsymmetric` | `:bidirectional` | no          |
+| Jacobian | mixed   | `:nonsymmetric` | `:bidirectional` | yes         |
 | Hessian  | -       | `:symmetric`    | `:column`        | yes         |
 | Hessian  | -       | `:symmetric`    | `:row`           | no          |
 """
@@ -223,6 +223,37 @@ function coloring(
     return TreeSetColoringResult(A, ag, color, tree_set, decompression_eltype)
 end
 
+function coloring(
+    A::AbstractMatrix,
+    ::ColoringProblem{:nonsymmetric,:bidirectional},
+    algo::GreedyColoringAlgorithm{decompression};
+    decompression_eltype::Type{R}=Float64,
+    symmetric_pattern::Bool=false,
+) where {decompression,R}
+    m, n = size(A)
+    T = eltype(A)
+    Aᵀ = if symmetric_pattern || A isa Union{Symmetric,Hermitian}
+        A
+    else
+        transpose(A)
+    end  # TODO: fuse with next step?
+    A_and_Aᵀ = [
+        spzeros(T, n, n) SparseMatrixCSC(Aᵀ)
+        SparseMatrixCSC(A) spzeros(T, m, m)
+    ]  # TODO: slow
+    ag = AdjacencyGraph(A_and_Aᵀ)
+    if decompression == :direct
+        color, star_set = star_coloring(ag, algo.order)
+        symmetric_result = StarSetColoringResult(A_and_Aᵀ, ag, color, star_set)
+    else
+        color, tree_set = acyclic_coloring(ag, algo.order)
+        symmetric_result = TreeSetColoringResult(
+            A_and_Aᵀ, ag, color, tree_set, decompression_eltype
+        )
+    end
+    return BicoloringResult(A, ag, symmetric_result, decompression_eltype)
+end
+
 ## ADTypes interface
 
 function ADTypes.column_coloring(A::AbstractMatrix, algo::GreedyColoringAlgorithm)
diff --git a/src/result.jl b/src/result.jl
index 46633f2..13a867d 100644
--- a/src/result.jl
+++ b/src/result.jl
@@ -55,6 +55,25 @@ Return a vector `group` such that for every color `c`, `group[c]` contains the i
 """
 function row_groups end
 
+"""
+    ncolors(result::AbstractColoringResult)
+
+Return the number of different colors used to color the matrix.
+
+For bidirectional partitions, this number is the sum of the number of row colors and the number of column colors.
+"""
+function ncolors(res::AbstractColoringResult{structure,:column}) where {structure}
+    return length(column_groups(res))
+end
+
+function ncolors(res::AbstractColoringResult{structure,:row}) where {structure}
+    return length(row_groups(res))
+end
+
+function ncolors(res::AbstractColoringResult{structure,:bidirectional}) where {structure}
+    return length(row_groups(res)) + length(column_groups(res))
+end
+
 """
     group_by_color(color::Vector{Int})
 
@@ -528,3 +547,107 @@ function LinearSystemColoringResult(
         T_factorization,
     )
 end
+
+## Bicoloring result
+
+"""
+    remap_colors(color::Vector{Int})
+
+Renumber the colors in `color` using their index in the vector `sort(unique(color))`, so that they are forced to go from `1` to some `cmax` contiguously.
+
+Return a tuple `(remapped_colors, color_to_ind)` such that `remapped_colors` is a vector containing the renumbered colors and `color_to_ind` is a dictionary giving the translation between old and new color numberings.
+
+For all vertex indices `i` we have:
+
+    remapped_color[i] = color_to_ind[color[i]]
+"""
+function remap_colors(color::Vector{Int})
+    color_to_ind = Dict(c => i for (i, c) in enumerate(sort(unique(color))))
+    remapped_colors = [color_to_ind[c] for c in color]
+    return remapped_colors, color_to_ind
+end
+
+"""
+$TYPEDEF
+
+Storage for the result of a bidirectional coloring with direct or substitution decompression, based on the symmetric coloring of a 2x2 block matrix.
+
+# Fields
+
+$TYPEDFIELDS
+
+# See also
+
+- [`AbstractColoringResult`](@ref)
+"""
+struct BicoloringResult{
+    M<:AbstractMatrix,
+    G<:AdjacencyGraph,
+    decompression,
+    V,
+    SR<:AbstractColoringResult{:symmetric,:column,decompression},
+    R,
+} <: AbstractColoringResult{:nonsymmetric,:bidirectional,decompression}
+    "matrix that was colored"
+    A::M
+    "adjacency graph that was used for coloring (constructed from the bipartite graph)"
+    abg::G
+    "one integer color for each column"
+    column_color::Vector{Int}
+    "one integer color for each row"
+    row_color::Vector{Int}
+    "color groups for columns"
+    column_group::V
+    "color groups for rows"
+    row_group::V
+    "result for the coloring of the symmetric 2x2 block matrix"
+    symmetric_result::SR
+    "column color to index"
+    col_color_ind::Dict{Int,Int}
+    "row color to index"
+    row_color_ind::Dict{Int,Int}
+    "combination of `Br` and `Bc` (almost a concatenation up to color remapping)"
+    Br_and_Bc::Matrix{R}
+    "CSC storage of `A_and_noAᵀ - `colptr`"
+    large_colptr::Vector{Int}
+    "CSC storage of `A_and_noAᵀ - `rowval`"
+    large_rowval::Vector{Int}
+end
+
+column_colors(result::BicoloringResult) = result.column_color
+column_groups(result::BicoloringResult) = result.column_group
+
+row_colors(result::BicoloringResult) = result.row_color
+row_groups(result::BicoloringResult) = result.row_group
+
+function BicoloringResult(
+    A::AbstractMatrix,
+    ag::AdjacencyGraph,
+    symmetric_result::AbstractColoringResult{:symmetric,:column},
+    decompression_eltype::Type{R},
+) where {R}
+    m, n = size(A)
+    symmetric_color = column_colors(symmetric_result)
+    column_color, col_color_ind = remap_colors(symmetric_color[1:n])
+    row_color, row_color_ind = remap_colors(symmetric_color[(n + 1):(n + m)])
+    column_group = group_by_color(column_color)
+    row_group = group_by_color(row_color)
+    Br_and_Bc = Matrix{R}(undef, n + m, maximum(column_colors(symmetric_result)))
+    large_colptr = copy(ag.S.colptr)
+    large_colptr[(n + 2):end] .= large_colptr[n + 1]  # last few columns are empty
+    large_rowval = ag.S.rowval[1:(end ÷ 2)]  # forget the second half of nonzeros
+    return BicoloringResult(
+        A,
+        ag,
+        column_color,
+        row_color,
+        column_group,
+        row_group,
+        symmetric_result,
+        col_color_ind,
+        row_color_ind,
+        Br_and_Bc,
+        large_colptr,
+        large_rowval,
+    )
+end
diff --git a/test/allocations.jl b/test/allocations.jl
index f0d3bbe..81272c9 100644
--- a/test/allocations.jl
+++ b/test/allocations.jl
@@ -28,49 +28,62 @@ function test_noallocs_sparse_decompression(
     result = coloring(
         A, ColoringProblem(; structure, partition), GreedyColoringAlgorithm(; decompression)
     )
-    B = compress(A, result)
 
-    @testset "Full decompression" begin
-        bench1_full = @be similar(A) decompress!(_, B, result) evals = 1
-        bench2_full = @be similar(Matrix(A)) decompress!(_, B, result) evals = 1
-        @test minimum(bench1_full).allocs == 0
-        @test minimum(bench2_full).allocs == 0
-    end
-    @testset "Single-color decompression" begin
-        if decompression == :direct
-            b = if partition == :column
-                B[:, 1]
-            else
-                B[1, :]
+    if partition == :bidirectional
+        Br, Bc = compress(A, result)
+        @testset "Full decompression" begin
+            bench1_full = @be similar(A) decompress!(_, Br, Bc, result) evals = 1
+            bench2_full = @be similar(Matrix(A)) decompress!(_, Br, Bc, result) evals = 1
+            @test minimum(bench1_full).allocs == 0
+            @test_broken minimum(bench2_full).allocs == 0
+        end
+    else
+        B = compress(A, result)
+        @testset "Full decompression" begin
+            bench1_full = @be similar(A) decompress!(_, B, result) evals = 1
+            bench2_full = @be similar(Matrix(A)) decompress!(_, B, result) evals = 1
+            @test minimum(bench1_full).allocs == 0
+            @test minimum(bench2_full).allocs == 0
+        end
+        @testset "Single-color decompression" begin
+            if decompression == :direct
+                b = if partition == :column
+                    B[:, 1]
+                else
+                    B[1, :]
+                end
+                bench1_singlecolor = @be similar(A) decompress_single_color!(
+                    _, b, 1, result
+                ) evals = 1
+                bench2_singlecolor = @be similar(Matrix(A)) decompress_single_color!(
+                    _, b, 1, result
+                ) evals = 1
+                @test minimum(bench1_singlecolor).allocs == 0
+                @test minimum(bench2_singlecolor).allocs == 0
             end
-            bench1_singlecolor = @be similar(A) decompress_single_color!(_, b, 1, result) evals =
-                1
-            bench2_singlecolor = @be similar(Matrix(A)) decompress_single_color!(
-                _, b, 1, result
-            ) evals = 1
-            @test minimum(bench1_singlecolor).allocs == 0
-            @test minimum(bench2_singlecolor).allocs == 0
         end
-    end
-    @testset "Triangle decompression" begin
-        if structure == :symmetric
-            bench1_triangle = @be similar(triu(A)) decompress!(_, B, result, :U) evals = 1
-            bench2_triangle = @be similar(Matrix(A)) decompress!(_, B, result, :U) evals = 1
-            @test minimum(bench1_triangle).allocs == 0
-            @test minimum(bench2_triangle).allocs == 0
+        @testset "Triangle decompression" begin
+            if structure == :symmetric
+                bench1_triangle = @be similar(triu(A)) decompress!(_, B, result, :U) evals =
+                    1
+                bench2_triangle = @be similar(Matrix(A)) decompress!(_, B, result, :U) evals =
+                    1
+                @test minimum(bench1_triangle).allocs == 0
+                @test minimum(bench2_triangle).allocs == 0
+            end
         end
-    end
-    @testset "Single-color triangle decompression" begin
-        if structure == :symmetric && decompression == :direct
-            b = B[:, 1]
-            bench1_singlecolor_triangle = @be similar(triu(A)) decompress_single_color!(
-                _, b, 1, result, :U
-            ) evals = 1
-            bench2_singlecolor_triangle = @be similar(Matrix(A)) decompress_single_color!(
-                _, b, 1, result, :U
-            ) evals = 1
-            @test minimum(bench1_singlecolor_triangle).allocs == 0
-            @test minimum(bench2_singlecolor_triangle).allocs == 0
+        @testset "Single-color triangle decompression" begin
+            if structure == :symmetric && decompression == :direct
+                b = B[:, 1]
+                bench1_singlecolor_triangle = @be similar(triu(A)) decompress_single_color!(
+                    _, b, 1, result, :U
+                ) evals = 1
+                bench2_singlecolor_triangle = @be similar(Matrix(A)) decompress_single_color!(
+                    _, b, 1, result, :U
+                ) evals = 1
+                @test minimum(bench1_singlecolor_triangle).allocs == 0
+                @test minimum(bench2_singlecolor_triangle).allocs == 0
+            end
         end
     end
 end
@@ -103,6 +116,8 @@ end
         (:nonsymmetric, :row, :direct),
         (:symmetric, :column, :direct),
         (:symmetric, :column, :substitution),
+        (:nonsymmetric, :bidirectional, :direct),
+        (:nonsymmetric, :bidirectional, :substitution),
     ]
         test_noallocs_sparse_decompression(1000; structure, partition, decompression)
     end
diff --git a/test/random.jl b/test/random.jl
index 82573b5..06e8565 100644
--- a/test/random.jl
+++ b/test/random.jl
@@ -68,3 +68,29 @@ end;
         test_coloring_decompression(A0, problem, algo)
     end
 end;
+
+@testset "Bicoloring & direct decompression" begin
+    problem = ColoringProblem(; structure=:nonsymmetric, partition=:bidirectional)
+    algo = GreedyColoringAlgorithm(RandomOrder(rng); decompression=:direct)
+    @testset "$((; m, n, p))" for (m, n, p) in asymmetric_params
+        A0 = sprand(rng, m, n, p)
+        test_bicoloring_decompression(A0, problem, algo)
+    end
+    @testset "$((; n, p))" for (n, p) in symmetric_params
+        A0 = sparse(Symmetric(sprand(rng, n, n, p)))
+        test_bicoloring_decompression(A0, problem, algo)
+    end
+end;
+
+@testset "Bicoloring & substitution decompression" begin
+    problem = ColoringProblem(; structure=:nonsymmetric, partition=:bidirectional)
+    algo = GreedyColoringAlgorithm(RandomOrder(rng); decompression=:substitution)
+    @testset "$((; m, n, p))" for (m, n, p) in asymmetric_params
+        A0 = sprand(rng, m, n, p)
+        test_bicoloring_decompression(A0, problem, algo)
+    end
+    @testset "$((; n, p))" for (n, p) in symmetric_params
+        A0 = sparse(Symmetric(sprand(rng, n, n, p)))
+        test_bicoloring_decompression(A0, problem, algo)
+    end
+end;
diff --git a/test/show_colors.jl b/test/show_colors.jl
index 102f204..b44278a 100644
--- a/test/show_colors.jl
+++ b/test/show_colors.jl
@@ -12,13 +12,15 @@ S = sparse([
 ]);
 algo = GreedyColoringAlgorithm(; decompression=:direct)
 
-@testset "$partition" for partition in (:column, :row)
+@testset "$partition" for partition in (:column, :row, :bidirectional)
     problem = ColoringProblem(; structure=:nonsymmetric, partition=partition)
     result = coloring(S, problem, algo)
 
-    img = show_colors(result)
-    @test size(img) == size(S)
-    @test img isa Matrix{<:Colorant}
+    if partition != :bidirectional
+        img = show_colors(result)
+        @test size(img) == size(S)
+        @test img isa Matrix{<:Colorant}
+    end
 
     h, w = size(S)
     scale = 3
@@ -31,15 +33,17 @@ algo = GreedyColoringAlgorithm(; decompression=:direct)
     @test size(img) == (h * (scale + pad) + pad, w * (scale + pad) + pad)
     @test img isa Matrix{<:Colorant}
 
-    @testset "color cycling" begin
-        colorscheme = [RGB(0, 0, 0), RGB(1, 1, 1)] # 2 colors, whereas S requires 3
-        img = @test_logs (:warn,) show_colors(result; colorscheme=colorscheme)
-        @test size(img) == size(S)
-        @test img isa Matrix{<:Colorant}
-
-        img = show_colors(result; colorscheme=colorscheme, warn=false)
-        @test size(img) == size(S)
-        @test img isa Matrix{<:Colorant}
+    if partition != :bidirectional
+        @testset "color cycling" begin
+            colorscheme = [RGB(0, 0, 0), RGB(1, 1, 1)] # 2 colors, whereas S requires 3
+            img = @test_logs (:warn,) show_colors(result; colorscheme)
+            @test size(img) == size(S)
+            @test img isa Matrix{<:Colorant}
+
+            img = show_colors(result; colorscheme, warn=false)
+            @test size(img) == size(S)
+            @test img isa Matrix{<:Colorant}
+        end
     end
 end
 
@@ -54,9 +58,4 @@ end
         @test_throws ArgumentError show_colors(result; pad=-1)
         @test_nowarn show_colors(result; pad=0)
     end
-    # @testset "Unsupported partitions" begin
-    #     problem = ColoringProblem(; structure=:nonsymmetric, partition=:bidirectional) # TODO: not implemented by SMC
-    #     result = coloring(S, problem, algo)
-    #     @test_throws ErrorException show_colors(result)
-    # end
 end
diff --git a/test/type_stability.jl b/test/type_stability.jl
index 81848ab..860a1a3 100644
--- a/test/type_stability.jl
+++ b/test/type_stability.jl
@@ -33,6 +33,8 @@ rng = StableRNG(63)
         (:nonsymmetric, :row, :direct),
         (:symmetric, :column, :direct),
         (:symmetric, :column, :substitution),
+        (:nonsymmetric, :bidirectional, :direct),
+        (:nonsymmetric, :bidirectional, :substitution),
     ]
         @test_opt target_modules = (SparseMatrixColorings,) coloring(
             A,
@@ -72,6 +74,8 @@ end;
         (:nonsymmetric, :row, :direct),
         (:symmetric, :column, :direct),
         (:symmetric, :column, :substitution),
+        (:nonsymmetric, :bidirectional, :direct),
+        (:nonsymmetric, :bidirectional, :substitution),
     ]
         @testset "A::$(typeof(A))" for A in matrix_versions(A0)
             result = coloring(
@@ -80,32 +84,43 @@ end;
                 GreedyColoringAlgorithm(; decompression);
                 decompression_eltype=eltype(A),
             )
-            B = compress(A, result)
-            @testset "Full decompression" begin
-                @test_opt compress(A, result)
-                @test_opt decompress(B, result) ≈ A0
-                @test_opt decompress!(respectful_similar(A), B, result)
-            end
-            @testset "Single-color decompression" begin
-                if decompression == :direct
-                    b = if partition == :column
-                        B[:, 1]
-                    else
-                        B[1, :]
+            if partition == :bidirectional
+                Br, Bc = compress(A, result)
+                @testset "Full decompression" begin
+                    @test_opt compress(A, result)
+                    @test_opt decompress(Br, Bc, result) ≈ A0
+                    @test_opt decompress!(respectful_similar(A), Br, Bc, result)
+                end
+            else
+                B = compress(A, result)
+                @testset "Full decompression" begin
+                    @test_opt compress(A, result)
+                    @test_opt decompress(B, result) ≈ A0
+                    @test_opt decompress!(respectful_similar(A), B, result)
+                end
+                @testset "Single-color decompression" begin
+                    if decompression == :direct
+                        b = if partition == :column
+                            B[:, 1]
+                        else
+                            B[1, :]
+                        end
+                        @test_opt decompress_single_color!(
+                            respectful_similar(A), b, 1, result
+                        )
                     end
-                    @test_opt decompress_single_color!(respectful_similar(A), b, 1, result)
                 end
-            end
-            @testset "Triangle decompression" begin
-                if structure == :symmetric
-                    @test_opt decompress!(respectful_similar(triu(A)), B, result, :U)
+                @testset "Triangle decompression" begin
+                    if structure == :symmetric
+                        @test_opt decompress!(respectful_similar(triu(A)), B, result, :U)
+                    end
                 end
-            end
-            @testset "Single-color triangle decompression" begin
-                if structure == :symmetric && decompression == :direct
-                    @test_opt decompress_single_color!(
-                        respectful_similar(triu(A)), B[:, 1], 1, result, :U
-                    )
+                @testset "Single-color triangle decompression" begin
+                    if structure == :symmetric && decompression == :direct
+                        @test_opt decompress_single_color!(
+                            respectful_similar(triu(A)), B[:, 1], 1, result, :U
+                        )
+                    end
                 end
             end
         end
diff --git a/test/utils.jl b/test/utils.jl
index 0cec537..de76624 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -26,7 +26,7 @@ function test_coloring_decompression(
 
         if structure == :nonsymmetric && issymmetric(A)
             result = coloring(
-                A, problem, algo; decompression_eltype=Float64, symmetric_pattern=true
+                A, problem, algo; decompression_eltype=Float32, symmetric_pattern=true
             )
         else
             result = coloring(A, problem, algo; decompression_eltype=Float64)
@@ -40,6 +40,12 @@ function test_coloring_decompression(
 
         B = compress(A, result)
 
+        if partition == :column
+            @test ncolors(result) == size(B, 2)
+        elseif partition == :row
+            @test ncolors(result) == size(B, 1)
+        end
+
         @testset "Reference" begin
             @test sparsity_pattern(result) === A  # identity of objects
             !isnothing(color0) && @test color == color0
@@ -142,6 +148,33 @@ function test_coloring_decompression(
     end
 end
 
+function test_bicoloring_decompression(
+    A0::AbstractMatrix,
+    problem::ColoringProblem{:nonsymmetric,:bidirectional},
+    algo::GreedyColoringAlgorithm{decompression};
+) where {decompression}
+    @testset "$(typeof(A))" for A in matrix_versions(A0)
+        yield()
+        if issymmetric(A)
+            result = coloring(
+                A, problem, algo; decompression_eltype=Float32, symmetric_pattern=true
+            )
+        else
+            result = coloring(A, problem, algo; decompression_eltype=Float64)
+        end
+        Br, Bc = compress(A, result)
+        @test size(Br, 1) == length(unique(row_colors(result)))
+        @test size(Bc, 2) == length(unique(column_colors(result)))
+        @test ncolors(result) == size(Br, 1) + size(Bc, 2)
+        @testset "Full decompression" begin
+            @test decompress(Br, Bc, result) ≈ A0
+            @test decompress(Br, Bc, result) ≈ A0  # check result wasn't modified
+            @test decompress!(respectful_similar(A), Br, Bc, result) ≈ A0
+            @test decompress!(respectful_similar(A), Br, Bc, result) ≈ A0
+        end
+    end
+end
+
 function test_structured_coloring_decompression(A::AbstractMatrix)
     column_problem = ColoringProblem(; structure=:nonsymmetric, partition=:column)
     row_problem = ColoringProblem(; structure=:nonsymmetric, partition=:row)