From 3ab0e0edd821c6f003b88583ea2fe39577612046 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Thu, 8 Oct 2026 19:34:57 -0400
Subject: [PATCH 01/12] SumSpace: component-wise sectorhash/sectorequal, fix
 axes(S, c)

The ElementarySpace fallbacks built and sorted a sector Set on every
structure-cache lookup keyed by a SumSpace HomSpace. The component-wise
key is finer than necessary but sound. axes(S, c) used an unimported
flatten on an untyped vector.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/vectorspaces/sumspace.jl  | 21 ++++++++++++++++++---
 test/vectorspaces/sumspace.jl | 14 ++++++++++++++
 2 files changed, 32 insertions(+), 3 deletions(-)

diff --git a/src/vectorspaces/sumspace.jl b/src/vectorspaces/sumspace.jl
index 1375e69..6cbff49 100644
--- a/src/vectorspaces/sumspace.jl
+++ b/src/vectorspaces/sumspace.jl
@@ -75,12 +75,12 @@ Base.axes(S::SumSpace) = Base.OneTo(dim(S))
 Base.axes(S::SumSpace, n::Int) = axes(S.spaces, n)
 function Base.axes(S::SumSpace, c::Sector)
     offset = 0
-    a = []
+    a = Int[]
     for s in S.spaces
-        a = push!(a, axes(s, c) .+ offset)
+        append!(a, axes(s, c) .+ offset)
         offset += dim(s)
     end
-    return collect(flatten(a))
+    return a
 end
 
 Base.hash(S::SumSpace, h::UInt) = hash(S.spaces, h)
@@ -134,6 +134,21 @@ function TensorKit._sectors(S::SumSpace, ::Type{I}) where {I}
     return sort!(collect(s))
 end
 
+function TensorKit.sectorequal(S₁::SumSpace, S₂::SumSpace)
+    isdual(S₁) == isdual(S₂) && length(S₁) == length(S₂) || return false
+    for (V₁, V₂) in zip(S₁.spaces, S₂.spaces)
+        TK.sectorequal(V₁, V₂) || return false
+    end
+    return true
+end
+function TensorKit.sectorhash(S::SumSpace, h::UInt)
+    h = hash(isdual(S), h)
+    for V in S.spaces
+        h = TK.sectorhash(V, h)
+    end
+    return h
+end
+
 TensorKit.dim(S::SumSpace, sector::Sector) = sum(v -> dim(v, sector), S.spaces; init = 0)
 # ambiguity fix:
 TensorKit.dim(S::SumSpace, ::Trivial) = sum(v -> dim(v, Trivial()), S.spaces; init = 0)
diff --git a/test/vectorspaces/sumspace.jl b/test/vectorspaces/sumspace.jl
index 0e9eb0c..b227b89 100644
--- a/test/vectorspaces/sumspace.jl
+++ b/test/vectorspaces/sumspace.jl
@@ -266,3 +266,17 @@ end
         @test issetequal(@constinferred(blocksectors(v ← v)), blocksectors(v))
     end
 end
+
+@testset "sectorhash and sectorequal" begin
+    V = Vect[U1Irrep](0 => 2, 1 => 2, -1 => 2) ⊞ Vect[U1Irrep](0 => 3, 1 => 1) ⊞ Vect[U1Irrep](2 => 1)
+    V′ = Vect[U1Irrep](0 => 1, 1 => 1, -1 => 3) ⊞ Vect[U1Irrep](0 => 1, 1 => 2) ⊞ Vect[U1Irrep](2 => 4)
+    @test TensorKit.sectorequal(V, V′)
+    @test TensorKit.sectorhash(V, UInt(1)) == TensorKit.sectorhash(V′, UInt(1))
+    @test !TensorKit.sectorequal(V, V')
+    @test !TensorKit.sectorequal(V, V ⊞ V)
+    @test axes(V, U1Irrep(1)) == [3, 4, 10]
+    for W in (V ⊗ V ← V, V ⊗ V' ← V ⊗ V, (V ⊞ V) ⊗ V' ← V)
+        @test blocksectors(W) == blocksectors(convert(TensorMapSpace, W))
+    end
+    @test blocksectors(V′ ⊗ V′' ← V′) == blocksectors(V ⊗ V' ← V)
+end

From 920711642ed6d72e8cca66cd190abd0851b0c0e0 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Thu, 8 Oct 2026 19:34:57 -0400
Subject: [PATCH 02/12] Cache sub-tensor block structure in blocks/foreachblock

blocks(t) now builds each sub-tensor's BlockIterator once and assembles
the mortar per sector; foreachblock and == for block tensors go through
it. Missing sparse entries get a zero block of the block size instead of
a full zero tensor, also in subblock. setindex! with a block array uses
a single view per block, and no longer calls getindex! on dense parents.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 .../abstractblocktensor/abstracttensormap.jl  | 118 ++++++++++++------
 test/abstracttensor/sparseblocktensor.jl      |  22 ++++
 2 files changed, 101 insertions(+), 39 deletions(-)

diff --git a/src/tensors/abstractblocktensor/abstracttensormap.jl b/src/tensors/abstractblocktensor/abstracttensormap.jl
index 3a68f2e..1ed8255 100644
--- a/src/tensors/abstractblocktensor/abstracttensormap.jl
+++ b/src/tensors/abstractblocktensor/abstracttensormap.jl
@@ -18,13 +18,15 @@ eachspace(t::AbstractBlockTensorMap) = SumSpaceIndices(space(t))
     numout(t) == length(f₁) && numin(t) == length(f₂) ||
         throw(DimensionMismatch("Invalid number of fusiontree legs for this tensor"))
 
-    subblocks = map(eachspace(t), parent(t)) do V, x
+    W = eachspace(t)
+    subblocks = map(eachindex(IndexCartesian(), t)) do I
+        V = W[I]
         sz = (dims(codomain(V), f₁.uncoupled)..., dims(domain(V), f₂.uncoupled)...)
-        if prod(sz) == 0
-            data = storagetype(t)(undef, 0)
+        if prod(sz) == 0 || (issparse(t) && !haskey(parent(t), I))
+            data = zerovector!(similar(storagetype(t), prod(sz)))
             return sreshape(StridedView(data), sz)
         else
-            return x[f₁, f₂]
+            return t[I][f₁, f₂]
         end
     end
 
@@ -34,10 +36,10 @@ end
         t::AbstractBlockTensorMap, v::AbstractBlockArray, f₁::FusionTree, f₂::FusionTree
     )
     for I in eachindex(t)
-        b = v[Block(I.I)]
-        if !isempty(b)
-            getindex!(t, I)[f₁, f₂] = v[Block(I.I)]
-        end
+        b = view(v, Block(I.I))
+        isempty(b) && continue
+        x = issparse(t) ? getindex!(t, I) : t[I]
+        x[f₁, f₂] = b
     end
     return t
 end
@@ -53,57 +55,95 @@ end
     return setindex!(t, v′, f₁, f₂)
 end
 
-function TensorKit.block(t::AbstractBlockTensorMap, c::Sector)
-    sectortype(t) == typeof(c) || throw(SectorMismatch())
-
-    rows = prod(TT.getindices(size(t), codomainind(t)))
-    cols = prod(TT.getindices(size(t), domainind(t)))
+TensorKit.block(t::AbstractBlockTensorMap, c::Sector) = _block(t, parent(t), c)
 
-    if rows == 0 || cols == 0
-        allblocks = Matrix{TK.blocktype(eltype(t))}(undef, rows, cols)
-
-        rowaxes = Int[]
-        if rows != 0
-            W′ = codomain(t) ← zero(spacetype(t))
-            for V in eachspace(W′)
-                push!(rowaxes, blockdim(codomain(V), c))
-            end
+function _block(t::AbstractBlockTensorMap, entries, c::Sector)
+    sectortype(t) == typeof(c) || throw(SectorMismatch())
+    rowdims = _blockdims(codomain(t), c)
+    coldims = _blockdims(domain(t), c)
+    allblocks = Matrix{_blockeltype(typeof(t))}(undef, length(rowdims), length(coldims))
+    for (k, I) in enumerate(eachindex(IndexCartesian(), t))
+        allblocks[k] = if !issparse(t) || haskey(entries, I)
+            _cachedblock(entries[I], c)
+        else
+            i, j = Tuple(CartesianIndices(allblocks)[k])
+            _zeroblock(eltype(t), rowdims[i], coldims[j], () -> eachspace(t)[I], c)
         end
+    end
+    return mortar(allblocks, rowdims, coldims)
+end
 
-        colaxes = Int[]
-        if cols != 0
-            W′ = zero(spacetype(t)) ← domain(t)
-            for V in eachspace(W′)
-                push!(colaxes, blockdim(domain(V), c))
-            end
+function _blockdims(P::ProductSumSpace{S, N}, c::Sector) where {S, N}
+    return vec(
+        map(CartesianIndices(map(length, P.spaces))) do I
+            return blockdim(ProductSpace{S, N}(map(getindex, P.spaces, Tuple(I))), c)
         end
+    )
+end
 
-        return mortar(allblocks, rowaxes, colaxes)
-    end
+_blockcache(x) = x
+_blockcache(x::AbstractTensorMap) = (b = TK.blocks(x); b isa TK.BlockIterator ? b : x)
+_cachedblock(x, c::Sector) = block(x, c)
+_cachedblock(b::TK.BlockIterator, c::Sector) = b[c]
 
-    allblocks = map(Base.Fix2(block, c), parent(t))
-    return mortar(reshape(allblocks, rows, cols))
+function _zeroblock(::Type{TT}, d₁::Int, d₂::Int, getspace, c::Sector) where {TT <: AbstractTensorMap}
+    TT <: TensorMap || return block(zerovector!(similar(TT, getspace())), c)
+    data = zerovector!(similar(storagetype(TT), d₁ * d₂))
+    return reshape(view(data, 1:(d₁ * d₂)), (d₁, d₂))
 end
 
 # TODO: this might get fixed once new tensormap is implemented
 TensorKit.blocksectors(t::AbstractBlockTensorMap) = blocksectors(space(t))
 TensorKit.hasblock(t::AbstractBlockTensorMap, c::Sector) = c in blocksectors(t)
 
