From 4c6fb62b2524768fef1210ddc10a51cd8c473b58 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 2 Jun 2026 12:35:10 -0400
Subject: [PATCH 01/31] Switch from Base.adjoint to Base.conj for sectors and
 axes

TensorKitSectors defines the dual representation via `Base.conj` on every sector type (and `dual(::Sector) = conj(a)` is an alias of it), and never defines `Base.adjoint` on sectors. The `Base.adjoint(::SectorRange) = dual` overload in GradedArrays drifted away from that convention, and the `U1(0)'`, `SU2(1//2)'` shorthand for "dual sector" overloaded `'` (adjoint) onto irrep labels, which is not a mathematically meaningful operation on labels.

Drops `Base.adjoint = dual` on `GradedOneTo`, `SectorOneTo`, and `SectorRange`. Adds `Base.conj = dual` on the same three. Migrates every `U1(0)'` / `SU2(1//2)'` test call site (and the in-source docstring example) to `conj(U1(0))` / `conj(SU2(1//2))`.

Pairs with NamedDimsArrays.jl PR #229's `Base.conj(::AbstractNamedUnitRange)` forwarder, so a named graded axis's `conj` propagates through the wrapper and flips sector arrows on the underlying axis. This is what `TensorAlgebra.similar_map` needs to compose with GradedArrays-backed prototypes.
---
 Project.toml                |  2 +-
 src/gradedoneto.jl          |  4 ++--
 src/sectoroneto.jl          |  2 +-
 src/sectorrange.jl          |  2 +-
 test/test_abelianarray.jl   | 10 +++++-----
 test/test_gradedoneto.jl    | 30 +++++++++++++++---------------
 test/test_sectorarray.jl    | 28 ++++++++++++++--------------
 test/test_sectoridentity.jl | 12 ++++++------
 test/test_sectormatrix.jl   |  6 +++---
 test/test_sectoroneto.jl    | 20 ++++++++++----------
 test/test_tensoralgebra.jl  |  4 ++--
 11 files changed, 60 insertions(+), 60 deletions(-)

diff --git a/Project.toml b/Project.toml
index e620873b..e71bc40c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "GradedArrays"
 uuid = "bc96ca6e-b7c8-4bb6-888e-c93f838762c2"
-version = "0.9.1"
+version = "0.9.2"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
 
 [workspace]
diff --git a/src/gradedoneto.jl b/src/gradedoneto.jl
index 0413ce1c..85ed4c33 100644
--- a/src/gradedoneto.jl
+++ b/src/gradedoneto.jl
@@ -125,7 +125,7 @@ function flip(g::GradedOneTo)
     return GradedOneTo(new_nondual, datalengths(g), !isdual(g))
 end
 flip_dual(g::GradedOneTo) = isdual(g) ? flip(g) : g
-Base.adjoint(g::GradedOneTo) = dual(g)
+Base.conj(g::GradedOneTo) = dual(g)
 
 to_gradedrange(g::GradedOneTo) = g
 
