From 9cf39b1331bcedcbaaed85007f0920f15dd8cedf Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Thu, 8 Oct 2026 17:02:51 -0400
Subject: [PATCH 1/5] Allocate sparse blocks through the allocator in
 `tensoralloc_add`/`tensoralloc_contract`

Sparse block tensors were allocated empty and filled lazily on the heap, while
`tensorfree!` released every block through the allocator. Allocate exactly the
blocks that will be written instead, and pre-insert the product blocks into a
sparse destination so that temporaries derived from it are complete.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/tensors/tensoroperations.jl | 73 +++++++++++++++++++++++++++++++++
 test/linalg/tensoroperations.jl | 52 +++++++++++++++++++++++
 2 files changed, 125 insertions(+)

diff --git a/src/tensors/tensoroperations.jl b/src/tensors/tensoroperations.jl
index a0342c2..3d8a1b2 100644
--- a/src/tensors/tensoroperations.jl
+++ b/src/tensors/tensoroperations.jl
@@ -105,6 +105,79 @@ function TO.tensoralloc(
     return C
 end
 
+# sparse results allocate exactly the blocks that will be written, so `tensorfree!` is consistent
+const BlockOrAdjoint = Union{AbstractBlockTensorMap, AdjointBlockTensorMap}
+
+_tensoralloc(ttype, structure, keys, istemp::Val, allocator) =
+    TO.tensoralloc(ttype, structure, istemp, allocator)
+function _tensoralloc(
+        ttype::Type{<:SparseBlockTensorMap}, structure, keys, istemp::Val, allocator
+    )
+    C = ttype(undef_blocks, structure)
+    Vs = eachspace(C)
+    for I in keys
+        C[I] = TO.tensoralloc(eltype(C), Vs[I], istemp, allocator)
+    end
+    return C
+end
+
+function contract_keys(A, pA::Index2Tuple, B, pB::Index2Tuple, pAB::Index2Tuple)
+    OB = NTuple{length(pB[2]), Int}
+    openB = Dict{NTuple{length(pB[1]), Int}, Vector{OB}}()
+    for IB in nonzero_keys(B)
+        push!(get!(Vector{OB}, openB, TT.getindices(IB.I, pB[1])), TT.getindices(IB.I, pB[2]))
+    end
+    p = TO.linearize(pAB)
+    keys = Set{CartesianIndex{length(p)}}()
+    for IA in nonzero_keys(A)
+        oA = TT.getindices(IA.I, pA[1])
+        for oB in get(openB, TT.getindices(IA.I, pA[2]), ())
+            push!(keys, CartesianIndex(TT.getindices((oA..., oB...), p)))
+        end
+    end
+    return keys
+end
+
+function TO.tensoralloc_add(
+        TC, A::BlockOrAdjoint, pA::Index2Tuple, conjA::Bool,
+        istemp::Val = Val(false), allocator = TO.DefaultAllocator()
+    )
+    ttype = TO.tensoradd_type(TC, A, pA, conjA)
+    structure = TO.tensoradd_structure(A, pA, conjA)
+    p = TO.linearize(pA)
+    keys = (CartesianIndex(TT.getindices(I.I, p)) for I in nonzero_keys(A))
+    return _tensoralloc(ttype, structure, keys, istemp, allocator)
+end
+
+function TO.tensoralloc_contract(
+        TC, A::BlockOrAdjoint, pA::Index2Tuple, conjA::Bool,
+        B::BlockOrAdjoint, pB::Index2Tuple, conjB::Bool,
+        pAB::Index2Tuple, istemp::Val = Val(false), allocator = TO.DefaultAllocator()
+    )
+    ttype = TO.tensorcontract_type(TC, A, pA, conjA, B, pB, conjB, pAB)
+    structure = TO.tensorcontract_structure(A, pA, conjA, B, pB, conjB, pAB)
+    keys = contract_keys(A, pA, B, pB, pAB)
+    return _tensoralloc(ttype, structure, keys, istemp, allocator)
+end
+
+# temporaries derived from `C` must already contain every block the contraction produces
+function TO.tensorcontract!(
+        C::SparseBlockTensorMap,
+        A::BlockOrAdjoint, pA::Index2Tuple, conjA::Bool,
+        B::BlockOrAdjoint, pB::Index2Tuple, conjB::Bool,
+        pAB::Index2Tuple, α::Number, β::Number, backend, allocator
+    )
+    for I in contract_keys(A, pA, B, pB, pAB)
+        getindex!(C, I)
+    end
+    return @invoke TO.tensorcontract!(
+        C::AbstractTensorMap,
+        A::AbstractTensorMap, pA::Index2Tuple, conjA::Bool,
+        B::AbstractTensorMap, pB::Index2Tuple, conjB::Bool,
+        pAB::Index2Tuple, α::Number, β::Number, backend::Any, allocator::Any
+    )
+end
+
 # tensorfree!
 # -----------
 function TO.tensorfree!(t::BlockTensorMap, allocator = TO.DefaultAllocator())
diff --git a/test/linalg/tensoroperations.jl b/test/linalg/tensoroperations.jl
index 67ad59a..a167a9d 100644
--- a/test/linalg/tensoroperations.jl
+++ b/test/linalg/tensoroperations.jl
@@ -150,3 +150,55 @@ end
         @test convert(TensorMap, F1) ≈ F2
     end
 end
+
+##
+
+struct TrackingAllocator
+    live::Set{UInt}
+    foreign::Base.RefValue{Int}
+end
+TrackingAllocator() = TrackingAllocator(Set{UInt}(), Ref(0))
+function TensorOperations.tensoralloc(
+        ::Type{A}, structure, ::Val, allocator::TrackingAllocator
+    ) where {A <: AbstractArray}
+    x = fill!(TensorOperations.tensoralloc(A, structure, Val(false)), NaN)
+    push!(allocator.live, objectid(x))
+    return x
+end
+function TensorOperations.tensorfree!(x::AbstractArray, allocator::TrackingAllocator)
+    if objectid(x) in allocator.live
+        delete!(allocator.live, objectid(x))
+    else
+        allocator.foreign[] += 1
+    end
+    return nothing
+end
+
+@testset "allocator only frees what it allocated ($(sectortype(a)))" for (a, b, c) in (
+        (SumSpace(ℂ^2, ℂ^1), SumSpace(ℂ^1, ℂ^2, ℂ^1), SumSpace(ℂ^2, ℂ^3)),
+        (
+            SumSpace(Vect[FermionParity](0 => 1, 1 => 1), Vect[FermionParity](0 => 2)),
+            SumSpace(Vect[FermionParity](1 => 2), Vect[FermionParity](0 => 1, 1 => 1)),
+            SumSpace(Vect[FermionParity](0 => 1), Vect[FermionParity](1 => 1)),
+        ),
+    )
+    A = sprand(Float64, a ⊗ b ← c, 0.4)
+    B = sprand(Float64, c ⊗ b' ← a, 0.4)
+    D = sprand(Float64, a ← b, 0.4)
+
+    allocator = TrackingAllocator()
+    @tensor allocator = allocator E[x; q] := A[x z; w] * B[w z; y] * D[y; q]
+    @tensor Eref[x; q] := A[x z; w] * B[w z; y] * D[y; q]
+    @test E ≈ Eref
+    @test allocator.foreign[] == 0
+
+    @tensor C[x y; p q] := A[x q; w] * B[w p; y]
+    C = spzeros(Float64, space(C))
+    C[1, 1, 1, 1] = randn(Float64, space(C[1, 1, 1, 1]))
+    Cref = copy(C)
+    allocator = TrackingAllocator()
+    @tensor allocator = allocator C[x y; p q] += A[x q; w] * B[w p; y]
+    @tensor Cref[x y; p q] += A[x q; w] * B[w p; y]
+    @test C ≈ Cref
+    @test allocator.foreign[] == 0
+end

From 1638bf41d84bae7f7e096cab64c5c4074a232ff0 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Thu, 8 Oct 2026 17:32:35 -0400
Subject: [PATCH 2/5] Route all sparse destination blocks through the allocator

- Contract into a sparse `C` through a temporary holding only the product blocks
  unless `C` is a valid BLAS destination, instead of copying all blocks of `C`.
- Allocate missing blocks of sparse destinations lazily through the operation's
  allocator (as non-temporaries), so that adding a differently-patterned tensor
  into an allocator-backed temporary stays consistent.
- Deduplicate the block keys in `tensoralloc_add`, which projects several source
  blocks onto the same output block for traces.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 .../abstractblocktensor/abstractarray.jl      |  2 +
 src/tensors/indexmanipulations.jl             | 20 ++++--
 src/tensors/sparseblocktensor.jl              | 10 +++
 src/tensors/tensoroperations.jl               | 26 ++++++--
 test/linalg/tensoroperations.jl               | 61 +++++++++++++------
 5 files changed, 89 insertions(+), 30 deletions(-)

diff --git a/src/tensors/abstractblocktensor/abstractarray.jl b/src/tensors/abstractblocktensor/abstractarray.jl
index 46ce188..e2581f2 100644
--- a/src/tensors/abstractblocktensor/abstractarray.jl
+++ b/src/tensors/abstractblocktensor/abstractarray.jl
@@ -88,6 +88,8 @@ end
     getindex!(parent(t), I...)
 @propagate_inbounds getindex!(t::AbstractBlockTensorMap, I::CartesianIndex{N}) where {N} =
     getindex!(parent(t), I)
+# missing blocks of a destination are obtained from the allocator of the operation
+@propagate_inbounds getindex!(t::AbstractBlockTensorMap, I::CartesianIndex, allocator) = t[I]
 
 # slicing getindex needs to correctly allocate output blocktensor:
 @propagate_inbounds Base.getindex(t::AbstractBlockTensorMap, indices::Vararg{SliceIndex}) =
diff --git a/src/tensors/indexmanipulations.jl b/src/tensors/indexmanipulations.jl
index 1194120..15f2243 100644
--- a/src/tensors/indexmanipulations.jl
+++ b/src/tensors/indexmanipulations.jl
@@ -24,7 +24,7 @@ end
     @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
+            getindex!(tdst, I′, allocator), v, p_lin, transformer, α, One(), backend, allocator
         )
     end
     return tdst
@@ -39,7 +39,7 @@ function TK.add_transform!(
     @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
+            getindex!(tdst, I′, allocator), v, p, transformer, α, One(), backend, allocator
         )
     end
     return tdst
@@ -88,7 +88,9 @@ for f in (:permute, :transpose)
             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.$f!(tdst[I′], v, p, α, One(), backend, allocator)
+                tdst[I′] = TK.$f!(
+                    getindex!(tdst, I′, allocator), v, p, α, One(), backend, allocator
+                )
             end
             return tdst
         end
@@ -101,7 +103,9 @@ for f in (:permute, :transpose)
             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)