-TensorKit.blocks(t::AbstractBlockTensorMap) = TK.BlockIterator(t, blocksectors(t))
-Base.@assume_effects :foldable function TensorKit.blocktype(::Type{TT}) where {TT <: AbstractBlockTensorMap}
+function TensorKit.blocks(t::AbstractBlockTensorMap)
+    entries = issparse(t) ? Dict(I => _blockcache(x) for (I, x) in nonzero_pairs(t)) :
+        map(_blockcache, parent(t))
+    return TK.BlockIterator(t, (blocksectors(t), entries))
+end
+
+Base.@assume_effects :foldable function _blockeltype(::Type{TT}) where {TT <: AbstractBlockTensorMap}
     T = scalartype(TT)
     B = TK.blocktype(eltype(TT))
-    (B <: AbstractMatrix{T}) || (B = AbstractMatrix{T}) # safeguard against type-instability
+    return B <: AbstractMatrix{T} ? B : AbstractMatrix{T} # safeguard against type-instability
+end
+function TensorKit.blocktype(::Type{TT}) where {TT <: AbstractBlockTensorMap}
     BS = NTuple{2, BlockedOneTo{Int, Vector{Int}}}
-    return BlockMatrix{T, Matrix{B}, BS}
+    return BlockMatrix{scalartype(TT), Matrix{_blockeltype(TT)}, BS}
 end
 
+Base.length(iter::TK.BlockIterator{<:AbstractBlockTensorMap}) = length(first(iter.structure))
+Base.isdone(iter::TK.BlockIterator{<:AbstractBlockTensorMap}, state...) =
+    Base.isdone(first(iter.structure), state...)
 function Base.iterate(iter::TK.BlockIterator{<:AbstractBlockTensorMap}, state...)
-    next = iterate(iter.structure, state...)
+    sectors, entries = iter.structure
+    next = iterate(sectors, state...)
     isnothing(next) && return next
     c, newstate = next
-    return c => block(iter.t, c), newstate
+    return c => _block(iter.t, entries, c), newstate
+end
+Base.getindex(iter::TK.BlockIterator{<:AbstractBlockTensorMap}, c::Sector) =
+    _block(iter.t, last(iter.structure), c)
+
+function Base.:(==)(t₁::AbstractBlockTensorMap, t₂::AbstractBlockTensorMap)
+    (codomain(t₁) == codomain(t₂) && domain(t₁) == domain(t₂)) || return false
+    b₁, b₂ = TK.blocks(t₁), TK.blocks(t₂)
+    return all(c -> b₁[c] == b₂[c], blocksectors(t₁))
+end
+
+function TensorKit.foreachblock(f, t::AbstractBlockTensorMap; scheduler = nothing)
+    foreach(TK.blocks(t)) do (c, b)
+        return f(c, (b,))
+    end
+    return nothing
+end
+function TensorKit.foreachblock(f, t::AbstractBlockTensorMap, ts...; scheduler = nothing)
+    tensors = (t, ts...)
+    iters = map(_blockcache, tensors)
+    foreach(union(blocksectors.(tensors)...)) do c
+        return f(c, map(Base.Fix2(_cachedblock, c), iters))
+    end
+    return nothing
 end
-Base.getindex(iter::TK.BlockIterator{<:AbstractBlockTensorMap}, c::Sector) = block(iter.t, c)
 
 TK.storagetype(::Type{TT}) where {TT <: AbstractBlockTensorMap} = storagetype(eltype(TT))
diff --git a/test/abstracttensor/sparseblocktensor.jl b/test/abstracttensor/sparseblocktensor.jl
index ae0d56c..4c32cb1 100644
--- a/test/abstracttensor/sparseblocktensor.jl
+++ b/test/abstracttensor/sparseblocktensor.jl
@@ -247,3 +247,25 @@ end
 
     @test HrA12array ≈ convert(TensorMap, HrA12)
 end