@@ -224,7 +224,7 @@ Non-dual inputs produce a non-dual axis; dual inputs produce a dual axis.
 
 ```julia
 gradedrange([U1(0) => 2, U1(1) => 3])     # non-dual
-gradedrange([U1(0)' => 2, U1(1)' => 3])   # dual
+gradedrange([conj(U1(0)) => 2, conj(U1(1)) => 3])   # dual
 ```
 """
 function gradedrange(
diff --git a/src/sectoroneto.jl b/src/sectoroneto.jl
index b74d0980..df97e03d 100644
--- a/src/sectoroneto.jl
+++ b/src/sectoroneto.jl
@@ -114,4 +114,4 @@ function Base.show(io::IO, r::SectorOneTo)
     return nothing
 end
 
-Base.adjoint(r::SectorOneTo) = dual(r)
+Base.conj(r::SectorOneTo) = dual(r)
diff --git a/src/sectorrange.jl b/src/sectorrange.jl
index 8a323461..5a236c20 100644
--- a/src/sectorrange.jl
+++ b/src/sectorrange.jl
@@ -105,7 +105,7 @@ nondual(r::SectorRange) = isdual(r) ? dual(r) : r
 fermionparity(c::SectorRange) = TKS.fermionparity(c.label)
 twist(c::SectorRange) = TKS.twist(c.label)
 
-Base.adjoint(r1::SectorRange) = dual(r1)
+Base.conj(r1::SectorRange) = dual(r1)
 
 # ===============================  Fusion rule interface  ==================================
 
diff --git a/test/test_abelianarray.jl b/test/test_abelianarray.jl
index 2ed234a6..8b7f44b0 100644
--- a/test/test_abelianarray.jl
+++ b/test/test_abelianarray.jl
@@ -50,12 +50,12 @@ using Test: @test, @test_throws, @testset
     end
 
     @testset "Block getindex returns correct sectors" begin
-        g1_dual = gradedrange([U1(0) => 2, U1(1) => 3])'
+        g1_dual = conj(gradedrange([U1(0) => 2, U1(1) => 3]))
         a = AbelianGradedArray{Float64}(undef, g1_dual, g2)
         a[Block(1, 1)] = ones(2, 1)
 
         blk = a[Block(1, 1)]
-        @test sectoraxes(blk) == (U1(0)', U1(0))
+        @test sectoraxes(blk) == (conj(U1(0)), U1(0))
     end
 
     @testset "Block getindex for unstored block errors" begin
@@ -98,8 +98,8 @@ using Test: @test, @test_throws, @testset
     end
 
     @testset "Dual axes" begin
-        g1_dual = gradedrange([U1(0) => 2, U1(1) => 3])'
-        g2_dual = gradedrange([U1(0) => 1, U1(-1) => 2])'
+        g1_dual = conj(gradedrange([U1(0) => 2, U1(1) => 3]))
+        g2_dual = conj(gradedrange([U1(0) => 1, U1(-1) => 2]))
         a = AbelianGradedArray{Float64}(undef, g1_dual, g2_dual)
 
         @test isdual(axes(a, 1)) == true
@@ -108,7 +108,7 @@ using Test: @test, @test_throws, @testset
 
         a[Block(1, 1)] = ones(2, 1)
         blk = a[Block(1, 1)]
-        @test sectoraxes(blk) == (U1(0)', U1(0)')
+        @test sectoraxes(blk) == (conj(U1(0)), conj(U1(0)))
     end
 
     @testset "similar" begin
diff --git a/test/test_gradedoneto.jl b/test/test_gradedoneto.jl
index 558e1b04..47073ef3 100644
--- a/test/test_gradedoneto.jl
+++ b/test/test_gradedoneto.jl
@@ -14,20 +14,20 @@ using Test: @test, @test_throws, @testset
     end
 
     @testset "gradedrange from dual SectorRange labels" begin
-        g = gradedrange([U1(0)' => 2, U1(1)' => 3])
+        g = gradedrange([conj(U1(0)) => 2, conj(U1(1)) => 3])
         @test g isa GradedOneTo{U1}
-        @test sectors(g) == [U1(0)', U1(1)']
+        @test sectors(g) == [conj(U1(0)), conj(U1(1))]
         @test datalengths(g) == [2, 3]
         @test isdual(g) == true
     end
 
     @testset "gradedrange mixed isdual throws" begin
-        @test_throws ArgumentError gradedrange([U1(0) => 2, U1(1)' => 3])
+        @test_throws ArgumentError gradedrange([U1(0) => 2, conj(U1(1)) => 3])
     end
 
-    @testset "dual via adjoint (U1)" begin
+    @testset "dual via conj (U1)" begin
         g = gradedrange([U1(0) => 2, U1(1) => 3])
-        gd = g'
+        gd = conj(g)
         @test isdual(gd) == true
         @test sectors(gd) == dual.(sectors(g))
         @test datalengths(gd) == datalengths(g)
@@ -36,7 +36,7 @@ using Test: @test, @test_throws, @testset
     @testset "double dual is identity" begin
         g = gradedrange([U1(0) => 2, U1(1) => 3])
         @test dual(dual(g)) == g
-        @test g'' == g
+        @test conj(conj(g)) == g
     end
 
     @testset "sectors accessor — non-dual" begin
@@ -45,8 +45,8 @@ using Test: @test, @test_throws, @testset
     end
 
     @testset "sectors accessor — dual" begin
-        g = gradedrange([U1(0) => 2, U1(1) => 3])'
-        @test sectors(g) == [U1(0)', U1(1)']
+        g = conj(gradedrange([U1(0) => 2, U1(1) => 3]))
+        @test sectors(g) == [conj(U1(0)), conj(U1(1))]
     end
 
     @testset "blocklength" begin
@@ -84,9 +84,9 @@ using Test: @test, @test_throws, @testset
     end
 
     @testset "flip on dual" begin
-        g = gradedrange([U1(1) => 3])'
+        g = conj(gradedrange([U1(1) => 3]))
         gf = flip(g)
-        @test sectors(gf) == [flip(U1(1)')]
+        @test sectors(gf) == [flip(conj(U1(1)))]
         @test isdual(gf) == false  # was dual, flip toggles
     end
 
@@ -94,7 +94,7 @@ using Test: @test, @test_throws, @testset
         g1 = gradedrange([U1(0) => 2, U1(1) => 3])
         g2 = gradedrange([U1(0) => 2, U1(1) => 3])
         g3 = gradedrange([U1(0) => 2, U1(1) => 4])
-        g4 = gradedrange([U1(0) => 2, U1(1) => 3])'
+        g4 = conj(gradedrange([U1(0) => 2, U1(1) => 3]))
 
         @test g1 == g2
         @test g1 != g3  # different multiplicity
@@ -125,7 +125,7 @@ using Test: @test, @test_throws, @testset
         @test !contains(str, "dual")
 
         # Dual axes factor the `dual` to the outside.
-        gd = g'
+        gd = conj(g)
         strd = sprint(show, gd)
         @test !endswith(strd, "'")
         @test startswith(strd, "dual(gradedrange(")
@@ -150,9 +150,9 @@ using Test: @test, @test_throws, @testset
         @test blocklength(g) == 2
         @test length(g) == 1 * 1 + 2 * 2  # 5
 
-        gd = g'
+        gd = conj(g)
         @test isdual(gd) == true
-        @test sectors(gd) == [SU2(0)', SU2(1 // 2)']
+        @test sectors(gd) == [conj(SU2(0)), conj(SU2(1 // 2))]
     end
 
     @testset "SU2 gradedrange from SectorRange" begin
@@ -180,7 +180,7 @@ using Test: @test, @test_throws, @testset
         g = gradedrange([U1(1) => 2, U1(0) => 3])
         @test tensor_product(g) == gradedrange([U1(0) => 3, U1(1) => 2])
 
-        gd = gradedrange([U1(1) => 2, U1(0) => 3])'
+        gd = conj(gradedrange([U1(1) => 2, U1(0) => 3]))
         tp_d = tensor_product(gd)
         @test !isdual(tp_d)
 
diff --git a/test/test_sectorarray.jl b/test/test_sectorarray.jl
index 0eb87432..2fd4d195 100644
--- a/test/test_sectorarray.jl
+++ b/test/test_sectorarray.jl
@@ -7,7 +7,7 @@ using Test: @test, @test_throws, @testset
 @testset "AbelianSectorArray" begin
     @testset "Construction from SectorRange tuples" begin
         data = [1.0 2.0; 3.0 4.0]
-        sa = AbelianSectorArray((U1(1), U1(-1)'), data)
+        sa = AbelianSectorArray((U1(1), conj(U1(-1))), data)
         @test sa isa AbelianSectorArray{Float64, 2, Matrix{Float64}, U1}
         @test sa isa AbstractArray{Float64, 2}
         @test !(sa isa GradedArrays.KroneckerArrays.AbstractKroneckerArray)
@@ -15,8 +15,8 @@ using Test: @test, @test_throws, @testset
 
     @testset "Construction with dual sectors" begin
         data = [1.0 2.0; 3.0 4.0]
-        sa = AbelianSectorArray((U1(1), U1(-1)'), data)
-        @test sectoraxes(sa) == (U1(1), U1(-1)')
+        sa = AbelianSectorArray((U1(1), conj(U1(-1))), data)
+        @test sectoraxes(sa) == (U1(1), conj(U1(-1)))
     end
 
     @testset "Undef constructor (SectorOneTo)" begin
@@ -31,11 +31,11 @@ using Test: @test, @test_throws, @testset
 
     @testset "Primitive accessors" begin
         data = ones(2, 3, 4)
-        sa = AbelianSectorArray((U1(1), U1(0)', U1(-1)), data)
+        sa = AbelianSectorArray((U1(1), conj(U1(0)), U1(-1)), data)
 
-        @test sectoraxes(sa) == (U1(1), U1(0)', U1(-1))
+        @test sectoraxes(sa) == (U1(1), conj(U1(0)), U1(-1))
         @test sectoraxes(sa, 1) == U1(1)
-        @test sectoraxes(sa, 2) == U1(0)'
+        @test sectoraxes(sa, 2) == conj(U1(0))
         @test sectoraxes(sa, 3) == U1(-1)
         @test isdual(sa, 1) == false
         @test isdual(sa, 2) == true
@@ -44,15 +44,15 @@ using Test: @test, @test_throws, @testset
 
     @testset "Derived accessors — sectoraxes" begin
         data = ones(2, 3)
-        sa = AbelianSectorArray((U1(1), U1(-1)'), data)
+        sa = AbelianSectorArray((U1(1), conj(U1(-1))), data)
         @test sectoraxes(sa, 1) == U1(1)
-        @test sectoraxes(sa, 2) == U1(-1)'
-        @test sectoraxes(sa) == (U1(1), U1(-1)')
+        @test sectoraxes(sa, 2) == conj(U1(-1))
+        @test sectoraxes(sa) == (U1(1), conj(U1(-1)))
     end
 
     @testset "sector(::AbelianSectorArray) returns AbelianSectorDelta" begin
         data = ones(2, 3)
-        sa = AbelianSectorArray((U1(1), U1(-1)'), data)
+        sa = AbelianSectorArray((U1(1), conj(U1(-1))), data)
         sd = sector(sa)
         @test sd isa AbelianSectorDelta{Float64, 2, U1}
         @test axes(sd) == sectoraxes(sa)
@@ -119,7 +119,7 @@ using Test: @test, @test_throws, @testset
 
     @testset "3D AbelianSectorArray" begin
         data = ones(2, 3, 4)
-        sa = AbelianSectorArray((U1(1), U1(0)', U1(-1)), data)
+        sa = AbelianSectorArray((U1(1), conj(U1(0)), U1(-1)), data)
         @test size(sa) == (2, 3, 4)
         @test ndims(sa) == 3
         @test sa[1, 2, 3] == 1.0
@@ -127,10 +127,10 @@ using Test: @test, @test_throws, @testset
 
     @testset "permutedims" begin
         data = [1.0 2.0 3.0; 4.0 5.0 6.0]
-        sa = AbelianSectorArray((U1(1), U1(0)'), data)
+        sa = AbelianSectorArray((U1(1), conj(U1(0))), data)
         sa_perm = permutedims(sa, (2, 1))
         @test size(sa_perm) == (3, 2)
-        @test sectoraxes(sa_perm) == (U1(0)', U1(1))
+        @test sectoraxes(sa_perm) == (conj(U1(0)), U1(1))
         @test sa_perm[1, 1] == 1.0
         @test sa_perm[1, 2] == 4.0
     end
@@ -141,7 +141,7 @@ using Test: @test, @test_throws, @testset
         b_data = [5.0 6.0; 7.0 8.0]
         c_data = zeros(2, 2)
         a = AbelianSectorArray((U1(0), U1(1)), a_data)
-        b = AbelianSectorArray((U1(1)', U1(0)), b_data)
+        b = AbelianSectorArray((conj(U1(1)), U1(0)), b_data)
         c = AbelianSectorArray((U1(0), U1(0)), c_data)
         mul!(c, a, b, 1.0, 0.0)
         @test data(c) ≈ a_data * b_data
diff --git a/test/test_sectoridentity.jl b/test/test_sectoridentity.jl
index 4567fe9a..56ec239a 100644
--- a/test/test_sectoridentity.jl
+++ b/test/test_sectoridentity.jl
@@ -11,8 +11,8 @@ using Test: @test, @testset
     end
 
     @testset "Construction from dual SectorRange (dual flag preserved)" begin
-        si = SectorIdentity{Float64}(U1(1)')
-        @test axes(si, 1) == U1(1)'
+        si = SectorIdentity{Float64}(conj(U1(1)))
+        @test axes(si, 1) == conj(U1(1))
         @test axes(si, 2) == U1(1)
     end
 
@@ -21,21 +21,21 @@ using Test: @test, @testset
         @test size(si) == (1, 1)
         ax1, ax2 = axes(si)
         @test ax1 == U1(0)
-        @test ax2 == U1(0)'
+        @test ax2 == conj(U1(0))
     end
 
     @testset "size and axes — SU2 j=1/2 (dim=2)" begin
         si = SectorIdentity{Float64}(SU2(1 // 2))
         @test size(si) == (2, 2)
         @test axes(si, 1) == SU2(1 // 2)
-        @test axes(si, 2) == SU2(1 // 2)'
+        @test axes(si, 2) == conj(SU2(1 // 2))
     end
 
     @testset "axes as sector ranges" begin
         si = SectorIdentity{Float64}(U1(1))
-        @test axes(si) == (U1(1), U1(1)')
+        @test axes(si) == (U1(1), conj(U1(1)))
         @test axes(si, 1) == U1(1)
-        @test axes(si, 2) == U1(1)'
+        @test axes(si, 2) == conj(U1(1))
     end
 
     @testset "sectortype" begin
diff --git a/test/test_sectormatrix.jl b/test/test_sectormatrix.jl
index 0020e7bd..c3f74ba7 100644
--- a/test/test_sectormatrix.jl
+++ b/test/test_sectormatrix.jl
@@ -23,9 +23,9 @@ using Test: @test, @test_throws, @testset
     @testset "sectoraxes" begin
         d = ones(2, 3)
         sm = SectorMatrix(U1(1), d)
-        @test sectoraxes(sm) == (U1(1), U1(1)')
+        @test sectoraxes(sm) == (U1(1), conj(U1(1)))
         @test sectoraxes(sm, 1) == U1(1)
-        @test sectoraxes(sm, 2) == U1(1)'
+        @test sectoraxes(sm, 2) == conj(U1(1))
     end
 
     @testset "sector returns SectorIdentity" begin
@@ -132,7 +132,7 @@ using Test: @test, @test_throws, @testset
         @test sm4 isa SectorMatrix
         @test data(sm4) ≈ 2.0 .* data(sm)
 
-        sa = AbelianSectorArray((U1(0), U1(0)'), zeros(2, 2))
+        sa = AbelianSectorArray((U1(0), conj(U1(0))), zeros(2, 2))
         sa .= sm
         @test data(sa) ≈ data(sm)
     end
diff --git a/test/test_sectoroneto.jl b/test/test_sectoroneto.jl
index e0a9b3d3..d682e9c8 100644
--- a/test/test_sectoroneto.jl
+++ b/test/test_sectoroneto.jl
@@ -19,9 +19,9 @@ using Test: @test, @testset
     end
 
     @testset "U1 dual sector accessor" begin
-        si = SectorOneTo(U1(1)', 3)
+        si = SectorOneTo(conj(U1(1)), 3)
         @test isdual(si) == true
-        @test sector(si) == U1(1)'
+        @test sector(si) == conj(U1(1))
     end
 
     @testset "SU2 construction and accessors" begin
@@ -33,8 +33,8 @@ using Test: @test, @testset
     end
 
     @testset "SU2 dual" begin
-        si = SectorOneTo(SU2(1)', 2)
-        @test sector(si) == SU2(1)'
+        si = SectorOneTo(conj(SU2(1)), 2)
+        @test sector(si) == conj(SU2(1))
     end
 
     @testset "Collection-like interface" begin
@@ -62,7 +62,7 @@ using Test: @test, @testset
     @testset "dual" begin
         si = SectorOneTo(U1(1), 3)
         sid = dual(si)
-        @test sector(sid) == U1(1)'
+        @test sector(sid) == conj(U1(1))
         @test datalength(sid) == 3
         @test isdual(sid) == true
         @test dual(sid) == si  # double dual is identity
@@ -87,7 +87,7 @@ using Test: @test, @testset
     @testset "equality" begin
         si1 = SectorOneTo(U1(1), 3)
         si2 = SectorOneTo(U1(1), 3)
-        si3 = SectorOneTo(U1(1)', 3)
+        si3 = SectorOneTo(conj(U1(1)), 3)
         si4 = SectorOneTo(U1(2), 3)
         si5 = SectorOneTo(U1(1), 4)
 
@@ -120,7 +120,7 @@ using Test: @test, @testset
         @test !contains(str, "dual")
 
         # Duality is factored to the outside: `dual(SectorOneTo(...))`.
-        sid = SectorOneTo(U1(1)', 3)
+        sid = SectorOneTo(conj(U1(1)), 3)
         strd = sprint(show, sid)
         @test startswith(strd, "dual(SectorOneTo(")
         @test endswith(strd, "))")
@@ -128,8 +128,8 @@ using Test: @test, @testset
     end
 
     @testset "dual sectors accessor for collection interface" begin
-        si = SectorOneTo(U1(1)', 3)
-        @test sectors(si) == [U1(1)']
+        si = SectorOneTo(conj(U1(1)), 3)
+        @test sectors(si) == [conj(U1(1))]
     end
 
     @testset "tensor_product (abelian)" begin
@@ -143,7 +143,7 @@ using Test: @test, @testset
 
         @test tensor_product(si1) == si1
 
-        si1d = SectorOneTo(U1(1)', 3)
+        si1d = SectorOneTo(conj(U1(1)), 3)
         tp1 = tensor_product(si1d)
         @test !isdual(tp1)
 
diff --git a/test/test_tensoralgebra.jl b/test/test_tensoralgebra.jl
index 65286eeb..84eaf7be 100644
--- a/test/test_tensoralgebra.jl
+++ b/test/test_tensoralgebra.jl
@@ -308,8 +308,8 @@ end
 @testset "unmatricize AbelianSectorMatrix with SectorOneTo axes" begin
     # Create a 3D AbelianSectorArray, matricize it, then unmatricize and verify roundtrip
     codomain_ax = SectorOneTo(U1(0), 2)
-    domain_ax1 = SectorOneTo(U1(0)', 3)
-    domain_ax2 = SectorOneTo(U1(1)', 4)
+    domain_ax1 = SectorOneTo(conj(U1(0)), 3)
+    domain_ax2 = SectorOneTo(conj(U1(1)), 4)
 
     data_3d = randn!(Array{Float64}(undef, 2, 3, 4))
     s = AbelianSectorArray(

From 3c935fc3049e4ff55f924b9f0224e7f449cf17b3 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 2 Jun 2026 16:59:07 -0400
Subject: [PATCH 02/31] Define block-aware `Base.conj`, `copy`, `norm`, `dot`,
 scalar ops, and `+`-mapreduce on `AbelianGradedArray`

The default `Base`/`LinearAlgebra` methods on `AbstractArray` iterate or scalar-index, both of which `AbstractGradedArray` rejects. This patch fills in the block-wise specializations needed for the common ops on `AbelianGradedArray`: `copy`, `conj` (which also flips axis duality, matching the existing `Base.conj` on the axis types), `norm`, `dot`, `:*`/`:/` with a scalar, and `mapreduce` over `+`/`add_sum` (so `sum` works without iterating). Unstored blocks are treated as implicit zero, which is correct for `+`-mapreduce, `norm`, and `dot`.

Also: `Base.fill!(::FusedGradedVector, val)` now fills every stored block with `val`. The previous `iszero(val) || throw(...)` was too strict, since each block of a fused vector is a per-sector segment whose entries are independently meaningful (`val * I` in the sector-fused basis is the canonical use). `Base.fill!(::AbelianGradedArray, val)` keeps the `iszero` guard, because filling stored blocks with a nonzero value produces a tensor that mixes the value on symmetry-allowed positions with implicit zeros on symmetry-forbidden ones, which is not a meaningful object.

And: scalar `TensorAlgebra.unmatricize(::SectorFusion, ::FusedGradedMatrix, (), ())` is implemented (was an explicit `error`). Reads the trivial-sector entry of `m` into a 0-D `Array{T}`, which is what the unfused codomain/domain axes describe.

Co-authored-by: Claude <noreply@anthropic.com>
---
 src/abeliangradedarray.jl | 70 +++++++++++++++++++++++++++++++++++++++
 src/fusedgradedvector.jl  |  8 ++---
 src/fusion.jl             | 11 ++++--
 3 files changed, 83 insertions(+), 6 deletions(-)

diff --git a/src/abeliangradedarray.jl b/src/abeliangradedarray.jl
index 6e450955..19f7b8c4 100644
--- a/src/abeliangradedarray.jl
+++ b/src/abeliangradedarray.jl
@@ -260,6 +260,25 @@ function Base.similar(a::AbelianGradedArray{T}) where {T}
     return similar(a, T)
 end
 
+# Block-wise copy; the default falls through to scalar `copyto!`.
+function Base.copy(a::AbelianGradedArray{T, N, D, S}) where {T, N, D, S}
+    return AbelianGradedArray{T, N, D, S}(
+        Dict{NTuple{N, Int}, D}(k => copy(v) for (k, v) in a.blockdata),
+        a.axes
+    )
+end
+
+# Conjugate element-wise *and* flip axis duality, mirroring `Base.conj` on the
+# axis types (`SectorRange`/`GradedOneTo`/`SectorOneTo`). Without the axis flip,
+# `conj(t)` would leave bra-layer tensors with the same duality as the ket, and
+# any contraction between them would silently pair non-dual against non-dual.
+function Base.conj(a::AbelianGradedArray{T, N, D, S}) where {T, N, D, S}
+    return AbelianGradedArray{T, N, D, S}(
+        Dict{NTuple{N, Int}, D}(k => conj(v) for (k, v) in a.blockdata),
+        map(conj, a.axes)
+    )
+end
+
 # ---------------------------------------------------------------------------
 #  sectortype
 # ---------------------------------------------------------------------------
@@ -311,6 +330,57 @@ function Base.fill!(a::AbelianGradedArray, v)
     return FI.zero!(a)
 end
 
+function LinearAlgebra.norm(a::AbelianGradedArray, p::Real = 2)
+    if p == Inf
+        isempty(a.blockdata) && return zero(float(real(eltype(a))))
+        return maximum(Base.Fix2(LinearAlgebra.norm, p), values(a.blockdata))
+    elseif p > 0
+        s = zero(float(real(eltype(a))))
+        for b in values(a.blockdata)
+            s += LinearAlgebra.norm(b, p)^p
+        end
+        return s^inv(p)
+    else
+        throw(ArgumentError("Norm with non-positive p ($p) is not defined"))
+    end
+end
+
+function LinearAlgebra.dot(a::AbelianGradedArray, b::AbelianGradedArray)
+    s = zero(LinearAlgebra.dot(zero(eltype(a)), zero(eltype(b))))
+    for (k, ablk) in pairs(a.blockdata)
+        haskey(b.blockdata, k) || continue
+        s += LinearAlgebra.dot(ablk, b.blockdata[k])
+    end
+    return s
+end
+
+# Reductions iterate scalar-by-scalar by default. Unstored blocks are zero, so
+# `+`-style reductions can skip them. Other ops require materializing — refuse
+# rather than silently dropping unstored contributions.
+function Base.mapreduce(
+        f, op::Union{typeof(+), typeof(Base.add_sum)},
+        a::AbelianGradedArray;
+        init = zero(eltype(a))
+    )
+    s = init
+    for b in values(a.blockdata)
+        s = op(s, mapreduce(f, op, b))
+    end
+    return s
+end
+
+# Block-wise scalar multiplication. The default `a .* x` / `a ./ x` paths broadcast
+# element-wise and scalar-index opaque storage.
+Base.:*(a::AbelianGradedArray, x::Number) = scale!(copy(a), x)
+Base.:*(x::Number, a::AbelianGradedArray) = a * x
+Base.:/(a::AbelianGradedArray, x::Number) = a * inv(x)
+
+# `LinearAlgebra.normalize` infers its result eltype via `typeof(first(a)/nrm)`,
+# which scalar-indexes opaque storage.
+function LinearAlgebra.normalize(a::AbelianGradedArray, p::Real = 2)
+    return a / LinearAlgebra.norm(a, p)
+end
+
 # ---------------------------------------------------------------------------
 #  twist!
 # ---------------------------------------------------------------------------
diff --git a/src/fusedgradedvector.jl b/src/fusedgradedvector.jl
index c3abc41c..95025c77 100644
--- a/src/fusedgradedvector.jl
+++ b/src/fusedgradedvector.jl
@@ -201,10 +201,10 @@ function FI.zero!(v::FusedGradedVector)
 end
 
 function Base.fill!(v::FusedGradedVector, val)
-    iszero(val) || throw(
-        ArgumentError("fill! with nonzero value is not supported for FusedGradedVector")
-    )
-    return FI.zero!(v)
+    for b in values(v.blocks)
+        fill!(b, val)
+    end
+    return v
 end
 
 # ========================  similar  ========================
diff --git a/src/fusion.jl b/src/fusion.jl
index 4b5b54f8..49154239 100644
--- a/src/fusion.jl
+++ b/src/fusion.jl
@@ -215,8 +215,15 @@ function TensorAlgebra.unmatricize(
         codomain_axes::Tuple{Vararg{GradedOneTo}},
         domain_axes::Tuple{Vararg{GradedOneTo}}
     )
-    isempty((codomain_axes..., domain_axes...)) &&
-        error("Scalar unmatricize not yet supported for FusedGradedMatrix")
+    if isempty(codomain_axes) && isempty(domain_axes)
+        # Scalar (rank-0) result: only the trivial-sector block contributes, and
+        # `cod ∩ dom = {trivial}` for a fused scalar means `m.blocks` is at most
+        # one 1×1 entry. Return a 0-D `Array` matching the eltype.
+        a = fill(zero(eltype(m)))
+        triv = trivial(sectortype(m))
+        haskey(m.blocks, triv) && (a[] = m.blocks[triv][1, 1])
+        return a
+    end
     # Compute the unfused 2D axes the N-D blocked axes produce, then their
     # sector-merged form (the axes matricize would have produced).
     unfused_axes = matricize_axes(BlockReshapeFusion(), m, codomain_axes, domain_axes)

From fd2009ec81af1feac19e85148a33ef57381112f1 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 2 Jun 2026 17:13:18 -0400
Subject: [PATCH 03/31] Implement `projectto!`, `Array`, and graded-axis
 `similar` disambiguators on `AbelianGradedArray`

`TensorAlgebra.projectto!(dest::AbelianGradedArray, src::AbstractArray)` copies the symmetry-allowed entries of `src` into `dest` block-by-block and zeroes anything else. Magnitude-blind by contract: any weight on forbidden blocks is dropped silently, leaving the tolerance check to `checked_projectto!`.

`Base.Array(a::AbelianGradedArray)` materializes the graded array into a dense `Array` by laying stored blocks into a zero-initialized destination. This is the densification primitive the default `checked_projectto!` reconstruct-compare check needs.

Two `similar(::AbstractArray|::StridedArray, ::Type{T}, ::Tuple{GradedOneTo, ...})` methods make `similar_map(prototype, T, codomain_axes, domain_axes)` route to an `AbelianGradedArray` allocation whenever the axes are graded, regardless of the prototype's storage type. Without these, dispatch is ambiguous against `BlockArrays`' `StridedArray`-specific overloads.

Co-authored-by: Claude <noreply@anthropic.com>
---
 src/abeliangradedarray.jl | 45 +++++++++++++++++++++++++++++++++++++++
 1 file changed, 45 insertions(+)

diff --git a/src/abeliangradedarray.jl b/src/abeliangradedarray.jl
index 19f7b8c4..2229be26 100644
--- a/src/abeliangradedarray.jl
+++ b/src/abeliangradedarray.jl
@@ -248,6 +248,23 @@ function Base.similar(
     D_N′ = isconcretetype(D_N) ? D_N : Array{T, N}
     return AbelianGradedArray{T, N, D_N′, S}(undef, axes)
 end
+
+# Allocate a graded array from any prototype when the requested axes are
+# `GradedOneTo`. The axes carry the structure; the prototype only fixes a
+# fallback datatype. Two overloads to resolve ambiguity with `BlockArrays`'
+# `StridedArray`-specific methods.
+function Base.similar(
+        ::AbstractArray, ::Type{T},
+        axes::Tuple{GradedOneTo{S}, Vararg{GradedOneTo{S}}}
+    ) where {T, S}
+    return AbelianGradedArray{T}(undef, axes)
+end
+function Base.similar(
+        ::StridedArray, ::Type{T},
+        axes::Tuple{GradedOneTo{S}, Vararg{GradedOneTo{S}}}
+    ) where {T, S}
+    return AbelianGradedArray{T}(undef, axes)
+end
 function Base.similar(
         a::AbstractGradedArray{T}, axes::Tuple{Vararg{GradedOneTo}}
     ) where {T}
@@ -330,6 +347,34 @@ function Base.fill!(a::AbelianGradedArray, v)
     return FI.zero!(a)
 end
 
+# Orthogonal projection of a dense source into the symmetry-allowed subspace.
+# Magnitude-blind: forbidden-block entries of `src` are dropped without inspection.
+# Use `TensorAlgebra.checked_projectto!` to verify the discarded weight is small.
+function TensorAlgebra.projectto!(dest::AbelianGradedArray, src::AbstractArray)
+    size(dest) == size(src) || throw(
+        DimensionMismatch(
+            "projectto!: dest has size $(size(dest)), src has size $(size(src))"
+        )
+    )
+    FI.zero!(dest)
+    for b in allowedblocks(axes(dest))
+        block_ranges = ntuple(d -> axes(dest, d)[Block(Int(Tuple(b)[d]))], ndims(dest))
+        view(dest, b) .= view(src, block_ranges...)
+    end
+    return dest
+end
+
+# Materialize the graded array into a dense `Array` for the default
+# `checked_projectto!`/`isapprox`-after path.
+function Base.Array(a::AbelianGradedArray{T, N}) where {T, N}
+    dest = zeros(T, size(a))
+    for bI in eachblockstoredindex(a)
+        block_ranges = ntuple(d -> axes(a, d)[Block(Int(Tuple(bI)[d]))], N)
+        copyto!(view(dest, block_ranges...), view(a, bI))
+    end
+    return dest
+end
+
 function LinearAlgebra.norm(a::AbelianGradedArray, p::Real = 2)
     if p == Inf
         isempty(a.blockdata) && return zero(float(real(eltype(a))))

From ab40c69610e2ac1ed668e30d78d7849cf0d973f9 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 2 Jun 2026 18:45:12 -0400
Subject: [PATCH 04/31] Block-aware methods on graded storage for the BP
 simple-update gauge path

- `LinearAlgebra.isdiag` on `AbelianGradedMatrix` and `FusedGradedMatrix`
- `Base.mapreduce` on `FusedGradedVector`
- `Base.map` on `FusedGradedVector`
- `MatrixAlgebraKit.diagonal` on `FusedGradedVector`
- Base matrix functions (`sqrt`, `exp`, `log`, ...) on `FusedGradedMatrix`
- `TensorAlgebra.checked_projectto!` on `AbelianGradedArray`

Co-authored-by: Claude <noreply@anthropic.com>
---
 src/abeliangradedarray.jl  | 23 +++++++++++++++++++++++
 src/abstractgradedarray.jl |  7 +++++++
 src/fusedgradedmatrix.jl   | 18 ++++++++++++++++++
 src/fusedgradedvector.jl   | 19 +++++++++++++++++++
 src/matrixalgebrakit.jl    | 10 ++++++++++
 5 files changed, 77 insertions(+)

diff --git a/src/abeliangradedarray.jl b/src/abeliangradedarray.jl
index 2229be26..6553d1a7 100644
--- a/src/abeliangradedarray.jl
+++ b/src/abeliangradedarray.jl
@@ -347,6 +347,16 @@ function Base.fill!(a::AbelianGradedArray, v)
     return FI.zero!(a)
 end
 
+# Block-aware diagonal check: block-diagonal (no off-diagonal stored blocks), and each
+# stored diagonal block is itself diagonal. Bypasses the generic scalar-indexing path.
+function LinearAlgebra.isdiag(A::AbelianGradedMatrix)
+    BlockSparseArrays.isblockdiagonal(A) || return false
+    for bI in eachblockstoredindex(A)
+        LinearAlgebra.isdiag(view(A, bI)) || return false
+    end
+    return true
+end
+
 # Orthogonal projection of a dense source into the symmetry-allowed subspace.
 # Magnitude-blind: forbidden-block entries of `src` are dropped without inspection.
 # Use `TensorAlgebra.checked_projectto!` to verify the discarded weight is small.
@@ -364,6 +374,19 @@ function TensorAlgebra.projectto!(dest::AbelianGradedArray, src::AbstractArray)
     return dest
 end
 
+# Compare via `Array(dest)` so the generic `isapprox(::AbstractArray, ::AbelianGradedArray)`
+# path doesn't fall back to a `src - dest` broadcast that scalar-indexes the block storage.
+function TensorAlgebra.checked_projectto!(
+        dest::AbelianGradedArray, src::AbstractArray;
+        atol::Real = 0,
+        rtol::Real = Base.rtoldefault(real(eltype(src)))
+    )
+    TensorAlgebra.projectto!(dest, src)
+    isapprox(src, Array(dest); atol, rtol) ||
+        throw(InexactError(:checked_projectto!, typeof(dest), src))
+    return dest
+end
+
 # Materialize the graded array into a dense `Array` for the default
 # `checked_projectto!`/`isapprox`-after path.
 function Base.Array(a::AbelianGradedArray{T, N}) where {T, N}
diff --git a/src/abstractgradedarray.jl b/src/abstractgradedarray.jl
index 0fb25471..4ff57ea9 100644
--- a/src/abstractgradedarray.jl
+++ b/src/abstractgradedarray.jl
@@ -15,6 +15,13 @@ function BlockSparseArrays.isblockdiagonal(A::AbstractGradedMatrix)
     return true
 end
 
+# Block-aware `LinearAlgebra.isdiag` for graded matrices. The generic
+# `LinearAlgebra.isdiag` would fall through to `_isbanded_impl`'s scalar-indexed
+# `iszero(view)` iteration, which throws on block storage. A graded matrix is
+# diagonal iff it is block-diagonal (no off-diagonal blocks stored) and each
+# stored block is itself diagonal — the latter checked block-by-block to avoid
+# materializing the whole matrix.
+
 # Scalar indexing is not supported for graded arrays.
 function Base.getindex(::AbstractGradedArray, ::Vararg{Int})
     return error(
diff --git a/src/fusedgradedmatrix.jl b/src/fusedgradedmatrix.jl
index 61fbb902..9b5324da 100644
--- a/src/fusedgradedmatrix.jl
+++ b/src/fusedgradedmatrix.jl
@@ -66,6 +66,24 @@ function FusedGradedMatrix(
     return FusedGradedMatrix{eltype(D), D, S}(codomain, domain, blocks)
 end
 
+# Block-diagonal by construction (one block per sector), so just check each block.
+LinearAlgebra.isdiag(A::FusedGradedMatrix) = all(LinearAlgebra.isdiag, A.blocks)
+
+# Block-diagonal by construction, so any matrix function `f(A) = blkdiag(f(blk_i))` for
+# each stored block — covers `sqrt`, `exp`, `log`, etc. Routes around the generic
+# `LinearAlgebra` impls that scalar-index for triangular / Hermitian detection.
+# Per-block result eltypes may differ (e.g. `sqrt(::Matrix{Float64})` returns
+# `Matrix{ComplexF64}` via Schur even when each block is real-PSD), so unify to the
+# `promote_type` of all returned blocks before reconstructing.
+for f in TensorAlgebra.MATRIX_FUNCTIONS
+    @eval function Base.$f(A::FusedGradedMatrix)
+        raw = map(Base.$f, A.blocks)
+        T = mapreduce(eltype, promote_type, raw)
+        blocks = map(b -> eltype(b) === T ? b : convert(AbstractMatrix{T}, b), raw)
+        return FusedGradedMatrix(A.codomain, A.domain, blocks)
+    end
+end
+
 """
     FusedGradedMatrix(sectors::Vector{S}, blocks::Vector{D})
 
