Construct fusion tree blocks from the shorter side, and group the cached trees in fusionblocks

docs/src/Changelog.md

@@ -42,6 +42,7 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 ### Performance
 
 - Recoupling matrices are stored such that they can leverage BLAS, also for mixed complex tensor with real recoupling coefficients.
+- Constructing a `FusionTreeBlock` only takes the fusion product of the shorter side, and `fusionblocks` groups the cached fusion trees of the space instead of constructing a block for every combination of sectors (internal)
 
 ## [0.17.2](https://github.com/QuantumKitHub/TensorKit.jl/compare/v0.17.1...v0.17.2) - 2026-09-20
 

src/fusiontrees/fusiontrees.jl

@@ -140,18 +140,35 @@ function FusionTreeBlock{I}(
     end
 
     # note: cs is unique due to ⊗ returning unique fusion results
-    cs = if N₁ == 0
-        sort!(collect(filter(isunit, ⊗(uncoupled[2]...))))
-    elseif N₂ == 0
-        sort!(collect(filter(isunit, ⊗(uncoupled[1]...))))
+    # a superset of the coupled sectors suffices, as there are no fusion trees for the others:
+    # the unit of the other side, or the fusion product of the shorter side, which avoids that of
+    # the longer side (`vec`: products of `ProductSector`s have a shape)
+    cs = if N₁ == 0 || N₂ == 0
+        # the only unit that can couple: on the left of the first and on the right of the last sector
+        us = N₁ == 0 ? uncoupled[2] : uncoupled[1]
+        leftunit(first(us)) == rightunit(last(us)) ? [leftunit(first(us))] : I[]
     else
-        sort!(collect(intersect(⊗(uncoupled[1]...), ⊗(uncoupled[2]...))))
+        sort!(vec(collect(⊗(uncoupled[N₁ <= N₂ ? 1 : 2]...))))
     end
 
-    for c in cs,
-            f₂ in fusiontrees(uncoupled[2], c, isdual[2]),
-            f₁ in fusiontrees(uncoupled[1], c, isdual[1])
-        push!(trees, (f₁, f₂))
+    for c in cs
+        # the shorter side has fusion trees for every candidate, the longer side not necessarily
+        # (with several units, `c` must also have the units at the ends of the longer side)
+        us = N₁ <= N₂ ? uncoupled[2] : uncoupled[1]
+        UnitStyle(I) isa SimpleUnit || N₁ == 0 || N₂ == 0 ||
+            (leftunit(c) == leftunit(first(us)) && rightunit(c) == rightunit(last(us))) || continue
+        if N₁ <= N₂
+            ftrees₂ = fusiontrees(uncoupled[2], c, isdual[2])
+            isempty(ftrees₂) && continue
+            ftrees₁ = fusiontrees(uncoupled[1], c, isdual[1])
+        else
+            ftrees₁ = fusiontrees(uncoupled[1], c, isdual[1])
+            isempty(ftrees₁) && continue
+            ftrees₂ = fusiontrees(uncoupled[2], c, isdual[2])
+        end
+        for f₂ in ftrees₂, f₁ in ftrees₁
+            push!(trees, (f₁, f₂))
+        end
     end
     return FusionTreeBlock(trees)
 end

src/spaces/homspace.jl

@@ -210,13 +210,20 @@ grouped by their uncoupled charges.
 function fusionblocks(W::HomSpace)
     I = sectortype(W)
     N₁, N₂ = numout(W), numin(W)
-    isdual_src = (map(isdual, codomain(W)), map(isdual, domain(W)))
-    fblocks = Vector{FusionTreeBlock{I, N₁, N₂, fusiontreetype(I, N₁, N₂)}}()
-    for dom_uncoupled_src in sectors(domain(W)), cod_uncoupled_src in sectors(codomain(W))
-        fs_src = FusionTreeBlock{I}((cod_uncoupled_src, dom_uncoupled_src), isdual_src)
-        isempty(fs_src) || push!(fblocks, fs_src)
+    F = fusiontreetype(I, N₁, N₂)
+    # group the (cached) fusion trees of `W` by their uncoupled sectors,
+    # in the order in which these first occur
+    blocks = FusionTreeBlock{I, N₁, N₂, F}[]
+    trees_of = Dict{Tuple{NTuple{N₁, I}, NTuple{N₂, I}}, Vector{F}}()
+    for (f₁, f₂) in fusiontrees(W)
+        trees = get!(trees_of, (f₁.uncoupled, f₂.uncoupled)) do
+            newtrees = F[]
+            push!(blocks, FusionTreeBlock(newtrees))
+            return newtrees
+        end
+        push!(trees, (f₁, f₂))
     end
-    return fblocks
+    return blocks
 end
 
 function diagonalblockstructure(W::HomSpace)

test/symmetries/doubletree.jl

@@ -235,6 +235,19 @@ using TensorKitSectors
         end
     end
 
+    UnitStyle(I) isa SimpleUnit && @testset "Double fusion tree: blocks of a HomSpace" begin
+        V = Vect[I](Dict(s => 1 for s in (out..., incoming)))
+        W = V ⊗ V' ← V ⊗ V' # never empty: the unit couples both sides
+        trees = collect(fusiontrees(W))
+        blocks = TK.fusionblocks(W)
+        # every block holds the trees of `W` with its uncoupled sectors, in the same order
+        for block in blocks
+            @test fusiontrees(block) == filter(f -> (f[1].uncoupled, f[2].uncoupled) == block.uncoupled, trees)
+        end
+        @test sum(length, blocks) == length(trees)
+        @test allunique(map(block -> block.uncoupled, blocks))
+    end
+
     @testset "Double fusion tree: planar trace" begin
         if FusionStyle(I) isa UniqueFusion
             f1, f1 = src