diff --git a/docs/src/dev.md b/docs/src/dev.md
index 2f253db..c33febe 100644
--- a/docs/src/dev.md
+++ b/docs/src/dev.md
@@ -26,6 +26,7 @@ SparseMatrixColorings.partial_distance2_coloring
 SparseMatrixColorings.star_coloring
 SparseMatrixColorings.acyclic_coloring
 SparseMatrixColorings.group_by_color
+SparseMatrixColorings.remap_colors
 SparseMatrixColorings.Forest
 SparseMatrixColorings.StarSet
 SparseMatrixColorings.TreeSet
@@ -39,8 +40,8 @@ SparseMatrixColorings.RowColoringResult
 SparseMatrixColorings.StarSetColoringResult
 SparseMatrixColorings.TreeSetColoringResult
 SparseMatrixColorings.LinearSystemColoringResult
-SparseMatrixColorings.BicoloringResult
-SparseMatrixColorings.remap_colors
+SparseMatrixColorings.StarSetBicoloringResult
+SparseMatrixColorings.TreeSetBicoloringResult
 ```
 
 ## Testing
diff --git a/src/decompression.jl b/src/decompression.jl
index 97aa090..68ff0de 100644
--- a/src/decompression.jl
+++ b/src/decompression.jl
@@ -515,12 +515,14 @@ end
 ## TreeSetColoringResult
 
 function decompress!(
-    A::AbstractMatrix, B::AbstractMatrix, result::TreeSetColoringResult, uplo::Symbol=:F
-)
+    A::AbstractMatrix{R},
+    B::AbstractMatrix{R},
+    result::TreeSetColoringResult,
+    uplo::Symbol=:F,
+) where {R<:Real}
     (; ag, color, reverse_bfs_orders, tree_edge_indices, nt, buffer) = result
     (; S) = ag
     uplo == :F && check_same_pattern(A, S)
-    R = eltype(A)
     fill!(A, zero(R))
 
     if eltype(buffer) == R
@@ -714,61 +716,184 @@ function decompress!(
     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 vectors `symmetric_to_row` and `symmetric_to_column` to map from symmetric colors to row and column colors.
-    =#
-    (; A, symmetric_to_column, symmetric_to_row) = 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 symmetric_to_row[c] > 0  # some rows were colored with the symmetric color c
-            copyto!(view(Br_and_Bc, 1:n, c), view(Br, symmetric_to_row[c], :))
-        end
-        if symmetric_to_column[c] > 0  # some columns were colored with the symmetric color c
-            copyto!(
-                view(Br_and_Bc, (n + 1):(n + m), c), view(Bc, :, symmetric_to_column[c])
-            )
+## StarSetBicoloringResult
+
+function decompress!(
+    A::AbstractMatrix,
+    Br::AbstractMatrix,
+    Bc::AbstractMatrix,
+    result::StarSetBicoloringResult,
+)
+    (; S, A_indices, compressed_indices, pos_Bc) = result
+    fill!(A, zero(eltype(A)))
+
+    ind_Bc = 1
+    ind_Br = nnz(S)
+    rvS = rowvals(S)
+    for j in axes(S, 2)
+        for k in nzrange(S, j)
+            i = rvS[k]
+            index_Bc = compressed_indices[ind_Bc]
+            if A_indices[ind_Bc] == k
+                A[i, j] = Bc[index_Bc]
+                if ind_Bc < pos_Bc
+                    ind_Bc += 1
+                end
+            else
+                index_Br = compressed_indices[ind_Br]
+                A[i, j] = Br[index_Br]
+                ind_Br -= 1
+            end
         end
     end
-    return Br_and_Bc
+    return A
 end
 
 function decompress!(
-    A::AbstractMatrix, Br::AbstractMatrix, Bc::AbstractMatrix, result::BicoloringResult
+    A::SparseMatrixCSC,
+    Br::AbstractMatrix,
+    Bc::AbstractMatrix,
+    result::StarSetBicoloringResult,
 )
+    (; A_indices, compressed_indices, pos_Bc) = result
+    nzA = nonzeros(A)
+    for k in 1:pos_Bc
+        nzA[A_indices[k]] = Bc[compressed_indices[k]]
+    end
+    for k in (pos_Bc + 1):length(nzA)
+        nzA[A_indices[k]] = Br[compressed_indices[k]]
+    end
+    return A
+end
+
+## TreeSetBicoloringResult
+
+function decompress!(
+    A::AbstractMatrix{R},
+    Br::AbstractMatrix{R},
+    Bc::AbstractMatrix{R},
+    result::TreeSetBicoloringResult,
+) where {R<:Real}
+    (;
+        symmetric_color,
+        symmetric_to_row,
+        symmetric_to_column,
+        reverse_bfs_orders,
+        tree_edge_indices,
+        nt,
+        buffer,
+    ) = result
+
     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
+    fill!(A, zero(R))
+
+    if eltype(buffer) == R
+        buffer_right_type = buffer
+    else
+        buffer_right_type = similar(buffer, R)
+    end
+
+    for k in 1:nt
+        # Positions of the first and last edges of the tree
+        first = tree_edge_indices[k]
+        last = tree_edge_indices[k + 1] - 1
+
+        # Reset the buffer to zero for all vertices in the tree (except the root)
+        for pos in first:last
+            (vertex, _) = reverse_bfs_orders[pos]
+            buffer_right_type[vertex] = zero(R)
+        end
+        # Reset the buffer to zero for the root vertex
+        (_, root) = reverse_bfs_orders[last]
+        buffer_right_type[root] = zero(R)
+
+        for pos in first:last
+            (i, j) = reverse_bfs_orders[pos]
+            cj = symmetric_color[j]
+            if in_triangle(i, j, :L)
+                val = Bc[i - n, symmetric_to_column[cj]] - buffer_right_type[i]
+                buffer_right_type[j] = buffer_right_type[j] + val
+                A[i - n, j] = val
+            else
+                val = Br[symmetric_to_row[cj], i] - buffer_right_type[i]
+                buffer_right_type[j] = buffer_right_type[j] + val
+                A[j - n, i] = val
+            end
+        end
+    end
     return A
 end
 
 function decompress!(
-    A::SparseMatrixCSC, Br::AbstractMatrix, Bc::AbstractMatrix, result::BicoloringResult
-)
-    (; large_colptr, large_rowval, symmetric_result) = result
+    A::SparseMatrixCSC{R},
+    Br::AbstractMatrix{R},
+    Bc::AbstractMatrix{R},
+    result::TreeSetBicoloringResult,
+) where {R<:Real}
+    (;
+        symmetric_color,
+        symmetric_to_column,
+        symmetric_to_row,
+        reverse_bfs_orders,
+        tree_edge_indices,
+        nt,
+        A_indices,
+        buffer,
+    ) = 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)
+    nzA = nonzeros(A)
+
+    if eltype(buffer) == R
+        buffer_right_type = buffer
+    else
+        buffer_right_type = similar(buffer, R)
+    end
+
+    # Recover the coefficients of A
+    counter = 0
+
+    # Recover the off-diagonal coefficients of A
+    for k in 1:nt
+        # Positions of the first and last edges of the tree
+        first = tree_edge_indices[k]
+        last = tree_edge_indices[k + 1] - 1
+
+        # Reset the buffer to zero for all vertices in the tree (except the root)
+        for pos in first:last
+            (vertex, _) = reverse_bfs_orders[pos]
+            buffer_right_type[vertex] = zero(R)
+        end
+        # Reset the buffer to zero for the root vertex
+        (_, root) = reverse_bfs_orders[last]
+        buffer_right_type[root] = zero(R)
+
+        for pos in first:last
+            (i, j) = reverse_bfs_orders[pos]
+            cj = symmetric_color[j]
+            counter += 1
+
+            #! format: off
+            # A[i,j] is in the lower triangular part of A
+            if in_triangle(i, j, :L)
+                val = Bc[i - n, symmetric_to_column[cj]] - buffer_right_type[i]
+                buffer_right_type[j] = buffer_right_type[j] + val
+
+                # A[i,j] is stored at index_ij = A_indices[counter] in A.nzval
+                nzind = A_indices[counter]
+                nzA[nzind] = val
+
+            # A[i,j] is in the upper triangular part of A
+            else
+                val = Br[symmetric_to_row[cj], i] - buffer_right_type[i]
+                buffer_right_type[j] = buffer_right_type[j] + val
+
+                # A[j,i] is stored at index_ji = A_indices[counter] in A.nzval
+                nzind = A_indices[counter]
+                nzA[nzind] = val
+            end
+            #! format: on
+        end
+    end
     return A
 end
diff --git a/src/graph.jl b/src/graph.jl
index 2333178..96b153f 100644
--- a/src/graph.jl
+++ b/src/graph.jl
@@ -100,7 +100,8 @@ end
 Return a [`SparsityPatternCSC`](@ref) corresponding to the matrix `[0 Aᵀ; A 0]`, with a minimum of allocations.
 """
 function bidirectional_pattern(A::AbstractMatrix; symmetric_pattern::Bool)
-    bidirectional_pattern(SparsityPatternCSC(SparseMatrixCSC(A)); symmetric_pattern)
+    S = SparsityPatternCSC(SparseMatrixCSC(A))
+    bidirectional_pattern(S; symmetric_pattern)
 end
 
 function bidirectional_pattern(S::SparsityPatternCSC{T}; symmetric_pattern::Bool) where {T}
@@ -172,7 +173,7 @@ function bidirectional_pattern(S::SparsityPatternCSC{T}; symmetric_pattern::Bool
 
     # Create the SparsityPatternCSC of the augmented adjacency matrix
     S_and_Sᵀ = SparsityPatternCSC{T}(p, p, colptr, rowval)
-    return S_and_Sᵀ, edge_to_index
+    return S, S_and_Sᵀ, edge_to_index
 end
 
 function build_edge_to_index(S::SparsityPatternCSC{T}) where {T}
diff --git a/src/interface.jl b/src/interface.jl
index 79fe71e..c27a623 100644
--- a/src/interface.jl
+++ b/src/interface.jl
@@ -282,7 +282,7 @@ function _coloring(
     algo::GreedyColoringAlgorithm{:substitution},
     decompression_eltype::Type{R},
     symmetric_pattern::Bool,
-) where {R}
+) where {R<:Real}
     ag = AdjacencyGraph(A; has_diagonal=true)
     vertices_in_order = vertices(ag, algo.order)
     color, tree_set = acyclic_coloring(ag, vertices_in_order, algo.postprocessing)
@@ -298,19 +298,18 @@ function _coloring(
     A::AbstractMatrix,
     ::ColoringProblem{:nonsymmetric,:bidirectional},
     algo::GreedyColoringAlgorithm{:direct},
-    decompression_eltype::Type{R},
+    decompression_eltype::Type,
     symmetric_pattern::Bool,
-) where {R}
-    A_and_Aᵀ, edge_to_index = bidirectional_pattern(A; symmetric_pattern)
-    ag = AdjacencyGraph(A_and_Aᵀ, edge_to_index; has_diagonal=false)
+)
+    S, S_and_Sᵀ, edge_to_index = bidirectional_pattern(A; symmetric_pattern)
+    ag = AdjacencyGraph(S_and_Sᵀ, edge_to_index; has_diagonal=false)
     vertices_in_order = vertices(ag, algo.order)
-    color, star_set = star_coloring(ag, vertices_in_order, algo.postprocessing)
+    symmetric_color, star_set = star_coloring(ag, vertices_in_order, algo.postprocessing)
     if speed_setting isa WithResult
-        symmetric_result = StarSetColoringResult(A_and_Aᵀ, ag, color, star_set)
-        return BicoloringResult(A, ag, symmetric_result, R)
+        return StarSetBicoloringResult(A, S, ag, symmetric_color, star_set)
     else
         row_color, column_color, _ = remap_colors(
-            eltype(ag), color, maximum(color), size(A)...
+            eltype(ag), symmetric_color, maximum(symmetric_color), size(A)...
         )
         return row_color, column_color
     end
@@ -324,16 +323,15 @@ function _coloring(
     decompression_eltype::Type{R},
     symmetric_pattern::Bool,
 ) where {R}
-    A_and_Aᵀ, edge_to_index = bidirectional_pattern(A; symmetric_pattern)
-    ag = AdjacencyGraph(A_and_Aᵀ, edge_to_index; has_diagonal=false)
+    S, S_and_Sᵀ, edge_to_index = bidirectional_pattern(A; symmetric_pattern)
+    ag = AdjacencyGraph(S_and_Sᵀ, edge_to_index; has_diagonal=false)
     vertices_in_order = vertices(ag, algo.order)
-    color, tree_set = acyclic_coloring(ag, vertices_in_order, algo.postprocessing)
+    symmetric_color, tree_set = acyclic_coloring(ag, vertices_in_order, algo.postprocessing)
     if speed_setting isa WithResult
-        symmetric_result = TreeSetColoringResult(A_and_Aᵀ, ag, color, tree_set, R)
-        return BicoloringResult(A, ag, symmetric_result, R)
+        return TreeSetBicoloringResult(A, S, ag, symmetric_color, tree_set, R)
     else
         row_color, column_color, _ = remap_colors(
-            eltype(ag), color, maximum(color), size(A)...
+            eltype(ag), symmetric_color, maximum(symmetric_color), size(A)...
         )
         return row_color, column_color
     end
diff --git a/src/result.jl b/src/result.jl
index ad2ba34..6cd5943 100644
--- a/src/result.jl
+++ b/src/result.jl
@@ -558,10 +558,20 @@ function remap_colors(
     return row_color, column_color, symmetric_to_row, symmetric_to_column
 end
 
+function column_colors(result::AbstractColoringResult{:nonsymmetric,:bidirectional})
+    return result.column_color
+end
+function column_groups(result::AbstractColoringResult{:nonsymmetric,:bidirectional})
+    return result.column_group
+end
+
+row_colors(result::AbstractColoringResult{:nonsymmetric,:bidirectional}) = result.row_color
+row_groups(result::AbstractColoringResult{:nonsymmetric,:bidirectional}) = result.row_group
+
 """
 $TYPEDEF
 
-Storage for the result of a bidirectional coloring with direct or substitution decompression, based on the symmetric coloring of a 2x2 block matrix.
+Storage for the result of a bidirectional coloring with direct decompression, based on the symmetric star coloring of a 2 x 2 block matrix.
 
 # Fields
 
@@ -571,19 +581,17 @@ $TYPEDFIELDS
 
 - [`AbstractColoringResult`](@ref)
 """
-struct BicoloringResult{
-    M<:AbstractMatrix,
-    T<:Integer,
-    G<:AdjacencyGraph{T},
-    decompression,
-    GT<:AbstractGroups{T},
-    SR<:AbstractColoringResult{:symmetric,:column,decompression},
-    R,
-} <: AbstractColoringResult{:nonsymmetric,:bidirectional,decompression}
+struct StarSetBicoloringResult{
+    M<:AbstractMatrix,T<:Integer,G<:AdjacencyGraph{T},GT<:AbstractGroups{T}
+} <: AbstractColoringResult{:nonsymmetric,:bidirectional,:direct}
     "matrix that was colored"
     A::M
+    "sparsity pattern of A"
+    S::SparsityPatternCSC{T}
     "augmented adjacency graph that was used for bicoloring"
     abg::G
+    "one integer color for each axis of the augmented matrix"
+    symmetric_color::Vector{T}
     "one integer color for each column"
     column_color::Vector{T}
     "one integer color for each row"
@@ -592,56 +600,200 @@ struct BicoloringResult{
     column_group::GT
     "color groups for rows"
     row_group::GT
-    "result for the coloring of the symmetric 2 x 2 block matrix"
-    symmetric_result::SR
     "maps symmetric colors to column colors"
     symmetric_to_column::Vector{T}
     "maps symmetric colors to row colors"
     symmetric_to_row::Vector{T}
-    "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{T}
-    "CSC storage of `A_and_noAᵀ - `rowval`"
-    large_rowval::Vector{T}
+    "indices of the nonzeros in A that can recovered with Br and Bc"
+    A_indices::Vector{T}
+    compressed_indices::Vector{T}
+    pos_Bc::T
 end
 