+
+@testset "sparse block access" begin
+    Vu = Vect[U1Irrep](0 => 2, 1 => 2, -1 => 2) ⊞ Vect[U1Irrep](0 => 3, 1 => 1) ⊞ Vect[U1Irrep](0 => 1, 2 => 1)
+    for W in (Vu ⊗ Vu ← Vu, Vu ⊗ Vu' ← Vu ⊗ Vu)
+        for t in (sprand(W, 0.3), spzeros(W))
+            td = BlockTensorMap(t)
+            @test length(TensorKit.blocks(t)) == length(blocksectors(t))
+            for (c, b) in TensorKit.blocks(t)
+                @test b isa TensorKit.blocktype(t)
+                @test @constinferred(block(t, c)) == b == block(td, c)
+                @test size.(b.blocks) == size.(block(td, c).blocks)
+            end
+            TensorKit.foreachblock(t, td) do c, (b, bd)
+                @test b == bd
+            end
+            for (f₁, f₂) in fusiontrees(t)
+                @test t[f₁, f₂] == td[f₁, f₂]
+            end
+            @test t == td
+        end
+    end
+end

From 4da243056ca9df0c296a9d6983836d8677a7ac00 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Thu, 8 Oct 2026 19:30:24 -0400
Subject: [PATCH 03/12] Speed up block write-back, norm, ==, and tr for
 AbstractBlockTensorMap

- add copy_blocks! (blockwise dense -> BlockMatrix copy, avoids scalar indexing) and use it for all dense write-backs
- truncate_domain!/truncate_codomain!: do not copy the destination that is fully overwritten
- norm: one AbstractBlockTensorMap method, no temporary; also pass p to the inner norms
- add == and tr for AbstractBlockTensorMap over nonzero keys
- fix undefined S in real/imag error message

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/auxiliary/blockarrays.jl             | 11 ++++++++
 src/linalg/factorizations.jl             | 12 ++++-----
 src/linalg/linalg.jl                     | 34 ++++++++++++++++++++----
 test/abstracttensor/blocktensor.jl       | 28 +++++++++++++++++++
 test/abstracttensor/sparseblocktensor.jl | 28 +++++++++++++++++++
 test/linalg/factorizations.jl            | 19 +++++++++++++
 6 files changed, 121 insertions(+), 11 deletions(-)

diff --git a/src/auxiliary/blockarrays.jl b/src/auxiliary/blockarrays.jl
index 97ffd24..a155eda 100644
--- a/src/auxiliary/blockarrays.jl
+++ b/src/auxiliary/blockarrays.jl
@@ -16,6 +16,17 @@ function copy_dense!(Adense, A)
     return Adense
 end
 
+function copy_blocks!(A::BlockMatrix, Adense)
+    for block_index in Iterators.product(blockaxes(A)...)
+        a = view(A, block_index...)
+        isempty(a) && continue
+        indices = getindex.(axes(A), block_index)
+        a .= view(Adense, indices...)
+    end
+    return A
+end
+copy_blocks!(A, Adense) = copyto!(A, Adense)
+
 const BlockBlasMat{T <: MAK.BlasFloat} = BlockMatrix{T}
 
 function MAK.zero!(A::BlockBlasMat)
diff --git a/src/linalg/factorizations.jl b/src/linalg/factorizations.jl
index ce887f5..317e4d4 100644
--- a/src/linalg/factorizations.jl
+++ b/src/linalg/factorizations.jl
@@ -28,7 +28,7 @@ for f! in (
             Fblocks′ = MAK.$f!(copy_dense!(similar_dense(tblock), tblock), alg)
             # deal with the case where the output is not in-place
             for (b′, b) in zip(Fblocks′, Fblocks)
-                b === b′ || copyto!(b, b′)
+                b === b′ || copy_blocks!(b, b′)
             end
             return nothing
         end
@@ -46,7 +46,7 @@ for f! in (
         TensorKit.foreachblock(t, N) do _, (tblock, Nblock)
             Nblock′ = MAK.$f!(copy_dense!(similar_dense(tblock), tblock), alg)
             # deal with the case where the output is not the same as the input
-            Nblock === Nblock′ || copyto!(Nblock, Nblock′)
+            Nblock === Nblock′ || copy_blocks!(Nblock, Nblock′)
             return nothing
         end
         return N
@@ -186,12 +186,12 @@ end
 function TensorKit.Factorizations.truncate_domain!(tdst::AbstractBlockTensorMap, tsrc::AbstractBlockTensorMap, inds)
     TensorKit.foreachblock(tdst, tsrc) do c, (dst_block, src_block)
         I = get(inds, c, nothing)
-        dst_dense = copy_dense!(similar_dense(dst_block), dst_block)
+        dst_dense = similar_dense(dst_block)
         src_dense = copy_dense!(similar_dense(src_block), src_block)
         @assert !isnothing(I)
         @views dst_dense .= src_dense[:, I]
         # deal with the case where the output is not in-place
-        dst_dense === dst_block || copyto!(dst_block, dst_dense)
+        dst_dense === dst_block || copy_blocks!(dst_block, dst_dense)
         return nothing
     end
     return tdst
@@ -199,12 +199,12 @@ end
 function TensorKit.Factorizations.truncate_codomain!(tdst::AbstractBlockTensorMap, tsrc::AbstractBlockTensorMap, inds)
     TensorKit.foreachblock(tdst, tsrc) do c, (dst_block, src_block)
         I = get(inds, c, nothing)
-        dst_dense = copy_dense!(similar_dense(dst_block), dst_block)
+        dst_dense = similar_dense(dst_block)
         src_dense = copy_dense!(similar_dense(src_block), src_block)
         @assert !isnothing(I)
         @views dst_dense .= src_dense[I, :]
         # deal with the case where the output is not in-place
-        dst_dense === dst_block || copyto!(dst_block, dst_dense)
+        dst_dense === dst_block || copy_blocks!(dst_block, dst_dense)
         return nothing
     end
     return tdst
diff --git a/src/linalg/linalg.jl b/src/linalg/linalg.jl
index 0abd9a3..cca8bd3 100644
--- a/src/linalg/linalg.jl
+++ b/src/linalg/linalg.jl
@@ -157,10 +157,34 @@ for (TC, TA, TB) in Iterators.product(
     end
 end
 
-function LinearAlgebra.norm(tA::BlockTensorMap, p::Real = 2)
+function LinearAlgebra.norm(tA::AbstractBlockTensorMap, p::Real = 2)
     vals = nonzero_values(tA)
     isempty(vals) && return norm(zero(scalartype(tA)), p)
-    return LinearAlgebra.norm(norm.(vals), p)
+    return LinearAlgebra.norm(Iterators.map(Base.Fix2(norm, p), vals), p)
+end
+
+function Base.:(==)(t1::AbstractBlockTensorMap, t2::AbstractBlockTensorMap)
+    space(t1) == space(t2) || return false
+    for I in nonzero_keys(t1)
+        v1 = t1[I]
+        (haskey(t2, I) ? v1 == t2[I] : iszero(norm(v1))) || return false
+    end
+    if issparse(t1)
+        for I in nonzero_keys(t2)
+            haskey(t1, I) || iszero(norm(t2[I])) || return false
+        end
+    end
+    return true
+end
+
+function LinearAlgebra.tr(t::AbstractBlockTensorMap)
+    domain(t) == codomain(t) ||
+        throw(SpaceMismatch("Trace of a tensor only exist when domain == codomain"))
+    s = zero(scalartype(t)) * zero(TK.dimscalartype(sectortype(t)))
+    for I in nonzero_keys(t)
+        TT.getindices(I.I, codomainind(t)) == TT.getindices(I.I, domainind(t)) && (s += tr(t[I]))
+    end
+    return s
 end
 
 for f in (:real, :imag)
@@ -173,7 +197,7 @@ for f in (:real, :imag)
             end
             return t′
         else
-            msg = "`$f` has not been implemented for `BlockTensorMap{$(S)}`."
+            msg = "`$f` has not been implemented for `BlockTensorMap{$(spacetype(t))}`."
             throw(ArgumentError(msg))
         end
     end
@@ -208,7 +232,7 @@ function LinearAlgebra.lmul!(D::DiagonalTensorMap, t::AbstractBlockTensorMap)
     domain(D) == codomain(t) || throw(SpaceMismatch())
     TensorKit.foreachblock(t, D) do c, (tblock, Dblock)
         tblock′ = lmul!(Dblock, copy_dense!(similar_dense(tblock), tblock))
-        tblock === tblock′ || copyto!(tblock, tblock′)
+        tblock === tblock′ || copy_blocks!(tblock, tblock′)
         return tblock
     end
     return t
@@ -218,7 +242,7 @@ function LinearAlgebra.rmul!(t::AbstractBlockTensorMap, D::DiagonalTensorMap)
     codomain(D) == domain(t) || throw(SpaceMismatch())
     TensorKit.foreachblock(t, D) do c, (tblock, Dblock)
         tblock′ = rmul!(copy_dense!(similar_dense(tblock), tblock), Dblock)
-        tblock === tblock′ || copyto!(tblock, tblock′)
+        tblock === tblock′ || copy_blocks!(tblock, tblock′)
         return tblock
     end
     return t
diff --git a/test/abstracttensor/blocktensor.jl b/test/abstracttensor/blocktensor.jl
index 2da67da..5658e3d 100644
--- a/test/abstracttensor/blocktensor.jl
+++ b/test/abstracttensor/blocktensor.jl
@@ -170,6 +170,34 @@ end
     end
 end
 
+@testset "norm, ==, tr: test via conversion" begin
+    W = V1 ⊗ V2 ← V1 ⊗ V2
+    for T in (Float64, ComplexF64)
+        t = rand(T, W)
+        t′ = convert(TensorMap, t)
+        for p in (1, 2, 3, Inf)
+            @test norm(t, p) ≈ norm(t′, p)
+        end
+        @test tr(t) ≈ tr(t′)
+        @test t == copy(t)
+        ts, td = SparseBlockTensorMap(t), BlockTensorMap(t)
+        @test t == ts == td
+        @test td == ts
+        @test td != 2 * ts
+        @test ts != 2 * td
+        @test norm(ts, 1) ≈ norm(td, 1)
+        @test t != 2 * t
+        t2 = copy(t)
+        t2[1] = zero(t2[1]) + t2[1]
+        @test t2 == t
+        @test (t2 + t) != t
+        @test t != rand(T, V1 ⊗ V2 ← V1)
+        @test iszero(norm(spzeros(T, W)))
+        @test spzeros(T, W) == zero(t)
+        @test iszero(tr(spzeros(T, W)))
+    end
+end
+
 @testset "Basic linear algebra: test via conversion" begin
     W = V1 ⊗ V2 ⊗ V3 ← V4 ⊗ V5
     for T in (Float32, ComplexF64)
diff --git a/test/abstracttensor/sparseblocktensor.jl b/test/abstracttensor/sparseblocktensor.jl
index 4c32cb1..4642d06 100644
--- a/test/abstracttensor/sparseblocktensor.jl
+++ b/test/abstracttensor/sparseblocktensor.jl
@@ -145,6 +145,34 @@ end
     end
 end
 
+@testset "norm, ==, tr: test via conversion" begin
+    W = V1 ⊗ V2 ← V1 ⊗ V2
+    for T in (Float64, ComplexF64)
+        t = sprand(T, W, 0.5)
+        t′ = convert(TensorMap, t)
+        for p in (1, 2, 3, Inf)
+            @test norm(t, p) ≈ norm(t′, p)
+        end
+        @test tr(t) ≈ tr(t′)
+        @test t == copy(t)
+        ts, td = SparseBlockTensorMap(t), BlockTensorMap(t)
+        @test t == ts == td
+        @test td == ts
+        @test td != 2 * ts
+        @test ts != 2 * td
+        @test norm(ts, 1) ≈ norm(td, 1)
+        @test t != 2 * t
+        t2 = copy(t)
+        t2[1] = zero(t2[1]) + t2[1]
+        @test t2 == t
+        @test (t2 + t) != t
+        @test t != rand(T, V1 ⊗ V2 ← V1)
+        @test iszero(norm(spzeros(T, W)))
+        @test spzeros(T, W) == zero(t)
+        @test iszero(tr(spzeros(T, W)))
+    end
+end
+
 @testset "Basic linear algebra: test via conversion" begin
     W = V1 ⊗ V2 ⊗ V3 ← V4 ⊗ V5
     for T in (Float32, ComplexF64)
diff --git a/test/linalg/factorizations.jl b/test/linalg/factorizations.jl
index dce9da3..9a1b457 100644
--- a/test/linalg/factorizations.jl
+++ b/test/linalg/factorizations.jl
@@ -526,3 +526,22 @@ end
 #         end
 #     end
 # end
+
+@testset "In-place write-back into block storage" begin
+    V1 = SumSpace(ℂ^2, ℂ^2)
+    V2 = SumSpace(ℂ^3, ℂ^1)
+    D = DiagonalTensorMap(rand(dim(V1)), V1)
+    D′ = DiagonalTensorMap(D.data, ℂ^4)
+    t = rand(V1 ← V1 ⊗ V2)
+    t0 = TensorMap(t)
+    @test TensorMap(lmul!(D, t)) ≈ D′ * t0
+    t = rand(V1 ⊗ V2 ← V1)
+    t0 = TensorMap(t)
+    @test TensorMap(rmul!(t, D)) ≈ t0 * D′
+
+    t = rand(V1 ← V1 ⊗ V2)
+    b = block(t, first(blocksectors(t)))
+    d = BlockTensorKit.copy_dense!(BlockTensorKit.similar_dense(b), b) .* 2
+    BlockTensorKit.copy_blocks!(b, d)
+    @test BlockTensorKit.copy_dense!(BlockTensorKit.similar_dense(b), b) == d
+end

From 6d9efe01eeeb3352373057c27b7ad38bf60e82d2 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Thu, 8 Oct 2026 19:36:17 -0400
Subject: [PATCH 04/12] Drop duplicate block-tensor `==` in favour of the
 entrywise one in linalg.jl

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/tensors/abstractblocktensor/abstracttensormap.jl | 6 ------
 1 file changed, 6 deletions(-)

diff --git a/src/tensors/abstractblocktensor/abstracttensormap.jl b/src/tensors/abstractblocktensor/abstracttensormap.jl
index 1ed8255..1cea560 100644
--- a/src/tensors/abstractblocktensor/abstracttensormap.jl
+++ b/src/tensors/abstractblocktensor/abstracttensormap.jl
@@ -125,12 +125,6 @@ end
 Base.getindex(iter::TK.BlockIterator{<:AbstractBlockTensorMap}, c::Sector) =
     _block(iter.t, last(iter.structure), c)
 
-function Base.:(==)(t₁::AbstractBlockTensorMap, t₂::AbstractBlockTensorMap)
-    (codomain(t₁) == codomain(t₂) && domain(t₁) == domain(t₂)) || return false
-    b₁, b₂ = TK.blocks(t₁), TK.blocks(t₂)
-    return all(c -> b₁[c] == b₂[c], blocksectors(t₁))
-end
-
 function TensorKit.foreachblock(f, t::AbstractBlockTensorMap; scheduler = nothing)
     foreach(TK.blocks(t)) do (c, b)
         return f(c, (b,))

From efa177ad7da0d1027428cb15e50ea42bc77ee58d Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Thu, 8 Oct 2026 19:39:24 -0400
Subject: [PATCH 05/12] Copy subblocks once per block in TensorMap <-> block
 tensor conversion
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

`_copy_subblocks!` looped over destination fusion trees and, for each, over
all source blocks calling `v[f₁, f₂]` (a structure cache lookup) even when the
slice was empty, rebuilding `blockedrange` axes per tree. The reverse
constructors `BlockTensorMap(t, space)` / `SparseBlockTensorMap(t, space)`
indexed the block tensor per fusion tree, building a mortar over all blocks.

Both directions now share one generic body (branching on `issparse`): loop
over the blocks once, iterate each block's `subblocks` once, and copy to/from
the matching slice of the full tensor's subblock, looked up through a single
`subblocks` iterator. Per-leg block offsets are precomputed per uncoupled
sector. Trivial sectors take the cheap `subblock` path to avoid a structure
lookup per block. Sparse targets keep the previous semantics of storing every
block with a nonempty space.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/tensors/abstractblocktensor/conversion.jl | 43 +++++++++++++++----
 src/tensors/blocktensor.jl                    | 11 +----
 src/tensors/sparseblocktensor.jl              |  6 +--
 test/abstracttensor/blocktensor.jl            | 15 +++++++
 test/abstracttensor/sparseblocktensor.jl      | 17 ++++++++
 5 files changed, 68 insertions(+), 24 deletions(-)

diff --git a/src/tensors/abstractblocktensor/conversion.jl b/src/tensors/abstractblocktensor/conversion.jl
index e7e742b..c55d0bc 100644
--- a/src/tensors/abstractblocktensor/conversion.jl
+++ b/src/tensors/abstractblocktensor/conversion.jl
@@ -1,16 +1,41 @@
 # Conversion
 # ----------
 
-function _copy_subblocks!(tdst, tsrc)
-    N₁, N₂ = numout(tsrc), numin(tsrc)
-    offsets = ntuple(i -> _sumspace_offsets(i <= N₁ ? codomain(tsrc)[i] : domain(tsrc)[i - N₁]), N₁ + N₂)
-    # a single `SubblockIterator` holds the structure: `tdst[f₁, f₂]` would look it up in TensorKit's locked cache every time
-    dstblocks = subblocks(tdst)
-    for (k, v) in nonzero_pairs(tsrc), ((f₁, f₂), src) in subblocks(v)
-        ranges = map(offsets, (f₁.uncoupled..., f₂.uncoupled...), Tuple(k)) do o, c, kᵢ
-            return (o[c][kᵢ] + 1):o[c][kᵢ + 1]
+function _subblock_ranges(offsets, (f₁, f₂), I::CartesianIndex)
+    return map(offsets, (f₁.uncoupled..., f₂.uncoupled...), Tuple(I)) do o, c, kᵢ
+        return (o[c][kᵢ] + 1):o[c][kᵢ + 1]
+    end
+end
+
+function _subblock_pairs(t::AbstractTensorMap)
+    sectortype(t) === Trivial || return subblocks(t)
+    f = TK.trivial_fusiontree(t)
+    return (f => subblock(t, f),)
+end
+function _subblock_getter(t::AbstractTensorMap)
+    return sectortype(t) === Trivial ? Base.Fix1(subblock, t) : Base.Fix1(getindex, subblocks(t))
+end
+
+function _copy_subblocks!(tdst::TensorMap, tsrc::AbstractBlockTensorMap)
+    offsets = map(_sumspace_offsets, eachspace(tsrc).sumspaces)
+    dst = _subblock_getter(tdst)
+    for (I, v) in nonzero_pairs(tsrc), (f, src) in _subblock_pairs(v)
+        copy!(view(dst(f), _subblock_ranges(offsets, f, I)...), src)
+    end
+    return tdst
+end
+
+function _copy_subblocks!(tdst::AbstractBlockTensorMap, tsrc::AbstractTensorMap)
+    offsets = map(_sumspace_offsets, eachspace(tdst).sumspaces)
+    src = _subblock_getter(tsrc)
+    for (I, V) in pairs(eachspace(tdst))
+        v = similar(eltype(tdst), V)
+        vblocks = _subblock_pairs(v)
+        issparse(tdst) && all(isempty ∘ last, vblocks) && continue
+        for (f, dst) in vblocks
+            copy!(dst, view(src(f), _subblock_ranges(offsets, f, I)...))
         end
-        copy!(dstblocks[(f₁, f₂)][ranges...], src)
+        @inbounds tdst[I] = v
     end
     return tdst
 end
diff --git a/src/tensors/blocktensor.jl b/src/tensors/blocktensor.jl
index 3d2b6d9..58146a9 100644
--- a/src/tensors/blocktensor.jl
+++ b/src/tensors/blocktensor.jl
@@ -112,16 +112,7 @@ end
 # AbstractTensorMap -> BlockTensorMap
 function BlockTensorMap(t::AbstractTensorMap, space::TensorMapSumSpace)
     TT = tensormaptype(spacetype(t), numout(t), numin(t), storagetype(t))
-    tdst = BlockTensorMap{TT}(undef, space)
-    @inbounds for (f₁, f₂) in fusiontrees(t)
-        dst = tdst[f₁, f₂]
-        src = t[f₁, f₂]
-        for block_index in Iterators.product(blockaxes(dst)...)
-            indices = getindex.(axes(dst), block_index)
-            dst[block_index...] .= @view src[indices...]
-        end
-    end
-    return tdst
+    return _copy_subblocks!(BlockTensorMap{TT}(undef_blocks, space), t)
 end
 
 # Convenience constructors
diff --git a/src/tensors/sparseblocktensor.jl b/src/tensors/sparseblocktensor.jl
index d919caf..4d15032 100644
--- a/src/tensors/sparseblocktensor.jl
+++ b/src/tensors/sparseblocktensor.jl
@@ -85,11 +85,7 @@ end
 # AbstractTensorMap -> SparseBlockTensorMap
 function SparseBlockTensorMap(t::AbstractTensorMap, space::TensorMapSumSpace)
     TT = tensormaptype(spacetype(t), numout(t), numin(t), storagetype(t))
-    tdst = SparseBlockTensorMap{TT}(undef, space)
-    for (f₁, f₂) in fusiontrees(tdst)
-        tdst[f₁, f₂] = t[f₁, f₂]
-    end
-    return tdst
+    return _copy_subblocks!(SparseBlockTensorMap{TT}(undef_blocks, space), t)
 end
 
 # Utility constructors
diff --git a/test/abstracttensor/blocktensor.jl b/test/abstracttensor/blocktensor.jl
index 5658e3d..13ebf7f 100644
--- a/test/abstracttensor/blocktensor.jl
+++ b/test/abstracttensor/blocktensor.jl
@@ -99,6 +99,21 @@ end
     @test t1 ≈ t1″
 end
 
+@testset "TensorMap conversion: sectors and empty blocks" begin
+    Vs = (
+        Vect[U1Irrep](0 => 2, 1 => 1) ⊞ Vect[U1Irrep](-1 => 1, 2 => 1) ⊞ Vect[U1Irrep](1 => 2),
+        Vect[SU2Irrep](0 => 1, 1 // 2 => 2) ⊞ Vect[SU2Irrep](1 => 1) ⊞ Vect[SU2Irrep](1 // 2 => 1, 3 // 2 => 1),
+        ℂ^0 ⊞ ℂ^2 ⊞ ℂ^3,
+    )
+    for V in Vs, W in (V ⊗ V' ← V ⊗ V, V ⊗ V ⊗ V' ← one(V)), T in scalartypes
+        t = randn(T, W)
+        t′ = @constinferred TensorMap(t)
+        @test norm(t) ≈ norm(t′)
+        @test BlockTensorMap(t′, W) ≈ t
+        @test SparseBlockTensorMap(t′, W) ≈ t
+    end
+end
+
 @testset "Adapt" begin
     W = V1 ⊗ V2 ⊗ V3 ← V4 ⊗ V5
     t1 = rand(Float32, W)
diff --git a/test/abstracttensor/sparseblocktensor.jl b/test/abstracttensor/sparseblocktensor.jl
index 4642d06..3f58892 100644
--- a/test/abstracttensor/sparseblocktensor.jl
+++ b/test/abstracttensor/sparseblocktensor.jl
@@ -91,6 +91,23 @@ end
     end
 end
 
+@testset "TensorMap conversion: sectors and empty blocks" begin
+    Vs = (
+        Vect[U1Irrep](0 => 2, 1 => 1) ⊞ Vect[U1Irrep](-1 => 1, 2 => 1) ⊞ Vect[U1Irrep](1 => 2),
+        Vect[SU2Irrep](0 => 1, 1 // 2 => 2) ⊞ Vect[SU2Irrep](1 => 1) ⊞ Vect[SU2Irrep](1 // 2 => 1, 3 // 2 => 1),
+        ℂ^0 ⊞ ℂ^2 ⊞ ℂ^3,
+    )
+    for V in Vs, W in (V ⊗ V' ← V ⊗ V, V ⊗ V ⊗ V' ← one(V)), T in scalartypes
+        t = sprand(T, W, 0.5)
+        t′ = @constinferred TensorMap(t)
+        @test norm(t) ≈ norm(t′)
+        t″ = SparseBlockTensorMap(t′, W)
+        @test t″ ≈ t
+        @test all(I -> dim(eachspace(t)[I]) > 0, nonzero_keys(t″))
+        @test TensorMap(BlockTensorMap(t)) ≈ t′
+    end
+end
+
 @testset "Adapt" begin
     W = V1 ⊗ V2 ⊗ V3 ← V4 ⊗ V5
     t1 = sprand(Float32, W, 0.5)

From 3d00a5b2f169baa69024959e6d575ea5e37bf150 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Thu, 8 Oct 2026 19:40:34 -0400
Subject: [PATCH 06/12] Write missing sparse entries directly instead of
 materializing zeros
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Blockwise add/scale/permute/braid/transpose/tensoradd!/trace_permute! into a
SparseBlockTensorMap used to read missing entries through getindex, which
allocates a zero-filled block that is then added into. Missing entries are now
allocated uninitialized (TO.tensoralloc, honouring the allocator) and the kernel
is called with β = Zero(); present entries receive β directly. Entries of the
destination that the (sparse) source does not touch are scaled up front, or
dropped for a sparse destination when β == 0, so preallocated destinations keep
their buffers.

VI.scale!(ty, tx, α) is now add!(ty, tx, α, Zero()), and VI.add allocates its
destination with similar instead of zerovector, since every entry is overwritten.

Remove the 8-argument TK.add_transform! methods: TensorKit 0.17.2's
add_transform! takes a conjsrc argument, so these were unreachable (and called
kernels that no longer exist). Merge the AdjointTensorMap-source permute!/
transpose!/braid! methods into the generic ones; the adjoint permute!/transpose!
referenced undefined p₁/p₂ and indexed with p instead of the linear permutation.
Fix SparseTensorArray numout/numin on types (referenced undefined A), and route
braid!(::TensorMap, ::BlockTensorMap) through braid! instead of deprecated add_braid!.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/auxiliary/sparsetensorarray.jl |   5 +-
 src/tensors/indexmanipulations.jl  | 107 ++++-------------------------
 src/tensors/tensoroperations.jl    |  14 ++--
 src/tensors/vectorinterface.jl     |  47 +++++++++----
 4 files changed, 58 insertions(+), 115 deletions(-)

diff --git a/src/auxiliary/sparsetensorarray.jl b/src/auxiliary/sparsetensorarray.jl
index 7831407..8dfc423 100644
--- a/src/auxiliary/sparsetensorarray.jl
+++ b/src/auxiliary/sparsetensorarray.jl
@@ -44,9 +44,9 @@ TensorKit.codomain(A::SparseTensorArray) = codomain(space(A))
 TensorKit.domain(A::SparseTensorArray) = domain(space(A))
 
 TensorKit.numout(A::SparseTensorArray) = numout(eltype(A))
-TensorKit.numout(::Type{T}) where {T <: SparseTensorArray} = numout(eltype(A))
+TensorKit.numout(::Type{T}) where {T <: SparseTensorArray} = numout(eltype(T))
 TensorKit.numin(A::SparseTensorArray) = numin(eltype(A))
-TensorKit.numin(::Type{T}) where {T <: SparseTensorArray} = numin(eltype(A))
+TensorKit.numin(::Type{T}) where {T <: SparseTensorArray} = numin(eltype(T))
 
 # AbstractArray interface
 # -----------------------
@@ -93,6 +93,7 @@ function Base.delete!(A::SparseTensorArray, I::Vararg{Int, N}) where {N}
 end
 Base.delete!(A::SparseTensorArray, I::CartesianIndex) = delete!(A.data, I)
 Base.empty!(A::SparseTensorArray) = empty!(A.data)
+Base.filter!(f, A::SparseTensorArray) = (filter!(f, A.data); A)
 function Base.haskey(A::SparseTensorArray, I::Vararg{Int, N}) where {N}
     return haskey(A.data, CartesianIndex(I))
 end
diff --git a/src/tensors/indexmanipulations.jl b/src/tensors/indexmanipulations.jl
index 1194120..905c567 100644
--- a/src/tensors/indexmanipulations.jl
+++ b/src/tensors/indexmanipulations.jl
@@ -1,65 +1,4 @@
-@propagate_inbounds function TK.add_transform!(
-        tdst::BlockTensorMap, tsrc::BlockTensorMap, p::Index2Tuple, transformer,
-        α::Number, β::Number, backend, allocator
-    )
-    @boundscheck TK.spacecheck_transform(permute, tdst, tsrc, p)
-
-    dstdata = parent(tdst)
-    srcdata = permutedims(StridedView(parent(tsrc)), (p[1]..., p[2]...))
-
-    @inbounds for I in eachindex(dstdata, srcdata)
-        dstdata[I] = TK.add_transform!(
-            dstdata[I], srcdata[I], p, transformer, α, β, backend, allocator
-        )
-    end
-    return tdst
-end
-@propagate_inbounds function TK.add_transform!(
-        tdst::AbstractBlockTensorMap, tsrc::AbstractBlockTensorMap,
-        p::Index2Tuple, transformer, α::Number, β::Number, backend, allocator
-    )
-    @boundscheck TK.spacecheck_transform(permute, tdst, tsrc, p)
-    scale!(tdst, β)
-    p_lin = (p[1]..., p[2]...)
-    @inbounds for (I, v) in nonzero_pairs(tsrc)
-        I′ = CartesianIndex(TT.getindices(I.I, p_lin))
-        tdst[I′] = TK.add_transform!(
-            tdst[I′], v, p_lin, transformer, α, One(), backend, allocator
-        )
-    end
-    return tdst
-end
-function TK.add_transform!(
-        tdst::AbstractBlockTensorMap, tsrc::AdjointTensorMap{T, S, N₁, N₂, TT},
-        p::Index2Tuple, transformer, α::Number, β::Number, backend, allocator
-    ) where {T, S, N₁, N₂, TT <: AbstractBlockTensorMap}
-    @boundscheck TK.spacecheck_transform(permute, tdst, tsrc, p)
-    scale!(tdst, β)
-    p_lin = (p[1]..., p[2]...)
-    @inbounds for (I, v) in nonzero_pairs(tsrc)
-        I′ = CartesianIndex(TT.getindices(I.I, p_lin))
-        tdst[I′] = TK.add_transform!(
-            tdst[I′], v, p, transformer, α, One(), backend, allocator
-        )
-    end
-    return tdst
-end
-function TK.add_transform!(
-        tdst::TensorMap, tsrc::BlockTensorMap, p::Index2Tuple, transformer,
-        α::Number, β::Number, backend, allocator
-    )
-    return TK.add_transform!(
-        tdst, only(tsrc), p, transformer, α, β, backend, allocator
-    )
-end
-function TK.add_transform!(
-        tdst::BlockTensorMap, tsrc::TensorMap, p::Index2Tuple, transformer,
-        α::Number, β::Number, backend, allocator
-    )
-    return TK.add_transform!(
-        only(tdst), tsrc, p, transformer, α, β, backend, allocator
-    )
-end
+const BlockTensorOrAdjoint = Union{AbstractBlockTensorMap, AdjointBlockTensorMap}
 
 # we need to capture the other functions earlier to enjoy the fast transformers...
 for f in (:permute, :transpose)
@@ -80,28 +19,17 @@ for f in (:permute, :transpose)
             return tdst
         end
         function TK.$f!(
-                tdst::AbstractBlockTensorMap, tsrc::AbstractBlockTensorMap,
+                tdst::AbstractBlockTensorMap, tsrc::BlockTensorOrAdjoint,
                 p::Index2Tuple, α::Number, β::Number, backend::AbstractBackend, allocator
             )
             @boundscheck TK.spacecheck_transform(TK.$f, tdst, tsrc, p)
-            scale!(tdst, β)
             p_lin = (p[1]..., p[2]...)
+            _scale_untouched!(tdst, tsrc, β, I -> CartesianIndex(TT.getindices(I.I, p_lin)))
             @inbounds for (I, v) in nonzero_pairs(tsrc)
                 I′ = CartesianIndex(TT.getindices(I.I, p_lin))
-                tdst[I′] = TK.$f!(tdst[I′], v, p, α, One(), backend, allocator)
-            end
-            return tdst
-        end
-        function TK.$f!(
-                tdst::AbstractBlockTensorMap, tsrc::AdjointTensorMap{T, S, N₁, N₂, TT},
-                p::Index2Tuple, α::Number, β::Number, backend::AbstractBackend, allocator
-            ) where {T, S, N₁, N₂, TT <: AbstractBlockTensorMap}
-            @boundscheck TK.spacecheck_transform(TK.$f, tdst, tsrc, p)
-            scale!(tdst, β)
-            p_lin = (p[1]..., p[2]...)
-            @inbounds for (I, v) in nonzero_pairs(tsrc)
-                I′ = CartesianIndex(TT.getindices(I.I, p))
-                tdst[I′] = TK.$f!(tdst[I′], v, (p₁, p₂), α, One(), backend, allocator)
+                _addblock!(tdst, I′, β, allocator) do c, γ
+                    TK.$f!(c, v, p, α, γ, backend, allocator)
+                end
             end
             return tdst
         end
@@ -138,29 +66,18 @@ end
     return tdst
 end
 @propagate_inbounds function TK.braid!(
-        tdst::AbstractBlockTensorMap, tsrc::AbstractBlockTensorMap,
+        tdst::AbstractBlockTensorMap, tsrc::BlockTensorOrAdjoint,
         p::Index2Tuple, levels::IndexTuple,
         α::Number, β::Number, backend::AbstractBackend, allocator
     )
     @boundscheck TK.spacecheck_transform(braid, tdst, tsrc, p, levels)
-    scale!(tdst, β)
-    p_lin = (p[1]..., p[2]...)
-    @inbounds for (I, v) in nonzero_pairs(tsrc)
-        I′ = CartesianIndex(TT.getindices(I.I, p_lin))
-        tdst[I′] = TK.braid!(tdst[I′], v, p, levels, α, One(), backend, allocator)
-    end
-    return tdst
-end
-function TK.braid!(
-        tdst::AbstractBlockTensorMap, tsrc::AdjointTensorMap{T, S, N₁, N₂, TT},
-        p::Index2Tuple, levels::IndexTuple, α::Number, β::Number, backend::AbstractBackend, allocator
-    ) where {T, S, N₁, N₂, TT <: AbstractBlockTensorMap}
-    @boundscheck TK.spacecheck_transform(braid, tdst, tsrc, p, levels)
-    scale!(tdst, β)
     p_lin = (p[1]..., p[2]...)
+    _scale_untouched!(tdst, tsrc, β, I -> CartesianIndex(TT.getindices(I.I, p_lin)))
     @inbounds for (I, v) in nonzero_pairs(tsrc)
         I′ = CartesianIndex(TT.getindices(I.I, p_lin))
-        tdst[I′] = TK.braid!(tdst[I′], v, p, levels, α, One(), backend, allocator)
+        _addblock!(tdst, I′, β, allocator) do c, γ
+            TK.braid!(c, v, p, levels, α, γ, backend, allocator)
+        end
     end
     return tdst
 end
@@ -169,7 +86,7 @@ function TK.braid!(
         p::Index2Tuple, levels::IndexTuple,
         α::Number, β::Number, backend::AbstractBackend, allocator
     )
-    return TK.add_braid!(tdst, only(tsrc), p, levels, α, β, backend, allocator)
+    return TK.braid!(tdst, only(tsrc), p, levels, α, β, backend, allocator)
 end
 function TK.braid!(
         tdst::BlockTensorMap, tsrc::TensorMap,
diff --git a/src/tensors/tensoroperations.jl b/src/tensors/tensoroperations.jl
index a0342c2..0ef7538 100644
--- a/src/tensors/tensoroperations.jl
+++ b/src/tensors/tensoroperations.jl
@@ -21,11 +21,13 @@ function TO.tensoradd!(
         C::AbstractBlockTensorMap, A::AbstractBlockTensorMap, pA::Index2Tuple, conjA::Bool,
         α::Number, β::Number, backend, allocator
     )
-    scale!(C, β)
     p_lin = (pA[1]..., pA[2]...)
+    _scale_untouched!(C, A, β, I -> CartesianIndex(TT.getindices(I.I, p_lin)))
     @inbounds for (I, v) in nonzero_pairs(A)
         I′ = CartesianIndex(TT.getindices(I.I, p_lin))
-        C[I′] = TO.tensoradd!(C[I′], v, pA, conjA, α, One(), backend, allocator)
+        _addblock!(C, I′, β, allocator) do c, γ
+            TO.tensoradd!(c, v, pA, conjA, α, γ, backend, allocator)
+        end
     end
     return C
 end
@@ -145,13 +147,13 @@ function TK.trace_permute!(
                     q₁ = $(q₁), q₂ = $(q₂)"))
     end
 
-    scale!(tdst, β)
+    _prescale!(tdst, β)
     @inbounds for (Isrc, vsrc) in nonzero_pairs(tsrc)
         TT.getindices(Isrc.I, q₁) == TT.getindices(Isrc.I, q₂) || continue
         Idst = CartesianIndex(TT.getindices(Isrc.I, (p₁..., p₂...)))
-        tdst[Idst] = TensorKit.trace_permute!(
-            tdst[Idst], vsrc, (p₁, p₂), (q₁, q₂), α, One(), backend
-        )
+        _addblock!(tdst, Idst) do c, γ
+            TensorKit.trace_permute!(c, vsrc, (p₁, p₂), (q₁, q₂), α, γ, backend)
+        end
     end
     return tdst
 end
diff --git a/src/tensors/vectorinterface.jl b/src/tensors/vectorinterface.jl
index f7fbe66..fd45d69 100644
--- a/src/tensors/vectorinterface.jl
+++ b/src/tensors/vectorinterface.jl
@@ -10,15 +10,8 @@ function VI.scale!(t::AbstractBlockTensorMap, α::Number)
     return t
 end
 
-function VI.scale!(ty::AbstractBlockTensorMap, tx::AbstractBlockTensorMap, α::Number)
-    space(ty) == space(tx) || throw(SpaceMismatch("$(space(ty)) ≠ $(space(tx))"))
-    # entries of ty that are structurally zero in tx have to be zeroed out
-    issparse(tx) && zerovector!(ty)
-    for (I, v) in nonzero_pairs(tx)
-        ty[I] = scale!!(ty[I], v, α)
-    end
-    return ty
-end
+VI.scale!(ty::AbstractBlockTensorMap, tx::AbstractBlockTensorMap, α::Number) =
+    add!(ty, tx, α, Zero())
 
 # add
 # ---
@@ -31,7 +24,7 @@ function VI.add(ty::AbstractBlockTensorMap, tx::AbstractBlockTensorMap, α::Numb
     # This would give issues for example with DiagonalTensorMap + TensorMap
     T = VectorInterface.promote_add(ty, tx, α, β)
     tdst = if typeof(ty) === typeof(tx)
-        zerovector(ty, T)
+        similar(ty, T)
     else
         M = TK.promote_storagetype(TK.similarstoragetype(ty, T), TK.similarstoragetype(tx, T))
         if issparse(ty) && issparse(tx)
@@ -46,13 +39,43 @@ end
 
 function VI.add!(ty::AbstractBlockTensorMap, tx::AbstractBlockTensorMap, α::Number, β::Number)
     space(ty) == space(tx) || throw(SpaceMismatch("$(space(ty)) ≠ $(space(tx))"))
-    isone(β) || scale!(ty, β)
+    _scale_untouched!(ty, tx, β)
     for (I, v) in nonzero_pairs(tx)
-        ty[I] = add!!(ty[I], v, α, One())
+        _addblock!((y, γ) -> add!!(y, v, α, γ), ty, I, β)
     end
     return ty
 end
 
+# blockwise updates
+# -----------------
+_newblock(t::AbstractBlockTensorMap, I, allocator = TO.DefaultAllocator()) =
+    TO.tensoralloc(eltype(t), eachspace(t)[I], Val(false), allocator)
+
+# `t[I] = f(t[I], β)`, where missing entries are allocated uninitialized and passed `Zero()`
+@propagate_inbounds function _addblock!(
+        f, t::AbstractBlockTensorMap, I::CartesianIndex, β::Number = One(), allocator = TO.DefaultAllocator()
+    )
+    t[I] = haskey(t, I) ? f(t[I], β) : f(_newblock(t, I, allocator), Zero())
+    return t
+end
+
+_prescale!(t::AbstractBlockTensorMap, β::Number) =
+    isone(β) ? t : iszero(β) ? zerovector!(t) : scale!(t, β)
+
+# scale the entries of `t` that no entry of `tsrc` is mapped onto by `Imap`
+function _scale_untouched!(t::AbstractBlockTensorMap, tsrc::AbstractTensorMap, β::Number, Imap = identity)
+    (issparse(tsrc) && !isone(β)) || return t
+    touched = Set(Imap(I) for I in nonzero_keys(tsrc))
+    if issparse(t) && iszero(β)
+        filter!(∈(touched) ∘ first, parent(t))
+    else
+        for (I, v) in nonzero_pairs(t)
+            I ∈ touched || scale!(v, β)
+        end
+    end
+    return t
+end
+
 # inner
 # -----
 function VI.inner(x::AbstractBlockTensorMap, y::AbstractBlockTensorMap)

From 4b774eb830c0cd7417a17ef786433377bf840f30 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Thu, 8 Oct 2026 19:42:01 -0400
Subject: [PATCH 07/12] Test permute!/braid!/transpose! into sparse and dense
 block tensors
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Cover sparse/dense destinations with missing and present entries,
β ∈ {Zero(), 0, One(), random}, adjoint sources, fermionic sectors, and
that β = 0 overwrites NaN-filled destinations.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 test/linalg/indexmanipulations.jl | 50 +++++++++++++++++++++++++++++++
 1 file changed, 50 insertions(+)

diff --git a/test/linalg/indexmanipulations.jl b/test/linalg/indexmanipulations.jl
index 8b13789..45e01eb 100644
--- a/test/linalg/indexmanipulations.jl
+++ b/test/linalg/indexmanipulations.jl
@@ -1 +1,51 @@
+using Test, TestExtras
+using TensorKit
+using BlockTensorKit
+using BlockTensorKit: sprand, nonzero_keys, nonzero_length
+using VectorInterface: Zero, One
 
+Vtr = (ℂ^2 ⊞ ℂ^1, (ℂ^1 ⊞ ℂ^2)', ℂ^2 ⊞ ℂ^3, ℂ^2 ⊞ ℂ^2 ⊞ ℂ^1, (ℂ^1 ⊞ ℂ^2)')
+F = Vect[FermionParity]
+Vf = (
+    F(0 => 1, 1 => 1) ⊞ F(0 => 2), F(1 => 1)' ⊞ F(0 => 1, 1 => 2)', F(0 => 2, 1 => 1),
+    F(0 => 1) ⊞ F(1 => 1) ⊞ F(0 => 1, 1 => 1), F(0 => 1, 1 => 1)' ⊞ F(1 => 2)',
+)
+
+maketensor(sparse, T, W) = sparse ? sprand(T, W, 0.4) : randn(T, W)
+dense(t) = convert(TensorMap, t)
+
+# `(p₁..., reverse(p₂)...)` is a cyclic rotation of the planar order of a 3 ← 2 tensor
+perms(f) = f === transpose! ? ((2, 3), (1, 4, 5)) : ((3, 1), (5, 2, 4))
+transform!(f, C, A, p, α, β) = f === braid! ? f(C, A, p, (1, 2, 3, 4, 5), α, β) : f(C, A, p, α, β)
+
+@testset "$f ($(sectortype(V[1])), C sparse = $Csparse, A sparse = $Asparse, adjoint = $adj)" for V in (Vtr, Vf),
+        f in (permute!, braid!, transpose!), Csparse in (false, true), Asparse in (false, true), adj in (false, true)
+    T = ComplexF64
+    W = V[1] ⊗ V[2] ⊗ V[3] ← V[4] ⊗ V[5]
+    A = maketensor(Asparse, T, adj ? W' : W)
+    A = adj ? A' : A
+    Ad = adj ? dense(A')' : dense(A)
+    p = perms(f)
+    Wdst = permute(space(A), p)
+    for β in (Zero(), zero(T), One(), randn(T))
+        α = randn(T)
+        C = maketensor(Csparse, T, Wdst)
+        Cd = transform!(f, dense(C), Ad, p, α, β)
+        @test transform!(f, C, A, p, α, β) === C
+        @test dense(C) ≈ Cd
+        if iszero(β) && Csparse && Asparse
+            @test nonzero_length(C) == nonzero_length(A)
+        end
+    end
+    C = scale!(maketensor(Csparse, T, Wdst), NaN)
+    @test dense(transform!(f, C, A, p, One(), Zero())) ≈ transform!(f, zerovector(dense(C)), Ad, p, One(), Zero())
+end
+
+@testset "permute of adjoint ($(sectortype(V[1])), sparse = $sparse)" for V in (Vtr, Vf), sparse in (false, true)
+    W = V[1] ⊗ V[2] ⊗ V[3] ← V[4] ⊗ V[5]
+    t = maketensor(sparse, ComplexF64, W)
+    p = ((3, 1), (5, 2, 4))
+    @test dense(permute(t', p)) ≈ permute(dense(t)', p)
+    C = similar(t, permute(space(t'), p))
+    @test dense(permute!(C, t', p)) ≈ permute(dense(t)', p)
+end

From 4ef2eb0563e6b817ece9941b1ec8f065a9fd7fa5 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Thu, 8 Oct 2026 19:54:22 -0400
Subject: [PATCH 08/12] Simplify `_zeroblock` signature and drop explanatory
 comments

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/tensors/abstractblocktensor/abstracttensormap.jl | 8 ++++----
 src/tensors/vectorinterface.jl                       | 2 --
 2 files changed, 4 insertions(+), 6 deletions(-)

diff --git a/src/tensors/abstractblocktensor/abstracttensormap.jl b/src/tensors/abstractblocktensor/abstracttensormap.jl
index 1cea560..8a64b1e 100644
--- a/src/tensors/abstractblocktensor/abstracttensormap.jl
+++ b/src/tensors/abstractblocktensor/abstracttensormap.jl
@@ -67,7 +67,7 @@ function _block(t::AbstractBlockTensorMap, entries, c::Sector)
             _cachedblock(entries[I], c)
         else
             i, j = Tuple(CartesianIndices(allblocks)[k])
-            _zeroblock(eltype(t), rowdims[i], coldims[j], () -> eachspace(t)[I], c)
+            _zeroblock(t, I, rowdims[i], coldims[j], c)
         end
     end
     return mortar(allblocks, rowdims, coldims)
@@ -86,9 +86,9 @@ _blockcache(x::AbstractTensorMap) = (b = TK.blocks(x); b isa TK.BlockIterator ?
 _cachedblock(x, c::Sector) = block(x, c)
 _cachedblock(b::TK.BlockIterator, c::Sector) = b[c]
 
-function _zeroblock(::Type{TT}, d₁::Int, d₂::Int, getspace, c::Sector) where {TT <: AbstractTensorMap}
-    TT <: TensorMap || return block(zerovector!(similar(TT, getspace())), c)
-    data = zerovector!(similar(storagetype(TT), d₁ * d₂))
+function _zeroblock(t::AbstractBlockTensorMap, I, d₁::Int, d₂::Int, c::Sector)
+    eltype(t) <: TensorMap || return block(zerovector!(similar(eltype(t), eachspace(t)[I])), c)
+    data = zerovector!(similar(storagetype(t), d₁ * d₂))
     return reshape(view(data, 1:(d₁ * d₂)), (d₁, d₂))
 end
 
diff --git a/src/tensors/vectorinterface.jl b/src/tensors/vectorinterface.jl
index fd45d69..1f56e19 100644
--- a/src/tensors/vectorinterface.jl
+++ b/src/tensors/vectorinterface.jl
@@ -51,7 +51,6 @@ end
 _newblock(t::AbstractBlockTensorMap, I, allocator = TO.DefaultAllocator()) =
     TO.tensoralloc(eltype(t), eachspace(t)[I], Val(false), allocator)
 
-# `t[I] = f(t[I], β)`, where missing entries are allocated uninitialized and passed `Zero()`
 @propagate_inbounds function _addblock!(
         f, t::AbstractBlockTensorMap, I::CartesianIndex, β::Number = One(), allocator = TO.DefaultAllocator()
     )
@@ -62,7 +61,6 @@ end
 _prescale!(t::AbstractBlockTensorMap, β::Number) =
     isone(β) ? t : iszero(β) ? zerovector!(t) : scale!(t, β)
 
-# scale the entries of `t` that no entry of `tsrc` is mapped onto by `Imap`
 function _scale_untouched!(t::AbstractBlockTensorMap, tsrc::AbstractTensorMap, β::Number, Imap = identity)
     (issparse(tsrc) && !isone(β)) || return t
     touched = Set(Imap(I) for I in nonzero_keys(tsrc))

From c4f505431a3a577348d26199435f0d6e957fcf23 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Fri, 9 Oct 2026 10:04:06 -0400
Subject: [PATCH 09/12] Test `copy_blocks!` write-back on JLArray storage

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 test/linalg/factorizations.jl | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/test/linalg/factorizations.jl b/test/linalg/factorizations.jl
index 9a1b457..4de3304 100644
--- a/test/linalg/factorizations.jl
+++ b/test/linalg/factorizations.jl
@@ -5,6 +5,8 @@ using BlockTensorKit
 using Random
 using Combinatorics
 using LinearAlgebra
+using Adapt
+using JLArrays
 
 Vtr = (
     SumSpace(ℂ^3),
@@ -544,4 +546,10 @@ end
     d = BlockTensorKit.copy_dense!(BlockTensorKit.similar_dense(b), b) .* 2
     BlockTensorKit.copy_blocks!(b, d)
     @test BlockTensorKit.copy_dense!(BlockTensorKit.similar_dense(b), b) == d
+
+    tg = adapt(JLVector{Float64}, t)
+    c = first(blocksectors(t))
+    BlockTensorKit.copy_blocks!(block(t, c), d .* 3)
+    BlockTensorKit.copy_blocks!(block(tg, c), JLArray(d .* 3))
+    @test adapt(Vector{Float64}, tg) ≈ t
 end

From 1d38cdf3469206b29bf4c5e87314a559b81bbe63 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Fri, 9 Oct 2026 10:14:36 -0400
Subject: [PATCH 10/12] Implement block and subblock access through cached
 `blocks`/`subblocks` iterators

`subblocks(t)` now caches each entry's subblock structure, mirroring `blocks(t)`;
`block`/`subblock` share the same per-entry accessors, and the conversion
code reuses them instead of its own getter.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 .../abstractblocktensor/abstracttensormap.jl  | 115 ++++++++++--------
 src/tensors/abstractblocktensor/conversion.jl |   7 +-
 test/abstracttensor/sparseblocktensor.jl      |   5 +
 3 files changed, 70 insertions(+), 57 deletions(-)

diff --git a/src/tensors/abstractblocktensor/abstracttensormap.jl b/src/tensors/abstractblocktensor/abstracttensormap.jl
index 8a64b1e..39a08dd 100644
--- a/src/tensors/abstractblocktensor/abstracttensormap.jl
+++ b/src/tensors/abstractblocktensor/abstracttensormap.jl
@@ -10,28 +10,6 @@ end
 
 eachspace(t::AbstractBlockTensorMap) = SumSpaceIndices(space(t))
 
-@inline function TensorKit.subblock(
-        t::AbstractBlockTensorMap, (f₁, f₂)::Tuple{FusionTree, FusionTree}
-    )
-    sectortype(t) === sectortype(f₁) === sectortype(f₂) ||
-        throw(SectorMismatch("Not a valid sectortype for this tensor"))
-    numout(t) == length(f₁) && numin(t) == length(f₂) ||
-        throw(DimensionMismatch("Invalid number of fusiontree legs for this tensor"))
-
-    W = eachspace(t)
-    subblocks = map(eachindex(IndexCartesian(), t)) do I
-        V = W[I]
-        sz = (dims(codomain(V), f₁.uncoupled)..., dims(domain(V), f₂.uncoupled)...)
-        if prod(sz) == 0 || (issparse(t) && !haskey(parent(t), I))
-            data = zerovector!(similar(storagetype(t), prod(sz)))
-            return sreshape(StridedView(data), sz)
-        else
-            return t[I][f₁, f₂]
-        end
-    end
-
-    return mortar(subblocks)
-end
 @inline function Base.setindex!(
         t::AbstractBlockTensorMap, v::AbstractBlockArray, f₁::FusionTree, f₂::FusionTree
     )
@@ -55,8 +33,53 @@ end
     return setindex!(t, v′, f₁, f₂)
 end
 
+# Blocks and subblocks
+# --------------------
+_entrycache(f, t::AbstractBlockTensorMap) =
+    issparse(t) ? Dict(I => f(x) for (I, x) in nonzero_pairs(t)) : map(f, parent(t))
+
+_blockcache(x) = x
+_blockcache(x::AbstractTensorMap) = (b = TK.blocks(x); b isa TK.BlockIterator ? b : x)
+_cachedblock(x, c::Sector) = block(x, c)
+_cachedblock(b::TK.BlockIterator, c::Sector) = b[c]
+
+_subblockcache(x) = x
+_subblockcache(x::AbstractTensorMap) = sectortype(x) === Trivial ? x : TK.subblocks(x)
+_cachedsubblock(x, f) = subblock(x, f)
+_cachedsubblock(s::TK.SubblockIterator, f) = s[f]
+
+TensorKit.blocks(t::AbstractBlockTensorMap) =
+    TK.BlockIterator(t, (blocksectors(t), _entrycache(_blockcache, t)))
 TensorKit.block(t::AbstractBlockTensorMap, c::Sector) = _block(t, parent(t), c)
 
+TensorKit.subblocks(t::AbstractBlockTensorMap) =
+    TK.SubblockIterator(t, (TK.fusiontrees(t), _entrycache(_subblockcache, t)))
+TensorKit.subblock(t::AbstractBlockTensorMap, f::Tuple{FusionTree, FusionTree}) = _subblock(t, parent(t), f)
+
+Base.length(iter::TK.BlockIterator{<:AbstractBlockTensorMap}) = length(first(iter.structure))
+Base.isdone(iter::TK.BlockIterator{<:AbstractBlockTensorMap}, state...) =
+    Base.isdone(first(iter.structure), state...)
+function Base.iterate(iter::TK.BlockIterator{<:AbstractBlockTensorMap}, state...)
+    next = iterate(first(iter.structure), state...)
+    isnothing(next) && return next
+    c, newstate = next
+    return c => iter[c], newstate
+end
+Base.getindex(iter::TK.BlockIterator{<:AbstractBlockTensorMap}, c::Sector) =
+    _block(iter.t, last(iter.structure), c)
+
+Base.length(iter::TK.SubblockIterator{<:AbstractBlockTensorMap}) = length(first(iter.structure))
+Base.isdone(iter::TK.SubblockIterator{<:AbstractBlockTensorMap}, i::Int = 1) = i > length(iter)
+function Base.iterate(iter::TK.SubblockIterator{<:AbstractBlockTensorMap}, i::Int = 1)
+    i > length(iter) && return nothing
+    f = TK.gettokenvalue(first(iter.structure), i)
+    return f => iter[f], i + 1
+end
+Base.getindex(iter::TK.SubblockIterator{<:AbstractBlockTensorMap}, i::Int) =
+    iter[TK.gettokenvalue(first(iter.structure), i)]
+Base.getindex(iter::TK.SubblockIterator{<:AbstractBlockTensorMap}, f::Tuple{FusionTree, FusionTree}) =
+    _subblock(iter.t, last(iter.structure), f)
+
 function _block(t::AbstractBlockTensorMap, entries, c::Sector)
     sectortype(t) == typeof(c) || throw(SectorMismatch())
     rowdims = _blockdims(codomain(t), c)
@@ -73,6 +96,24 @@ function _block(t::AbstractBlockTensorMap, entries, c::Sector)
     return mortar(allblocks, rowdims, coldims)
 end
 
+function _subblock(t::AbstractBlockTensorMap, entries, (f₁, f₂)::Tuple{FusionTree, FusionTree})
+    sectortype(t) === sectortype(f₁) === sectortype(f₂) ||
+        throw(SectorMismatch("Not a valid sectortype for this tensor"))
+    numout(t) == length(f₁) && numin(t) == length(f₂) ||
+        throw(DimensionMismatch("Invalid number of fusiontree legs for this tensor"))
+    W = eachspace(t)
+    subblocks = map(eachindex(IndexCartesian(), t)) do I
+        V = W[I]
+        sz = (dims(codomain(V), f₁.uncoupled)..., dims(domain(V), f₂.uncoupled)...)
+        if prod(sz) == 0 || (issparse(t) && !haskey(entries, I))
+            return sreshape(StridedView(zerovector!(similar(storagetype(t), prod(sz)))), sz)
+        else
+            return _cachedsubblock(entries[I], (f₁, f₂))
+        end
+    end
+    return mortar(subblocks)
+end
+
 function _blockdims(P::ProductSumSpace{S, N}, c::Sector) where {S, N}
     return vec(
         map(CartesianIndices(map(length, P.spaces))) do I
@@ -81,11 +122,6 @@ function _blockdims(P::ProductSumSpace{S, N}, c::Sector) where {S, N}
     )
 end
 
-_blockcache(x) = x
-_blockcache(x::AbstractTensorMap) = (b = TK.blocks(x); b isa TK.BlockIterator ? b : x)
-_cachedblock(x, c::Sector) = block(x, c)
-_cachedblock(b::TK.BlockIterator, c::Sector) = b[c]
-
 function _zeroblock(t::AbstractBlockTensorMap, I, d₁::Int, d₂::Int, c::Sector)
     eltype(t) <: TensorMap || return block(zerovector!(similar(eltype(t), eachspace(t)[I])), c)
     data = zerovector!(similar(storagetype(t), d₁ * d₂))
@@ -96,12 +132,6 @@ end
 TensorKit.blocksectors(t::AbstractBlockTensorMap) = blocksectors(space(t))
 TensorKit.hasblock(t::AbstractBlockTensorMap, c::Sector) = c in blocksectors(t)
 
-function TensorKit.blocks(t::AbstractBlockTensorMap)
-    entries = issparse(t) ? Dict(I => _blockcache(x) for (I, x) in nonzero_pairs(t)) :
-        map(_blockcache, parent(t))
-    return TK.BlockIterator(t, (blocksectors(t), entries))
-end
-
 Base.@assume_effects :foldable function _blockeltype(::Type{TT}) where {TT <: AbstractBlockTensorMap}
     T = scalartype(TT)
     B = TK.blocktype(eltype(TT))
@@ -112,25 +142,6 @@ function TensorKit.blocktype(::Type{TT}) where {TT <: AbstractBlockTensorMap}
     return BlockMatrix{scalartype(TT), Matrix{_blockeltype(TT)}, BS}
 end
 
-Base.length(iter::TK.BlockIterator{<:AbstractBlockTensorMap}) = length(first(iter.structure))
-Base.isdone(iter::TK.BlockIterator{<:AbstractBlockTensorMap}, state...) =
-    Base.isdone(first(iter.structure), state...)
-function Base.iterate(iter::TK.BlockIterator{<:AbstractBlockTensorMap}, state...)
-    sectors, entries = iter.structure
-    next = iterate(sectors, state...)
-    isnothing(next) && return next
-    c, newstate = next
-    return c => _block(iter.t, entries, c), newstate
-end
-Base.getindex(iter::TK.BlockIterator{<:AbstractBlockTensorMap}, c::Sector) =
-    _block(iter.t, last(iter.structure), c)
-
-function TensorKit.foreachblock(f, t::AbstractBlockTensorMap; scheduler = nothing)
-    foreach(TK.blocks(t)) do (c, b)
-        return f(c, (b,))
-    end
-    return nothing
-end
 function TensorKit.foreachblock(f, t::AbstractBlockTensorMap, ts...; scheduler = nothing)
     tensors = (t, ts...)
     iters = map(_blockcache, tensors)
diff --git a/src/tensors/abstractblocktensor/conversion.jl b/src/tensors/abstractblocktensor/conversion.jl
index c55d0bc..581215b 100644
--- a/src/tensors/abstractblocktensor/conversion.jl
+++ b/src/tensors/abstractblocktensor/conversion.jl
@@ -12,13 +12,10 @@ function _subblock_pairs(t::AbstractTensorMap)
     f = TK.trivial_fusiontree(t)
     return (f => subblock(t, f),)
 end
-function _subblock_getter(t::AbstractTensorMap)
-    return sectortype(t) === Trivial ? Base.Fix1(subblock, t) : Base.Fix1(getindex, subblocks(t))
-end
 
 function _copy_subblocks!(tdst::TensorMap, tsrc::AbstractBlockTensorMap)
     offsets = map(_sumspace_offsets, eachspace(tsrc).sumspaces)
-    dst = _subblock_getter(tdst)
+    dst = Base.Fix1(_cachedsubblock, _subblockcache(tdst))
     for (I, v) in nonzero_pairs(tsrc), (f, src) in _subblock_pairs(v)
         copy!(view(dst(f), _subblock_ranges(offsets, f, I)...), src)
     end
@@ -27,7 +24,7 @@ end
 
 function _copy_subblocks!(tdst::AbstractBlockTensorMap, tsrc::AbstractTensorMap)
     offsets = map(_sumspace_offsets, eachspace(tdst).sumspaces)
-    src = _subblock_getter(tsrc)
+    src = Base.Fix1(_cachedsubblock, _subblockcache(tsrc))
     for (I, V) in pairs(eachspace(tdst))
         v = similar(eltype(tdst), V)
         vblocks = _subblock_pairs(v)
diff --git a/test/abstracttensor/sparseblocktensor.jl b/test/abstracttensor/sparseblocktensor.jl
index 3f58892..2ae01d5 100644
--- a/test/abstracttensor/sparseblocktensor.jl
+++ b/test/abstracttensor/sparseblocktensor.jl
@@ -310,6 +310,11 @@ end
             for (f₁, f₂) in fusiontrees(t)
                 @test t[f₁, f₂] == td[f₁, f₂]
             end
+            sbs = TensorKit.subblocks(t)
+            @test length(sbs) == length(fusiontrees(t))
+            for (i, (f, b)) in enumerate(sbs)
+                @test b == sbs[i] == TensorKit.subblock(td, f)
+            end
             @test t == td
         end
     end

From b2f3f694bd8a4d41b09ad2423f0fce281c9ba968 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Fri, 9 Oct 2026 13:20:01 -0400
Subject: [PATCH 11/12] Make `norm, ==, tr` sparse test independent of an empty
 `sprand` draw

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 test/abstracttensor/sparseblocktensor.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/abstracttensor/sparseblocktensor.jl b/test/abstracttensor/sparseblocktensor.jl
index 2ae01d5..ad415d0 100644
--- a/test/abstracttensor/sparseblocktensor.jl
+++ b/test/abstracttensor/sparseblocktensor.jl
@@ -166,6 +166,7 @@ end
     W = V1 ⊗ V2 ← V1 ⊗ V2
     for T in (Float64, ComplexF64)
         t = sprand(T, W, 0.5)
+        t[1] = randn(T, space(t[1]))
         t′ = convert(TensorMap, t)
         for p in (1, 2, 3, Inf)
             @test norm(t, p) ≈ norm(t′, p)

From 6fda5855231bdf89f21f1add74cabbbb36d77264 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Fri, 9 Oct 2026 14:14:32 -0400
Subject: [PATCH 12/12] Make SumSpace `sectorequal`/`sectorhash` match
 TensorKit's set semantics
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The component-wise versions were sound as a cache key but finer than TensorKit's
"same set of sectors" meaning, reporting e.g. `⊞(A, B)` and `⊞(B, A)` as unequal.
Equality now checks mutual sector containment via `hassector`, and the hash
combines the min and max sector hash, both set functions, without allocating.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/vectorspaces/sumspace.jl  | 16 +++++++---------
 test/vectorspaces/sumspace.jl | 17 ++++++++++++++++-
 2 files changed, 23 insertions(+), 10 deletions(-)

diff --git a/src/vectorspaces/sumspace.jl b/src/vectorspaces/sumspace.jl
index 6cbff49..a859809 100644
--- a/src/vectorspaces/sumspace.jl
+++ b/src/vectorspaces/sumspace.jl
@@ -134,19 +134,17 @@ function TensorKit._sectors(S::SumSpace, ::Type{I}) where {I}
     return sort!(collect(s))
 end
 
+_sectorsubset(S₁::SumSpace, S₂::SumSpace) = all(V -> all(s -> hassector(S₂, s), sectors(V)), S₁.spaces)
 function TensorKit.sectorequal(S₁::SumSpace, S₂::SumSpace)
-    isdual(S₁) == isdual(S₂) && length(S₁) == length(S₂) || return false
-    for (V₁, V₂) in zip(S₁.spaces, S₂.spaces)
-        TK.sectorequal(V₁, V₂) || return false
-    end
-    return true
+    return isdual(S₁) == isdual(S₂) && _sectorsubset(S₁, S₂) && _sectorsubset(S₂, S₁)
 end
 function TensorKit.sectorhash(S::SumSpace, h::UInt)
-    h = hash(isdual(S), h)
-    for V in S.spaces
-        h = TK.sectorhash(V, h)
+    lo, hi = typemax(UInt), typemin(UInt)
+    for V in S.spaces, s in sectors(V)
+        x = hash(s)
+        lo, hi = min(lo, x), max(hi, x)
     end
-    return h
+    return hash(lo, hash(hi, hash(isdual(S), h)))
 end
 
 TensorKit.dim(S::SumSpace, sector::Sector) = sum(v -> dim(v, sector), S.spaces; init = 0)
diff --git a/test/vectorspaces/sumspace.jl b/test/vectorspaces/sumspace.jl
index b227b89..64c1fc6 100644
--- a/test/vectorspaces/sumspace.jl
+++ b/test/vectorspaces/sumspace.jl
@@ -273,7 +273,22 @@ end
     @test TensorKit.sectorequal(V, V′)
     @test TensorKit.sectorhash(V, UInt(1)) == TensorKit.sectorhash(V′, UInt(1))
     @test !TensorKit.sectorequal(V, V')
-    @test !TensorKit.sectorequal(V, V ⊞ V)
+    for V″ in (V ⊞ V, V[[3, 1, 2]], ⊞(⊕(V)), ⊞(Vect[U1Irrep](0 => 1, 1 => 1, -1 => 1, 2 => 1)))
+        @test TensorKit.sectorequal(V, V″)
+        @test TensorKit.sectorhash(V, UInt(1)) == TensorKit.sectorhash(V″, UInt(1))
+    end
+    @test !TensorKit.sectorequal(V, V[1:2])
+
+    fallback_equal(a, b) = invoke(TensorKit.sectorequal, Tuple{ElementarySpace, ElementarySpace}, a, b)
+    randcomp() = Vect[U1Irrep](c => rand(0:1) for c in -2:2)
+    for _ in 1:500
+        dual = rand(Bool)
+        a, b = (SumSpace([randcomp() for _ in 1:rand(1:4)]) for _ in 1:2)
+        a, b = dual ? (a', b') : (a, b)
+        @test TensorKit.sectorequal(a, b) == fallback_equal(a, b)
+        TensorKit.sectorequal(a, b) &&
+            @test TensorKit.sectorhash(a, UInt(1)) == TensorKit.sectorhash(b, UInt(1))
+    end
     @test axes(V, U1Irrep(1)) == [3, 4, 10]
     for W in (V ⊗ V ← V, V ⊗ V' ← V ⊗ V, (V ⊞ V) ⊗ V' ← V)
         @test blocksectors(W) == blocksectors(convert(TensorMapSpace, W))