From f6aef5fa2d63a068278a6eb8fc4b37390103fe53 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Tue, 10 Mar 2026 14:31:41 -0400
Subject: [PATCH 01/19] Another attempt to fix matricize of GradedArray

---
 Project.toml         |  2 +-
 src/tensoralgebra.jl | 46 +++++++++++++++++++++++++++++++++++++++++---
 2 files changed, 44 insertions(+), 4 deletions(-)

diff --git a/Project.toml b/Project.toml
index d2234b0d..49b74ce3 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "GradedArrays"
 uuid = "bc96ca6e-b7c8-4bb6-888e-c93f838762c2"
-version = "0.6.19"
+version = "0.6.20"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
 
 [workspace]
diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl
index 5ec64a41..da5aa6e5 100644
--- a/src/tensoralgebra.jl
+++ b/src/tensoralgebra.jl
@@ -136,9 +136,49 @@ end
 
 # Sort the blocks by sector and then merge the common sectors.
 function sectormergesort(a::AbstractArray)
-    # TODO: fix this, no clue why broken and no clue how to fix
-    return a
-
     I = sectormergesortperm.(axes(a))
     return a[I...]
 end
+
+using BlockArrays: AbstractBlockVector, Block
+function Base.getindex(
+        a::GradedArray{<:Any, N},
+        I::Vararg{AbstractBlockVector{<:Block{1}}, N}
+    ) where {N}
+    axes_dest = ntuple(d -> only(axes(axes(a, d)[I[d]])), N)
+    a_dest = BlockSparseArray{eltype(a)}(undef, axes_dest)
+
+    grouped_blocks = ntuple(N) do d
+        return map(blocks(I[d])) do bI
+            return map(k -> I[d][k], only(bI.indices))
+        end
+    end
+    grouped_ranges = ntuple(N) do d
+        return map(grouped_blocks[d]) do bs
+            lengths = map(b -> length(axes(a, d)[b]), bs)
+            starts = cumsum(vcat(1, lengths[1:(end - 1)]))
+            return map((s, n) -> s:(s + n - 1), starts, lengths)
+        end
+    end
+
+    for I_dest in CartesianIndices(map(gs -> Base.OneTo(length(gs)), grouped_blocks))
+        src_blocks = ntuple(d -> grouped_blocks[d][I_dest[d]], N)
+        src_ranges = ntuple(d -> grouped_ranges[d][I_dest[d]], N)
+        block_dims = ntuple(d -> sum(length.(src_ranges[d])), N)
+        block_data = zeros(eltype(a), block_dims)
+        block_stored = false
+
+        src_group_indices = Iterators.product(map(Base.OneTo ∘ length, src_blocks)...)
+        for I_src in src_group_indices
+            bI_src = Block(ntuple(d -> src_blocks[d][I_src[d]], N))
+            !isstored(a, bI_src) && continue
+            block_stored = true
+            rI_dest = ntuple(d -> src_ranges[d][I_src[d]], N)
+            block_data[rI_dest...] = a[bI_src]
+        end
+
+        block_stored || continue
+        a_dest[Block(Tuple(I_dest))] = block_data
+    end
+    return a_dest
+end

From 904d0fd1ba6cb05724971fd8385f425cb4f1e051 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Tue, 10 Mar 2026 14:50:30 -0400
Subject: [PATCH 02/19] Refactor

---
 src/tensoralgebra.jl | 77 ++++++++++++++++++++++++++++++--------------
 1 file changed, 53 insertions(+), 24 deletions(-)

diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl
index da5aa6e5..05cf627e 100644
--- a/src/tensoralgebra.jl
+++ b/src/tensoralgebra.jl
@@ -146,39 +146,68 @@ function Base.getindex(
         I::Vararg{AbstractBlockVector{<:Block{1}}, N}
     ) where {N}
     axes_dest = ntuple(d -> only(axes(axes(a, d)[I[d]])), N)
-    a_dest = BlockSparseArray{eltype(a)}(undef, axes_dest)
+    # Preserve `blocktype(a)` (for example GPU blocks) by using `similar`.
+    a_dest = similar(a, eltype(a), axes_dest)
 
+    # Group selected source blocks per destination block index along each dimension.
     grouped_blocks = ntuple(N) do d
-        return map(blocks(I[d])) do bI
-            return map(k -> I[d][k], only(bI.indices))
+        map(blocks(I[d])) do bI
+            return collect(bI)
         end
     end
-    grouped_ranges = ntuple(N) do d
-        return map(grouped_blocks[d]) do bs
-            lengths = map(b -> length(axes(a, d)[b]), bs)
-            starts = cumsum(vcat(1, lengths[1:(end - 1)]))
-            return map((s, n) -> s:(s + n - 1), starts, lengths)
+
+    # Map source block index -> destination block indices containing it.
+    source_to_dest = ntuple(N) do d
+        block_map = Dict{Int, Vector{Int}}()
+        for j in eachindex(grouped_blocks[d]), b in grouped_blocks[d][j]
+            push!(get!(block_map, Int(b), Int[]), j)
         end
+        return block_map
     end
 
-    for I_dest in CartesianIndices(map(gs -> Base.OneTo(length(gs)), grouped_blocks))
-        src_blocks = ntuple(d -> grouped_blocks[d][I_dest[d]], N)
-        src_ranges = ntuple(d -> grouped_ranges[d][I_dest[d]], N)
-        block_dims = ntuple(d -> sum(length.(src_ranges[d])), N)
-        block_data = zeros(eltype(a), block_dims)
-        block_stored = false
-
-        src_group_indices = Iterators.product(map(Base.OneTo ∘ length, src_blocks)...)
-        for I_src in src_group_indices
-            bI_src = Block(ntuple(d -> src_blocks[d][I_src[d]], N))
-            !isstored(a, bI_src) && continue
-            block_stored = true
-            rI_dest = ntuple(d -> src_ranges[d][I_src[d]], N)
-            block_data[rI_dest...] = a[bI_src]
+    # Only process destination blocks that can receive data from stored source blocks.
+    dest_blocks = Set{NTuple{N, Int}}()
+    for bI_src in eachblockstoredindex(a)
+        dest_choices = map(1:N) do d
+            return get(source_to_dest[d], Int(Tuple(bI_src)[d]), Int[])
+        end
+        any(isempty, dest_choices) && continue
+        for bI_dest in Iterators.product(dest_choices...)
+            push!(dest_blocks, ntuple(d -> bI_dest[d], N))
         end
+    end
 
-        block_stored || continue
-        a_dest[Block(Tuple(I_dest))] = block_data
+    for bI_dest in dest_blocks
+        b = Block(bI_dest)
+        bs = ntuple(d -> grouped_blocks[d][bI_dest[d]], N)
+        view(a_dest, b) .= view(a, bs...)
     end
     return a_dest
 end
+
+function Base.copyto!(C::SectorArray, A::SectorArray)
+    axes(C) == axes(A) || throw(DimensionMismatch())
+    copyto!(C.data, A.data)
+    return C
+end
+function Base.copyto!(
+        C::SectorArray,
+        A::SubArray{<:Any, <:Any, <:BlockSparseArrays.AbstractBlockSparseArray}
+    )
+    size(C) == size(A) || throw(DimensionMismatch())
+    copyto!(C.data, Array(A))
+    return C
+end
+function Base.copyto!(C::SectorArray, A::AbstractArray)
+    size(C) == size(A) || throw(DimensionMismatch())
+    copyto!(C.data, A)
+    return C
+end
+
+# Resolve ambiguity for copying from block-sliced GradedArray views into a block view.
+function Base.copyto!(
+        a_dest::BlockSparseArrays.BlockView{<:Any, <:Any, <:GradedArray},
+        a_src::SubArray{<:Any, <:Any, <:GradedArray}
+    )
+    return copyto!(a_dest, Array(a_src))
+end