-column_colors(result::BicoloringResult) = result.column_color
-column_groups(result::BicoloringResult) = result.column_group
+function StarSetBicoloringResult(
+    A::AbstractMatrix,
+    S::SparsityPatternCSC{T},
+    ag::AdjacencyGraph{T},
+    symmetric_color::Vector{<:Integer},
+    star_set::StarSet{<:Integer},
+) where {T<:Integer}
+    m, n = size(A)
+    (; star, hub) = star_set
+    num_sym_colors = maximum(symmetric_color)
+    row_color, column_color, symmetric_to_row, symmetric_to_column = remap_colors(
+        T, symmetric_color, num_sym_colors, m, n
+    )
+    column_group = group_by_color(T, column_color)
+    row_group = group_by_color(T, row_color)
+    num_row_colors = length(row_group)
 
-row_colors(result::BicoloringResult) = result.row_color
-row_groups(result::BicoloringResult) = result.row_group
+    rv = rowvals(S)
+    nnzA = nnz(S)
+    A_indices = Vector{T}(undef, nnzA)
+    compressed_indices = Vector{T}(undef, nnzA)
+    pos_Bc = 0
+    pos_Br = nnzA + 1
+    for j in 1:n
+        for k in nzrange(S, j)
+            i = rv[k]
+            s = star[k]
+            h = abs(hub[s])
+            # Assign the non-hub vertex (spoke) to the correct position in spokes
+            if j == h
+                # j is the hub and i is the spoke
+                c = symmetric_color[j]
+
+                # A[i, j] = Bc[i, symmetric_to_column[c]]
+                pos_Bc += 1
+                A_indices[pos_Bc] = k
+                compressed_indices[pos_Bc] = (symmetric_to_column[c] - 1) * m + i
+            else  # i + n == h
+                # (i + n) is the hub and j is the spoke
+                c = symmetric_color[i + n]
+
+                # A[i, j] = Br[symmetric_to_row[c], j]
+                pos_Br -= 1
+                A_indices[pos_Br] = k
+                compressed_indices[pos_Br] = (j - 1) * num_row_colors + symmetric_to_row[c]
+            end
+        end
+    end
 