diff --git a/src/fusedgradedvector.jl b/src/fusedgradedvector.jl
index 95025c77..549edfa1 100644
--- a/src/fusedgradedvector.jl
+++ b/src/fusedgradedvector.jl
@@ -168,6 +168,25 @@ Base.size(v::FusedGradedVector) = map(length, axes(v))
 Base.eltype(::Type{FusedGradedVector{T}}) where {T} = T
 Base.eltype(::Type{<:FusedGradedVector{T}}) where {T} = T
 
+# Block-wise `mapreduce`: reduce each block locally (so GPU blocks stay on the device for
+# their reduction kernel) and combine per-block scalars on the CPU. Routes
+# `maximum(abs, v; init=…)`, `sum`, `LinearAlgebra.norm`, etc. without ever falling
+# through to `getindex(v, ::Int)`. `init` (and any other kwargs) flow only to the outer
+# cross-block fold; double-applying `init` per block would be wrong for non-idempotent
+# reductions (e.g. `sum(v; init=10)`).
+function Base.mapreduce(f, op, v::FusedGradedVector; kwargs...)
+    return mapfoldl(b -> mapreduce(f, op, b), op, values(v.blocks); kwargs...)
+end
+
+# Block-wise `map`: returns a `FusedGradedVector` with the same axis and `f` applied to
+# each stored block, instead of falling through to `collect_similar` which would
+# allocate an `AbelianGradedVector` and scalar-setindex! into it. Each per-block `map`
+# dispatches to the storage backend's `map` (e.g. GPU kernel for `CuVector` blocks).
+function Base.map(f, v::FusedGradedVector)
+    blocks = dictionary(s => map(f, b) for (s, b) in pairs(v.blocks))
+    return FusedGradedVector(v.axis, blocks)
+end
+
 # ========================  Block indexing (primitive)  ========================
 
 function Base.view(v::FusedGradedVector, I::Block{1})
diff --git a/src/matrixalgebrakit.jl b/src/matrixalgebrakit.jl
index d8a89258..6850dc6d 100644
--- a/src/matrixalgebrakit.jl
+++ b/src/matrixalgebrakit.jl
@@ -293,6 +293,16 @@ function MAK.diagview(m::FusedGradedMatrix)
     return FusedGradedVector(diag_axis, diag_blocks)
 end
 
+# Inverse of `diagview`: wrap a FusedGradedVector as a block-diagonal FusedGradedMatrix
+# whose inner blocks are `Diagonal`. Keeps the structured form through
+# `pow_diag_safe(D) = MAK.diagonal(map(f, MAK.diagview(D)))`, so the next
+# `V * MAK.diagonal(...)` stays in the block-diagonal multiplication path instead of
+# falling through to LinearAlgebra's scalar-indexing `Diagonal*Matrix` impl.
+function MAK.diagonal(v::FusedGradedVector)
+    diag_blocks = map(Diagonal, v.blocks)
+    return FusedGradedMatrix(v.axis, v.axis, diag_blocks)
+end
+
 # Count how many elements are kept for a given index specification and block size
 _count_kept(::Colon, n) = n
 _count_kept(ind::AbstractVector{Bool}, _) = count(ind)

From ea95880b5e7c4a78802e8d8dea7e7330fc91fa44 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 2 Jun 2026 19:59:59 -0400
Subject: [PATCH 05/31] AbelianGradedMatrix alias, block-wise dot on
 FusedGradedMatrix

- Add `const AbelianGradedMatrix = AbelianGradedArray{T,2,D,S}` for the block-aware `LinearAlgebra.isdiag` dispatch that uses it (precompile was failing without the alias defined at parse time).
- Define `LinearAlgebra.dot(::FusedGradedMatrix, ::FusedGradedMatrix)` block-wise as the sum of `length(c) * dot(A.blocks[c], B.blocks[c])` over shared sectors. The `length(c)` factor matches the per-sector multiplicity convention in `LinearAlgebra.norm`, so `dot(A, A) == norm(A, 2)^2`.
---
 src/abeliangradedarray.jl |  2 ++
 src/fusedgradedmatrix.jl  | 13 +++++++++++++
 2 files changed, 15 insertions(+)

diff --git a/src/abeliangradedarray.jl b/src/abeliangradedarray.jl
index 6553d1a7..d529f80a 100644
--- a/src/abeliangradedarray.jl
+++ b/src/abeliangradedarray.jl
@@ -18,6 +18,8 @@ struct AbelianGradedArray{T, N, D <: AbstractArray{T, N}, S <: SectorRange} <:
     axes::NTuple{N, GradedOneTo{S}}
 end
 
+const AbelianGradedMatrix{T, D, S} = AbelianGradedArray{T, 2, D, S}
+
 # ---------------------------------------------------------------------------
 #  Constructors
 # ---------------------------------------------------------------------------