From 6c780e1d450576c31c6ce203acb0c93b4e495c12 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Tue, 10 Mar 2026 16:50:31 -0400
Subject: [PATCH 03/19] Refactor

---
 src/tensoralgebra.jl | 84 +++++++++++++++++++-------------------------
 1 file changed, 36 insertions(+), 48 deletions(-)

diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl
index 05cf627e..72fbb8fe 100644
--- a/src/tensoralgebra.jl
+++ b/src/tensoralgebra.jl
@@ -145,9 +145,8 @@ function Base.getindex(
         a::GradedArray{<:Any, N},
         I::Vararg{AbstractBlockVector{<:Block{1}}, N}
     ) where {N}
-    axes_dest = ntuple(d -> only(axes(axes(a, d)[I[d]])), N)
-    # Preserve `blocktype(a)` (for example GPU blocks) by using `similar`.
-    a_dest = similar(a, eltype(a), axes_dest)
+    axes_dest = ntuple(d -> only(axes(axes(a, d)[I[d]])), Val(N))
+    a_dest = similar(a, axes_dest)
 
     # Group selected source blocks per destination block index along each dimension.
     grouped_blocks = ntuple(N) do d
@@ -156,58 +155,47 @@ function Base.getindex(
         end
     end
 
-    # Map source block index -> destination block indices containing it.
+    # Map source block index -> (destination block index, destination subrange).
     source_to_dest = ntuple(N) do d
-        block_map = Dict{Int, Vector{Int}}()
-        for j in eachindex(grouped_blocks[d]), b in grouped_blocks[d][j]
-            push!(get!(block_map, Int(b), Int[]), j)
+        block_map = Dict{Int, Tuple{Int, UnitRange{Int}}}()
+        for (j, src_blocks) in pairs(grouped_blocks[d])
+            offset = 1
+            for b in src_blocks
+                b_int = Int(b)
+                haskey(block_map, b_int) &&
+                    throw(
+                    ArgumentError(
+                        "Source block appears in multiple destination groups."
+                    )
+                )
+                len_b = length(axes(a, d)[b])
+                block_map[b_int] = (j, offset:(offset + len_b - 1))
+                offset += len_b
+            end
         end
         return block_map
     end
 
-    # Only process destination blocks that can receive data from stored source blocks.
-    dest_blocks = Set{NTuple{N, Int}}()
+    # Populate destination blocks by placing each stored source block into the
+    # corresponding destination subblock.
     for bI_src in eachblockstoredindex(a)
-        dest_choices = map(1:N) do d
-            return get(source_to_dest[d], Int(Tuple(bI_src)[d]), Int[])
+        bI_dest = Vector{Int}(undef, N)
+        rI_dest = Vector{UnitRange{Int}}(undef, N)
+        valid_dest = true
+        for d in 1:N
+            dst = get(source_to_dest[d], Int(Tuple(bI_src)[d]), nothing)
+            if isnothing(dst)
+                valid_dest = false
+                break
+            end
+            j, r = dst
+            bI_dest[d] = j
+            rI_dest[d] = r
         end
-        any(isempty, dest_choices) && continue
-        for bI_dest in Iterators.product(dest_choices...)
-            push!(dest_blocks, ntuple(d -> bI_dest[d], N))
-        end
-    end
-
-    for bI_dest in dest_blocks
-        b = Block(bI_dest)
-        bs = ntuple(d -> grouped_blocks[d][bI_dest[d]], N)
-        view(a_dest, b) .= view(a, bs...)
+        valid_dest || continue
+        b_dest = Block(Tuple(bI_dest))
+        a_dest_b = @view!(a_dest[b_dest])
+        copyto!(@view(a_dest_b[Tuple(rI_dest)...]), a[bI_src])
     end
     return a_dest
 end
-
-function Base.copyto!(C::SectorArray, A::SectorArray)
-    axes(C) == axes(A) || throw(DimensionMismatch())
-    copyto!(C.data, A.data)
-    return C
-end
-function Base.copyto!(
-        C::SectorArray,
-        A::SubArray{<:Any, <:Any, <:BlockSparseArrays.AbstractBlockSparseArray}
-    )
-    size(C) == size(A) || throw(DimensionMismatch())
-    copyto!(C.data, Array(A))
-    return C
-end
-function Base.copyto!(C::SectorArray, A::AbstractArray)
-    size(C) == size(A) || throw(DimensionMismatch())
-    copyto!(C.data, A)
-    return C
-end
-
-# Resolve ambiguity for copying from block-sliced GradedArray views into a block view.
-function Base.copyto!(
-        a_dest::BlockSparseArrays.BlockView{<:Any, <:Any, <:GradedArray},
-        a_src::SubArray{<:Any, <:Any, <:GradedArray}
-    )
-    return copyto!(a_dest, Array(a_src))
-end

From 8f7b58dae433d5d16952b5263fe17c3db6c7423f Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Tue, 10 Mar 2026 17:09:27 -0400
Subject: [PATCH 04/19] Simplify

---
 src/tensoralgebra.jl | 60 +++++++++++---------------------------------
 1 file changed, 15 insertions(+), 45 deletions(-)

diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl
index 72fbb8fe..f0e3afe6 100644
--- a/src/tensoralgebra.jl
+++ b/src/tensoralgebra.jl
@@ -147,55 +147,25 @@ function Base.getindex(
     ) where {N}
     axes_dest = ntuple(d -> only(axes(axes(a, d)[I[d]])), Val(N))
     a_dest = similar(a, axes_dest)
-
-    # Group selected source blocks per destination block index along each dimension.
-    grouped_blocks = ntuple(N) do d
-        map(blocks(I[d])) do bI
-            return collect(bI)
-        end
-    end
-
-    # Map source block index -> (destination block index, destination subrange).
-    source_to_dest = ntuple(N) do d
-        block_map = Dict{Int, Tuple{Int, UnitRange{Int}}}()
-        for (j, src_blocks) in pairs(grouped_blocks[d])
-            offset = 1
-            for b in src_blocks
-                b_int = Int(b)
-                haskey(block_map, b_int) &&
-                    throw(
-                    ArgumentError(
-                        "Source block appears in multiple destination groups."
-                    )
-                )
-                len_b = length(axes(a, d)[b])
-                block_map[b_int] = (j, offset:(offset + len_b - 1))
-                offset += len_b
+    ax = axes(a)
+    # Map source block index -> (dest block index, subrange); 0 means not selected.
+    src_to_dest = ntuple(N) do d
+        result = fill((0, 0:0), length(blocks(ax[d])))
+        for j in 1:length(blocks(I[d]))
+            grp = I[d][Block(j)]
+            grp_start = first(ax[d][first(grp)])
+            for b in grp
+                result[Int(b)] = (j, ax[d][b] .- (grp_start - 1))
             end
         end
-        return block_map
+        return result
     end
-
-    # Populate destination blocks by placing each stored source block into the
-    # corresponding destination subblock.
     for bI_src in eachblockstoredindex(a)
-        bI_dest = Vector{Int}(undef, N)
-        rI_dest = Vector{UnitRange{Int}}(undef, N)
-        valid_dest = true
-        for d in 1:N
-            dst = get(source_to_dest[d], Int(Tuple(bI_src)[d]), nothing)
-            if isnothing(dst)
-                valid_dest = false
-                break
-            end
-            j, r = dst
-            bI_dest[d] = j
-            rI_dest[d] = r
-        end
-        valid_dest || continue
-        b_dest = Block(Tuple(bI_dest))
-        a_dest_b = @view!(a_dest[b_dest])
-        copyto!(@view(a_dest_b[Tuple(rI_dest)...]), a[bI_src])
+        src_tuple = Tuple(bI_src)
+        info = ntuple(d -> src_to_dest[d][Int(src_tuple[d])], Val(N))
+        any(iszero ∘ first, info) && continue
+        a_dest_b = @view!(a_dest[Block(map(first, info))])
+        copyto!(@view(a_dest_b[map(last, info)...]), a[bI_src])
     end
     return a_dest
 end

From c4e30ef964f18d3d03a677bf161cd87dfa93e308 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Tue, 10 Mar 2026 17:25:19 -0400
Subject: [PATCH 05/19] More simplification

---
 src/tensoralgebra.jl | 26 ++++++++++++++------------
 1 file changed, 14 insertions(+), 12 deletions(-)

diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl
index f0e3afe6..70f824ed 100644
--- a/src/tensoralgebra.jl
+++ b/src/tensoralgebra.jl
@@ -140,32 +140,34 @@ function sectormergesort(a::AbstractArray)
     return a[I...]
 end
 
-using BlockArrays: AbstractBlockVector, Block
+using BlockArrays: AbstractBlockVector, Block, BlockIndexRange
 function Base.getindex(
         a::GradedArray{<:Any, N},
         I::Vararg{AbstractBlockVector{<:Block{1}}, N}
     ) where {N}
-    axes_dest = ntuple(d -> only(axes(axes(a, d)[I[d]])), Val(N))
-    a_dest = similar(a, axes_dest)
+    ax_dest = ntuple(d -> only(axes(axes(a, d)[I[d]])), Val(N))
+    a_dest = similar(a, ax_dest)
     ax = axes(a)
-    # Map source block index -> (dest block index, subrange); 0 means not selected.
-    src_to_dest = ntuple(N) do d
-        result = fill((0, 0:0), length(blocks(ax[d])))
+    # Map source Block -> BlockIndexRange encoding dest block + subrange within it
+    src_to_dest = ntuple(Val(N)) do d
+        dict = Dict{Block{1}, BlockIndexRange{1}}()
         for j in 1:length(blocks(I[d]))
             grp = I[d][Block(j)]
             grp_start = first(ax[d][first(grp)])
             for b in grp
-                result[Int(b)] = (j, ax[d][b] .- (grp_start - 1))
+                dict[b] = Block(j)[ax[d][b] .- (grp_start - 1)]
             end
         end
-        return result
+        return dict
     end
     for bI_src in eachblockstoredindex(a)
         src_tuple = Tuple(bI_src)
-        info = ntuple(d -> src_to_dest[d][Int(src_tuple[d])], Val(N))
-        any(iszero ∘ first, info) && continue
-        a_dest_b = @view!(a_dest[Block(map(first, info))])
-        copyto!(@view(a_dest_b[map(last, info)...]), a[bI_src])
+        all(d -> haskey(src_to_dest[d], src_tuple[d]), 1:N) || continue
+        dest_info = ntuple(d -> src_to_dest[d][src_tuple[d]], Val(N))
+        dest_b = Block(map(di -> only(Tuple(di.block)), dest_info))
+        a_dest_b = @view!(a_dest[dest_b])
+        dest_r = map(di -> only(di.indices), dest_info)
+        copyto!(@view(a_dest_b[dest_r...]), a[bI_src])
     end
     return a_dest
 end

From a9488122329bbfe6cc2b11555346072b3e2647f3 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Tue, 10 Mar 2026 18:09:39 -0400
Subject: [PATCH 06/19] Fix offset issue

---
 src/tensoralgebra.jl | 13 +++++++++----
 1 file changed, 9 insertions(+), 4 deletions(-)

diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl
index 70f824ed..3b93e7b3 100644
--- a/src/tensoralgebra.jl
+++ b/src/tensoralgebra.jl
@@ -140,7 +140,7 @@ function sectormergesort(a::AbstractArray)
     return a[I...]
 end
 
-using BlockArrays: AbstractBlockVector, Block, BlockIndexRange
+using BlockArrays: AbstractBlockVector, Block
 function Base.getindex(
         a::GradedArray{<:Any, N},
         I::Vararg{AbstractBlockVector{<:Block{1}}, N}
@@ -150,12 +150,17 @@ function Base.getindex(
     ax = axes(a)
     # Map source Block -> BlockIndexRange encoding dest block + subrange within it
     src_to_dest = ntuple(Val(N)) do d
-        dict = Dict{Block{1}, BlockIndexRange{1}}()
+        key_T = eltype(I[d])
+        range_T = typeof(firstindex(ax[d]):lastindex(ax[d]))
+        bix_T = Base.promote_op(getindex, key_T, range_T)
+        dict = Dict{key_T, bix_T}()
         for j in 1:length(blocks(I[d]))
             grp = I[d][Block(j)]
-            grp_start = first(ax[d][first(grp)])
+            offset = 1
             for b in grp
-                dict[b] = Block(j)[ax[d][b] .- (grp_start - 1)]
+                len = length(ax[d][b])
+                dict[b] = Block(j)[offset:(offset + len - 1)]
+                offset += len
             end
         end
         return dict

From ec31e6d339d7655c960870287656682324f2cd46 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Tue, 10 Mar 2026 19:33:50 -0400
Subject: [PATCH 07/19] Simplify more, fix sector range printing

---
 src/sectorarray.jl   |  7 +++----
 src/tensoralgebra.jl | 21 ++++++++++-----------
 2 files changed, 13 insertions(+), 15 deletions(-)

diff --git a/src/sectorarray.jl b/src/sectorarray.jl
index 8a19832a..6acfa1a8 100644
--- a/src/sectorarray.jl
+++ b/src/sectorarray.jl
@@ -86,10 +86,9 @@ function Base.show(
         io::IO, g::SectorUnitRange{I, RB, R}
     ) where {I <: SectorRange, RB <: AbstractUnitRange{Int}, R <: AbstractUnitRange{Int}}
     a, b = kroneckerfactors(g)
-    if b isa Base.OneTo
-        print(io, "sectorrange(", a, ", ", unproduct(g), ")")
-    else
-        print(io, "sectorrange(", a, " => ", b, ", ", unproduct(g), ")")
+    print(io, "sectorrange(", a, ", ", b, ")")
+    if !isone(first(g))
+        print(io, " .+ ", first(g) - 1)
     end
     return nothing
 end
diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl
index 3b93e7b3..2ea0905f 100644
--- a/src/tensoralgebra.jl
+++ b/src/tensoralgebra.jl
@@ -150,24 +150,23 @@ function Base.getindex(
     ax = axes(a)
     # Map source Block -> BlockIndexRange encoding dest block + subrange within it
     src_to_dest = ntuple(Val(N)) do d
-        key_T = eltype(I[d])
-        range_T = typeof(firstindex(ax[d]):lastindex(ax[d]))
-        bix_T = Base.promote_op(getindex, key_T, range_T)
-        dict = Dict{key_T, bix_T}()
-        for j in 1:length(blocks(I[d]))
-            grp = I[d][Block(j)]
+        key_type = eltype(I[d])
+        range_type = UnitRange{Int}
+        val_type = Base.promote_op(getindex, key_type, range_type)
+        dict = Dict{key_type, val_type}()
+        for j in eachindex(blocks(I[d]))
+            sub_blocks = I[d][Block(j)]
             offset = 1
-            for b in grp
-                len = length(ax[d][b])
-                dict[b] = Block(j)[offset:(offset + len - 1)]
-                offset += len
+            for b in sub_blocks
+                r = Base.OneTo(length(ax[d][b])) .+ (offset - 1)
+                dict[b] = Block(j)[r]
+                offset += length(r)
             end
         end
         return dict
     end
     for bI_src in eachblockstoredindex(a)
         src_tuple = Tuple(bI_src)
-        all(d -> haskey(src_to_dest[d], src_tuple[d]), 1:N) || continue
         dest_info = ntuple(d -> src_to_dest[d][src_tuple[d]], Val(N))
         dest_b = Block(map(di -> only(Tuple(di.block)), dest_info))
         a_dest_b = @view!(a_dest[dest_b])

From 91877732bbe39f9c011667fc138ed1ad5f4fd361 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Tue, 10 Mar 2026 19:38:38 -0400
Subject: [PATCH 08/19] Fix print tests

---
 test/test_show.jl | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/test/test_show.jl b/test/test_show.jl
index 254c48fb..1830745d 100644
--- a/test/test_show.jl
+++ b/test/test_show.jl
@@ -38,17 +38,17 @@ end
     @test sprint(show, g1) == "GradedUnitRange[$x => 2, $y => 3, $z => 2]"
     @test sprint(show, MIME("text/plain"), g1) ==
         "GradedUnitRange{$U1}\n" *
-        "sectorrange($x, 1:2)\n" *
-        "sectorrange($y, 3:5)\n" *
-        "sectorrange($z, 6:7)"
+        "sectorrange($x, Base.OneTo(2))\n" *
+        "sectorrange($y, Base.OneTo(3)) .+ 2\n" *
+        "sectorrange($z, Base.OneTo(2)) .+ 5"
 
     g1d = dual(g1)
     @test sprint(show, g1d) == "GradedUnitRange[$x' => 2, $y' => 3, $z' => 2]"
     @test sprint(show, MIME("text/plain"), g1d) ==
         "GradedUnitRange{$U1}\n" *
-        "sectorrange($x', 1:2)\n" *
-        "sectorrange($y', 3:5)\n" *
-        "sectorrange($z', 6:7)"
+        "sectorrange($x', Base.OneTo(2))\n" *
+        "sectorrange($y', Base.OneTo(3)) .+ 2\n" *
+        "sectorrange($z', Base.OneTo(2)) .+ 5"
 end
 
 @testset "show GradedArray" begin

From 0de6f93b7fdf3609a52db09cf96dcd324aeb3693 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Wed, 11 Mar 2026 08:43:14 -0400
Subject: [PATCH 09/19] Rename offset to start

---
 src/tensoralgebra.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl
index 2ea0905f..134223e3 100644
--- a/src/tensoralgebra.jl
+++ b/src/tensoralgebra.jl
@@ -156,11 +156,11 @@ function Base.getindex(
         dict = Dict{key_type, val_type}()
         for j in eachindex(blocks(I[d]))
             sub_blocks = I[d][Block(j)]
-            offset = 1
+            start = 1
             for b in sub_blocks
-                r = Base.OneTo(length(ax[d][b])) .+ (offset - 1)
+                r = Base.OneTo(length(ax[d][b])) .+ (start - 1)
                 dict[b] = Block(j)[r]
-                offset += length(r)
+                start += length(r)
             end
         end
         return dict

From b63f9b515b00a155f404dfed74b986203586fea2 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Wed, 11 Mar 2026 09:16:47 -0400
Subject: [PATCH 10/19] Start fixing unmatricize

---
 src/tensoralgebra.jl | 92 ++++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 88 insertions(+), 4 deletions(-)

diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl
index 134223e3..e1f1046a 100644
--- a/src/tensoralgebra.jl
+++ b/src/tensoralgebra.jl
@@ -1,7 +1,8 @@
-using BlockArrays: blocks, eachblockaxes1
+using BlockArrays: BlockIndexRange, blocks, eachblockaxes1
 using BlockSparseArrays: BlockSparseArray, blockrange, blockreshape
-using GradedArrays: GradedArray, GradedUnitRange, SectorRange, flip, invblockperm,
-    sectormergesortperm, sectorsortperm, trivial, unmerged_tensor_product, ×
+using GradedArrays: GradedArray, GradedUnitRange, SectorRange, flip, gradedrange,
+    invblockperm, sectormergesortperm, sectors, sectorsortperm, trivial,
+    unmerged_tensor_product, ×
 using TensorAlgebra: TensorAlgebra, AbstractBlockPermutation, BlockedTuple, FusionStyle,
     ReshapeFusion, matricize, matricize_axes, tensor_product_axis, trivialbiperm,
     tuplemortar, unmatricize
@@ -140,7 +141,90 @@ function sectormergesort(a::AbstractArray)
     return a[I...]
 end
 
-using BlockArrays: AbstractBlockVector, Block
+using BlockArrays: AbstractBlockVector, Block, BlockVector
+
+# Splitting: each I[d][k] = Block(b)[r] means dest block k comes from source block b
+# at subrange r. This is the inverse of the merging getindex below.
+function Base.getindex(
+        a::GradedArray{<:Any, N},
+        I::Vararg{AbstractVector{<:BlockIndexRange{1}}, N}
+    ) where {N}
+    ax = axes(a)
+    ax_dest = ntuple(Val(N)) do d
+        return gradedrange(
+            [
+                sectors(ax[d])[only(Tuple(I[d][k].block))] => length(only(I[d][k].indices))
+                    for k in eachindex(I[d])
+            ]
+        )
+    end
+    a_dest = similar(a, ax_dest)
+    # Map source block b → list of (dest block k, src subrange r, dest subrange 1:length(r))
+    src_to_dests = ntuple(Val(N)) do d
+        dict = Dict{Block{1}, Vector{Tuple{Int, UnitRange{Int}, Base.OneTo{Int}}}}()
+        for k in eachindex(I[d])
+            bir = I[d][k]
+            b = Block(only(Tuple(bir.block)))
+            r = only(bir.indices)
+            push!(
+                get!(dict, b, Tuple{Int, UnitRange{Int}, Base.OneTo{Int}}[]),
+                (k, r, Base.OneTo(length(r)))
+            )
+        end
+        return dict
+    end
+    for bI_src in eachblockstoredindex(a)
+        src_tuple = Tuple(bI_src)
+        all(d -> haskey(src_to_dests[d], src_tuple[d]), 1:N) || continue
+        dest_refs = ntuple(d -> src_to_dests[d][src_tuple[d]], Val(N))
+        for combo in Iterators.product(dest_refs...)
+            dest_b = Block(ntuple(d -> combo[d][1], Val(N)))
+            a_dest_b = @view!(a_dest[dest_b])
+            src_r = ntuple(d -> combo[d][2], Val(N))
+            dest_r = ntuple(d -> combo[d][3], Val(N))
+            copyto!(@view(a_dest_b[dest_r...]), @view(a[bI_src][src_r...]))
+        end
+    end
+    return a_dest
+end
+
+# GradedUnitRange index: compute BlockIndexRange per block by linear scan,
+# then delegate to the splitting getindex above.
+# Assumes blocks of each I[d] subdivide blocks of the corresponding axis of a
+# (as is the case in unmatricize, where I[d] is derived from the unmerged fused axis).
+function Base.getindex(
+        a::GradedArray{<:Any, N},
+        I::Vararg{GradedUnitRange, N}
+    ) where {N}
+    J = map(axes(a), I) do src_ax, tgt_ax
+        n = length(sectors(tgt_ax))
+        J_d = Vector{BlockIndexRange{1}}(undef, n)
+        j = 1
+        offset = 0
+        for k in 1:n
+            size_k = length(tgt_ax[Block(k)])
+            J_d[k] = Block(j)[(offset + 1):(offset + size_k)]
+            offset += size_k
+            if offset == length(src_ax[Block(j)])
+                j += 1
+                offset = 0
+            end
+        end
+        return J_d
+    end
+    return a[J...]
+end
+
+# Vector{Block{1}} index: block permutation with no merging.
+# Wraps as a singleton-group BlockVector and delegates to the merging getindex below.
+function Base.getindex(
+        a::GradedArray{<:Any, N},
+        I::Vararg{Vector{<:Block{1}}, N}
+    ) where {N}
+    return a[map(v -> BlockVector(v, fill(1, length(v))), I)...]
+end
+
+# Merging: each I[d] groups source blocks into destination blocks.
 function Base.getindex(
         a::GradedArray{<:Any, N},
         I::Vararg{AbstractBlockVector{<:Block{1}}, N}

From afa4563faaf8653b5f28ce0399758a7405a37b3a Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Wed, 11 Mar 2026 09:57:29 -0400
Subject: [PATCH 11/19] Fix

---
 src/tensoralgebra.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl
index e1f1046a..9c2fe44e 100644
--- a/src/tensoralgebra.jl
+++ b/src/tensoralgebra.jl
@@ -153,7 +153,7 @@ function Base.getindex(
     ax_dest = ntuple(Val(N)) do d
         return gradedrange(
             [
-                sectors(ax[d])[only(Tuple(I[d][k].block))] => length(only(I[d][k].indices))
+                sectors(ax[d])[Int(I[d][k].block)] => length(only(I[d][k].indices))
                     for k in eachindex(I[d])
             ]
         )
@@ -164,7 +164,7 @@ function Base.getindex(
         dict = Dict{Block{1}, Vector{Tuple{Int, UnitRange{Int}, Base.OneTo{Int}}}}()
         for k in eachindex(I[d])
             bir = I[d][k]
-            b = Block(only(Tuple(bir.block)))
+            b = Block(Int(bir.block))
             r = only(bir.indices)
             push!(
                 get!(dict, b, Tuple{Int, UnitRange{Int}, Base.OneTo{Int}}[]),

From 3f4b4abe86a15dbc8a3b5ffd1f0d6220e2d7ed31 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Wed, 11 Mar 2026 10:04:29 -0400
Subject: [PATCH 12/19] Reenable fixed tests

---
 test/test_gradedarray.jl      | 2 +-
 test/test_tensoralgebraext.jl | 8 ++++----
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/test/test_gradedarray.jl b/test/test_gradedarray.jl
index 7d820bbe..70ab3abf 100644
--- a/test/test_gradedarray.jl
+++ b/test/test_gradedarray.jl
@@ -368,7 +368,7 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
             I = [Block(1)[1:1]]
             @test_broken size(b[I, :]) == (1, 4)
             @test_broken size(b[:, I]) == (4, 1)
-            @test_broken size(b[I, I]) == (1, 1)
+            @test size(b[I, I]) == (1, 1)
         end
     end
     @testset "Matrix multiplication" begin
diff --git a/test/test_tensoralgebraext.jl b/test/test_tensoralgebraext.jl
index 88402667..cd506b54 100644
--- a/test/test_tensoralgebraext.jl
+++ b/test/test_tensoralgebraext.jl
@@ -6,7 +6,7 @@ using GradedArrays: GradedArray, GradedMatrix, SU2, SectorDelta, U1, dual, flip,
 using Random: randn!
 using TensorAlgebra:
     FusionStyle, contract, matricize, tensor_product_axis, trivial_axis, unmatricize
-using Test: @test, @testset
+using Test: @test, @test_broken, @testset
 
 function randn_blockdiagonal(elt::Type, axes::Tuple)
     a = BlockSparseArray{elt}(undef, axes)
@@ -65,10 +65,10 @@ end
     @test unmatricize(m, (U1(1), U1(1)), (U1(-2), U1(-1))) isa SectorDelta
 end
 
-broken = true
+const contract_broken = true
 
 const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
-broken || @testset "`contract` `GradedArray` (eltype=$elt)" for elt in elts
+@testset "`contract` `GradedArray` (eltype=$elt)" for elt in elts
     @testset "matricize" begin
         d1 = gradedrange([U1(0) => 1, U1(1) => 1])
         d2 = gradedrange([U1(0) => 1, U1(1) => 1])
@@ -115,7 +115,7 @@ broken || @testset "`contract` `GradedArray` (eltype=$elt)" for elt in elts
         @test a == unmatricize(m, (), (d1, d2, dual(d1), dual(d2)))
     end
 
-    @testset "contract with U(1)" begin
+    contract_broken || @testset "contract with U(1)" begin
         d = gradedrange([U1(0) => 2, U1(1) => 3])
         a1 = randn_blockdiagonal(elt, (d, d, dual(d), dual(d)))
         a2 = randn_blockdiagonal(elt, (d, d, dual(d), dual(d)))

From fbd3c0e5422331fd3762291cba79acbf0cf9b52d Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Wed, 11 Mar 2026 10:39:03 -0400
Subject: [PATCH 13/19] Simplify unmatricize

---
 src/tensoralgebra.jl | 81 ++++++++++++++++++--------------------------
 1 file changed, 33 insertions(+), 48 deletions(-)

diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl
index 9c2fe44e..691a8292 100644
--- a/src/tensoralgebra.jl
+++ b/src/tensoralgebra.jl
@@ -119,20 +119,41 @@ function TensorAlgebra.unmatricize(
         return a
     end
 
-    # First, fuse axes to get `sectormergesortperm`.
-    # Then unpermute the blocks.
     fused_axes = matricize_axes(BlockReshapeFusion(), m, codomain_axes, domain_axes)
-
     blockperms = sectorsortperm.(fused_axes)
-    sorted_axes = map((r, I) -> only(axes(r[I])), fused_axes, blockperms)
 
-    # TODO: This is doing extra copies of the blocks,
-    # use `@view a[axes_prod...]` instead.
-    # That will require implementing some reindexing logic
-    # for this combination of slicing.
-    m_unblocked = m[sorted_axes...]
-    m_blockpermed = m_unblocked[invblockperm.(blockperms)...]
-    return unmatricize(FusionStyle(BlockSparseArray), m_blockpermed, blocked_axes)
+    # Build a BlockIndexRange index for each axis that composes the inverse block
+    # permutation with the splitting of m's merged blocks into fused sub-blocks.
+    # Each J[d][k] = Block(j)[r] means: dest block k comes from block j of m at subrange r.
+    J = map(fused_axes, blockperms, axes(m)) do fused_ax, blockperm, m_ax
+        n = length(blockperm)
+        # Linear scan: map sorted block k' (= fused block blockperm[k']) to (j, r) in m_ax.
+        # Requires that blocks of fused_ax subdivide blocks of m_ax (commensurate).
+        scan = Vector{Tuple{Int, UnitRange{Int}}}(undef, n)
+        j = 1
+        offset = 0
+        for k′ in 1:n
+            size_k′ = length(fused_ax[blockperm[k′]])
+            m_block_size = length(m_ax[Block(j)])
+            offset + size_k′ ≤ m_block_size ||
+                throw(ArgumentError("fused_ax blocks do not subdivide m_ax blocks"))
+            scan[k′] = (j, (offset + 1):(offset + size_k′))
+            offset += size_k′
+            if offset == m_block_size
+                j += 1
+                offset = 0
+            end
+        end
+        # Compose with inverse permutation: dest block k comes from scan[iperm[k]].
+        iperm = invblockperm(blockperm)
+        return [
+            let (j_k, r_k) = scan[Int(iperm[k])]
+                    Block(j_k)[r_k]
+            end for k in 1:n
+        ]
+    end
+
+    return unmatricize(FusionStyle(BlockSparseArray), m[J...], blocked_axes)
 end
 
 # Sort the blocks by sector and then merge the common sectors.
@@ -141,7 +162,7 @@ function sectormergesort(a::AbstractArray)
     return a[I...]
 end
 
-using BlockArrays: AbstractBlockVector, Block, BlockVector
+using BlockArrays: AbstractBlockVector, Block
 
 # Splitting: each I[d][k] = Block(b)[r] means dest block k comes from source block b
 # at subrange r. This is the inverse of the merging getindex below.
@@ -188,42 +209,6 @@ function Base.getindex(
     return a_dest
 end
 
-# GradedUnitRange index: compute BlockIndexRange per block by linear scan,
-# then delegate to the splitting getindex above.
-# Assumes blocks of each I[d] subdivide blocks of the corresponding axis of a
-# (as is the case in unmatricize, where I[d] is derived from the unmerged fused axis).
-function Base.getindex(
-        a::GradedArray{<:Any, N},
-        I::Vararg{GradedUnitRange, N}
-    ) where {N}
-    J = map(axes(a), I) do src_ax, tgt_ax
-        n = length(sectors(tgt_ax))
-        J_d = Vector{BlockIndexRange{1}}(undef, n)
-        j = 1
-        offset = 0
-        for k in 1:n
-            size_k = length(tgt_ax[Block(k)])
-            J_d[k] = Block(j)[(offset + 1):(offset + size_k)]
-            offset += size_k
-            if offset == length(src_ax[Block(j)])
-                j += 1
-                offset = 0
-            end
-        end
-        return J_d
-    end
-    return a[J...]
-end
-
-# Vector{Block{1}} index: block permutation with no merging.
-# Wraps as a singleton-group BlockVector and delegates to the merging getindex below.
-function Base.getindex(
-        a::GradedArray{<:Any, N},
-        I::Vararg{Vector{<:Block{1}}, N}
-    ) where {N}
-    return a[map(v -> BlockVector(v, fill(1, length(v))), I)...]
-end
-
 # Merging: each I[d] groups source blocks into destination blocks.
 function Base.getindex(
         a::GradedArray{<:Any, N},

From bf03c83b14080ad930d98aefb8843e6c075ea7ac Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Wed, 11 Mar 2026 10:57:42 -0400
Subject: [PATCH 14/19] Refactor and simplify

---
 src/tensoralgebra.jl | 68 ++++++++++++++++++--------------------------
 1 file changed, 27 insertions(+), 41 deletions(-)

diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl
index 691a8292..91901c5a 100644
--- a/src/tensoralgebra.jl
+++ b/src/tensoralgebra.jl
@@ -121,38 +121,7 @@ function TensorAlgebra.unmatricize(
 
     fused_axes = matricize_axes(BlockReshapeFusion(), m, codomain_axes, domain_axes)
     blockperms = sectorsortperm.(fused_axes)
-
-    # Build a BlockIndexRange index for each axis that composes the inverse block
-    # permutation with the splitting of m's merged blocks into fused sub-blocks.
-    # Each J[d][k] = Block(j)[r] means: dest block k comes from block j of m at subrange r.
-    J = map(fused_axes, blockperms, axes(m)) do fused_ax, blockperm, m_ax
-        n = length(blockperm)
-        # Linear scan: map sorted block k' (= fused block blockperm[k']) to (j, r) in m_ax.
-        # Requires that blocks of fused_ax subdivide blocks of m_ax (commensurate).
-        scan = Vector{Tuple{Int, UnitRange{Int}}}(undef, n)
-        j = 1
-        offset = 0
-        for k′ in 1:n
-            size_k′ = length(fused_ax[blockperm[k′]])
-            m_block_size = length(m_ax[Block(j)])
-            offset + size_k′ ≤ m_block_size ||
-                throw(ArgumentError("fused_ax blocks do not subdivide m_ax blocks"))
-            scan[k′] = (j, (offset + 1):(offset + size_k′))
-            offset += size_k′
-            if offset == m_block_size
-                j += 1
-                offset = 0
-            end
-        end
-        # Compose with inverse permutation: dest block k comes from scan[iperm[k]].
-        iperm = invblockperm(blockperm)
-        return [
-            let (j_k, r_k) = scan[Int(iperm[k])]
-                    Block(j_k)[r_k]
-            end for k in 1:n
-        ]
-    end
-
+    J = map(sectorunmatricize_index, fused_axes, blockperms, axes(m))
     return unmatricize(FusionStyle(BlockSparseArray), m[J...], blocked_axes)
 end
 
@@ -162,6 +131,31 @@ function sectormergesort(a::AbstractArray)
     return a[I...]
 end
 
+# Build a Vector{BlockIndexRange{1}} for one axis of unmatricize, composing the inverse
+# block permutation with the splitting of m_ax's merged blocks into fused sub-blocks.
+# Requires that blocks of fused_ax subdivide blocks of m_ax (commensurate by construction
+# in unmatricize, since fused_ax is derived from the block-reshape of m's axes).
+function sectorunmatricize_index(fused_ax, blockperm, m_ax)
+    n = length(blockperm)
+    scan = Vector{Tuple{Int, UnitRange{Int}}}(undef, n)
+    j = 1
+    offset = 0
+    for k′ in 1:n
+        size_k′ = length(fused_ax[blockperm[k′]])
+        m_block_size = length(m_ax[Block(j)])
+        offset + size_k′ ≤ m_block_size ||
+            throw(ArgumentError("fused_ax blocks do not subdivide m_ax blocks"))
+        scan[k′] = (j, (offset + 1):(offset + size_k′))
+        offset += size_k′
+        if offset == m_block_size
+            j += 1
+            offset = 0
+        end
+    end
+    iperm = invblockperm(blockperm)
+    return [let (j_k, r_k) = scan[Int(iperm[k])]; Block(j_k)[r_k]; end for k in 1:n]
+end
+
 using BlockArrays: AbstractBlockVector, Block
 
 # Splitting: each I[d][k] = Block(b)[r] means dest block k comes from source block b
@@ -170,15 +164,7 @@ function Base.getindex(
         a::GradedArray{<:Any, N},
         I::Vararg{AbstractVector{<:BlockIndexRange{1}}, N}
     ) where {N}
-    ax = axes(a)
-    ax_dest = ntuple(Val(N)) do d
-        return gradedrange(
-            [
-                sectors(ax[d])[Int(I[d][k].block)] => length(only(I[d][k].indices))
-                    for k in eachindex(I[d])
-            ]
-        )
-    end
+    ax_dest = ntuple(d -> only(axes(axes(a, d)[I[d]])), Val(N))
     a_dest = similar(a, ax_dest)
     # Map source block b → list of (dest block k, src subrange r, dest subrange 1:length(r))
     src_to_dests = ntuple(Val(N)) do d

From ada21f9f6d3d654d4fa4c9324189cdf767ab4328 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Wed, 11 Mar 2026 11:23:54 -0400
Subject: [PATCH 15/19] Small tweaks

---
 src/tensoralgebra.jl | 30 ++++++++++++++++++++++++------
 1 file changed, 24 insertions(+), 6 deletions(-)

diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl
index 91901c5a..2f5aeb92 100644
--- a/src/tensoralgebra.jl
+++ b/src/tensoralgebra.jl
@@ -158,24 +158,42 @@ end
 
 using BlockArrays: AbstractBlockVector, Block
 
+function checkindices(
+        a::GradedArray{<:Any, N},
+        I::NTuple{N, AbstractVector{<:BlockIndexRange{1}}}
+    ) where {N}
+    for d in 1:N
+        nblocks_d = length(axes(a, d))
+        for bir in I[d]
+            Int(bir.block) ≤ nblocks_d ||
+                throw(BlockBoundsError(a, ntuple(i -> i == d ? bir : I[i][1], Val(N))))
+        end
+    end
+end
+
 # Splitting: each I[d][k] = Block(b)[r] means dest block k comes from source block b
 # at subrange r. This is the inverse of the merging getindex below.
 function Base.getindex(
         a::GradedArray{<:Any, N},
         I::Vararg{AbstractVector{<:BlockIndexRange{1}}, N}
     ) where {N}
+    checkindices(a, I)
     ax_dest = ntuple(d -> only(axes(axes(a, d)[I[d]])), Val(N))
     a_dest = similar(a, ax_dest)
-    # Map source block b → list of (dest block k, src subrange r, dest subrange 1:length(r))
+    # Map source block b → list of (dest BlockIndexRange, src subrange).
+    # Stored blocks of a not referenced by I are skipped (partial block selection).
     src_to_dests = ntuple(Val(N)) do d
-        dict = Dict{Block{1}, Vector{Tuple{Int, UnitRange{Int}, Base.OneTo{Int}}}}()
+        key_type = typeof(Block(1))
+        dest_bir_type = Base.promote_op(getindex, key_type, Base.OneTo{Int})
+        val_type = Tuple{dest_bir_type, UnitRange{Int}}
+        dict = Dict{key_type, Vector{val_type}}()
         for k in eachindex(I[d])
             bir = I[d][k]
             b = Block(Int(bir.block))
             r = only(bir.indices)
             push!(
-                get!(dict, b, Tuple{Int, UnitRange{Int}, Base.OneTo{Int}}[]),
-                (k, r, Base.OneTo(length(r)))
+                get!(dict, b, val_type[]),
+                (Block(k)[Base.OneTo(length(r))], r)
             )
         end
         return dict
@@ -185,10 +203,10 @@ function Base.getindex(
         all(d -> haskey(src_to_dests[d], src_tuple[d]), 1:N) || continue
         dest_refs = ntuple(d -> src_to_dests[d][src_tuple[d]], Val(N))
         for combo in Iterators.product(dest_refs...)
-            dest_b = Block(ntuple(d -> combo[d][1], Val(N)))
+            dest_b = Block(ntuple(d -> only(Tuple(combo[d][1].block)), Val(N)))
             a_dest_b = @view!(a_dest[dest_b])
             src_r = ntuple(d -> combo[d][2], Val(N))
-            dest_r = ntuple(d -> combo[d][3], Val(N))
+            dest_r = ntuple(d -> only(combo[d][1].indices), Val(N))
             copyto!(@view(a_dest_b[dest_r...]), @view(a[bI_src][src_r...]))
         end
     end

From ce1ce2fb280c214f1b067c971bef9cf3ba3bec89 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Wed, 11 Mar 2026 11:39:47 -0400
Subject: [PATCH 16/19] Style

---
 src/tensoralgebra.jl | 22 +++++++++++++---------
 1 file changed, 13 insertions(+), 9 deletions(-)

diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl
index 2f5aeb92..f4a076f2 100644
--- a/src/tensoralgebra.jl
+++ b/src/tensoralgebra.jl
@@ -153,14 +153,17 @@ function sectorunmatricize_index(fused_ax, blockperm, m_ax)
         end
     end
     iperm = invblockperm(blockperm)
-    return [let (j_k, r_k) = scan[Int(iperm[k])]; Block(j_k)[r_k]; end for k in 1:n]
+    return [
+        let (j_k, r_k) = scan[Int(iperm[k])]
+                Block(j_k)[r_k]
+        end for k in 1:n
+    ]
 end
 
 using BlockArrays: AbstractBlockVector, Block
 
 function checkindices(
-        a::GradedArray{<:Any, N},
-        I::NTuple{N, AbstractVector{<:BlockIndexRange{1}}}
+        a::GradedArray{<:Any, N}, I::NTuple{N, AbstractVector{<:BlockIndexRange{1}}}
     ) where {N}
     for d in 1:N
         nblocks_d = length(axes(a, d))
@@ -169,13 +172,13 @@ function checkindices(
                 throw(BlockBoundsError(a, ntuple(i -> i == d ? bir : I[i][1], Val(N))))
         end
     end
+    return nothing
 end
 
 # Splitting: each I[d][k] = Block(b)[r] means dest block k comes from source block b
 # at subrange r. This is the inverse of the merging getindex below.
 function Base.getindex(
-        a::GradedArray{<:Any, N},
-        I::Vararg{AbstractVector{<:BlockIndexRange{1}}, N}
+        a::GradedArray{<:Any, N}, I::Vararg{AbstractVector{<:BlockIndexRange{1}}, N}
     ) where {N}
     checkindices(a, I)
     ax_dest = ntuple(d -> only(axes(axes(a, d)[I[d]])), Val(N))
@@ -203,11 +206,13 @@ function Base.getindex(
         all(d -> haskey(src_to_dests[d], src_tuple[d]), 1:N) || continue
         dest_refs = ntuple(d -> src_to_dests[d][src_tuple[d]], Val(N))
         for combo in Iterators.product(dest_refs...)
+            src_r = ntuple(d -> combo[d][2], Val(N))
+            src_data = @view(a[bI_src][src_r...])
+            iszero(src_data) && continue
             dest_b = Block(ntuple(d -> only(Tuple(combo[d][1].block)), Val(N)))
             a_dest_b = @view!(a_dest[dest_b])
-            src_r = ntuple(d -> combo[d][2], Val(N))
             dest_r = ntuple(d -> only(combo[d][1].indices), Val(N))
-            copyto!(@view(a_dest_b[dest_r...]), @view(a[bI_src][src_r...]))
+            copyto!(@view(a_dest_b[dest_r...]), src_data)
         end
     end
     return a_dest
@@ -215,8 +220,7 @@ end
 
 # Merging: each I[d] groups source blocks into destination blocks.
 function Base.getindex(
-        a::GradedArray{<:Any, N},
-        I::Vararg{AbstractBlockVector{<:Block{1}}, N}
+        a::GradedArray{<:Any, N}, I::Vararg{AbstractBlockVector{<:Block{1}}, N}
     ) where {N}
     ax_dest = ntuple(d -> only(axes(axes(a, d)[I[d]])), Val(N))
     a_dest = similar(a, ax_dest)

From 8b09cf3711360e66b575b5df6e46386b42428a67 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Wed, 11 Mar 2026 11:43:14 -0400
Subject: [PATCH 17/19] Style improvements

---
 src/tensoralgebra.jl | 21 +++++++++------------
 1 file changed, 9 insertions(+), 12 deletions(-)

diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl
index f4a076f2..55a64e82 100644
--- a/src/tensoralgebra.jl
+++ b/src/tensoralgebra.jl
@@ -137,27 +137,24 @@ end
 # in unmatricize, since fused_ax is derived from the block-reshape of m's axes).
 function sectorunmatricize_index(fused_ax, blockperm, m_ax)
     n = length(blockperm)
-    scan = Vector{Tuple{Int, UnitRange{Int}}}(undef, n)
+    bir_type = Base.promote_op(getindex, Block{1, Int}, UnitRange{Int})
+    J = Vector{bir_type}(undef, n)
     j = 1
     offset = 0
     for k′ in 1:n
-        size_k′ = length(fused_ax[blockperm[k′]])
+        k = Int(blockperm[k′])
+        size_k = length(fused_ax[Block(k)])
         m_block_size = length(m_ax[Block(j)])
-        offset + size_k′ ≤ m_block_size ||
+        offset + size_k ≤ m_block_size ||
             throw(ArgumentError("fused_ax blocks do not subdivide m_ax blocks"))
-        scan[k′] = (j, (offset + 1):(offset + size_k′))
-        offset += size_k′
+        J[k] = Block(j)[(offset + 1):(offset + size_k)]
+        offset += size_k
         if offset == m_block_size
             j += 1
             offset = 0
         end
     end
-    iperm = invblockperm(blockperm)
-    return [
-        let (j_k, r_k) = scan[Int(iperm[k])]
-                Block(j_k)[r_k]
-        end for k in 1:n
-    ]
+    return J
 end
 
 using BlockArrays: AbstractBlockVector, Block
@@ -186,7 +183,7 @@ function Base.getindex(
     # Map source block b → list of (dest BlockIndexRange, src subrange).
     # Stored blocks of a not referenced by I are skipped (partial block selection).
     src_to_dests = ntuple(Val(N)) do d
-        key_type = typeof(Block(1))
+        key_type = Block{1, Int}
         dest_bir_type = Base.promote_op(getindex, key_type, Base.OneTo{Int})
         val_type = Tuple{dest_bir_type, UnitRange{Int}}
         dict = Dict{key_type, Vector{val_type}}()

From 68dcc2a6f86fcd0dbd980ff3758aaba1e16e1da4 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Wed, 11 Mar 2026 11:52:28 -0400
Subject: [PATCH 18/19] Rename

---
 src/tensoralgebra.jl | 22 +++++++++++-----------
 1 file changed, 11 insertions(+), 11 deletions(-)

diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl
index 55a64e82..abeb9ffa 100644
--- a/src/tensoralgebra.jl
+++ b/src/tensoralgebra.jl
@@ -121,7 +121,7 @@ function TensorAlgebra.unmatricize(
 
     fused_axes = matricize_axes(BlockReshapeFusion(), m, codomain_axes, domain_axes)
     blockperms = sectorsortperm.(fused_axes)
-    J = map(sectorunmatricize_index, fused_axes, blockperms, axes(m))
+    J = map(invblockmergeperm, fused_axes, blockperms, axes(m))
     return unmatricize(FusionStyle(BlockSparseArray), m[J...], blocked_axes)
 end
 
@@ -131,11 +131,11 @@ function sectormergesort(a::AbstractArray)
     return a[I...]
 end
 
-# Build a Vector{BlockIndexRange{1}} for one axis of unmatricize, composing the inverse
-# block permutation with the splitting of m_ax's merged blocks into fused sub-blocks.
-# Requires that blocks of fused_ax subdivide blocks of m_ax (commensurate by construction
-# in unmatricize, since fused_ax is derived from the block-reshape of m's axes).
-function sectorunmatricize_index(fused_ax, blockperm, m_ax)
+# Returns a Vector{BlockIndexRange{1}} mapping each block of fine_ax (in original order)
+# to its position (block + subrange) within the merged axis merged_ax, given the block
+# permutation blockperm used to sort and merge fine_ax into merged_ax.
+# Requires that blocks of fine_ax subdivide blocks of merged_ax.
+function invblockmergeperm(fine_ax, blockperm, merged_ax)
     n = length(blockperm)
     bir_type = Base.promote_op(getindex, Block{1, Int}, UnitRange{Int})
     J = Vector{bir_type}(undef, n)
@@ -143,13 +143,13 @@ function sectorunmatricize_index(fused_ax, blockperm, m_ax)
     offset = 0
     for k′ in 1:n
         k = Int(blockperm[k′])
-        size_k = length(fused_ax[Block(k)])
-        m_block_size = length(m_ax[Block(j)])
-        offset + size_k ≤ m_block_size ||
-            throw(ArgumentError("fused_ax blocks do not subdivide m_ax blocks"))
+        size_k = length(fine_ax[Block(k)])
+        merged_block_size = length(merged_ax[Block(j)])
+        offset + size_k ≤ merged_block_size ||
+            throw(ArgumentError("fine_ax blocks do not subdivide merged_ax blocks"))
         J[k] = Block(j)[(offset + 1):(offset + size_k)]
         offset += size_k
-        if offset == m_block_size
+        if offset == merged_block_size
             j += 1
             offset = 0
         end

From 4eef0ed1e7b737747c63fc46cfd043400c5a8d96 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Wed, 11 Mar 2026 11:54:52 -0400
Subject: [PATCH 19/19] Style

---
 src/tensoralgebra.jl | 7 ++-----
 1 file changed, 2 insertions(+), 5 deletions(-)

diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl
index abeb9ffa..81b15590 100644
--- a/src/tensoralgebra.jl
+++ b/src/tensoralgebra.jl
@@ -191,10 +191,7 @@ function Base.getindex(
             bir = I[d][k]
             b = Block(Int(bir.block))
             r = only(bir.indices)
-            push!(
-                get!(dict, b, val_type[]),
-                (Block(k)[Base.OneTo(length(r))], r)
-            )
+            push!(get!(dict, b, val_type[]), (Block(k)[Base.axes1(r)], r))
         end
         return dict
     end
@@ -219,7 +216,7 @@ end
 function Base.getindex(
         a::GradedArray{<:Any, N}, I::Vararg{AbstractBlockVector{<:Block{1}}, N}
     ) where {N}
-    ax_dest = ntuple(d -> only(axes(axes(a, d)[I[d]])), Val(N))
+    ax_dest = ntuple(d -> Base.axes1(axes(a, d)[I[d]]), Val(N))
     a_dest = similar(a, ax_dest)
     ax = axes(a)
     # Map source Block -> BlockIndexRange encoding dest block + subrange within it