-function BicoloringResult(
+    return StarSetBicoloringResult(
+        A,
+        S,
+        ag,
+        symmetric_color,
+        column_color,
+        row_color,
+        column_group,
+        row_group,
+        symmetric_to_column,
+        symmetric_to_row,
+        A_indices,
+        compressed_indices,
+        T(pos_Bc),
+    )
+end
+
+"""
+$TYPEDEF
+
+Storage for the result of a bidirectional coloring with decompression by substitution, based on the symmetric acyclic coloring of a 2 x 2 block matrix.
+
+# Fields
+
+$TYPEDFIELDS
+
+# See also
+
+- [`AbstractColoringResult`](@ref)
+"""
+struct TreeSetBicoloringResult{
+    M<:AbstractMatrix,T<:Integer,G<:AdjacencyGraph{T},GT<:AbstractGroups{T},R
+} <: AbstractColoringResult{:nonsymmetric,:bidirectional,:substitution}
+    "matrix that was colored"
+    A::M
+    "sparsity pattern of A"
+    S::SparsityPatternCSC{T}
+    "augmented adjacency graph that was used for bicoloring"
+    abg::G
+    "one integer color for each axis of the augmented matrix"
+    symmetric_color::Vector{T}
+    "one integer color for each column"
+    column_color::Vector{T}
+    "one integer color for each row"
+    row_color::Vector{T}
+    "color groups for columns"
+    column_group::GT
+    "color groups for rows"
+    row_group::GT
+    "maps symmetric colors to column colors"
+    symmetric_to_column::Vector{T}
+    "maps symmetric colors to row colors"
+    symmetric_to_row::Vector{T}
+    "indices of the nonzeros in A that can recovered with Br and Bc"
+    A_indices::Vector{T}
+    "reverse BFS orders of the trees"
+    reverse_bfs_orders::Vector{Tuple{T,T}}
+    tree_edge_indices::Vector{T}
+    nt::T
+    "buffer needed during acyclic decompression"
+    buffer::Vector{R}
+end
+
+function TreeSetBicoloringResult(
     A::AbstractMatrix,
+    S::SparsityPatternCSC{T},
     ag::AdjacencyGraph{T},
-    symmetric_result::AbstractColoringResult{:symmetric,:column},
+    symmetric_color::Vector{<:Integer},
+    tree_set::TreeSet{<:Integer},
     decompression_eltype::Type{R},
 ) where {T,R}
+    (; reverse_bfs_orders, tree_edge_indices, nt) = tree_set
+
     m, n = size(A)
-    symmetric_color = column_colors(symmetric_result)
     num_sym_colors = maximum(symmetric_color)
     row_color, column_color, symmetric_to_row, symmetric_to_column = remap_colors(
         T, symmetric_color, num_sym_colors, m, n
     )
     column_group = group_by_color(T, column_color)
     row_group = group_by_color(T, row_color)
-    Br_and_Bc = Matrix{R}(undef, n + m, num_sym_colors)
-    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(
+
+    rv = rowvals(S)
+    nnzA = nnz(S)
+    A_indices = Vector{T}(undef, nnzA)
+
+    index = 0
+    for k in 1:nt
+        # Positions of the edges for each tree
+        first = tree_edge_indices[k]
+        last = tree_edge_indices[k + 1] - 1
+
+        for pos in first:last
+            (leaf, neighbor) = reverse_bfs_orders[pos]
+            i = leaf
+            j = neighbor
+            index += 1
+
+            #! format: off
+            # S[i,j] is in the lower triangular part of S
+            if in_triangle(i, j, :L)
+                # S[i,j] is stored at index_ij = (S.colptr[j] + offset) in S.nzval
+                col_j = view(rv, nzrange(S, j))
+                offset = searchsortedfirst(col_j, i-n)::Int - 1
+                A_indices[index] = S.colptr[j] + offset
+
+            # S[i,j] is in the upper triangular part of S
+            else
+                # S[j,i] is stored at index_ji = (S.colptr[i] + offset) in S.nzval
+                col_i = view(rv, nzrange(S, i))
+                offset = searchsortedfirst(col_i, j-n)::Int - 1
+                A_indices[index] = S.colptr[i] + offset
+            end
+            #! format: on
+        end
+    end
+
+    # buffer holds the sum of edge values for subtrees in a tree.
+    # For each vertex i, buffer[i] is the sum of edge values in the subtree rooted at i.
+    buffer = Vector{R}(undef, n + m)
+
+    return TreeSetBicoloringResult(
         A,
+        S,
         ag,
+        symmetric_color,
         column_color,
         row_color,
         column_group,
         row_group,
-        symmetric_result,
         symmetric_to_column,
         symmetric_to_row,
-        Br_and_Bc,
-        large_colptr,
-        large_rowval,
+        A_indices,
+        reverse_bfs_orders,
+        tree_edge_indices,
+        nt,
+        buffer,
     )
 end
diff --git a/test/allocations.jl b/test/allocations.jl
index caed10f..a2ca80a 100644
--- a/test/allocations.jl
+++ b/test/allocations.jl
@@ -35,7 +35,7 @@ function test_noallocs_sparse_decompression(
             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
+            @test minimum(bench2_full).allocs == 0
         end
     else
         B = compress(A, result)
diff --git a/test/graph.jl b/test/graph.jl
index 9c8f66c..e65673b 100644
--- a/test/graph.jl
+++ b/test/graph.jl
@@ -35,7 +35,9 @@ using Test
                 p = 0.05 * rand()
                 A = sprand(Bool, m, n, p)
                 A_and_Aᵀ = [spzeros(Bool, n, n) transpose(A); A spzeros(Bool, m, m)]
-                S_and_Sᵀ, edge_to_index = bidirectional_pattern(A; symmetric_pattern=false)
+                S, S_and_Sᵀ, edge_to_index = bidirectional_pattern(
+                    A; symmetric_pattern=false
+                )
                 @test S_and_Sᵀ.colptr == A_and_Aᵀ.colptr
                 @test S_and_Sᵀ.rowval == A_and_Aᵀ.rowval
                 M = SparseMatrixCSC(
@@ -50,7 +52,9 @@ using Test
                 p = 0.05 * rand()
                 A = sparse(Symmetric(sprand(Bool, m, m, p)))
                 A_and_Aᵀ = [spzeros(Bool, m, m) transpose(A); A spzeros(Bool, m, m)]
-                S_and_Sᵀ, edge_to_index = bidirectional_pattern(A; symmetric_pattern=true)
+                S, S_and_Sᵀ, edge_to_index = bidirectional_pattern(
+                    A; symmetric_pattern=true
+                )
                 @test S_and_Sᵀ.colptr == A_and_Aᵀ.colptr
                 @test S_and_Sᵀ.rowval == A_and_Aᵀ.rowval
                 M = SparseMatrixCSC(