+                tdst[I′] = TK.$f!(
+                    getindex!(tdst, I′, allocator), v, (p₁, p₂), α, One(), backend, allocator
+                )
             end
             return tdst
         end
@@ -147,7 +151,9 @@ end
     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)
+        tdst[I′] = TK.braid!(
+            getindex!(tdst, I′, allocator), v, p, levels, α, One(), backend, allocator
+        )
     end
     return tdst
 end
@@ -160,7 +166,9 @@ function TK.braid!(
     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)
+        tdst[I′] = TK.braid!(
+            getindex!(tdst, I′, allocator), v, p, levels, α, One(), backend, allocator
+        )
     end
     return tdst
 end
diff --git a/src/tensors/sparseblocktensor.jl b/src/tensors/sparseblocktensor.jl
index d919caf..f3739d2 100644
--- a/src/tensors/sparseblocktensor.jl
+++ b/src/tensors/sparseblocktensor.jl
@@ -189,6 +189,16 @@ function Base.delete!(t::SparseBlockTensorMap{TT}, I::Vararg{Int, N}) where {TT,
     return delete!(t, CartesianIndex(I...))
 end
 
+# missing blocks of a destination are obtained from the allocator of the operation
+@propagate_inbounds getindex!(
+    t::SparseBlockTensorMap, I::CartesianIndex, ::TO.DefaultAllocator
+) = getindex!(t, I)
+@propagate_inbounds function getindex!(t::SparseBlockTensorMap, I::CartesianIndex, allocator)
+    haskey(t, I) && return t[I]
+    tI = TO.tensoralloc(eltype(t), eachspace(t)[I], Val(false), allocator)
+    return t[I] = zerovector!(tI)
+end
+
 # Show
 # ----
 function Base.showarg(io::IO, t::SparseBlockTensorMap, toplevel::Bool)
diff --git a/src/tensors/tensoroperations.jl b/src/tensors/tensoroperations.jl
index 3d8a1b2..98d9c28 100644
--- a/src/tensors/tensoroperations.jl
+++ b/src/tensors/tensoroperations.jl
@@ -25,7 +25,9 @@ function TO.tensoradd!(
     p_lin = (pA[1]..., pA[2]...)
     @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)
+        C[I′] = TO.tensoradd!(
+            getindex!(C, I′, allocator), v, pA, conjA, α, One(), backend, allocator
+        )
     end
     return C
 end
@@ -116,7 +118,7 @@ function _tensoralloc(
     C = ttype(undef_blocks, structure)
     Vs = eachspace(C)
     for I in keys
-        C[I] = TO.tensoralloc(eltype(C), Vs[I], istemp, allocator)
+        haskey(C, I) || (C[I] = TO.tensoralloc(eltype(C), Vs[I], istemp, allocator))
     end
     return C
 end
@@ -160,16 +162,30 @@ function TO.tensoralloc_contract(
     return _tensoralloc(ttype, structure, keys, istemp, allocator)
 end
 
-# temporaries derived from `C` must already contain every block the contraction produces
+# contract directly into `C` when it is a valid BLAS destination, otherwise go through a
+# temporary holding only the product blocks instead of a copy of all blocks of `C`
 function TO.tensorcontract!(
         C::SparseBlockTensorMap,
         A::BlockOrAdjoint, pA::Index2Tuple, conjA::Bool,
         B::BlockOrAdjoint, pB::Index2Tuple, conjB::Bool,
         pAB::Index2Tuple, α::Number, β::Number, backend, allocator
     )
-    for I in contract_keys(A, pA, B, pB, pAB)
-        getindex!(C, I)
+    ipAB = TO.oindABinC(pAB, pA, pB)
+    if TO.isblasdestination(C, ipAB) || TO.isblasdestination(C, reverse(ipAB))
+        foreach(I -> getindex!(C, I, allocator), contract_keys(A, pA, B, pB, pAB))
+        return _tensorcontract!(C, A, pA, conjA, B, pB, conjB, pAB, α, β, backend, allocator)
     end
+    N₁ = length(pA[1])
+    pAB′ = (ntuple(identity, N₁), ntuple(i -> N₁ + i, length(pB[2])))
+    AB = TO.tensoralloc_contract(
+        scalartype(C), A, pA, conjA, B, pB, conjB, pAB′, Val(true), allocator
+    )
+    _tensorcontract!(AB, A, pA, conjA, B, pB, conjB, pAB′, One(), Zero(), backend, allocator)
+    TO.tensoradd!(C, AB, pAB, false, α, β, backend, allocator)
+    TO.tensorfree!(AB, allocator)
+    return C
+end
+function _tensorcontract!(C, A, pA, conjA, B, pB, conjB, pAB, α, β, backend, allocator)
     return @invoke TO.tensorcontract!(
         C::AbstractTensorMap,
         A::AbstractTensorMap, pA::Index2Tuple, conjA::Bool,
diff --git a/test/linalg/tensoroperations.jl b/test/linalg/tensoroperations.jl
index a167a9d..502d4ef 100644
--- a/test/linalg/tensoroperations.jl
+++ b/test/linalg/tensoroperations.jl
@@ -154,15 +154,17 @@ end
 ##
 
 struct TrackingAllocator
-    live::Set{UInt}
-    foreign::Base.RefValue{Int}
+    live::Set{UInt} # handed out and not freed yet
+    temps::Set{UInt} # handed out as temporaries
+    foreign::Base.RefValue{Int} # freed arrays that were not handed out by the allocator
 end
-TrackingAllocator() = TrackingAllocator(Set{UInt}(), Ref(0))
+TrackingAllocator() = TrackingAllocator(Set{UInt}(), Set{UInt}(), Ref(0))
 function TensorOperations.tensoralloc(
-        ::Type{A}, structure, ::Val, allocator::TrackingAllocator
-    ) where {A <: AbstractArray}
+        ::Type{A}, structure, ::Val{istemp}, allocator::TrackingAllocator
+    ) where {A <: AbstractArray, istemp}
     x = fill!(TensorOperations.tensoralloc(A, structure, Val(false)), NaN)
     push!(allocator.live, objectid(x))
+    istemp && push!(allocator.temps, objectid(x))
     return x
 end
 function TensorOperations.tensorfree!(x::AbstractArray, allocator::TrackingAllocator)
@@ -173,6 +175,12 @@ function TensorOperations.tensorfree!(x::AbstractArray, allocator::TrackingAlloc
     end
     return nothing
 end
+function check_allocator(f)
+    allocator = TrackingAllocator()
+    @test f(allocator) ≈ f(TensorOperations.DefaultAllocator())
+    @test allocator.foreign[] == 0
+    return @test isdisjoint(allocator.live, allocator.temps)
+end
 
 @testset "allocator only frees what it allocated ($(sectortype(a)))" for (a, b, c) in (
         (SumSpace(ℂ^2, ℂ^1), SumSpace(ℂ^1, ℂ^2, ℂ^1), SumSpace(ℂ^2, ℂ^3)),
@@ -185,20 +193,35 @@ end
     A = sprand(Float64, a ⊗ b ← c, 0.4)
     B = sprand(Float64, c ⊗ b' ← a, 0.4)
     D = sprand(Float64, a ← b, 0.4)
+    E = sprand(Float64, a ← b, 0.4)
+    X = sprand(Float64, a ⊗ b ← a ⊗ b, 0.5)
 
-    allocator = TrackingAllocator()
-    @tensor allocator = allocator E[x; q] := A[x z; w] * B[w z; y] * D[y; q]
-    @tensor Eref[x; q] := A[x z; w] * B[w z; y] * D[y; q]
-    @test E ≈ Eref
-    @test allocator.foreign[] == 0
+    check_allocator() do al
+        return @tensor allocator = al R[x; q] := A[x z; w] * B[w z; y] * D[y; q]
+    end
+    check_allocator() do al # temporary receiving blocks outside of its initial pattern
+        return @tensor allocator = al R[x; q] := A[x z; w] * B[w z; y] * (D[y; q] + 2 * E[y; q])
+    end
+    check_allocator() do al # several blocks of `X` are traced into the same output block
+        return @tensor allocator = al R[x; q] := X[x z; y z] * D[y; q]
+    end
+    check_allocator() do al
+        return @planar allocator = al R[x; q] := A[x z; w] * B[w z; y] * D[y; q]
+    end
 
-    @tensor C[x y; p q] := A[x q; w] * B[w p; y]
-    C = spzeros(Float64, space(C))
-    C[1, 1, 1, 1] = randn(Float64, space(C[1, 1, 1, 1]))
-    Cref = copy(C)
-    allocator = TrackingAllocator()
-    @tensor allocator = allocator C[x y; p q] += A[x q; w] * B[w p; y]
-    @tensor Cref[x y; p q] += A[x q; w] * B[w p; y]
-    @test C ≈ Cref
-    @test allocator.foreign[] == 0
+    # permuted and BLAS-compatible destinations, without and with blocks outside of the product
+    @tensor C1[x y; p q] := A[x q; w] * B[w p; y]
+    @tensor C2[x q; p y] := A[x q; w] * B[w p; y]
+    for C in (spzeros(Float64, space(C1)), sprand(Float64, space(C1), 1.0))
+        check_allocator() do al
+            C′ = copy(C)
+            return @tensor allocator = al C′[x y; p q] += A[x q; w] * B[w p; y]
+        end
+    end
+    for C in (spzeros(Float64, space(C2)), sprand(Float64, space(C2), 1.0))
+        check_allocator() do al
+            C′ = copy(C)
+            return @tensor allocator = al C′[x q; p y] += A[x q; w] * B[w p; y]
+        end
+    end
 end

From cbc716a6b1e6f4d14978e280539bdaecaf7c2479 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Thu, 8 Oct 2026 18:14:36 -0400
Subject: [PATCH 3/5] Track temporaries with an `istemp` flag on
 `SparseBlockTensorMap`

`TO.tensoralloc` returns an empty sparse tensor that remembers whether it is a
temporary. Blocks added later are requested from the allocator with the matching
`istemp`, and `tensorfree!` (and `zerovector!`) only release blocks of temporaries,
so `tensorfree!` never receives arrays that were not allocated as temporaries.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/tensors/sparseblocktensor.jl | 26 ++++++++++++++++++--------
 src/tensors/tensoroperations.jl  | 20 ++++++++++++--------
 src/tensors/vectorinterface.jl   |  6 +++++-
 test/linalg/tensoroperations.jl  | 12 +++++-------
 4 files changed, 40 insertions(+), 24 deletions(-)

diff --git a/src/tensors/sparseblocktensor.jl b/src/tensors/sparseblocktensor.jl
index f3739d2..cc903bd 100644
--- a/src/tensors/sparseblocktensor.jl
+++ b/src/tensors/sparseblocktensor.jl
@@ -2,19 +2,22 @@
     struct SparseBlockTensorMap{TT<:AbstractTensorMap{E,S,N₁,N₂}} <: AbstractBlockTensorMap{E,S,N₁,N₂}
 
 Sparse `SparseBlockTensorMap` type that stores tensors of type `TT` in a sparse dictionary.
+`istemp` marks temporaries of TensorOperations, whose blocks are obtained from and released
+to an allocator.
 """
 struct SparseBlockTensorMap{TT <: AbstractTensorMap, E, S, N₁, N₂, N} <:
     AbstractBlockTensorMap{E, S, N₁, N₂}
     data::Dict{CartesianIndex{N}, TT}
     space::TensorMapSumSpace{S, N₁, N₂}
+    istemp::Bool
 
     # uninitialized constructor
     function SparseBlockTensorMap{TT, E, S, N₁, N₂, N}(
-            ::UndefBlocksInitializer, space::TensorMapSumSpace{S, N₁, N₂}
+            ::UndefBlocksInitializer, space::TensorMapSumSpace{S, N₁, N₂}; istemp::Bool = false
         ) where {E, S, N₁, N₂, N, TT <: AbstractTensorMap{E, S, N₁, N₂}}
         @assert N₁ + N₂ == N "SparseBlockTensorMap: data has wrong number of dimensions"
         data = Dict{CartesianIndex{N}, TT}()
-        return new{TT, E, S, N₁, N₂, N}(data, space)
+        return new{TT, E, S, N₁, N₂, N}(data, space, istemp)
     end
 
     # constructor from data
@@ -22,7 +25,7 @@ struct SparseBlockTensorMap{TT <: AbstractTensorMap, E, S, N₁, N₂, N} <:
             data::Dict{CartesianIndex{N}, TT}, space::TensorMapSumSpace{S, N₁, N₂}
         ) where {E, S, N₁, N₂, N, TT <: AbstractTensorMap{E, S, N₁, N₂}}
         @assert N₁ + N₂ == N "SparseBlockTensorMap: data has wrong number of dimensions"
-        return new{TT, E, S, N₁, N₂, N}(data, space)
+        return new{TT, E, S, N₁, N₂, N}(data, space, false)
     end
 end
 
@@ -34,10 +37,11 @@ end
 
 # uninitialized constructor
 function SparseBlockTensorMap{TT}(
-        ::Union{UndefBlocksInitializer, UndefInitializer}, space::TensorMapSumSpace{S, N₁, N₂}
+        ::Union{UndefBlocksInitializer, UndefInitializer}, space::TensorMapSumSpace{S, N₁, N₂};
+        kwargs...
     ) where {E, S, N₁, N₂, TT <: AbstractTensorMap{E, S, N₁, N₂}}
     N = N₁ + N₂
-    return SparseBlockTensorMap{TT, E, S, N₁, N₂, N}(undef_blocks, space)
+    return SparseBlockTensorMap{TT, E, S, N₁, N₂, N}(undef_blocks, space; kwargs...)
 end
 
 # constructor from data
@@ -194,9 +198,15 @@ end
     t::SparseBlockTensorMap, I::CartesianIndex, ::TO.DefaultAllocator
 ) = getindex!(t, I)
 @propagate_inbounds function getindex!(t::SparseBlockTensorMap, I::CartesianIndex, allocator)
-    haskey(t, I) && return t[I]
-    tI = TO.tensoralloc(eltype(t), eachspace(t)[I], Val(false), allocator)
-    return t[I] = zerovector!(tI)
+    return haskey(t, I) ? t[I] : zerovector!(allocblock!(t, I, allocator))
+end
+function allocblock!(t::SparseBlockTensorMap, I::CartesianIndex, allocator)
+    V = eachspace(t)[I]
+    return t[I] = if t.istemp
+        TO.tensoralloc(eltype(t), V, Val(true), allocator)
+    else
+        TO.tensoralloc(eltype(t), V, Val(false), allocator)
+    end
 end
 
 # Show
diff --git a/src/tensors/tensoroperations.jl b/src/tensors/tensoroperations.jl
index 98d9c28..32d99e9 100644
--- a/src/tensors/tensoroperations.jl
+++ b/src/tensors/tensoroperations.jl
@@ -101,13 +101,20 @@ function TO.tensoralloc(
         ::Type{BT}, structure::TensorMapSumSpace, istemp::Val, allocator = TO.DefaultAllocator()
     ) where {BT <: AbstractBlockTensorMap}
     C = BT(undef_blocks, structure)
-    issparse(C) && return C # don't fill up sparse blocks
     blockallocator(V) = TO.tensoralloc(eltype(C), V, istemp, allocator)
     map!(blockallocator, parent(C), eachspace(C))
     return C
 end
 
-# sparse results allocate exactly the blocks that will be written, so `tensorfree!` is consistent
+# sparse tensors start out empty, and only remember whether their blocks are temporaries
+function TO.tensoralloc(
+        ::Type{BT}, structure::TensorMapSumSpace, ::Val{istemp},
+        allocator = TO.DefaultAllocator()
+    ) where {BT <: SparseBlockTensorMap, istemp}
+    return BT(undef_blocks, structure; istemp)
+end
+
+# sparse results allocate the blocks that will be written up front
 const BlockOrAdjoint = Union{AbstractBlockTensorMap, AdjointBlockTensorMap}
 
 _tensoralloc(ttype, structure, keys, istemp::Val, allocator) =
@@ -115,11 +122,8 @@ _tensoralloc(ttype, structure, keys, istemp::Val, allocator) =
 function _tensoralloc(
         ttype::Type{<:SparseBlockTensorMap}, structure, keys, istemp::Val, allocator
     )
-    C = ttype(undef_blocks, structure)
-    Vs = eachspace(C)
-    for I in keys
-        haskey(C, I) || (C[I] = TO.tensoralloc(eltype(C), Vs[I], istemp, allocator))
-    end
+    C = TO.tensoralloc(ttype, structure, istemp, allocator)
+    foreach(I -> haskey(C, I) || allocblock!(C, I, allocator), keys)
     return C
 end
 
@@ -201,7 +205,7 @@ function TO.tensorfree!(t::BlockTensorMap, allocator = TO.DefaultAllocator())
     return nothing
 end
 function TO.tensorfree!(t::SparseBlockTensorMap, allocator = TO.DefaultAllocator())
-    foreach(Base.Fix2(TO.tensorfree!, allocator), nonzero_values(t))
+    t.istemp && foreach(Base.Fix2(TO.tensorfree!, allocator), nonzero_values(t))
     return nothing
 end
 
diff --git a/src/tensors/vectorinterface.jl b/src/tensors/vectorinterface.jl
index f7fbe66..fdbfd26 100644
--- a/src/tensors/vectorinterface.jl
+++ b/src/tensors/vectorinterface.jl
@@ -1,7 +1,11 @@
 # zerovector
 # ----------
 VI.zerovector!(t::AbstractBlockTensorMap) = (zerovector!(parent(t)); t)
-VI.zerovector!(t::SparseBlockTensorMap) = (empty!(t.data); t)
+# blocks of temporaries belong to an allocator, and are kept until `tensorfree!`
+function VI.zerovector!(t::SparseBlockTensorMap)
+    t.istemp ? foreach(zerovector!, nonzero_values(t)) : empty!(t.data)
+    return t
+end
 
 # scale
 # -----
diff --git a/test/linalg/tensoroperations.jl b/test/linalg/tensoroperations.jl
index 502d4ef..7c0b02c 100644
--- a/test/linalg/tensoroperations.jl
+++ b/test/linalg/tensoroperations.jl
@@ -154,17 +154,15 @@ end
 ##
 
 struct TrackingAllocator
-    live::Set{UInt} # handed out and not freed yet
-    temps::Set{UInt} # handed out as temporaries
-    foreign::Base.RefValue{Int} # freed arrays that were not handed out by the allocator
+    live::Set{UInt} # temporaries that were handed out and not freed yet
+    foreign::Base.RefValue{Int} # freed arrays that were not handed out as temporaries
 end
-TrackingAllocator() = TrackingAllocator(Set{UInt}(), Set{UInt}(), Ref(0))
+TrackingAllocator() = TrackingAllocator(Set{UInt}(), Ref(0))
 function TensorOperations.tensoralloc(
         ::Type{A}, structure, ::Val{istemp}, allocator::TrackingAllocator
     ) where {A <: AbstractArray, istemp}
     x = fill!(TensorOperations.tensoralloc(A, structure, Val(false)), NaN)
-    push!(allocator.live, objectid(x))
-    istemp && push!(allocator.temps, objectid(x))
+    istemp && push!(allocator.live, objectid(x))
     return x
 end
 function TensorOperations.tensorfree!(x::AbstractArray, allocator::TrackingAllocator)
@@ -179,7 +177,7 @@ function check_allocator(f)
     allocator = TrackingAllocator()
     @test f(allocator) ≈ f(TensorOperations.DefaultAllocator())
     @test allocator.foreign[] == 0
-    return @test isdisjoint(allocator.live, allocator.temps)
+    return @test isempty(allocator.live)
 end
 
 @testset "allocator only frees what it allocated ($(sectortype(a)))" for (a, b, c) in (

From 5b0a55927e9b4bab894913b038b8f841a841c230 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Fri, 9 Oct 2026 10:17:41 -0400
Subject: [PATCH 4/5] Store allocated blocks as the type the allocator hands
 out

`TO.tensoralloc` for block tensors kept the block type of the inputs, so blocks
from allocators with a different storage type (e.g. `PtrArray` temporaries of
`ManualAllocator`) were silently converted on insertion: the copy was passed to
`tensorfree!`, and the original was never released. Infer the block type from
the allocator instead, for both dense and sparse block tensors.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/tensors/tensoroperations.jl | 24 +++++++++++++++++-------
 test/linalg/tensoroperations.jl | 10 ++++++++--
 2 files changed, 25 insertions(+), 9 deletions(-)

diff --git a/src/tensors/tensoroperations.jl b/src/tensors/tensoroperations.jl
index 32d99e9..d0a4213 100644
--- a/src/tensors/tensoroperations.jl
+++ b/src/tensors/tensoroperations.jl
@@ -97,21 +97,31 @@ function similarblocktype(::Type{A}, ::Type{TT}) where {A, TT}
     return Core.Compiler.return_type(similar, Tuple{A, Type{TT}, NTuple{numind(TT), Int}})
 end
 
+# blocks are stored as the type the allocator hands out, e.g. `PtrArray`-backed temporaries
+function allocated_blocktype(::Type{TT}, structure, istemp::Val, allocator) where {TT}
+    V = eltype(SumSpaceIndices(structure))
+    TT′ = Base.promote_op(TO.tensoralloc, Type{TT}, V, typeof(istemp), typeof(allocator))
+    TA = AbstractTensorMap{scalartype(TT), spacetype(TT), numout(TT), numin(TT)}
+    return typeintersect(TT′, TA)
+end
+
 function TO.tensoralloc(
         ::Type{BT}, structure::TensorMapSumSpace, istemp::Val, allocator = TO.DefaultAllocator()
     ) where {BT <: AbstractBlockTensorMap}
-    C = BT(undef_blocks, structure)
-    blockallocator(V) = TO.tensoralloc(eltype(C), V, istemp, allocator)
-    map!(blockallocator, parent(C), eachspace(C))
+    TT = eltype(BT)
+    C = BlockTensorMap{allocated_blocktype(TT, structure, istemp, allocator)}(
+        undef_blocks, structure
+    )
+    map!(V -> TO.tensoralloc(TT, V, istemp, allocator), parent(C), eachspace(C))
     return C
 end
 
 # sparse tensors start out empty, and only remember whether their blocks are temporaries
 function TO.tensoralloc(
-        ::Type{BT}, structure::TensorMapSumSpace, ::Val{istemp},
-        allocator = TO.DefaultAllocator()
-    ) where {BT <: SparseBlockTensorMap, istemp}
-    return BT(undef_blocks, structure; istemp)
+        ::Type{BT}, structure::TensorMapSumSpace, istemp::Val, allocator = TO.DefaultAllocator()
+    ) where {BT <: SparseBlockTensorMap}
+    TT = allocated_blocktype(eltype(BT), structure, istemp, allocator)
+    return SparseBlockTensorMap{TT}(undef_blocks, structure; istemp = istemp === Val(true))
 end
 
 # sparse results allocate the blocks that will be written up front
diff --git a/test/linalg/tensoroperations.jl b/test/linalg/tensoroperations.jl
index 7c0b02c..fcfb53d 100644
--- a/test/linalg/tensoroperations.jl
+++ b/test/linalg/tensoroperations.jl
@@ -153,6 +153,7 @@ end
 
 ##
 
+# temporaries are `PtrArray`s, such that converting them into other storage is detected
 struct TrackingAllocator
     live::Set{UInt} # temporaries that were handed out and not freed yet
     foreign::Base.RefValue{Int} # freed arrays that were not handed out as temporaries
@@ -161,13 +162,18 @@ TrackingAllocator() = TrackingAllocator(Set{UInt}(), Ref(0))
 function TensorOperations.tensoralloc(
         ::Type{A}, structure, ::Val{istemp}, allocator::TrackingAllocator
     ) where {A <: AbstractArray, istemp}
-    x = fill!(TensorOperations.tensoralloc(A, structure, Val(false)), NaN)
-    istemp && push!(allocator.live, objectid(x))
+    istemp || return TensorOperations.tensoralloc(A, structure, Val(false))
+    x = fill!(
+        TensorOperations.tensoralloc(A, structure, Val(true), TensorOperations.ManualAllocator()),
+        NaN
+    )
+    push!(allocator.live, objectid(x))
     return x
 end
 function TensorOperations.tensorfree!(x::AbstractArray, allocator::TrackingAllocator)
     if objectid(x) in allocator.live
         delete!(allocator.live, objectid(x))
+        TensorOperations.tensorfree!(x, TensorOperations.ManualAllocator())
     else
         allocator.foreign[] += 1
     end

From 5875e322a894356b07dd42f893d89bb0b64dcf15 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Fri, 9 Oct 2026 14:07:26 -0400
Subject: [PATCH 5/5] Let sparse temporaries keep their allocator

Replace the `istemp` flag by the allocator itself, so that kernels which do not
receive an allocator (`mul!` from TensorKit's `blas_contract!`/`planarcontract!`,
`trace_permute!` from `tensortrace!`) can still obtain missing blocks of a
temporary through it. Allocating the blocks up front in `tensoralloc_add` and
`tensoralloc_contract` is now an optimization rather than required for correctness.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/linalg/linalg.jl             |  2 ++
 src/tensors/blocktensor.jl       |  1 +
 src/tensors/sparseblocktensor.jl | 26 +++++++++++++++++---------
 src/tensors/tensoroperations.jl  | 10 ++++++----
 src/tensors/vectorinterface.jl   |  2 +-
 test/linalg/tensoroperations.jl  | 12 +++++++++++-
 6 files changed, 38 insertions(+), 15 deletions(-)

diff --git a/src/linalg/linalg.jl b/src/linalg/linalg.jl
index 0abd9a3..15ec1cf 100644
--- a/src/linalg/linalg.jl
+++ b/src/linalg/linalg.jl
@@ -118,6 +118,8 @@ end
     )
     return if haskey(C, I)
         C[I] = _mul!!(C[I], A, B, α, β)
+    elseif istemp(C)
+        C[I] = _mul!!(getindex!(C, I), A, B, α, β)
     else
         C[I] = _mul!!(nothing, A, B, α, β)
     end
diff --git a/src/tensors/blocktensor.jl b/src/tensors/blocktensor.jl
index 3d2b6d9..244f162 100644
--- a/src/tensors/blocktensor.jl
+++ b/src/tensors/blocktensor.jl
@@ -219,3 +219,4 @@ end
 # -------
 Base.haskey(t::BlockTensorMap, I::CartesianIndex) = checkbounds(Bool, t.data, I)
 Base.haskey(t::BlockTensorMap, i::Int) = checkbounds(Bool, t.data, i)
+@propagate_inbounds getindex!(t::BlockTensorMap, I::CartesianIndex{N}) where {N} = t[I]
diff --git a/src/tensors/sparseblocktensor.jl b/src/tensors/sparseblocktensor.jl
index cc903bd..393c642 100644
--- a/src/tensors/sparseblocktensor.jl
+++ b/src/tensors/sparseblocktensor.jl
@@ -2,22 +2,22 @@
     struct SparseBlockTensorMap{TT<:AbstractTensorMap{E,S,N₁,N₂}} <: AbstractBlockTensorMap{E,S,N₁,N₂}
 
 Sparse `SparseBlockTensorMap` type that stores tensors of type `TT` in a sparse dictionary.
-`istemp` marks temporaries of TensorOperations, whose blocks are obtained from and released
-to an allocator.
+Temporaries of TensorOperations keep their `allocator` (otherwise `nothing`), from which their
+blocks are obtained and to which they are released.
 """
 struct SparseBlockTensorMap{TT <: AbstractTensorMap, E, S, N₁, N₂, N} <:
     AbstractBlockTensorMap{E, S, N₁, N₂}
     data::Dict{CartesianIndex{N}, TT}
     space::TensorMapSumSpace{S, N₁, N₂}
-    istemp::Bool
+    allocator::Any
 
     # uninitialized constructor
     function SparseBlockTensorMap{TT, E, S, N₁, N₂, N}(
-            ::UndefBlocksInitializer, space::TensorMapSumSpace{S, N₁, N₂}; istemp::Bool = false
+            ::UndefBlocksInitializer, space::TensorMapSumSpace{S, N₁, N₂}; allocator = nothing
         ) where {E, S, N₁, N₂, N, TT <: AbstractTensorMap{E, S, N₁, N₂}}
         @assert N₁ + N₂ == N "SparseBlockTensorMap: data has wrong number of dimensions"
         data = Dict{CartesianIndex{N}, TT}()
-        return new{TT, E, S, N₁, N₂, N}(data, space, istemp)
+        return new{TT, E, S, N₁, N₂, N}(data, space, allocator)
     end
 
     # constructor from data
@@ -25,7 +25,7 @@ struct SparseBlockTensorMap{TT <: AbstractTensorMap, E, S, N₁, N₂, N} <:
             data::Dict{CartesianIndex{N}, TT}, space::TensorMapSumSpace{S, N₁, N₂}
         ) where {E, S, N₁, N₂, N, TT <: AbstractTensorMap{E, S, N₁, N₂}}
         @assert N₁ + N₂ == N "SparseBlockTensorMap: data has wrong number of dimensions"
-        return new{TT, E, S, N₁, N₂, N}(data, space, false)
+        return new{TT, E, S, N₁, N₂, N}(data, space, nothing)
     end
 end
 
@@ -193,7 +193,15 @@ function Base.delete!(t::SparseBlockTensorMap{TT}, I::Vararg{Int, N}) where {TT,
     return delete!(t, CartesianIndex(I...))
 end
 
-# missing blocks of a destination are obtained from the allocator of the operation
+# missing blocks of temporaries are obtained from their own allocator, those of other
+# destinations from the allocator of the operation, if any
+istemp(t::AbstractBlockTensorMap) = false
+istemp(t::SparseBlockTensorMap) = !isnothing(t.allocator)
+
+@propagate_inbounds function getindex!(t::SparseBlockTensorMap, I::CartesianIndex{N}) where {N}
+    istemp(t) || return getindex!(parent(t), I)
+    return haskey(t, I) ? t[I] : zerovector!(allocblock!(t, I, t.allocator))
+end
 @propagate_inbounds getindex!(
     t::SparseBlockTensorMap, I::CartesianIndex, ::TO.DefaultAllocator
 ) = getindex!(t, I)
@@ -202,8 +210,8 @@ end
 end
 function allocblock!(t::SparseBlockTensorMap, I::CartesianIndex, allocator)
     V = eachspace(t)[I]
-    return t[I] = if t.istemp
-        TO.tensoralloc(eltype(t), V, Val(true), allocator)
+    return t[I] = if istemp(t)
+        TO.tensoralloc(eltype(t), V, Val(true), t.allocator)
     else
         TO.tensoralloc(eltype(t), V, Val(false), allocator)
     end
diff --git a/src/tensors/tensoroperations.jl b/src/tensors/tensoroperations.jl
index d0a4213..e6d14bd 100644
--- a/src/tensors/tensoroperations.jl
+++ b/src/tensors/tensoroperations.jl
@@ -116,12 +116,14 @@ function TO.tensoralloc(
     return C
 end
 
-# sparse tensors start out empty, and only remember whether their blocks are temporaries
+# sparse tensors start out empty, and temporaries remember their allocator
 function TO.tensoralloc(
         ::Type{BT}, structure::TensorMapSumSpace, istemp::Val, allocator = TO.DefaultAllocator()
     ) where {BT <: SparseBlockTensorMap}
     TT = allocated_blocktype(eltype(BT), structure, istemp, allocator)
-    return SparseBlockTensorMap{TT}(undef_blocks, structure; istemp = istemp === Val(true))
+    return SparseBlockTensorMap{TT}(
+        undef_blocks, structure; allocator = istemp === Val(true) ? allocator : nothing
+    )
 end
 
 # sparse results allocate the blocks that will be written up front
@@ -215,7 +217,7 @@ function TO.tensorfree!(t::BlockTensorMap, allocator = TO.DefaultAllocator())
     return nothing
 end
 function TO.tensorfree!(t::SparseBlockTensorMap, allocator = TO.DefaultAllocator())
-    t.istemp && foreach(Base.Fix2(TO.tensorfree!, allocator), nonzero_values(t))
+    istemp(t) && foreach(Base.Fix2(TO.tensorfree!, allocator), nonzero_values(t))
     return nothing
 end
 
@@ -253,7 +255,7 @@ function TK.trace_permute!(
         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
+            getindex!(tdst, Idst), vsrc, (p₁, p₂), (q₁, q₂), α, One(), backend
         )
     end
     return tdst
diff --git a/src/tensors/vectorinterface.jl b/src/tensors/vectorinterface.jl
index fdbfd26..32bf017 100644
--- a/src/tensors/vectorinterface.jl
+++ b/src/tensors/vectorinterface.jl
@@ -3,7 +3,7 @@
 VI.zerovector!(t::AbstractBlockTensorMap) = (zerovector!(parent(t)); t)
 # blocks of temporaries belong to an allocator, and are kept until `tensorfree!`
 function VI.zerovector!(t::SparseBlockTensorMap)
-    t.istemp ? foreach(zerovector!, nonzero_values(t)) : empty!(t.data)
+    istemp(t) ? foreach(zerovector!, nonzero_values(t)) : empty!(t.data)
     return t
 end
 
diff --git a/test/linalg/tensoroperations.jl b/test/linalg/tensoroperations.jl
index fcfb53d..4e972ca 100644
--- a/test/linalg/tensoroperations.jl
+++ b/test/linalg/tensoroperations.jl
@@ -3,7 +3,7 @@ using BlockTensorKit
 using TensorKit
 using TensorOperations
 using Random
-using LinearAlgebra: norm
+using LinearAlgebra: norm, mul!
 
 ##
 Vtr = (
@@ -212,6 +212,16 @@ end
     check_allocator() do al
         return @planar allocator = al R[x; q] := A[x z; w] * B[w z; y] * D[y; q]
     end
+    check_allocator() do al # kernels without allocator argument fill an empty temporary
+        T = TensorOperations.tensoralloc(typeof(D), a ← a, Val(true), al)
+        mul!(T, D, D')
+        TensorOperations.tensortrace!(
+            T, X, ((1,), (3,)), ((2,), (4,)), false, 1, 1, TensorOperations.DefaultBackend(), al
+        )
+        n = norm(T)
+        TensorOperations.tensorfree!(T, al)
+        return n
+    end
 
     # permuted and BLAS-compatible destinations, without and with blocks outside of the product
     @tensor C1[x y; p q] := A[x q; w] * B[w p; y]