diff --git a/src/fusedgradedmatrix.jl b/src/fusedgradedmatrix.jl
index 9b5324da..090f2d78 100644
--- a/src/fusedgradedmatrix.jl
+++ b/src/fusedgradedmatrix.jl
@@ -265,6 +265,19 @@ LinearAlgebra.istriu(A::FusedGradedMatrix) = all(LinearAlgebra.istriu, values(A.
 LinearAlgebra.istril(A::FusedGradedMatrix) = all(LinearAlgebra.istril, values(A.blocks))
 LinearAlgebra.isposdef(A::FusedGradedMatrix) = all(LinearAlgebra.isposdef, values(A.blocks))
 
+function LinearAlgebra.dot(A::FusedGradedMatrix, B::FusedGradedMatrix)
+    A.codomain == B.codomain ||
+        throw(DimensionMismatch("dot: codomain mismatch"))
+    A.domain == B.domain ||
+        throw(DimensionMismatch("dot: domain mismatch"))
+    T = promote_type(eltype(A), eltype(B))
+    s = zero(T)
+    for c in intersect(keys(A.blocks), keys(B.blocks))
+        s += length(c) * LinearAlgebra.dot(A.blocks[c], B.blocks[c])
+    end
+    return s
+end
+
 # ========================  similar  ========================
 
 function Base.similar(m::FusedGradedMatrix, ::Type{T}) where {T}

From d983a8b3f079b17caff27ba3418e23ca7130cbdb Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Wed, 3 Jun 2026 10:16:44 -0400
Subject: [PATCH 06/31] Disambiguate N=1 block view/getindex/setindex! on
 graded vectors

For an `AbstractGradedArray{T, N}`, the splatted method `view(::AbstractGradedArray{T, N}, ::Vararg{Block{1}, N})` collapses to a single `Block{1}` argument when N=1, which is ambiguous with the single-`Block{N}` method. Add explicit `Block{1}`-only methods on `AbstractGradedArray{T, 1}` and `AbelianGradedArray{T, 1}` for `view`, `getindex`, and `setindex!` to resolve the dispatch.

Co-authored-by: Claude <noreply@anthropic.com>
---
 src/abeliangradedarray.jl  | 9 +++++++++
 src/abstractgradedarray.jl | 8 ++++++++
 2 files changed, 17 insertions(+)

diff --git a/src/abeliangradedarray.jl b/src/abeliangradedarray.jl
index d529f80a..590453e6 100644
--- a/src/abeliangradedarray.jl
+++ b/src/abeliangradedarray.jl
@@ -92,6 +92,15 @@ function Base.view(a::AbelianGradedArray{T, N}, I::Block{N}) where {T, N}
     return AbelianSectorArray(sects, a.blockdata[bk])
 end
 
+# Disambiguate against `view(::AbstractGradedArray{T, N}, ::Vararg{Block{1}, N})` for
+# N=1, where the splatted form collapses to a single Block{1} argument.
+function Base.view(a::AbelianGradedArray{T, 1}, I::Block{1}) where {T}
+    bk = (Int(I),)
+    haskey(a.blockdata, bk) || error("Block $bk is not stored.")
+    sects = (sectors(axes(a, 1))[bk[1]],)
+    return AbelianSectorArray(sects, a.blockdata[bk])
+end
+
 # ---------------------------------------------------------------------------
 #  blocks — lazy view delegating to view (following BlockArrays convention)
 # ---------------------------------------------------------------------------
diff --git a/src/abstractgradedarray.jl b/src/abstractgradedarray.jl
index 4ff57ea9..2387d407 100644
--- a/src/abstractgradedarray.jl
+++ b/src/abstractgradedarray.jl
@@ -46,6 +46,8 @@ end
 function Base.view(a::AbstractGradedArray{T, N}, I::Vararg{Block{1}, N}) where {T, N}
     return view(a, Block(Int.(I)))
 end
+# Disambiguate against subtype-specific `view(::ConcreteGradedVector, ::Block{1})` methods.
+Base.view(a::AbstractGradedArray{T, 1}, I::Block{1}) where {T} = view(a, Block((Int(I),)))
 
 function Base.getindex(a::AbstractGradedArray{T, N}, I::Block{N}) where {T, N}
     return copy(view(a, I))
@@ -53,6 +55,8 @@ end
 function Base.getindex(a::AbstractGradedArray{T, N}, I::Vararg{Block{1}, N}) where {T, N}
     return a[Block(Int.(I))]
 end
+# Disambiguate the N=1 case: route through the `Block{N}` method to avoid recursion.
+Base.getindex(a::AbstractGradedArray{T, 1}, I::Block{1}) where {T} = copy(view(a, I))
 
 function Base.setindex!(a::AbstractGradedArray{<:Any, N}, value, I::Block{N}) where {N}
     return setindex!(a, value, Tuple(I)...)
@@ -63,6 +67,10 @@ function Base.setindex!(
     view(a, I...) .= value
     return a
 end
+function Base.setindex!(a::AbstractGradedArray{<:Any, 1}, value, I::Block{1})
+    view(a, I) .= value
+    return a
+end
 
 # ---------------------------------------------------------------------------
 #  Data indexing — raw block data without sector wrappers

From 41a65946c2e366ce15aac64014f2a85021d33189 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Wed, 3 Jun 2026 10:17:04 -0400
Subject: [PATCH 07/31] AbelianGradedArray scale!, zero!, fill! via direct
 blockdata iteration

Iterate `values(a.blockdata)` directly instead of going through `eachblockstoredindex(a)` + `view`, which routed through the AbstractGradedArray-generic path and was the wrong shape for the new top-level array ops. `fill!` now actually fills with any value block-wise rather than rejecting nonzero arguments.

Co-authored-by: Claude <noreply@anthropic.com>
---
 src/abeliangradedarray.jl | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/src/abeliangradedarray.jl b/src/abeliangradedarray.jl
index 590453e6..9b25fb87 100644
--- a/src/abeliangradedarray.jl
+++ b/src/abeliangradedarray.jl
@@ -337,25 +337,25 @@ end
 # ---------------------------------------------------------------------------
 
 scale!(a::AbstractArray, β::Number) = (a .*= β; a)
-function scale!(a::AbstractGradedArray, β::Number)
-    for bI in eachblockstoredindex(a)
-        scale!(view(a, bI), β)
+function scale!(a::AbelianGradedArray, β::Number)
+    for b in values(a.blockdata)
+        scale!(b, β)
     end
     return a
 end
 
-function FI.zero!(a::AbstractGradedArray)
-    for bI in eachblockstoredindex(a)
-        FI.zero!(view(a, bI))
+function FI.zero!(a::AbelianGradedArray)
+    for b in values(a.blockdata)
+        FI.zero!(b)
     end
     return a
 end
 
 function Base.fill!(a::AbelianGradedArray, v)
-    iszero(v) || throw(
-        ArgumentError("fill! with nonzero value is not supported for AbelianGradedArray")
-    )
-    return FI.zero!(a)
+    for b in values(a.blockdata)
+        fill!(b, v)
+    end
+    return a
 end
 
 # Block-aware diagonal check: block-diagonal (no off-diagonal stored blocks), and each

From 8ae503ad672e277a87477d7cb9916bf504597d0a Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Wed, 3 Jun 2026 10:17:19 -0400
Subject: [PATCH 08/31] Block-wise + and - on FusedGradedMatrix

`FusedGradedMatrix` inherits `BroadcastStyle = GradedStyle{2}()` from `AbstractGradedArray`, but the generic `bipermutedimsopadd!(::AbstractGradedArray, ...)` path on which the broadcast machinery depends assumes the cartesian `(row_block, col_block)` block layout of `AbelianGradedArray`. On a `FusedGradedMatrix`, whose blocks live in a `Dictionary{Sector, Matrix}` keyed at coupled sectors, that path writes to the wrong slots and silently produces wrong results. `A - B` on the residual of an SVD reconstruction was the symptom that surfaced this.

Define `+` and `-` directly block-wise on the union of stored sectors, so a sector present on one side and absent on the other is treated as a zero block. The proper fix is to make broadcasting itself work on `FusedGradedMatrix` (either a per-method specialization of `bipermutedimsopadd!` or a dedicated `FusedGradedStyle`), tracked as a follow-up.

Co-authored-by: Claude <noreply@anthropic.com>
---
 src/fusedgradedmatrix.jl | 45 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 45 insertions(+)

diff --git a/src/fusedgradedmatrix.jl b/src/fusedgradedmatrix.jl
index 090f2d78..aa1688ef 100644
--- a/src/fusedgradedmatrix.jl
+++ b/src/fusedgradedmatrix.jl
@@ -238,6 +238,51 @@ function Base.:(*)(A::FusedGradedMatrix, B::FusedGradedMatrix)
     return mul!(C, A, B)
 end
 
+# ========================  Block-wise +, - ========================
+#
+# FusedGradedMatrix has no `BroadcastStyle`, so the AbstractArray fallback for `+`/`-`
+# (`broadcast_preserving_zero_d`) ends up trying scalar indexing on a sector-block
+# structure — which silently produces wrong results. Define the operations directly
+# block-by-block, taking the union of sectors so an absent block on one side is
+# treated as zero.
+
+function _check_add_axes(A::FusedGradedMatrix, B::FusedGradedMatrix, op)
+    A.codomain == B.codomain ||
+        throw(DimensionMismatch("$(op): codomain mismatch"))
+    A.domain == B.domain ||
+        throw(DimensionMismatch("$(op): domain mismatch"))
+    return nothing
+end
+
+function _block_combine(op, A::FusedGradedMatrix, B::FusedGradedMatrix)
+    T = promote_type(eltype(A), eltype(B))
+    DA = datatype(BlockSparseArrays.blocktype(A))
+    DB = datatype(BlockSparseArrays.blocktype(B))
+    D = Base.promote_op(op, DA, DB)
+    D′ = isconcretetype(D) ? D : Matrix{T}
+    S = sectortype(typeof(A))
+    sectors = sort!(collect(union(keys(A.blocks), keys(B.blocks))))
+    blocks = Dictionary{S, D′}()
+    for c in sectors
+        rows = get(A.codomain, c, get(B.codomain, c, 0))
+        cols = get(A.domain, c, get(B.domain, c, 0))
+        a = get(() -> zeros(eltype(A), rows, cols), A.blocks, c)
+        b = get(() -> zeros(eltype(B), rows, cols), B.blocks, c)
+        insert!(blocks, c, op(a, b))
+    end
+    return FusedGradedMatrix(A.codomain, A.domain, blocks)
+end
+
+function Base.:(+)(A::FusedGradedMatrix, B::FusedGradedMatrix)
+    _check_add_axes(A, B, :+)
+    return _block_combine(+, A, B)
+end
+
+function Base.:(-)(A::FusedGradedMatrix, B::FusedGradedMatrix)
+    _check_add_axes(A, B, :-)
+    return _block_combine(-, A, B)
+end
+
 # ======================== LinearAlgebra ======================
 
 function Base.adjoint(A::FusedGradedMatrix)

From 9606a28958fb33101ee4095b31983cc6923fb4e0 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Wed, 3 Jun 2026 10:18:16 -0400
Subject: [PATCH 09/31] Drop fully truncated coupled sectors from
 FusedGradedMatrix SVD truncate
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

`MAK.truncate(svd_trunc!, ::NTuple{3, FusedGradedMatrix}, ...)` previously kept every coupled sector in the bond, including those whose singular values were all truncated to zero. Downstream gate applications that fused against this bond then produced row/col block indices outside the destination's `allowedblocks` and failed deep in `unmatricize` as `Block (i, j, k) is not stored`.

Drop sectors with zero retained columns from the bond axis only. For U the row (codomain) axis is the input codomain and keeps its full sector set, and only the column (bond) axis shrinks. Analogously for Vᴴ. For S, both sides are the bond and shrink together. Matches the `truncate_space` convention used by `TensorMap` factorizations.

Co-authored-by: Claude <noreply@anthropic.com>
---
 src/matrixalgebrakit.jl | 58 ++++++++++++++++++++++++++++++++++-------
 1 file changed, 48 insertions(+), 10 deletions(-)

diff --git a/src/matrixalgebrakit.jl b/src/matrixalgebrakit.jl
index 6850dc6d..badadd0d 100644
--- a/src/matrixalgebrakit.jl
+++ b/src/matrixalgebrakit.jl
@@ -422,8 +422,13 @@ function MAK.findtruncated_svd(v::FusedGradedVector, strategy::MAK.TruncationInt
     ]
 end
 
-# truncate for FusedGradedMatrix: build reduced-dimension output, drop empty sectors
-# TODO: how do we handle discarde sectors while keeping U and V square in the sectors?
+# truncate for FusedGradedMatrix: build reduced-dimension output, dropping fully
+# truncated sectors from the bond side only. For U the row (codomain) axis is the
+# input codomain and must keep its full sector set, and only the column (domain)
+# axis shrinks. Analogously for Vᴴ. For S, both sides are the bond and shrink
+# together. Sectors whose singular values are all truncated to zero are dropped
+# entirely from the bond, matching the `truncate_space` convention used in
+# `TensorMap` factorizations.
 function MAK.truncate(
         ::typeof(MAK.svd_trunc!),
         (U, S, Vᴴ)::NTuple{3, FusedGradedMatrix},
@@ -431,14 +436,47 @@ function MAK.truncate(
     )
     sv = MAK.diagview(S)
     inds = MAK.findtruncated_svd(sv, strategy)
-    sectors = collect(keys(U.blocks))
-    U_blocks = [U.blocks[s][:, inds[i]] for (i, s) in enumerate(sectors)]
-    S_blocks =
-        [Diagonal(MAK.diagview(S.blocks[s])[inds[i]]) for (i, s) in enumerate(sectors)]
-    Vᴴ_blocks = [Vᴴ.blocks[s][inds[i], :] for (i, s) in enumerate(sectors)]
-    Ũ = FusedGradedMatrix(sectors, U_blocks)
-    S̃ = FusedGradedMatrix(sectors, S_blocks)
-    Ṽᴴ = FusedGradedMatrix(sectors, Vᴴ_blocks)
+    sectors_all = collect(keys(U.blocks))
+
+    # Slice every sector's blocks first. `inds[i]` may be `Colon()` (notrunc) or a
+    # `Vector{Int}` (rank/tol/error truncations), so check emptiness via the resulting
+    # column count rather than `isempty(inds[i])`.
+    U_blocks_all =
+        [U.blocks[sectors_all[i]][:, inds[i]] for i in eachindex(inds)]
+    S_blocks_all = [
+        Diagonal(MAK.diagview(S.blocks[sectors_all[i]])[inds[i]])
+            for i in eachindex(inds)
+    ]
+    Vᴴ_blocks_all =
+        [Vᴴ.blocks[sectors_all[i]][inds[i], :] for i in eachindex(inds)]
+
+    keep = [i for i in eachindex(inds) if size(U_blocks_all[i], 2) > 0]
+    sectors_kept = sectors_all[keep]
+    bond_dims = [size(U_blocks_all[i], 2) for i in keep]
+
+    # U: rows = input codomain (full), cols = bond (shrunk).
+    U_cod = U.codomain
+    U_dom = Dictionary{eltype(sectors_kept), Int}(sectors_kept, bond_dims)
+    U_blks = Dictionary{eltype(sectors_kept), eltype(typeof(U.blocks))}(
+        sectors_kept, U_blocks_all[keep]
+    )
+    Ũ = FusedGradedMatrix(U_cod, U_dom, U_blks)
+
+    # S: both sides are the bond (shrunk).
+    S_side = Dictionary{eltype(sectors_kept), Int}(sectors_kept, bond_dims)
+    S_blks = Dictionary{eltype(sectors_kept), eltype(typeof(S.blocks))}(
+        sectors_kept, S_blocks_all[keep]
+    )
+    S̃ = FusedGradedMatrix(S_side, S_side, S_blks)
+
+    # Vᴴ: rows = bond (shrunk), cols = input domain (full).
+    Vᴴ_cod = Dictionary{eltype(sectors_kept), Int}(sectors_kept, bond_dims)
+    Vᴴ_dom = Vᴴ.domain
+    Vᴴ_blks = Dictionary{eltype(sectors_kept), eltype(typeof(Vᴴ.blocks))}(
+        sectors_kept, Vᴴ_blocks_all[keep]
+    )
+    Ṽᴴ = FusedGradedMatrix(Vᴴ_cod, Vᴴ_dom, Vᴴ_blks)
+
     return (Ũ, S̃, Ṽᴴ), inds
 end
 

From 966b5e29a65bfa815a1b9a08f603f47454d5a459 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Wed, 3 Jun 2026 11:32:59 -0400
Subject: [PATCH 10/31] Convert `import` to `using` for module-as-name imports

Match the project convention of `using Module: Module [as Alias]` over `import Module [as Alias]`. No behavior change.

Co-authored-by: Claude <noreply@anthropic.com>
---
 src/GradedArrays.jl     | 2 +-
 src/matrixalgebrakit.jl | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/GradedArrays.jl b/src/GradedArrays.jl
index 51e30dd9..6e5129f4 100644
--- a/src/GradedArrays.jl
+++ b/src/GradedArrays.jl
@@ -22,13 +22,13 @@ export dual, flip, gradedrange, isdual,
 
 # imports
 # -------
-import FunctionImplementations as FI
 using BlockArrays: BlockArrays, AbstractBlockArray, AbstractBlockVector,
     AbstractBlockedUnitRange, Block, BlockIndexRange, BlockVector, BlockedOneTo,
     blockedrange, blocklasts, blocklength, blocklengths, blocks, eachblockaxes1
 using BlockSparseArrays: BlockSparseArrays, blockdiagindices, blockstoredlength,
     eachblockaxis, eachblockstoredindex, mortar_axis
 using Dictionaries: Dictionaries, Dictionary, dictionary, gettoken, gettokenvalue
+using FunctionImplementations: FunctionImplementations as FI
 using KroneckerArrays: KroneckerArrays, kroneckerfactors, ×, ⊗
 using LinearAlgebra: LinearAlgebra, Adjoint, Diagonal, mul!
 using SparseArraysBase: SparseArraysBase
diff --git a/src/matrixalgebrakit.jl b/src/matrixalgebrakit.jl
index badadd0d..b32d90f7 100644
--- a/src/matrixalgebrakit.jl
+++ b/src/matrixalgebrakit.jl
@@ -1,4 +1,4 @@
-import MatrixAlgebraKit as MAK
+using MatrixAlgebraKit: MatrixAlgebraKit as MAK
 
 struct GradedBlockAlgorithm{A <: MAK.AbstractAlgorithm} <: MAK.AbstractAlgorithm
     alg::A

From ca41b4c7cd68c3a281d06dfe67dd3bf52978a0d3 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Wed, 3 Jun 2026 11:41:09 -0400
Subject: [PATCH 11/31] Pin TensorAlgebra to mf/gram-eigh-balanced via
 `[sources]`

GradedArrays now references `TensorAlgebra.projectto!`, added in https://github.com/ITensor/TensorAlgebra.jl/pull/177. Without the pin, CI resolves against the registered TensorAlgebra 0.9.4 which does not have `projectto!`, and precompilation fails. Drop the pin once that branch registers.

Co-authored-by: Claude <noreply@anthropic.com>
---
 Project.toml | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/Project.toml b/Project.toml
index e71bc40c..0b726297 100644
--- a/Project.toml
+++ b/Project.toml
@@ -29,6 +29,10 @@ TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138"
 NamedDimsArrays = "60cbd0c0-df58-4cb7-918c-6f5607b73fde"
 SUNRepresentations = "1a50b95c-7aac-476d-a9ce-2bfc675fc617"
 
+[sources.TensorAlgebra]
+rev = "mf/gram-eigh-balanced"
+url = "https://github.com/ITensor/TensorAlgebra.jl"
+
 [extensions]
 GradedArraysNamedDimsArraysExt = "NamedDimsArrays"
 GradedArraysSUNRepresentationsExt = "SUNRepresentations"

From f27e1a1116f6806eddfa6186c110519173a4eae8 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Wed, 3 Jun 2026 11:53:54 -0400
Subject: [PATCH 12/31] Block-wise `Base.copyto!` on `AbelianGradedArray`
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The default `copyto!(::AbstractArray, ::AbstractArray)` fallback uses scalar indexing, which isn't supported on `AbstractGradedArray`. With this specialization, an `AbelianGradedArray → AbelianGradedArray` copy goes through the block dictionary directly.

Surfaces in broadcasts that collapse to the bare source array, e.g. `conj.(a)` on a `Real`-eltype `AbelianGradedArray`: `tryflattenlinear` returns the source itself (since `conj` is a no-op on reals), and the subsequent `copyto!(dest, src)` then needed this method.

Co-authored-by: Claude <noreply@anthropic.com>
---
 src/abeliangradedarray.jl | 17 +++++++++++++++++
 1 file changed, 17 insertions(+)

diff --git a/src/abeliangradedarray.jl b/src/abeliangradedarray.jl
index 9b25fb87..aa03f0ac 100644
--- a/src/abeliangradedarray.jl
+++ b/src/abeliangradedarray.jl
@@ -296,6 +296,23 @@ function Base.copy(a::AbelianGradedArray{T, N, D, S}) where {T, N, D, S}
     )
 end
 
+function Base.copyto!(
+        dest::AbelianGradedArray{<:Any, N},
+        src::AbelianGradedArray{<:Any, N}
+    ) where {N}
+    axes(dest) == axes(src) ||
+        throw(
+        DimensionMismatch(
+            "copyto! axes mismatch: dest $(axes(dest)), src $(axes(src))"
+        )
+    )
+    empty!(dest.blockdata)
+    for (k, v) in src.blockdata
+        dest.blockdata[k] = copy(v)
+    end
+    return dest
+end
+
 # Conjugate element-wise *and* flip axis duality, mirroring `Base.conj` on the
 # axis types (`SectorRange`/`GradedOneTo`/`SectorOneTo`). Without the axis flip,
 # `conj(t)` would leave bra-layer tensors with the same duality as the ket, and

From 336ef0f975554afe36f12479e5f1d312c69c696a Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Wed, 3 Jun 2026 12:48:08 -0400
Subject: [PATCH 13/31] Loud-error broadcasting on FusedGradedMatrix

Reject all broadcasts on `FusedGradedMatrix` with an explicit error via a new `FusedGradedStyle`. The generic `GradedStyle` path lowers to `bipermutedimsopadd!` over a cartesian block layout, but `FusedGradedMatrix` stores blocks keyed by coupled sector, so the result was silently wrong. The error points at `+`/`-` or block-wise ops as the supported alternatives until a sector-keyed broadcast path lands.

Add block-wise `Base.:*(::FusedGradedMatrix, ::Number)`, the reverse, and `Base.:/` so the common scalar-multiply cases keep working without going through broadcast.

Co-authored-by: Claude <noreply@anthropic.com>
---
 src/broadcast.jl           | 34 ++++++++++++++++++++++++++++++++++
 src/fusedgradedmatrix.jl   | 10 ++++++++++
 test/test_tensoralgebra.jl | 13 +++++++++----
 3 files changed, 53 insertions(+), 4 deletions(-)

diff --git a/src/broadcast.jl b/src/broadcast.jl
index eee2be78..8863eea3 100644
--- a/src/broadcast.jl
+++ b/src/broadcast.jl
@@ -88,3 +88,37 @@ function Base.copyto!(dest::AbstractGradedArray, bc::BC.Broadcasted{<:GradedStyl
         throw(ArgumentError("AbstractGradedArray broadcasting requires linear operations"))
     return copyto!(dest, lb)
 end
+
+# ====================  FusedGradedMatrix broadcasting  ====================
+#
+# `FusedGradedMatrix` stores blocks keyed by coupled sector rather than by
+# cartesian block index, so the `GradedStyle` path (which lowers to
+# `bipermutedimsopadd!` over cartesian block storage) silently produces wrong
+# results. Until a proper sector-keyed broadcast path lands, FGM broadcasting
+# is an explicit error; use `Base.:(+)` / `Base.:(-)` or block-wise operations
+# instead.
+
+struct FusedGradedStyle <: BC.AbstractArrayStyle{2} end
+FusedGradedStyle(::Val{N}) where {N} = FusedGradedStyle()
+
+BC.BroadcastStyle(::Type{<:FusedGradedMatrix}) = FusedGradedStyle()
+BC.BroadcastStyle(s::FusedGradedStyle, ::BC.DefaultArrayStyle{0}) = s
+BC.BroadcastStyle(::BC.DefaultArrayStyle{0}, s::FusedGradedStyle) = s
+BC.BroadcastStyle(s::FusedGradedStyle, ::FusedGradedStyle) = s
+
+function Base.copy(::BC.Broadcasted{FusedGradedStyle})
+    return throw(
+        ArgumentError(
+            "Broadcasting on `FusedGradedMatrix` is not supported; use `+`/`-` " *
+                "or operate block-wise via `blocks(A)` instead."
+        )
+    )
+end
+function Base.copyto!(::AbstractArray, ::BC.Broadcasted{FusedGradedStyle})
+    return throw(
+        ArgumentError(
+            "Broadcasting on `FusedGradedMatrix` is not supported; use `+`/`-` " *
+                "or operate block-wise via `blocks(A)` instead."
+        )
+    )
+end
diff --git a/src/fusedgradedmatrix.jl b/src/fusedgradedmatrix.jl
index aa1688ef..e4dbef11 100644
--- a/src/fusedgradedmatrix.jl
+++ b/src/fusedgradedmatrix.jl
@@ -283,6 +283,16 @@ function Base.:(-)(A::FusedGradedMatrix, B::FusedGradedMatrix)
     return _block_combine(-, A, B)
 end
 
+function Base.:(*)(A::FusedGradedMatrix, x::Number)
+    new_blocks = map(b -> b * x, A.blocks)
+    return FusedGradedMatrix(A.codomain, A.domain, new_blocks)
+end
+Base.:(*)(x::Number, A::FusedGradedMatrix) = A * x
+function Base.:(/)(A::FusedGradedMatrix, x::Number)
+    new_blocks = map(b -> b / x, A.blocks)
+    return FusedGradedMatrix(A.codomain, A.domain, new_blocks)
+end
+
 # ======================== LinearAlgebra ======================
 
 function Base.adjoint(A::FusedGradedMatrix)
diff --git a/test/test_tensoralgebra.jl b/test/test_tensoralgebra.jl
index 84eaf7be..cdcc8bad 100644
--- a/test/test_tensoralgebra.jl
+++ b/test/test_tensoralgebra.jl
@@ -383,7 +383,7 @@ end
     @test all(iszero, data(c[Block(2, 2)]))
 end
 
-@testset "FusedGradedMatrix broadcasting" begin
+@testset "FusedGradedMatrix block-wise arithmetic" begin
     m = FusedGradedMatrix([U1(0), U1(1)], [[1.0 2.0; 3.0 4.0], [5.0 0.0; 0.0 6.0]])
 
     m2 = 3 * m
@@ -394,9 +394,14 @@ end
     s = m + n
     @test data(s[Block(1, 1)]) == [2.0 3.0; 4.0 5.0]
     @test data(s[Block(2, 2)]) == [6.0 1.0; 1.0 7.0]
+end
 
+@testset "FusedGradedMatrix broadcasting errors" begin
+    m = FusedGradedMatrix([U1(0), U1(1)], [[1.0 2.0; 3.0 4.0], [5.0 0.0; 0.0 6.0]])
+    n = FusedGradedMatrix([U1(0), U1(1)], [ones(2, 2), ones(2, 2)])
+    @test_throws ArgumentError m .+ n
+    @test_throws ArgumentError 3 .* m
+    @test_throws ArgumentError conj.(m)
     c = similar(m, Float64)
-    c .= 3 .* m .+ 2 .* n
-    @test data(c[Block(1, 1)]) == [5.0 8.0; 11.0 14.0]
-    @test data(c[Block(2, 2)]) == [17.0 2.0; 2.0 20.0]
+    @test_throws ArgumentError c .= 3 .* m .+ 2 .* n
 end

From 9923935ba096cab38d5ef92e281f2c97c8a703ef Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Wed, 3 Jun 2026 17:26:02 -0400
Subject: [PATCH 14/31] Overload `TensorAlgebra.trivialrange` for `GradedOneTo`

Overloads `TensorAlgebra.trivialrange(::Type{<:GradedOneTo})` to return a charge-0 one-dimensional graded range, the identity for `tensor_product_axis` on graded ranges, by delegating to the existing `trivial(R)` method.

Co-authored-by: Claude <noreply@anthropic.com>
---
 src/gradedoneto.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/gradedoneto.jl b/src/gradedoneto.jl
index 85ed4c33..133a7ba5 100644
--- a/src/gradedoneto.jl
+++ b/src/gradedoneto.jl
@@ -55,6 +55,8 @@ function trivial(::Type{GradedOneTo{S}}) where {S}
 end
 trivial(g::GradedOneTo) = trivial(typeof(g))
 
+TensorAlgebra.trivialrange(R::Type{<:GradedOneTo}) = trivial(R)
+
 """
     gradedrange(xs::AbstractVector{<:Pair})
 

From 0cdf49902cdb298f8573ceb67be3941dbd2c440f Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Wed, 3 Jun 2026 18:03:22 -0400
Subject: [PATCH 15/31] Block-aware `Random.rand!` and `Random.randn!` on
 `AbelianGradedArray`

Dispatch to each block's `rand!` / `randn!`, alongside the existing block-aware `fill!`, `zero!`, and `scale!`. The generic `AbstractArray` fallbacks go through scalar indexing, which errors on graded arrays.

Co-authored-by: Claude <noreply@anthropic.com>
---
 src/GradedArrays.jl       |  1 +
 src/abeliangradedarray.jl | 15 +++++++++++++++
 2 files changed, 16 insertions(+)

diff --git a/src/GradedArrays.jl b/src/GradedArrays.jl
index 6e5129f4..bb5acade 100644
--- a/src/GradedArrays.jl
+++ b/src/GradedArrays.jl
@@ -31,6 +31,7 @@ using Dictionaries: Dictionaries, Dictionary, dictionary, gettoken, gettokenvalu
 using FunctionImplementations: FunctionImplementations as FI
 using KroneckerArrays: KroneckerArrays, kroneckerfactors, ×, ⊗
 using LinearAlgebra: LinearAlgebra, Adjoint, Diagonal, mul!
+using Random: Random, AbstractRNG
 using SparseArraysBase: SparseArraysBase
 using TensorAlgebra: TensorAlgebra, BlockedTuple, FusionStyle, bipermutedimsopadd!,
     check_input, matricize, matricize_axes, permutedimsadd!, tensor_product_axis,
diff --git a/src/abeliangradedarray.jl b/src/abeliangradedarray.jl
index aa03f0ac..90b5e846 100644
--- a/src/abeliangradedarray.jl
+++ b/src/abeliangradedarray.jl
@@ -375,6 +375,21 @@ function Base.fill!(a::AbelianGradedArray, v)
     return a
 end
 
+# Block-aware random fills: dispatch to the underlying block's `rand!`/`randn!`,
+# bypassing the generic `AbstractArray` fallbacks that go through scalar indexing.
+function Random.rand!(rng::AbstractRNG, a::AbelianGradedArray)
+    for b in values(a.blockdata)
+        Random.rand!(rng, b)
+    end
+    return a
+end
+function Random.randn!(rng::AbstractRNG, a::AbelianGradedArray)
+    for b in values(a.blockdata)
+        Random.randn!(rng, b)
+    end
+    return a
+end
+
 # Block-aware diagonal check: block-diagonal (no off-diagonal stored blocks), and each
 # stored diagonal block is itself diagonal. Bypasses the generic scalar-indexing path.
 function LinearAlgebra.isdiag(A::AbelianGradedMatrix)

From 00ddc6b099dc6028a87bfa33a70c0810084a93ef Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Wed, 3 Jun 2026 19:36:11 -0400
Subject: [PATCH 16/31] Update fill! test on AbelianGradedArray to match the
 implementation

Commit 41a65946c2e3 changed `Base.fill!(::AbelianGradedArray, val)` to fill stored blocks block-wise with any value rather than throwing on nonzero. The test still expected an `ArgumentError`.

Co-authored-by: Claude <noreply@anthropic.com>
---
 test/test_abelianarray.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/test/test_abelianarray.jl b/test/test_abelianarray.jl
index 8b7f44b0..64ae6656 100644
--- a/test/test_abelianarray.jl
+++ b/test/test_abelianarray.jl
@@ -205,8 +205,9 @@ using Test: @test, @test_throws, @testset
         @test !isempty(a.blockdata)
         @test all(iszero, a.blockdata[(1, 1)])
 
-        # fill! with nonzero errors
-        @test_throws ArgumentError fill!(a, 1.0)
+        # fill! fills stored blocks block-wise with any value
+        fill!(a, 1.0)
+        @test all(==(1.0), a.blockdata[(1, 1)])
 
         # zero! zeros stored blocks in place (blocks stay allocated)
         a[Block(1, 1)] = AbelianSectorArray((U1(0), dual(U1(0))), ones(2, 2))

From 1badce4051cc007e8ff892b2d757d66bd398c1e4 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Wed, 3 Jun 2026 20:40:44 -0400
Subject: [PATCH 17/31] Allow empty GradedOneTo through mortar_axis, sqrt, and
 equality

`mortar_axis(::AbstractVector{SectorOneTo{S}})` and `mortar_axis(::AbstractVector{GradedOneTo{S}})` now return `GradedOneTo(S[], Int[])` on empty input. The signatures are parameterized on the sector type `S` so it can be recovered from the eltype.

`sqrt(::FusedGradedMatrix)` (and the rest of the `MATRIX_FUNCTIONS` loop) passes `init = eltype(A)` to the eltype-promoting `mapreduce`, so the reduction works when the block list is empty.

`Base.isequal(::GradedOneTo, ::GradedOneTo)` and `Base.hash` treat two empty graded ranges as equal regardless of `isdual`, since with no sectors there is nothing to flip.

These come up in a U(1) gate-application QR whose codomain index list is empty, producing an empty bond axis that flows through `tensor_product_gradedrange`, `unmatricize`, and `check_input(*, FGM, FGM)`.
---
 src/fusedgradedmatrix.jl |  2 +-
 src/gradedoneto.jl       | 22 +++++++++++-----------
 2 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/src/fusedgradedmatrix.jl b/src/fusedgradedmatrix.jl
index e4dbef11..1ee72f63 100644
--- a/src/fusedgradedmatrix.jl
+++ b/src/fusedgradedmatrix.jl
@@ -78,7 +78,7 @@ LinearAlgebra.isdiag(A::FusedGradedMatrix) = all(LinearAlgebra.isdiag, A.blocks)
 for f in TensorAlgebra.MATRIX_FUNCTIONS
     @eval function Base.$f(A::FusedGradedMatrix)
         raw = map(Base.$f, A.blocks)
-        T = mapreduce(eltype, promote_type, raw)
+        T = mapreduce(eltype, promote_type, raw; init = eltype(A))
         blocks = map(b -> eltype(b) === T ? b : convert(AbstractMatrix{T}, b), raw)
         return FusedGradedMatrix(A.codomain, A.domain, blocks)
     end
diff --git a/src/gradedoneto.jl b/src/gradedoneto.jl
index 133a7ba5..6b78ad83 100644
--- a/src/gradedoneto.jl
+++ b/src/gradedoneto.jl
@@ -82,25 +82,21 @@ end
 eachdataaxis(g::GradedOneTo) = data.(eachblockaxis(g))
 eachsectoraxis(g::GradedOneTo) = sector.(eachblockaxis(g))
 
-function BlockSparseArrays.mortar_axis(axs::AbstractVector{<:SectorOneTo})
-    isempty(axs) && throw(
-        ArgumentError("Cannot create GradedOneTo from empty vector without type info")
-    )
+function BlockSparseArrays.mortar_axis(axs::AbstractVector{SectorOneTo{S}}) where {S}
+    isempty(axs) && return GradedOneTo(S[], Int[])
     allequal(isdual, axs) ||
         throw(ArgumentError("Cannot combine sectors with different arrows"))
     d = isdual(first(axs))
     # Store non-dual sectors; apply isdual via dual() if needed
-    ss = [d ? dual(sector(r)) : sector(r) for r in axs]
-    ms = [datalength(r) for r in axs]
+    ss = S[d ? dual(sector(r)) : sector(r) for r in axs]
+    ms = Int[datalength(r) for r in axs]
     g = GradedOneTo(ss, ms)
     return d ? dual(g) : g
 end
 
 # Non-abelian fusion: flatten GradedOneTo elements into a single GradedOneTo
-function BlockSparseArrays.mortar_axis(axs::AbstractVector{<:GradedOneTo})
-    isempty(axs) && throw(
-        ArgumentError("Cannot create GradedOneTo from empty vector without type info")
-    )
+function BlockSparseArrays.mortar_axis(axs::AbstractVector{GradedOneTo{S}}) where {S}
+    isempty(axs) && return GradedOneTo(S[], Int[])
     return mortar_axis(mapreduce(eachblockaxis, vcat, axs))
 end
 
@@ -169,14 +165,18 @@ function Base.checkindex(::Type{Bool}, g::GradedOneTo, i::Int)
     return 1 <= i <= length(g)
 end
 
-# Equality and hashing
+# Equality and hashing. Empty graded ranges have no sectors and the `isdual`
+# flag has no observable effect, so any two empty graded ranges of the same
+# sector type compare equal.
 function Base.isequal(a::GradedOneTo, b::GradedOneTo)
+    isempty(a.nondual_sectors) && isempty(b.nondual_sectors) && return true
     return isequal(a.nondual_sectors, b.nondual_sectors) &&
         isequal(datalengths(a), datalengths(b)) &&
         isequal(isdual(a), isdual(b))
 end
 Base.:(==)(a::GradedOneTo, b::GradedOneTo) = isequal(a, b)
 function Base.hash(g::GradedOneTo, h::UInt)
+    isempty(g.nondual_sectors) && return hash(GradedOneTo, h)
     return hash(g.nondual_sectors, hash(datalengths(g), hash(isdual(g), h)))
 end
 

From a727126effa13d8ca1b4b7e27da6c4b9dc62efa1 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Wed, 3 Jun 2026 21:43:39 -0400
Subject: [PATCH 18/31] Check contracted axes are duals of each other in
 bosonic graded contract

`contractopadd!(::Matricize{SectorFusion}, ...)` now verifies that every contracted-axis pair satisfies `dual(axes(a1)[i]) == axes(a2)[j]`. A same-`isdual` pairing (or a sector-content mismatch) was previously aligned block-wise without complaint, producing a number that did not correspond to the intended contraction.

Skipped when the sector's `BraidingStyle` is fermionic, since the supertrace formalism intentionally pairs same-`isdual` axes together with `contraction_twist!` to pick up the right phase.
---
 src/tensoralgebra.jl | 24 ++++++++++++++++++++++++
 1 file changed, 24 insertions(+)

diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl
index c776b3c2..24b6da84 100644
--- a/src/tensoralgebra.jl
+++ b/src/tensoralgebra.jl
@@ -217,6 +217,29 @@ function contraction_twist!(a::AbstractArray, ndims_codomain::Int)
     return twist!(a, (i for i in 1:ndims_codomain if isdual(a, i)))
 end
 
+# For every contracted axis pair, require `dual(axes(a1)[i]) == axes(a2)[j]`.
+# A same-`isdual` pairing (or a sector-content mismatch) means the two ends do not
+# agree on what is being summed, and the graded contract code would otherwise
+# align blocks index-wise and produce a result that does not correspond to the
+# intended contraction.
+#
+# Skipped for fermionic braiding: the supertrace formalism intentionally uses
+# same-duality pairings together with `contraction_twist!` to pick up the right
+# fermion signs, so flagging those would be wrong.
+function check_contracted_axes_dual(a1, perm1_domain, a2, perm2_codomain)
+    TKS.BraidingStyle(sectortype(a1)) isa TKS.Bosonic || return nothing
+    for (i, j) in zip(perm1_domain, perm2_codomain)
+        ax1 = axes(a1)[i]
+        ax2 = axes(a2)[j]
+        dual(ax1) == ax2 || throw(
+            ArgumentError(
+                "Contracted axes do not match: `axes(a1)[$i] = $ax1` and `axes(a2)[$j] = $ax2` are not duals of each other"
+            )
+        )
+    end
+    return nothing
+end
+
 #=
 This is an overload that follows the standard TensorAlgebra implementation,
 with the single exception of inserting a `contraction_twist!` call before the matricization.
@@ -238,6 +261,7 @@ function TensorAlgebra.contractopadd!(
         a1, perm1_codomain, perm1_domain,
         a2, perm2_codomain, perm2_domain
     )
+    check_contracted_axes_dual(a1, perm1_domain, a2, perm2_codomain)
 
     a1_mat = TensorAlgebra.matricizeop(
         algorithm.fusion_style,

From 00b782d8690f8dabb0df1f180f16333388ca7aa6 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Wed, 3 Jun 2026 21:51:28 -0400
Subject: [PATCH 19/31] Refactor dual-orientation check as a
 check_input(contract) overload

Same behavior as the previous commit's freestanding helper, but routed through `TensorAlgebra.check_input(::typeof(contract), ::AbstractGradedArray, ...)` so it follows the existing GradedArrays convention for graded-specific input validation. `check_input(contract!)` already delegates to `check_input(contract)`, so the explicit call in `contractopadd!` is no longer needed.
---
 src/tensoralgebra.jl | 12 ++++++++++--
 1 file changed, 10 insertions(+), 2 deletions(-)

diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl
index 24b6da84..4bdc0532 100644
--- a/src/tensoralgebra.jl
+++ b/src/tensoralgebra.jl
@@ -226,7 +226,16 @@ end
 # Skipped for fermionic braiding: the supertrace formalism intentionally uses
 # same-duality pairings together with `contraction_twist!` to pick up the right
 # fermion signs, so flagging those would be wrong.
-function check_contracted_axes_dual(a1, perm1_domain, a2, perm2_codomain)
+function TensorAlgebra.check_input(
+        f::typeof(TensorAlgebra.contract),
+        a1::AbstractGradedArray, perm1_codomain, perm1_domain,
+        a2::AbstractGradedArray, perm2_codomain, perm2_domain
+    )
+    @invoke TensorAlgebra.check_input(
+        f,
+        a1::AbstractArray, perm1_codomain, perm1_domain,
+        a2::AbstractArray, perm2_codomain, perm2_domain
+    )
     TKS.BraidingStyle(sectortype(a1)) isa TKS.Bosonic || return nothing
     for (i, j) in zip(perm1_domain, perm2_codomain)
         ax1 = axes(a1)[i]
@@ -261,7 +270,6 @@ function TensorAlgebra.contractopadd!(
         a1, perm1_codomain, perm1_domain,
         a2, perm2_codomain, perm2_domain
     )
-    check_contracted_axes_dual(a1, perm1_domain, a2, perm2_codomain)
 
     a1_mat = TensorAlgebra.matricizeop(
         algorithm.fusion_style,

From 0db5da00875333ccbb4e7f8608665cb0cb6e59ba Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Wed, 3 Jun 2026 21:57:47 -0400
Subject: [PATCH 20/31] Refine contracted-axis check to validate sector content
 for fermionic too

Factor the per-axis-pair logic into `axes_match_for_contraction(ax1, ax2)` and call it from `check_input(::typeof(contract), ::AbstractGradedArray, ...)`. The helper always requires sector content to match modulo `isdual`. For bosonic braiding it additionally requires opposite `isdual`, while for fermionic braiding `ax1 == ax2` is also accepted to leave room for the supertrace formalism's same-`isdual` pairings.
---
 src/tensoralgebra.jl | 25 +++++++++++++------------
 1 file changed, 13 insertions(+), 12 deletions(-)

diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl
index 4bdc0532..982d350f 100644
--- a/src/tensoralgebra.jl
+++ b/src/tensoralgebra.jl
@@ -217,15 +217,17 @@ function contraction_twist!(a::AbstractArray, ndims_codomain::Int)
     return twist!(a, (i for i in 1:ndims_codomain if isdual(a, i)))
 end
 
-# For every contracted axis pair, require `dual(axes(a1)[i]) == axes(a2)[j]`.
-# A same-`isdual` pairing (or a sector-content mismatch) means the two ends do not
-# agree on what is being summed, and the graded contract code would otherwise
-# align blocks index-wise and produce a result that does not correspond to the
-# intended contraction.
-#
-# Skipped for fermionic braiding: the supertrace formalism intentionally uses
-# same-duality pairings together with `contraction_twist!` to pick up the right
-# fermion signs, so flagging those would be wrong.
+# True iff `ax1` and `ax2` are valid as partners on a contracted axis. Always
+# requires sector content to match modulo the `isdual` flag (`dual(ax1) == ax2`
+# is the canonical non-dual ↔ dual pairing). For fermionic braiding `ax1 == ax2`
+# is also accepted, since the supertrace formalism uses `contraction_twist!`
+# to pick up the right phase on same-`isdual` pairings.
+function axes_match_for_contraction(ax1, ax2)
+    dual(ax1) == ax2 && return true
+    TKS.BraidingStyle(sectortype(typeof(ax1))) isa TKS.Bosonic && return false
+    return ax1 == ax2
+end
+
 function TensorAlgebra.check_input(
         f::typeof(TensorAlgebra.contract),
         a1::AbstractGradedArray, perm1_codomain, perm1_domain,
@@ -236,13 +238,12 @@ function TensorAlgebra.check_input(
         a1::AbstractArray, perm1_codomain, perm1_domain,
         a2::AbstractArray, perm2_codomain, perm2_domain
     )
-    TKS.BraidingStyle(sectortype(a1)) isa TKS.Bosonic || return nothing
     for (i, j) in zip(perm1_domain, perm2_codomain)
         ax1 = axes(a1)[i]
         ax2 = axes(a2)[j]
-        dual(ax1) == ax2 || throw(
+        axes_match_for_contraction(ax1, ax2) || throw(
             ArgumentError(
-                "Contracted axes do not match: `axes(a1)[$i] = $ax1` and `axes(a2)[$j] = $ax2` are not duals of each other"
+                "Contracted axes do not match: `axes(a1)[$i] = $ax1` and `axes(a2)[$j] = $ax2`"
             )
         )
     end

From ac1ffce31cf02ea4d0ae24275831ccbe902aa043 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Wed, 3 Jun 2026 22:00:53 -0400
Subject: [PATCH 21/31] Use `axes(a, i)` instead of `axes(a)[i]` in the
 contracted-axis check

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

diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl
index 982d350f..42803a65 100644
--- a/src/tensoralgebra.jl
+++ b/src/tensoralgebra.jl
@@ -239,11 +239,11 @@ function TensorAlgebra.check_input(
         a2::AbstractArray, perm2_codomain, perm2_domain
     )
     for (i, j) in zip(perm1_domain, perm2_codomain)
-        ax1 = axes(a1)[i]
-        ax2 = axes(a2)[j]
+        ax1 = axes(a1, i)
+        ax2 = axes(a2, j)
         axes_match_for_contraction(ax1, ax2) || throw(
             ArgumentError(
-                "Contracted axes do not match: `axes(a1)[$i] = $ax1` and `axes(a2)[$j] = $ax2`"
+                "Contracted axes do not match: `axes(a1, $i) = $ax1` and `axes(a2, $j) = $ax2`"
             )
         )
     end

From 62b8b79edb9e2b9e6c93f04424b4f139d16c6bb9 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Thu, 4 Jun 2026 14:52:32 -0400
Subject: [PATCH 22/31] Regression tests for the TA factorization axes-conj
 fixes
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Adds graded-array round-trip tests for `TensorAlgebra.svd` and
`gram_eigh_full_with_pinv`. The matrix-level tests in TensorAlgebra missed the
unmatricize-axis bugs the TA branch fixed (`axes_S` taken from U/Vᴴ rather
than S itself, and `axes_Y` not conjugated), because the bugs only surface
when the codomain and domain axes are conj-related rather than identical,
which only the graded path exercises. Both round-trips are routed through
`contract`. The natural `U * S * Vᴴ` and `X * X'` chains fall into the generic
scalar-indexing matmul on `AbstractGradedMatrix` and are `@test_broken`
pending a block-aware `*`.

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
---
 test/test_tensoralgebra.jl | 45 +++++++++++++++++++++++++++++++++++++-
 1 file changed, 44 insertions(+), 1 deletion(-)

diff --git a/test/test_tensoralgebra.jl b/test/test_tensoralgebra.jl
index cdcc8bad..e4fdccd7 100644
--- a/test/test_tensoralgebra.jl
+++ b/test/test_tensoralgebra.jl
@@ -8,7 +8,7 @@ using GradedArrays: AbelianGradedArray, AbelianGradedMatrix, AbelianSectorArray,
 using Random: randn!
 using TensorAlgebra:
     TensorAlgebra, FusionStyle, contract, linearbroadcasted, matricize, unmatricize
-using Test: @test, @test_throws, @testset
+using Test: @test, @test_broken, @test_throws, @testset
 
 @testset "AbelianSectorArray linear broadcasting" begin
     s = AbelianSectorArray(
@@ -405,3 +405,46 @@ end
     c = similar(m, Float64)
     @test_throws ArgumentError c .= 3 .* m .+ 2 .* n
 end
+
+# Regression coverage for the TensorAlgebra-level unmatricize-axis bugs that
+# previously surfaced through this stack:
+#  (a) `svd!!` assigned `axes_S = (axes(U, 2), axes(Vᴴ, 1))`, but for graded
+#      operators U's rank axis and S's first axis are conj-related rather than
+#      identical, so the unmatricized `S` had the conj of its own axes and
+#      `U * S * Vᴴ` failed to contract.
+#  (b) `gram_eigh_full_with_pinv` assigned `axes_Y = (axes(Y, 1), axes_codomain)`,
+#      but Y's domain contracts with X's codomain and must carry `conj(axes_codomain)`
+#      on a graded backend; the previous form wrote forbidden-charge blocks into
+#      the unmatricize buffer (`Block (i, j) is not stored`).
+@testset "TA.svd round-trip on AbelianGradedArray (axes_S regression)" begin
+    s = gradedrange([U1(0) => 2, U1(1) => 3, U1(2) => 2])
+    A = AbelianGradedArray{Float64}(undef, s, dual(s))
+    randn!(A)
+    U, S, Vᴴ = TensorAlgebra.svd(A, (1,), (2,))
+    # The natural `U * S * Vᴴ` form falls into LinearAlgebra's `_tri_matmul`,
+    # which scalar-indexes on `AbstractGradedArray`. `contract` is the
+    # block-wise route, but the chain form should also work once a block-aware
+    # matmul lands on `AbstractGradedMatrix`.
+    US = contract((:a, :r), U, (:a, :i), S, (:i, :r))
+    USV = contract((:a, :b), US, (:a, :r), Vᴴ, (:r, :b))
+    @test A ≈ USV
+    @test_broken A ≈ U * S * Vᴴ
+end
+
+@testset "TA.gram_eigh_full_with_pinv on AbelianGradedMatrix (axes_Y regression)" begin
+    s = gradedrange([U1(0) => 2, U1(1) => 3, U1(2) => 2])
+    B = AbelianGradedArray{Float64}(undef, s, dual(s))
+    randn!(B)
+    # PSD by construction. Build `A = B * B'` block-wise via `contract`
+    # so we stay on the graded matmul path; the natural `*` form is broken
+    # against the same scalar-indexing path as the SVD round-trip above.
+    A = contract((:a, :b), B, (:a, :r), conj(B), (:b, :r))
+    @test_broken A ≈ B * B'
+    X, Y = TensorAlgebra.gram_eigh_full_with_pinv(A, (1,), (2,))
+    # X · conj(X) ≈ A on the rank subspace.
+    @test A ≈ contract((:a, :b), X, (:a, :r), conj(X), (:b, :r))
+    @test_broken A ≈ X * X'
+    # Y is a left inverse of X on the rank subspace.
+    YX = contract((:r, :s), Y, (:r, :a), X, (:a, :s))
+    @test YX ≈ TensorAlgebra.one(YX, (:r, :s), (:r,), (:s,))
+end

From b488d9373bf1475e1a8f7716133ddf45b3b490f8 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Thu, 4 Jun 2026 14:52:47 -0400
Subject: [PATCH 23/31] Block-aware `rand`/`randn` constructors and `iszero` on
 `AbelianGradedArray`

Adds `Base.rand`/`Base.randn(rng, T, ax)` constructors for
`Tuple{GradedOneTo, ...}` axes, a 3-arg
`Random.rand!(rng, ::AbelianGradedArray, ::Random.Sampler)`, and `Base.iszero`.
The constructors allocate an `AbelianGradedArray{T}` with the right block
structure and fill via the in-place methods. The 3-arg `rand!` replaces the
prior 2-arg form so that the `rand!(A)` / `rand!(A, X)` / `rand!(rng, A, X)`
shims, which all converge on the 3-arg `Sampler` path, reach the block-aware
fill. The generic `AbstractArray` paths for all four scalar-index and throw.

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
---
 src/abeliangradedarray.jl | 32 ++++++++++++++++++++++++++++++--
 test/test_abelianarray.jl | 33 +++++++++++++++++++++++++++++++++
 2 files changed, 63 insertions(+), 2 deletions(-)

diff --git a/src/abeliangradedarray.jl b/src/abeliangradedarray.jl
index 90b5e846..029ce7d2 100644
--- a/src/abeliangradedarray.jl
+++ b/src/abeliangradedarray.jl
@@ -375,11 +375,22 @@ function Base.fill!(a::AbelianGradedArray, v)
     return a
 end
 
+# Block-aware iszero: non-stored blocks are implicitly zero, so an
+# `AbelianGradedArray` is zero iff every stored block is zero. The generic
+# `Base.iszero(::AbstractArray) = all(iszero, x)` path iterates elements,
+# which throws on the no-scalar-indexing guard.
+function Base.iszero(a::AbelianGradedArray)
+    return all(iszero, values(a.blockdata))
+end
+
 # Block-aware random fills: dispatch to the underlying block's `rand!`/`randn!`,
 # bypassing the generic `AbstractArray` fallbacks that go through scalar indexing.
-function Random.rand!(rng::AbstractRNG, a::AbelianGradedArray)
+# The 3-arg `Random.rand!(rng, A, sp::Sampler)` form is what Random's `rand!(A)`
+# / `rand!(A, X)` / `rand!(rng, A, X)` shims ultimately call, so overriding it
+# covers every entry point.
+function Random.rand!(rng::AbstractRNG, a::AbelianGradedArray, sp::Random.Sampler)
     for b in values(a.blockdata)
-        Random.rand!(rng, b)
+        Random.rand!(rng, b, sp)
     end
     return a
 end
@@ -390,6 +401,23 @@ function Random.randn!(rng::AbstractRNG, a::AbelianGradedArray)
     return a
 end
 
+# Constructors `rand(rng, T, axes)` / `randn(rng, T, axes)` for graded axes:
+# allocate an `AbelianGradedArray` from the axes, then fill via the block-aware
+# in-place methods above. The generic `Base.rand`/`randn` fallbacks build a
+# `Matrix` from `length.(axes)`, which loses the graded structure.
+function Base.rand(
+        rng::AbstractRNG, ::Type{T},
+        ax::Tuple{GradedOneTo, Vararg{GradedOneTo}}
+    ) where {T}
+    return Random.rand!(rng, AbelianGradedArray{T}(undef, ax))
+end
+function Base.randn(
+        rng::AbstractRNG, ::Type{T},
+        ax::Tuple{GradedOneTo, Vararg{GradedOneTo}}
+    ) where {T}
+    return Random.randn!(rng, AbelianGradedArray{T}(undef, ax))
+end
+
 # Block-aware diagonal check: block-diagonal (no off-diagonal stored blocks), and each
 # stored diagonal block is itself diagonal. Bypasses the generic scalar-indexing path.
 function LinearAlgebra.isdiag(A::AbelianGradedMatrix)
diff --git a/test/test_abelianarray.jl b/test/test_abelianarray.jl
index 64ae6656..cf7a62b8 100644
--- a/test/test_abelianarray.jl
+++ b/test/test_abelianarray.jl
@@ -5,6 +5,7 @@ using GradedArrays: GradedArrays, AbelianGradedArray, AbelianSectorArray,
     AbstractGradedArray, FusedGradedMatrix, GradedOneTo, SU2, SectorRange, U1, data,
     datalengths, dual, gradedrange, isdual, sectoraxes, sectors, sectortype
 using LinearAlgebra: LinearAlgebra
+using Random: Random
 using TensorKitSectors: TensorKitSectors as TKS
 using Test: @test, @test_throws, @testset
 
@@ -373,3 +374,35 @@ end
     @test size(m_undef.blocks[U1(0)]) == (2, 4)
     @test size(m_undef.blocks[U1(1)]) == (3, 5)
 end
+
+@testset "Block-aware random fills and iszero on AbelianGradedArray" begin
+    g1 = gradedrange([U1(0) => 2, U1(1) => 3])
+    rng = Random.Xoshiro(42)
+
+    # In-place rand!/randn! fill each stored block via the underlying
+    # block's method, no scalar indexing. Both the no-rng and rng-explicit
+    # entry points must work.
+    a = AbelianGradedArray{Float64}(undef, g1, dual(g1))
+    fill!(a, 0)
+    @test iszero(a)
+    Random.randn!(rng, a)
+    @test !iszero(a)
+    fill!(a, 0)
+    Random.randn!(a)
+    @test !iszero(a)
+    Random.rand!(rng, a)
+    @test !iszero(a)
+    fill!(a, 0)
+    Random.rand!(a)
+    @test !iszero(a)
+
+    # Constructor form: `rand(T, axes)` / `randn(T, axes)` for graded axes
+    # builds an `AbelianGradedArray` with the right block structure.
+    r = randn(rng, Float64, (g1, dual(g1)))
+    @test r isa AbelianGradedArray{Float64, 2}
+    @test axes(r) == (g1, dual(g1))
+    @test !iszero(r)
+    u = rand(rng, Float64, (g1, dual(g1)))
+    @test u isa AbelianGradedArray{Float64, 2}
+    @test !iszero(u)
+end

From 5c51d0dbf6d48abb939422d942008c46597bc363 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Fri, 5 Jun 2026 23:00:22 -0400
Subject: [PATCH 24/31] Define `FI.permuteddims` eagerly on the abstract array
 types

Add `FI.permuteddims(a::T, perm) = permutedims(a, perm)` on `AbstractGradedArray`, `AbstractSectorArray`, and `AbstractSectorDelta`. The default `FI.permuteddims` returns a `Base.PermutedDimsArray` view, which falls through GradedArrays' block-aware `bipermutedimsopadd!` overloads (the wrapper itself is not a graded array) and reaches TensorAlgebra's generic implementation. That generic implementation broadcasts back into the graded destination, which lowers through the `GradedStyle` pipeline again and infinite-recurses. NamedDimsArrays' broadcast alignment uses `FI.permuteddims` to align same-name-different-order operands, so any such broadcast on an `AbstractGradedArray` would hang without an eager override.

Defining the override at the abstract types covers every concrete GradedArrays array type at once and lets the previous narrower per-subtype overrides go away. Handling lazy permutations systematically, and the future of the `FI.permuteddims` interface itself, remain as follow-ups.
---
 src/abeliansectorarray.jl  | 9 ---------
 src/abeliansectordelta.jl  | 3 ---
 src/abstractgradedarray.jl | 5 +++++
 src/abstractsectorarray.jl | 5 +++++
 src/abstractsectordelta.jl | 5 +++++
 5 files changed, 15 insertions(+), 12 deletions(-)

diff --git a/src/abeliansectorarray.jl b/src/abeliansectorarray.jl
index b68d9c50..3ea7fb47 100644
--- a/src/abeliansectorarray.jl
+++ b/src/abeliansectorarray.jl
@@ -104,15 +104,6 @@ function Base.permutedims!(y::AbelianSectorArray, x::AbelianSectorArray, perm)
     TensorAlgebra.permutedimsopadd!(y, identity, x, perm, true, false)
     return y
 end
-# Eager fallback for non-bosonic braiding: a lazy view of a fermionic array
-# would skip the permutation phase under scalar indexing.
-function FI.permuteddims(x::AbelianSectorArray, perm)
-    return if TKS.BraidingStyle(sectortype(x)) isa TKS.Bosonic
-        AbelianSectorArray(permutedims(sector(x), perm), FI.permuteddims(data(x), perm))
-    else
-        permutedims(x, perm)
-    end
-end
 
 # ========================  mul!  ========================
 
diff --git a/src/abeliansectordelta.jl b/src/abeliansectordelta.jl
index 8936f81f..d7a7164c 100644
--- a/src/abeliansectordelta.jl
+++ b/src/abeliansectordelta.jl
@@ -38,9 +38,6 @@ function Base.permutedims(x::AbelianSectorDelta, perm)
     new_sectors = ntuple(n -> x.sectors[perm[n]], Val(ndims(x)))
     return AbelianSectorDelta{eltype(x)}(new_sectors)
 end
-function FI.permuteddims(x::AbelianSectorDelta, perm)
-    return permutedims(x, perm)
-end
 
 # ========================  adjoint / broadcasting  ========================
 
diff --git a/src/abstractgradedarray.jl b/src/abstractgradedarray.jl
index 2387d407..1439e89e 100644
--- a/src/abstractgradedarray.jl
+++ b/src/abstractgradedarray.jl
@@ -7,6 +7,11 @@ Concrete subtypes include [`AbelianGradedArray`](@ref) and [`FusedGradedMatrix`]
 abstract type AbstractGradedArray{T, N} <: AbstractArray{T, N} end
 const AbstractGradedMatrix{T} = AbstractGradedArray{T, 2}
 
+# Used by NamedDimsArrays broadcast alignment. Eager for simplicity for now,
+# pending the follow-ups on lazy permutations and the `FI.permuteddims`
+# interface itself.
+FI.permuteddims(a::AbstractGradedArray, perm) = permutedims(a, perm)
+
 function BlockSparseArrays.isblockdiagonal(A::AbstractGradedMatrix)
     for bI in eachblockstoredindex(A)
         row, col = Tuple(bI)
diff --git a/src/abstractsectorarray.jl b/src/abstractsectorarray.jl
index 7959c86e..e04c1541 100644
--- a/src/abstractsectorarray.jl
+++ b/src/abstractsectorarray.jl
@@ -9,6 +9,11 @@ Concrete subtypes:
 """
 abstract type AbstractSectorArray{T, N} <: AbstractArray{T, N} end
 
+# Used by NamedDimsArrays broadcast alignment. Eager for simplicity for now,
+# pending the follow-ups on lazy permutations and the `FI.permuteddims`
+# interface itself.
+FI.permuteddims(a::AbstractSectorArray, perm) = permutedims(a, perm)
+
 """
     data(sa::AbstractSectorArray)
 
diff --git a/src/abstractsectordelta.jl b/src/abstractsectordelta.jl
index 2dacba80..18aa3acf 100644
--- a/src/abstractsectordelta.jl
+++ b/src/abstractsectordelta.jl
@@ -9,5 +9,10 @@ Concrete subtypes:
 """
 abstract type AbstractSectorDelta{T, N} <: AbstractArray{T, N} end
 
+# Used by NamedDimsArrays broadcast alignment. Eager for simplicity for now,
+# pending the follow-ups on lazy permutations and the `FI.permuteddims`
+# interface itself.
+FI.permuteddims(a::AbstractSectorDelta, perm) = permutedims(a, perm)
+
 Base.copy(A::AbstractSectorDelta) = A
 Base.size(A::AbstractSectorDelta) = length.(axes(A))

From 4f07f35b03875de9ce38246b8afb88d2b76e441f Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Mon, 15 Jun 2026 13:38:20 -0400
Subject: [PATCH 25/31] Add tests for the block-aware graded array methods

Cover behaviors that previously had no direct tests: the `contract` input check that rejects contracted-axis pairs with mismatched duality, `projectto!` / `checked_projectto!` on `AbelianGradedArray`, the `conj` axis-duality flip, `isdiag` on `AbelianGradedMatrix`, and dropping fully truncated sectors from the bond in truncated SVD.
---
 test/test_abelianarray.jl   | 61 +++++++++++++++++++++++++++++++++++++
 test/test_factorizations.jl | 15 +++++++++
 test/test_tensoralgebra.jl  | 18 +++++++++++
 3 files changed, 94 insertions(+)

diff --git a/test/test_abelianarray.jl b/test/test_abelianarray.jl
index cf7a62b8..f81fee46 100644
--- a/test/test_abelianarray.jl
+++ b/test/test_abelianarray.jl
@@ -6,6 +6,7 @@ using GradedArrays: GradedArrays, AbelianGradedArray, AbelianSectorArray,
     datalengths, dual, gradedrange, isdual, sectoraxes, sectors, sectortype
 using LinearAlgebra: LinearAlgebra
 using Random: Random
+using TensorAlgebra: TensorAlgebra
 using TensorKitSectors: TensorKitSectors as TKS
 using Test: @test, @test_throws, @testset
 
@@ -406,3 +407,63 @@ end
     @test u isa AbelianGradedArray{Float64, 2}
     @test !iszero(u)
 end
+
+@testset "conj flips axis duality" begin
+    g = gradedrange([U1(0) => 2, U1(1) => 3])
+    a = AbelianGradedArray{ComplexF64}(undef, g, dual(g))
+    a[Block(1, 1)] = AbelianSectorArray((U1(0), dual(U1(0))), randn(ComplexF64, 2, 2))
+    a[Block(2, 2)] = AbelianSectorArray((U1(1), dual(U1(1))), randn(ComplexF64, 3, 3))
+
+    ca = conj(a)
+    @test ca isa AbelianGradedArray
+    # Each axis flips its duality, mirroring `conj` on the axis types.
+    @test isdual(axes(ca, 1)) == !isdual(axes(a, 1))
+    @test isdual(axes(ca, 2)) == !isdual(axes(a, 2))
+    @test axes(ca, 1) == conj(axes(a, 1))
+    @test axes(ca, 2) == conj(axes(a, 2))
+    # The data is conjugated element-wise.
+    @test Array(ca) ≈ conj(Array(a))
+end
+
+@testset "isdiag on AbelianGradedMatrix" begin
+    g = gradedrange([U1(0) => 2, U1(1) => 2])
+
+    # Block-diagonal with each stored block diagonal.
+    a = AbelianGradedArray{Float64}(undef, g, dual(g))
+    a[Block(1, 1)] = AbelianSectorArray((U1(0), dual(U1(0))), [1.0 0.0; 0.0 2.0])
+    a[Block(2, 2)] = AbelianSectorArray((U1(1), dual(U1(1))), [3.0 0.0; 0.0 4.0])
+    @test LinearAlgebra.isdiag(a)
+
+    # A non-diagonal stored block breaks it.
+    b = AbelianGradedArray{Float64}(undef, g, dual(g))
+    b[Block(1, 1)] = AbelianSectorArray((U1(0), dual(U1(0))), [1.0 5.0; 0.0 2.0])
+    b[Block(2, 2)] = AbelianSectorArray((U1(1), dual(U1(1))), [3.0 0.0; 0.0 4.0])
+    @test !LinearAlgebra.isdiag(b)
+end
+
+@testset "projectto! / checked_projectto!" begin
+    g = gradedrange([U1(0) => 2, U1(1) => 3])
+    dest = AbelianGradedArray{Float64}(undef, g, dual(g))
+
+    # Project a dense source: the symmetry-allowed (block-diagonal) regions are
+    # copied, the forbidden off-diagonal regions are dropped.
+    src = randn(5, 5)
+    @test TensorAlgebra.projectto!(dest, src) === dest
+    @test data(dest[Block(1, 1)]) ≈ src[1:2, 1:2]
+    @test data(dest[Block(2, 2)]) ≈ src[3:5, 3:5]
+    proj = Array(dest)
+    @test iszero(proj[1:2, 3:5])
+    @test iszero(proj[3:5, 1:2])
+
+    # `checked_projectto!` accepts a source already in the allowed subspace.
+    src_allowed = Array(dest)
+    dest_ok = AbelianGradedArray{Float64}(undef, g, dual(g))
+    @test TensorAlgebra.checked_projectto!(dest_ok, src_allowed) === dest_ok
+    @test Array(dest_ok) ≈ src_allowed
+
+    # ... and rejects one carrying significant forbidden-block weight.
+    src_bad = copy(src_allowed)
+    src_bad[1, 5] += 10.0
+    dest_bad = AbelianGradedArray{Float64}(undef, g, dual(g))
+    @test_throws InexactError TensorAlgebra.checked_projectto!(dest_bad, src_bad)
+end
diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl
index 7822d9e0..2c89715a 100644
--- a/test/test_factorizations.jl
+++ b/test/test_factorizations.jl
@@ -360,6 +360,21 @@ end
             @test U isa FusedGradedMatrix
             @test sum(size(b, 2) for b in values(U.blocks)) <= 3
         end
+
+        @testset "drops fully truncated sectors from the bond" begin
+            # U1(0) carries singular values of order 1, U1(1) only of order 1e-3,
+            # so a tolerance between the two scales removes U1(1) from the bond
+            # entirely (not just shrinks it).
+            A = FusedGradedMatrix(
+                [U1(0), U1(1)], [Matrix(1.0I, 2, 2), 1.0e-3 * Matrix(1.0I, 2, 2)]
+            )
+            U, S, Vᴴ, ε = MAK.svd_trunc(A; trunc = trunctol(; atol = 1.0e-2))
+            @test collect(keys(U.blocks)) == [U1(0)]
+            @test collect(keys(S.blocks)) == [U1(0)]
+            @test collect(keys(Vᴴ.blocks)) == [U1(0)]
+            # The dropped sector's weight shows up as the truncation error.
+            @test ε ≈ norm(1.0e-3 * Matrix(1.0I, 2, 2)) atol = precision(eltype(A))
+        end
     end
 
     # -----------------------------------------------------------------------
diff --git a/test/test_tensoralgebra.jl b/test/test_tensoralgebra.jl
index e4fdccd7..85375771 100644
--- a/test/test_tensoralgebra.jl
+++ b/test/test_tensoralgebra.jl
@@ -448,3 +448,21 @@ end
     YX = contract((:r, :s), Y, (:r, :a), X, (:a, :s))
     @test YX ≈ TensorAlgebra.one(YX, (:r, :s), (:r,), (:s,))
 end
+
+@testset "contract rejects mismatched contracted-axis duality (bosonic)" begin
+    g = gradedrange([U1(0) => 2, U1(1) => 3])
+    a = AbelianGradedArray{Float64}(undef, g, dual(g))
+    randn!(a)
+
+    # The contracted leg of `a` is `dual(g)`; here `b`'s contracted leg is also
+    # `dual(g)`, which is neither the canonical dual pairing nor (for bosonic
+    # U1) an accepted same-`isdual` pair, so the contraction is rejected.
+    b = AbelianGradedArray{Float64}(undef, dual(g), dual(g))
+    @test_throws ArgumentError contract(a, (1, -1), b, (-1, 2))
+
+    # Sanity: the canonically dual-paired contraction is accepted.
+    b_ok = AbelianGradedArray{Float64}(undef, g, dual(g))
+    randn!(b_ok)
+    result, = contract(a, (1, -1), b_ok, (-1, 2))
+    @test result isa AbelianGradedArray
+end

From d60be2c70bd05857dccaa7d05355f05a5691c9ad Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Mon, 15 Jun 2026 13:38:31 -0400
Subject: [PATCH 26/31] Share the AbelianGradedArray view via view_abelian

The `Block{N}` view method and its `Block{1}` disambiguator duplicated the `haskey`/sectors/wrap logic, so both now delegate to a shared `view_abelian(a, bk)`.
---
 src/abeliangradedarray.jl | 12 +++++-------
 1 file changed, 5 insertions(+), 7 deletions(-)

diff --git a/src/abeliangradedarray.jl b/src/abeliangradedarray.jl
index 029ce7d2..f9dd5c2d 100644
--- a/src/abeliangradedarray.jl
+++ b/src/abeliangradedarray.jl
@@ -85,21 +85,19 @@ datatype(::Type{<:AbelianGradedArray{T, N, D, S}}) where {T, N, D, S} = D
 #  view (primitive): returns AbelianSectorArray sharing data with blockdata
 # ---------------------------------------------------------------------------
 
-function Base.view(a::AbelianGradedArray{T, N}, I::Block{N}) where {T, N}
+# Shared implementation: build the `AbelianSectorArray` view for a stored block.
+function view_abelian(a::AbelianGradedArray{T, N}, I::Block{N}) where {T, N}
     bk = Int.(Tuple(I))
     haskey(a.blockdata, bk) || error("Block $bk is not stored.")
     sects = ntuple(d -> sectors(axes(a, d))[bk[d]], Val(N))
     return AbelianSectorArray(sects, a.blockdata[bk])
 end
 
+Base.view(a::AbelianGradedArray{T, N}, I::Block{N}) where {T, N} = view_abelian(a, I)
+
 # Disambiguate against `view(::AbstractGradedArray{T, N}, ::Vararg{Block{1}, N})` for
 # N=1, where the splatted form collapses to a single Block{1} argument.
-function Base.view(a::AbelianGradedArray{T, 1}, I::Block{1}) where {T}
-    bk = (Int(I),)
-    haskey(a.blockdata, bk) || error("Block $bk is not stored.")
-    sects = (sectors(axes(a, 1))[bk[1]],)
-    return AbelianSectorArray(sects, a.blockdata[bk])
-end
+Base.view(a::AbelianGradedArray{T, 1}, I::Block{1}) where {T} = view_abelian(a, I)
 
 # ---------------------------------------------------------------------------
 #  blocks — lazy view delegating to view (following BlockArrays convention)

From 1cc2ad7a88e46c1773100f2618ba6841306bf46c Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Mon, 15 Jun 2026 14:32:25 -0400
Subject: [PATCH 27/31] Consolidate scalar, fill, and copy operations on
 AbstractGradedArray

Define `scale!`, `zero!`, and `fill!` once on `AbstractGradedArray` via the `eachblockstoredindex`/`view` block interface, replacing the per-type versions on `AbelianGradedArray`, `FusedGradedMatrix`, and `FusedGradedVector`. `fill!` is now uniformly permissive: it fills the stored, symmetry-allowed blocks with any value.

`AbelianGradedArray` scalar `*` and `/` route through the Base `AbstractArray`-scalar broadcasting, so the dedicated overrides are removed. `FusedGradedMatrix` keeps explicit scalar methods because its broadcasting is disabled, with a TODO to drop them once that is supported.

`copyto!` between `AbelianGradedArray`s with matching axes copies each block into the existing destination buffer instead of clearing and reallocating, since matching axes imply matching allocated blocks.

`checked_projectto!` forwards its keyword arguments to `isapprox`.
---
 src/abeliangradedarray.jl  | 46 ++++++++------------------------------
 src/abstractgradedarray.jl | 32 ++++++++++++++++++++++++++
 src/fusedgradedmatrix.jl   | 21 +++++------------
 src/fusedgradedvector.jl   | 16 -------------
 4 files changed, 46 insertions(+), 69 deletions(-)

diff --git a/src/abeliangradedarray.jl b/src/abeliangradedarray.jl
index f9dd5c2d..fc5d6e0c 100644
--- a/src/abeliangradedarray.jl
+++ b/src/abeliangradedarray.jl
@@ -304,9 +304,10 @@ function Base.copyto!(
             "copyto! axes mismatch: dest $(axes(dest)), src $(axes(src))"
         )
     )
-    empty!(dest.blockdata)
+    # Matching axes mean matching allowed-block keys (every allowed block is
+    # allocated), so copy each block into the existing destination buffer.
     for (k, v) in src.blockdata
-        dest.blockdata[k] = copy(v)
+        copyto!(dest.blockdata[k], v)
     end
     return dest
 end
@@ -347,32 +348,6 @@ function Base.permutedims!(
     return y
 end
 
-# ---------------------------------------------------------------------------
-#  fill! / zero! / scale!
-# ---------------------------------------------------------------------------
-
-scale!(a::AbstractArray, β::Number) = (a .*= β; a)
-function scale!(a::AbelianGradedArray, β::Number)
-    for b in values(a.blockdata)
-        scale!(b, β)
-    end
-    return a
-end
-
-function FI.zero!(a::AbelianGradedArray)
-    for b in values(a.blockdata)
-        FI.zero!(b)
-    end
-    return a
-end
-
-function Base.fill!(a::AbelianGradedArray, v)
-    for b in values(a.blockdata)
-        fill!(b, v)
-    end
-    return a
-end
-
 # Block-aware iszero: non-stored blocks are implicitly zero, so an
 # `AbelianGradedArray` is zero iff every stored block is zero. The generic
 # `Base.iszero(::AbstractArray) = all(iszero, x)` path iterates elements,
@@ -445,13 +420,12 @@ end
 
 # Compare via `Array(dest)` so the generic `isapprox(::AbstractArray, ::AbelianGradedArray)`
 # path doesn't fall back to a `src - dest` broadcast that scalar-indexes the block storage.
+# `kwargs` (e.g. `atol`/`rtol`) are forwarded to `isapprox`, matching the generic verb.
 function TensorAlgebra.checked_projectto!(
-        dest::AbelianGradedArray, src::AbstractArray;
-        atol::Real = 0,
-        rtol::Real = Base.rtoldefault(real(eltype(src)))
+        dest::AbelianGradedArray, src::AbstractArray; kwargs...
     )
     TensorAlgebra.projectto!(dest, src)
-    isapprox(src, Array(dest); atol, rtol) ||
+    isapprox(src, Array(dest); kwargs...) ||
         throw(InexactError(:checked_projectto!, typeof(dest), src))
     return dest
 end
@@ -506,11 +480,9 @@ function Base.mapreduce(
     return s
 end
 
-# Block-wise scalar multiplication. The default `a .* x` / `a ./ x` paths broadcast
-# element-wise and scalar-index opaque storage.
-Base.:*(a::AbelianGradedArray, x::Number) = scale!(copy(a), x)
-Base.:*(x::Number, a::AbelianGradedArray) = a * x
-Base.:/(a::AbelianGradedArray, x::Number) = a * inv(x)
+# Scalar `*` / `/` are inherited from Base's `AbstractArray`-scalar methods, which
+# forward to broadcasting (`a .* x` / `a ./ x`). `AbelianGradedArray` supports the
+# linear-broadcast path, so no dedicated overrides are needed here.
 
 # `LinearAlgebra.normalize` infers its result eltype via `typeof(first(a)/nrm)`,
 # which scalar-indexes opaque storage.
diff --git a/src/abstractgradedarray.jl b/src/abstractgradedarray.jl
index 1439e89e..d54144b0 100644
--- a/src/abstractgradedarray.jl
+++ b/src/abstractgradedarray.jl
@@ -105,6 +105,38 @@ end
 datatype(a::AbstractGradedArray) = datatype(typeof(a))
 sectortype(a::AbstractGradedArray) = sectortype(typeof(a))
 
+# ---------------------------------------------------------------------------
+#  fill! / zero! / scale! — block-wise over the stored blocks
+#
+#  Defined once via the `eachblockstoredindex`/`view` interface every
+#  `AbstractGradedArray` implements, so every concrete subtype is covered.
+#  These only touch stored (symmetry-allowed) blocks, so a nonzero `fill!`
+#  value leaves the forbidden positions at zero.
+# ---------------------------------------------------------------------------
+
+scale!(a::AbstractArray, β::Number) = (a .*= β; a)
+
+function scale!(a::AbstractGradedArray, β::Number)
+    for bI in eachblockstoredindex(a)
+        scale!(view(a, bI), β)
+    end
+    return a
+end
+
+function FI.zero!(a::AbstractGradedArray)
+    for bI in eachblockstoredindex(a)
+        FI.zero!(view(a, bI))
+    end
+    return a
+end
+
+function Base.fill!(a::AbstractGradedArray, v)
+    for bI in eachblockstoredindex(a)
+        fill!(view(a, bI), v)
+    end
+    return a
+end
+
 # ---------------------------------------------------------------------------
 #  Display — convert to BlockSparseArray for printing
 # ---------------------------------------------------------------------------
diff --git a/src/fusedgradedmatrix.jl b/src/fusedgradedmatrix.jl
index 1ee72f63..aa1d41a4 100644
--- a/src/fusedgradedmatrix.jl
+++ b/src/fusedgradedmatrix.jl
@@ -171,22 +171,6 @@ function BlockArrays.blocks(m::FusedGradedMatrix)
     return [view(m, I) for I in eachblockstoredindex(m)]
 end
 
-# ========================  fill! / zero!  ========================
-
-function FI.zero!(m::FusedGradedMatrix)
-    for b in values(m.blocks)
-        fill!(b, zero(eltype(b)))
-    end
-    return m
-end
-
-function Base.fill!(m::FusedGradedMatrix, v)
-    iszero(v) || throw(
-        ArgumentError("fill! with nonzero value is not supported for FusedGradedMatrix")
-    )
-    return FI.zero!(m)
-end
-
 # ========================  mul!  ========================
 
 function TensorAlgebra.check_input(::typeof(*), A::FusedGradedMatrix, B::FusedGradedMatrix)
@@ -283,6 +267,11 @@ function Base.:(-)(A::FusedGradedMatrix, B::FusedGradedMatrix)
     return _block_combine(-, A, B)
 end
 
+# TODO: these explicit scalar-op methods exist only because broadcasting is
+# disabled for `FusedGradedMatrix`. Once structure-preserving broadcasting is
+# supported, drop them and let Base's `AbstractArray`-scalar `*` / `/` forward to
+# broadcasting (as `AbelianGradedArray` already does). See the
+# `fusedgradedmatrix_broadcasting` project.
 function Base.:(*)(A::FusedGradedMatrix, x::Number)
     new_blocks = map(b -> b * x, A.blocks)
     return FusedGradedMatrix(A.codomain, A.domain, new_blocks)
diff --git a/src/fusedgradedvector.jl b/src/fusedgradedvector.jl
index 549edfa1..7f193102 100644
--- a/src/fusedgradedvector.jl
+++ b/src/fusedgradedvector.jl
@@ -210,22 +210,6 @@ function BlockArrays.blocks(v::FusedGradedVector)
     return [view(v, I) for I in eachblockstoredindex(v)]
 end
 
-# ========================  fill! / zero!  ========================
-
-function FI.zero!(v::FusedGradedVector)
-    for b in values(v.blocks)
-        fill!(b, zero(eltype(v)))
-    end
-    return v
-end
-
-function Base.fill!(v::FusedGradedVector, val)
-    for b in values(v.blocks)
-        fill!(b, val)
-    end
-    return v
-end
-
 # ========================  similar  ========================
 
 function Base.similar(v::FusedGradedVector, ::Type{T}) where {T}

From b1e25484cfd7ff63551b99115b1788b93f373026 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Mon, 15 Jun 2026 15:56:29 -0400
Subject: [PATCH 28/31] Tidy graded-array dot, reductions, and view dispatch

Check axes equality in `dot` on `AbelianGradedArray` and drop the redundant per-block `haskey` guard, since matching axes imply matching allowed-block keys. Replace the `+`-only `mapreduce` method with a direct `sum`, the reduction actually needed, which avoids the footgun of skipping forbidden (zero) blocks for reductions where `f(0)` is nonzero.

Remove the unused `view(::AbstractGradedArray{T,1}, ::Block{1})` method, which rebuilt the same `Block{1}` and re-dispatched. It is shadowed by the concrete per-type `Block{1}` views and would otherwise recurse. Rename `axes_match_for_contraction` to `are_axes_contractible`.

Fix the `FusedGradedVector` docstring to say the stored blocks match the axis sectors exactly, which the constructor enforces. Add `dot` and `sum` tests, trim the unmatricize regression comment, and note that block-aware matmul and adjoint on `AbelianGradedMatrix` are still needed.
---
 src/abeliangradedarray.jl  | 19 ++++++++-----------
 src/abstractgradedarray.jl |  2 --
 src/fusedgradedvector.jl   |  4 ++--
 src/tensoralgebra.jl       |  4 ++--
 test/test_abelianarray.jl  | 22 ++++++++++++++++++++++
 test/test_tensoralgebra.jl | 17 +++++++----------
 6 files changed, 41 insertions(+), 27 deletions(-)

diff --git a/src/abeliangradedarray.jl b/src/abeliangradedarray.jl
index fc5d6e0c..64091940 100644
--- a/src/abeliangradedarray.jl
+++ b/src/abeliangradedarray.jl
@@ -457,25 +457,22 @@ function LinearAlgebra.norm(a::AbelianGradedArray, p::Real = 2)
 end
 
 function LinearAlgebra.dot(a::AbelianGradedArray, b::AbelianGradedArray)
+    axes(a) == axes(b) ||
+        throw(DimensionMismatch("dot axes mismatch: a $(axes(a)), b $(axes(b))"))
+    # Matching axes mean matching allowed-block keys, so each `a` block has a
+    # counterpart in `b`.
     s = zero(LinearAlgebra.dot(zero(eltype(a)), zero(eltype(b))))
     for (k, ablk) in pairs(a.blockdata)
-        haskey(b.blockdata, k) || continue
         s += LinearAlgebra.dot(ablk, b.blockdata[k])
     end
     return s
 end
 
-# Reductions iterate scalar-by-scalar by default. Unstored blocks are zero, so
-# `+`-style reductions can skip them. Other ops require materializing — refuse
-# rather than silently dropping unstored contributions.
-function Base.mapreduce(
-        f, op::Union{typeof(+), typeof(Base.add_sum)},
-        a::AbelianGradedArray;
-        init = zero(eltype(a))
-    )
-    s = init
+# Forbidden blocks are zero, so the total is the sum over the stored blocks.
+function Base.sum(a::AbelianGradedArray)
+    s = zero(eltype(a))
     for b in values(a.blockdata)
-        s = op(s, mapreduce(f, op, b))
+        s += sum(b)
     end
     return s
 end
diff --git a/src/abstractgradedarray.jl b/src/abstractgradedarray.jl
index d54144b0..1f2a8aeb 100644
--- a/src/abstractgradedarray.jl
+++ b/src/abstractgradedarray.jl
@@ -51,8 +51,6 @@ end
 function Base.view(a::AbstractGradedArray{T, N}, I::Vararg{Block{1}, N}) where {T, N}
     return view(a, Block(Int.(I)))
 end
-# Disambiguate against subtype-specific `view(::ConcreteGradedVector, ::Block{1})` methods.
-Base.view(a::AbstractGradedArray{T, 1}, I::Block{1}) where {T} = view(a, Block((Int(I),)))
 
 function Base.getindex(a::AbstractGradedArray{T, N}, I::Block{N}) where {T, N}
     return copy(view(a, I))
diff --git a/src/fusedgradedvector.jl b/src/fusedgradedvector.jl
index 7f193102..b3f6296f 100644
--- a/src/fusedgradedvector.jl
+++ b/src/fusedgradedvector.jl
@@ -82,8 +82,8 @@ Fields:
 
   - `axis::Dictionary{S,Int}` — axis layout, mapping each sector to its block
     size. Keys are sorted and unique. Stored non-dual (codomain convention).
-  - `blocks::Dictionary{S,D}` — stored data blocks, keyed by sector. Keys are
-    a subset of `keys(axis)` and `length(blocks[s]) == axis[s]`.
+  - `blocks::Dictionary{S,D}` — stored data blocks, keyed by sector. Keys match
+    `keys(axis)` exactly and `length(blocks[s]) == axis[s]`.
 """
 struct FusedGradedVector{T, D <: AbstractVector{T}, S <: SectorRange} <:
     AbstractGradedArray{T, 1}
diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl
index 42803a65..ae12aa13 100644
--- a/src/tensoralgebra.jl
+++ b/src/tensoralgebra.jl
@@ -222,7 +222,7 @@ end
 # is the canonical non-dual ↔ dual pairing). For fermionic braiding `ax1 == ax2`
 # is also accepted, since the supertrace formalism uses `contraction_twist!`
 # to pick up the right phase on same-`isdual` pairings.
-function axes_match_for_contraction(ax1, ax2)
+function are_axes_contractible(ax1, ax2)
     dual(ax1) == ax2 && return true
     TKS.BraidingStyle(sectortype(typeof(ax1))) isa TKS.Bosonic && return false
     return ax1 == ax2
@@ -241,7 +241,7 @@ function TensorAlgebra.check_input(
     for (i, j) in zip(perm1_domain, perm2_codomain)
         ax1 = axes(a1, i)
         ax2 = axes(a2, j)
-        axes_match_for_contraction(ax1, ax2) || throw(
+        are_axes_contractible(ax1, ax2) || throw(
             ArgumentError(
                 "Contracted axes do not match: `axes(a1, $i) = $ax1` and `axes(a2, $j) = $ax2`"
             )
diff --git a/test/test_abelianarray.jl b/test/test_abelianarray.jl
index f81fee46..85d9ef7c 100644
--- a/test/test_abelianarray.jl
+++ b/test/test_abelianarray.jl
@@ -467,3 +467,25 @@ end
     dest_bad = AbelianGradedArray{Float64}(undef, g, dual(g))
     @test_throws InexactError TensorAlgebra.checked_projectto!(dest_bad, src_bad)
 end
+
+@testset "dot" begin
+    g = gradedrange([U1(0) => 2, U1(1) => 3])
+    a = AbelianGradedArray{ComplexF64}(undef, g, dual(g))
+    b = AbelianGradedArray{ComplexF64}(undef, g, dual(g))
+    Random.randn!(a)
+    Random.randn!(b)
+    @test LinearAlgebra.dot(a, b) ≈ LinearAlgebra.dot(Array(a), Array(b))
+
+    # Mismatched axes are rejected.
+    h = gradedrange([U1(0) => 1, U1(1) => 2])
+    c = AbelianGradedArray{ComplexF64}(undef, h, dual(h))
+    Random.randn!(c)
+    @test_throws DimensionMismatch LinearAlgebra.dot(a, c)
+end
+
+@testset "sum" begin
+    g = gradedrange([U1(0) => 2, U1(1) => 3])
+    a = AbelianGradedArray{ComplexF64}(undef, g, dual(g))
+    Random.randn!(a)
+    @test sum(a) ≈ sum(Array(a))
+end
diff --git a/test/test_tensoralgebra.jl b/test/test_tensoralgebra.jl
index 85375771..1070c500 100644
--- a/test/test_tensoralgebra.jl
+++ b/test/test_tensoralgebra.jl
@@ -406,16 +406,9 @@ end
     @test_throws ArgumentError c .= 3 .* m .+ 2 .* n
 end
 
-# Regression coverage for the TensorAlgebra-level unmatricize-axis bugs that
-# previously surfaced through this stack:
-#  (a) `svd!!` assigned `axes_S = (axes(U, 2), axes(Vᴴ, 1))`, but for graded
-#      operators U's rank axis and S's first axis are conj-related rather than
-#      identical, so the unmatricized `S` had the conj of its own axes and
-#      `U * S * Vᴴ` failed to contract.
-#  (b) `gram_eigh_full_with_pinv` assigned `axes_Y = (axes(Y, 1), axes_codomain)`,
-#      but Y's domain contracts with X's codomain and must carry `conj(axes_codomain)`
-#      on a graded backend; the previous form wrote forbidden-charge blocks into
-#      the unmatricize buffer (`Block (i, j) is not stored`).
+# Regression coverage for TensorAlgebra-level unmatricize-axis bugs on graded
+# operators: a factor's reconstructed axes must respect the conj/dual pairing
+# between contracted bonds rather than reuse the factor's own axes.
 @testset "TA.svd round-trip on AbelianGradedArray (axes_S regression)" begin
     s = gradedrange([U1(0) => 2, U1(1) => 3, U1(2) => 2])
     A = AbelianGradedArray{Float64}(undef, s, dual(s))
@@ -443,6 +436,10 @@ end
     X, Y = TensorAlgebra.gram_eigh_full_with_pinv(A, (1,), (2,))
     # X · conj(X) ≈ A on the rank subspace.
     @test A ≈ contract((:a, :b), X, (:a, :r), conj(X), (:b, :r))
+    # Matmul (`*`) on `AbelianGradedMatrix` is unimplemented, so any `X * Y`
+    # falls through to LinearAlgebra's scalar-indexing path and throws. The
+    # adjoint forms (`X * X'`, `B * B'`) additionally need a block-aware
+    # `adjoint`. Both will pass once those land.
     @test_broken A ≈ X * X'
     # Y is a left inverse of X on the rank subspace.
     YX = contract((:r, :s), Y, (:r, :a), X, (:a, :s))

From fc4444a16624fc471c99ba813d652506cb1e4790 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Mon, 15 Jun 2026 16:14:27 -0400
Subject: [PATCH 29/31] Funnel FusedGradedMatrix arithmetic through a
 broadcasting stand-in

Route the block-wise `+`, `-`, `*`, and `/` on `FusedGradedMatrix` through a single `_broadcast_fusedgradedmatrix` helper, renamed from `_block_combine`. The binary form for `+` and `-` checks axis agreement and combines blocks over the union of sectors, and the single-array form applies a function to each stored block for the scalar ops. The name and a comment make clear that these methods stand in for the broadcasting currently disabled on `FusedGradedMatrix`, and are superseded once structure-preserving broadcasting lands.
---
 src/fusedgradedmatrix.jl | 45 +++++++++++++++++++---------------------
 1 file changed, 21 insertions(+), 24 deletions(-)

diff --git a/src/fusedgradedmatrix.jl b/src/fusedgradedmatrix.jl
index aa1d41a4..acbd79d9 100644
--- a/src/fusedgradedmatrix.jl
+++ b/src/fusedgradedmatrix.jl
@@ -226,19 +226,15 @@ end
 #
 # FusedGradedMatrix has no `BroadcastStyle`, so the AbstractArray fallback for `+`/`-`
 # (`broadcast_preserving_zero_d`) ends up trying scalar indexing on a sector-block
-# structure — which silently produces wrong results. Define the operations directly
-# block-by-block, taking the union of sectors so an absent block on one side is
-# treated as zero.
-
-function _check_add_axes(A::FusedGradedMatrix, B::FusedGradedMatrix, op)
-    A.codomain == B.codomain ||
-        throw(DimensionMismatch("$(op): codomain mismatch"))
-    A.domain == B.domain ||
-        throw(DimensionMismatch("$(op): domain mismatch"))
-    return nothing
-end
-
-function _block_combine(op, A::FusedGradedMatrix, B::FusedGradedMatrix)
+# structure, which silently produces wrong results. As a stand-in, define the
+# operations directly block-by-block, taking the union of sectors so an absent block
+# on one side is treated as zero. Once structure-preserving broadcasting is
+# implemented for FusedGradedMatrix, `_broadcast_fusedgradedmatrix` (and the
+# scalar-op methods below) are superseded by it.
+
+function _broadcast_fusedgradedmatrix(op, A::FusedGradedMatrix, B::FusedGradedMatrix)
+    axes(A) == axes(B) ||
+        throw(DimensionMismatch("axes mismatch: A $(axes(A)), B $(axes(B))"))
     T = promote_type(eltype(A), eltype(B))
     DA = datatype(BlockSparseArrays.blocktype(A))
     DB = datatype(BlockSparseArrays.blocktype(B))
@@ -257,29 +253,30 @@ function _block_combine(op, A::FusedGradedMatrix, B::FusedGradedMatrix)
     return FusedGradedMatrix(A.codomain, A.domain, blocks)
 end
 
-function Base.:(+)(A::FusedGradedMatrix, B::FusedGradedMatrix)
-    _check_add_axes(A, B, :+)
-    return _block_combine(+, A, B)
+# Single-array form: apply `f` to each stored block. Absent blocks stay absent
+# (zero), which is correct for `f` with `f(0) == 0` such as scalar `*` / `/`.
+function _broadcast_fusedgradedmatrix(f, A::FusedGradedMatrix)
+    blocks = map(f, A.blocks)
+    return FusedGradedMatrix(A.codomain, A.domain, blocks)
 end
 
+function Base.:(+)(A::FusedGradedMatrix, B::FusedGradedMatrix)
+    return _broadcast_fusedgradedmatrix(+, A, B)
+end
 function Base.:(-)(A::FusedGradedMatrix, B::FusedGradedMatrix)
-    _check_add_axes(A, B, :-)
-    return _block_combine(-, A, B)
+    return _broadcast_fusedgradedmatrix(-, A, B)
 end
 
 # TODO: these explicit scalar-op methods exist only because broadcasting is
 # disabled for `FusedGradedMatrix`. Once structure-preserving broadcasting is
 # supported, drop them and let Base's `AbstractArray`-scalar `*` / `/` forward to
-# broadcasting (as `AbelianGradedArray` already does). See the
-# `fusedgradedmatrix_broadcasting` project.
+# broadcasting (as `AbelianGradedArray` already does).
 function Base.:(*)(A::FusedGradedMatrix, x::Number)
-    new_blocks = map(b -> b * x, A.blocks)
-    return FusedGradedMatrix(A.codomain, A.domain, new_blocks)
+    return _broadcast_fusedgradedmatrix(Base.Fix2(*, x), A)
 end
 Base.:(*)(x::Number, A::FusedGradedMatrix) = A * x
 function Base.:(/)(A::FusedGradedMatrix, x::Number)
-    new_blocks = map(b -> b / x, A.blocks)
-    return FusedGradedMatrix(A.codomain, A.domain, new_blocks)
+    return _broadcast_fusedgradedmatrix(Base.Fix2(/, x), A)
 end
 
 # ======================== LinearAlgebra ======================

From 8095b913e249026c72841a8d9967add1fe95a71d Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Mon, 15 Jun 2026 17:19:14 -0400
Subject: [PATCH 30/31] Require canonical dual pairing for contracted graded
 axes

Contraction now requires every contracted axis pair to be a canonical dual pair (`dual(ax1) == ax2`) for all sector types, fermionic included, dropping the earlier relaxation that also accepted same-`isdual` fermionic pairings. A contraction always pairs a space with its dual.

Two fermionic contraction tests had inadvertently relied on same-`isdual` pairings. They now test the fermionic twist through proper canonical contractions, where the two odd contracted indices appear in swapped order so a single transposition introduces the expected -1.
---
 src/tensoralgebra.jl   | 15 +++------------
 test/test_fermionic.jl | 15 +++++++++------
 2 files changed, 12 insertions(+), 18 deletions(-)

diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl
index ae12aa13..b1fb48dd 100644
--- a/src/tensoralgebra.jl
+++ b/src/tensoralgebra.jl
@@ -217,17 +217,6 @@ function contraction_twist!(a::AbstractArray, ndims_codomain::Int)
     return twist!(a, (i for i in 1:ndims_codomain if isdual(a, i)))
 end
 
-# True iff `ax1` and `ax2` are valid as partners on a contracted axis. Always
-# requires sector content to match modulo the `isdual` flag (`dual(ax1) == ax2`
-# is the canonical non-dual ↔ dual pairing). For fermionic braiding `ax1 == ax2`
-# is also accepted, since the supertrace formalism uses `contraction_twist!`
-# to pick up the right phase on same-`isdual` pairings.
-function are_axes_contractible(ax1, ax2)
-    dual(ax1) == ax2 && return true
-    TKS.BraidingStyle(sectortype(typeof(ax1))) isa TKS.Bosonic && return false
-    return ax1 == ax2
-end
-
 function TensorAlgebra.check_input(
         f::typeof(TensorAlgebra.contract),
         a1::AbstractGradedArray, perm1_codomain, perm1_domain,
@@ -238,10 +227,12 @@ function TensorAlgebra.check_input(
         a1::AbstractArray, perm1_codomain, perm1_domain,
         a2::AbstractArray, perm2_codomain, perm2_domain
     )
+    # Contracted axes must be a canonical dual pair (`dual(ax1) == ax2`), so a
+    # contraction always pairs a space with its dual, for every sector type.
     for (i, j) in zip(perm1_domain, perm2_codomain)
         ax1 = axes(a1, i)
         ax2 = axes(a2, j)
-        are_axes_contractible(ax1, ax2) || throw(
+        dual(ax1) == ax2 || throw(
             ArgumentError(
                 "Contracted axes do not match: `axes(a1, $i) = $ax1` and `axes(a2, $j) = $ax2`"
             )
diff --git a/test/test_fermionic.jl b/test/test_fermionic.jl
index e950aa86..e03853ba 100644
--- a/test/test_fermionic.jl
+++ b/test/test_fermionic.jl
@@ -197,13 +197,15 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
     end
 
     @testset "matrix-matrix contraction B: total phase -1" begin
-        # labels (1,-1,-2,2) × (2,1,-3,-4): P2=-1, T1=-1, T2=+1, U=-1 → total=-1
+        # Canonical contraction (a1 codomain ↔ a2 domain). The two odd contracted
+        # indices appear in swapped order on a2, one transposition, so the
+        # contraction twists by -1 relative to the dense result.
         a1 = randn_blockdiagonal(elt, (r_odd, r_odd, dual(r_odd), dual(r_odd)))
         a2 = randn_blockdiagonal(elt, (r_odd, r_odd, dual(r_odd), dual(r_odd)))
         a1_dense = to_dense(a1)
         a2_dense = to_dense(a2)
-        a_dest, _ = contract(a1, (1, -1, -2, 2), a2, (2, 1, -3, -4))
-        a_dest_dense, _ = contract(a1_dense, (1, -1, -2, 2), a2_dense, (2, 1, -3, -4))
+        a_dest, _ = contract(a1, (1, 2, -1, -2), a2, (-3, -4, 2, 1))
+        a_dest_dense, _ = contract(a1_dense, (1, 2, -1, -2), a2_dense, (-3, -4, 2, 1))
         @test to_dense(a_dest) ≈ -1 * a_dest_dense
     end
 
@@ -219,13 +221,14 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
     end
 
     @testset "matrix-matrix contraction G: total phase -1" begin
-        # labels (1,2,-1,-2) × (2,-3,1,-4): P1=+1, T1=+1, P2=+1, T2=-1, U=+1 → total=-1
+        # Canonical contraction with the contracted indices on a1's domain (dual)
+        # side and a2's codomain side. Same single-transposition twist as case B.
         a1 = randn_blockdiagonal(elt, (r_odd, r_odd, dual(r_odd), dual(r_odd)))
         a2 = randn_blockdiagonal(elt, (r_odd, r_odd, dual(r_odd), dual(r_odd)))
         a1_dense = to_dense(a1)
         a2_dense = to_dense(a2)
-        a_dest, _ = contract(a1, (1, 2, -1, -2), a2, (2, -3, 1, -4))
-        a_dest_dense, _ = contract(a1_dense, (1, 2, -1, -2), a2_dense, (2, -3, 1, -4))
+        a_dest, _ = contract(a1, (-1, -2, 1, 2), a2, (2, 1, -3, -4))
+        a_dest_dense, _ = contract(a1_dense, (-1, -2, 1, 2), a2_dense, (2, 1, -3, -4))
         @test to_dense(a_dest) ≈ -1 * a_dest_dense
     end
 end

From a93193c314cc31b7a605db7bf9019b0ef8a648dd Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Mon, 15 Jun 2026 19:05:55 -0400
Subject: [PATCH 31/31] Drop TensorAlgebra [sources] pin, require 0.9.5

TensorAlgebra 0.9.5 is registered with the projectto! and trivialrange primitives this PR overloads, so resolve it from the registry and set it as the compat floor.
---
 Project.toml      | 6 +-----
 test/Project.toml | 2 +-
 2 files changed, 2 insertions(+), 6 deletions(-)

diff --git a/Project.toml b/Project.toml
index 0b726297..2bc7c289 100644
--- a/Project.toml
+++ b/Project.toml
@@ -29,10 +29,6 @@ TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138"
 NamedDimsArrays = "60cbd0c0-df58-4cb7-918c-6f5607b73fde"
 SUNRepresentations = "1a50b95c-7aac-476d-a9ce-2bfc675fc617"
 
-[sources.TensorAlgebra]
-rev = "mf/gram-eigh-balanced"
-url = "https://github.com/ITensor/TensorAlgebra.jl"
-
 [extensions]
 GradedArraysNamedDimsArraysExt = "NamedDimsArrays"
 GradedArraysSUNRepresentationsExt = "SUNRepresentations"
@@ -54,7 +50,7 @@ Random = "1.10"
 SUNRepresentations = "0.3, 0.4"
 SparseArraysBase = "0.9"
 SplitApplyCombine = "1.2.3"
-TensorAlgebra = "0.8, 0.9"
+TensorAlgebra = "0.9.5"
 TensorKitSectors = "0.3"
 TypeParameterAccessors = "0.4"
 julia = "1.10"
diff --git a/test/Project.toml b/test/Project.toml
index 30cd7ab5..9bde4db4 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -38,7 +38,7 @@ SUNRepresentations = "0.3, 0.4"
 SafeTestsets = "0.1"
 SparseArraysBase = "0.9"
 Suppressor = "0.2.8"
-TensorAlgebra = "0.8, 0.9"
+TensorAlgebra = "0.9.5"
 TensorKitSectors = "0.3"
 Test = "1.10"
 TestExtras = "0.3.1"