From 1ea4f1b9c87801b448fda7f7a57e553c2e79f544 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Sat, 11 Apr 2026 13:03:33 -0400 Subject: [PATCH 01/12] Add Data{N} type, undef constructors, and data/block indexing - Add `Data{N}` indexing type with `Block`/`Data` conversions - Add `SectorMatrix` undef constructors using `similar(D, axes)` - Add `FusedGradedMatrix` undef constructors (BlockedOneTo innermost) - Add `Data` view/getindex/setindex! on FusedGradedMatrix and AbelianGradedArray - Add `Block` setindex! on FusedGradedMatrix for SectorMatrix values - Restructure SectorMatrix accessors: sector/data primitive, sectoraxes/axes derived - Use `similar(D, axes)` instead of `D(undef, ...)` in all undef constructors - Use explicit type params in convert/similar, infer elsewhere - Add `SectorOneTo(::SectorRange, ::Base.OneTo)` constructor - Use `blockdiagindices` in FusedGradedMatrix blocks accessor Co-Authored-By: Claude Opus 4.6 (1M context) --- Project.toml | 2 +- src/GradedArrays.jl | 11 ++-- src/abeliangradedarray.jl | 24 ++++++++- src/abeliansectorarray.jl | 12 ++--- src/data.jl | 13 +++++ src/fusedgradedmatrix.jl | 89 ++++++++++++++++++++++++++++++- src/sectormatrix.jl | 43 +++++++++++---- src/sectoroneto.jl | 1 + test/test_abelianarray.jl | 52 +++++++++++++++++- test/test_data.jl | 108 ++++++++++++++++++++++++++++++++++++++ test/test_exports.jl | 1 + test/test_sectormatrix.jl | 24 ++++++++- 12 files changed, 351 insertions(+), 29 deletions(-) create mode 100644 src/data.jl create mode 100644 test/test_data.jl diff --git a/Project.toml b/Project.toml index 2a420a93..261be7df 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "GradedArrays" uuid = "bc96ca6e-b7c8-4bb6-888e-c93f838762c2" -version = "0.8.1" +version = "0.8.2" authors = ["ITensor developers and contributors"] [workspace] diff --git a/src/GradedArrays.jl b/src/GradedArrays.jl index 6f273a79..2cbd5489 100644 --- a/src/GradedArrays.jl +++ b/src/GradedArrays.jl @@ -11,19 +11,21 @@ export AbstractSectorArray, export AbstractGradedArray, AbstractGradedMatrix export AbelianGradedArray, AbelianGradedVector, AbelianGradedMatrix export FusedGradedMatrix -export gradedrange export dual, flip, gradedrange, isdual, data, dataaxes, dataaxes1, datalength, datalengths, sector, sectoraxes, sectoraxes1, sectorlength, sectorlengths, - sectors, sectortype + sectors, sectortype, + Data # imports # ------- import FunctionImplementations as FI using BlockArrays: BlockArrays, AbstractBlockArray, AbstractBlockVector, Block, - BlockIndexRange, BlockVector, blocklength, blocklengths, blocks, eachblockaxes1 -using BlockSparseArrays: BlockSparseArrays, eachblockaxis, eachblockstoredindex, mortar_axis + BlockIndexRange, BlockVector, BlockedOneTo, blockedrange, blocklength, blocklengths, + blocks, eachblockaxes1 +using BlockSparseArrays: + BlockSparseArrays, blockdiagindices, eachblockaxis, eachblockstoredindex, mortar_axis using KroneckerArrays: KroneckerArrays, kroneckerfactors, ×, ⊗ using LinearAlgebra: LinearAlgebra, Adjoint, mul! using TensorAlgebra: TensorAlgebra, BlockedTuple, FusionStyle, matricize, matricize_axes, @@ -33,6 +35,7 @@ using TensorKitSectors: TensorKitSectors as TKS using TypeParameterAccessors: type_parameters, unspecify_type_parameters include("sectorrange.jl") +include("data.jl") include("sectoroneto.jl") include("gradedoneto.jl") include("abstractsectordelta.jl") diff --git a/src/abeliangradedarray.jl b/src/abeliangradedarray.jl index 9d6cd673..f7835b6d 100644 --- a/src/abeliangradedarray.jl +++ b/src/abeliangradedarray.jl @@ -30,7 +30,10 @@ function AbelianGradedArray{T, N, D, S}( bks = allowedblocks(axs) blockdata = Dict{NTuple{N, Int}, D}( Int.(Tuple(bk)) => - D(undef, ntuple(d -> blocklengths(axs[d])[Int(Tuple(bk)[d])], Val(N))) + similar( + D, + ntuple(d -> Base.OneTo(blocklengths(axs[d])[Int(Tuple(bk)[d])]), Val(N)) + ) for bk in bks ) return AbelianGradedArray{T, N, D, S}(blockdata, axs) @@ -134,6 +137,25 @@ function Base.setindex!( return a end +# --------------------------------------------------------------------------- +# Data indexing — raw block data without sector wrappers +# --------------------------------------------------------------------------- + +function Base.view(a::AbelianGradedArray{T, N}, I::Data{N}) where {T, N} + return data(view(a, Block(I))) +end + +function Base.getindex(a::AbelianGradedArray{T, N}, I::Data{N}) where {T, N} + return copy(view(a, I)) +end + +function Base.setindex!( + a::AbelianGradedArray{<:Any, N}, value::AbstractArray{<:Any, N}, I::Data{N} + ) where {N} + view(a, I) .= value + return a +end + # --------------------------------------------------------------------------- # Splitting getindex: each I[d][k] = Block(b)[r] means dest block k comes # from source block b at subrange r. Inverse of the merging getindex. diff --git a/src/abeliansectorarray.jl b/src/abeliansectorarray.jl index 18e9dc45..e48726b4 100644 --- a/src/abeliansectorarray.jl +++ b/src/abeliansectorarray.jl @@ -55,9 +55,7 @@ datatype(::Type{AbelianSectorArray{T, N, A, S}}) where {T, N, A, S} = A # AbstractArray interface # ----------------------- function Base.axes(sa::AbelianSectorArray) - sa_axes = sectoraxes(sa) - da_axes = dataaxes(sa) - return ntuple(d -> SectorOneTo(sa_axes[d], length(da_axes[d])), Val(ndims(sa))) + return map(SectorOneTo, sectoraxes(sa), dataaxes(sa)) end Base.copy(A::AbelianSectorArray) = AbelianSectorArray(sector(A), copy(data(A))) @@ -65,13 +63,11 @@ Base.copy(A::AbelianSectorArray) = AbelianSectorArray(sector(A), copy(data(A))) # similar for AbelianSectorArray with SectorOneTo axes. # Delegates to similar on the data array for the data dimensions. function Base.similar( - a::AbelianSectorArray, + ::AbelianSectorArray, ::Type{T}, axes::Tuple{SectorOneTo, Vararg{SectorOneTo}} ) where {T} - sects = sector.(axes) - data_ax = data.(axes) - return AbelianSectorArray(sects, similar(data(a), T, data_ax)) + return AbelianSectorArray{T}(undef, axes) end function Base.fill!(A::AbelianSectorArray, v) @@ -88,7 +84,7 @@ function Base.convert( x::AbelianSectorArray{T₂, N, B, S} )::AbelianSectorArray{T₁, N, A, S} where {T₁, T₂, N, A, B, S} A === B && return x - return AbelianSectorArray(sector(x), convert(A, data(x))) + return AbelianSectorArray{T₁, N, A, S}(x.sectors, convert(A, data(x))) end # ======================== permutedims ======================== diff --git a/src/data.jl b/src/data.jl new file mode 100644 index 00000000..103f9fd2 --- /dev/null +++ b/src/data.jl @@ -0,0 +1,13 @@ +""" + Data{N} + +Block-data indexing type analogous to `BlockArrays.Block{N}`. Indexing a graded +array with `Data(i, j, ...)` accesses the raw data array for that block, without +sector metadata wrappers. +""" +struct Data{N} + n::NTuple{N, Int} +end +Data(n::Vararg{Int, N}) where {N} = Data{N}(n) +BlockArrays.Block(I::Data) = Block(I.n) +Data(I::Block) = Data(Int.(Tuple(I))) diff --git a/src/fusedgradedmatrix.jl b/src/fusedgradedmatrix.jl index a402ff56..7f95ae97 100644 --- a/src/fusedgradedmatrix.jl +++ b/src/fusedgradedmatrix.jl @@ -35,6 +35,57 @@ struct FusedGradedMatrix{T, D <: AbstractMatrix{T}, S <: SectorRange} <: end end +# ======================== undef constructors ======================== + +function _check_undef_args(sectors, axes::Tuple{BlockedOneTo, BlockedOneTo}) + n = length(sectors) + length(blocklengths(axes[1])) == n && length(blocklengths(axes[2])) == n || + throw( + ArgumentError( + "sectors length ($n) must match block count in axes" + ) + ) + return nothing +end + +# Innermost: fully parameterized, takes (BlockedOneTo, BlockedOneTo) axes. +function FusedGradedMatrix{T, D, S}( + ::UndefInitializer, + sectors::Vector{S}, + axes::Tuple{BlockedOneTo, BlockedOneTo} + ) where {T, D <: AbstractMatrix{T}, S <: SectorRange} + _check_undef_args(sectors, axes) + cod_bls = blocklengths(axes[1]) + dom_bls = blocklengths(axes[2]) + blks = [ + similar(D, (Base.OneTo(cod_bls[i]), Base.OneTo(dom_bls[i]))) + for i in eachindex(sectors) + ] + return FusedGradedMatrix(sectors, blks) +end + +# Convenience: default D = Matrix{T}. +function FusedGradedMatrix{T}( + ::UndefInitializer, + sectors::Vector{S}, + axes::Tuple{BlockedOneTo, BlockedOneTo} + ) where {T, S <: SectorRange} + return FusedGradedMatrix{T, Matrix{T}, S}(undef, sectors, axes) +end + +# Vector{Int} convenience: wraps into BlockedOneTo and delegates. +function FusedGradedMatrix{T}( + ::UndefInitializer, + sectors::Vector{<:SectorRange}, + codomain_blocklengths::Vector{Int}, + domain_blocklengths::Vector{Int} + ) where {T} + return FusedGradedMatrix{T}( + undef, sectors, + (blockedrange(codomain_blocklengths), blockedrange(domain_blocklengths)) + ) +end + # ======================== Accessors ======================== BlockArrays.blocklength(m::FusedGradedMatrix) = length(m.sectors) @@ -73,6 +124,40 @@ function Base.getindex(m::FusedGradedMatrix, i::Block{1}, j::Block{1}) return m[Block(Int(i), Int(j))] end +function Base.setindex!(m::FusedGradedMatrix, value::SectorMatrix, I::Block{2}) + i, j = Int.(Tuple(I)) + i == j || + error("Off-diagonal access not supported for block-diagonal FusedGradedMatrix") + value.sector == m.sectors[i] || + error("Sector mismatch: block has $(m.sectors[i]) but value has $(value.sector)") + size(data(value)) == size(m.blocks[i]) || + throw( + DimensionMismatch( + "data size $(size(data(value))) does not match block size $(size(m.blocks[i]))" + ) + ) + copyto!(m.blocks[i], data(value)) + return m +end +function Base.setindex!(m::FusedGradedMatrix, value::SectorMatrix, i::Block{1}, j::Block{1}) + return setindex!(m, value, Block(Int(i), Int(j))) +end + +# ======================== Data indexing ======================== + +function Base.view(m::FusedGradedMatrix, I::Data{2}) + return data(view(m, Block(I))) +end + +function Base.getindex(m::FusedGradedMatrix, I::Data{2}) + return copy(view(m, I)) +end + +function Base.setindex!(m::FusedGradedMatrix, value::AbstractMatrix, I::Data{2}) + view(m, I) .= value + return m +end + # ======================== eachblockstoredindex ======================== function BlockSparseArrays.eachblockstoredindex(m::FusedGradedMatrix) @@ -83,8 +168,8 @@ end using LinearAlgebra: Diagonal function BlockArrays.blocks(m::FusedGradedMatrix) - sector_blocks = [SectorMatrix(s, b) for (s, b) in zip(m.sectors, m.blocks)] - return Diagonal(sector_blocks) + diagblocks = map(b -> view(m, b), blockdiagindices(m)) + return Diagonal(collect(diagblocks)) end # ======================== fill! / zero! ======================== diff --git a/src/sectormatrix.jl b/src/sectormatrix.jl index 469ac201..1da6796b 100644 --- a/src/sectormatrix.jl +++ b/src/sectormatrix.jl @@ -16,26 +16,45 @@ struct SectorMatrix{T, D <: AbstractMatrix{T}, S <: SectorRange} <: data::D end -# ---- accessors ---- +# ---- undef constructors ---- + +# Innermost: fully parameterized, takes AbstractUnitRange axes. +function SectorMatrix{T, D, S}( + ::UndefInitializer, sector::S, r1::AbstractUnitRange, r2::AbstractUnitRange + ) where {T, D <: AbstractMatrix{T}, S <: SectorRange} + return SectorMatrix{T, D, S}(sector, similar(D, (r1, r2))) +end + +# Convenience: default D = Matrix{T}. +function SectorMatrix{T}( + ::UndefInitializer, sector::S, r1::AbstractUnitRange, r2::AbstractUnitRange + ) where {T, S <: SectorRange} + return SectorMatrix{T, Matrix{T}, S}(undef, sector, r1, r2) +end -function sectoraxes(sm::SectorMatrix) - return (sm.sector, dual(sm.sector)) +# Int convenience: maps to Base.OneTo. +function SectorMatrix{T}( + ::UndefInitializer, sector::SectorRange, m::Int, n::Int + ) where {T} + return SectorMatrix{T}(undef, sector, Base.OneTo(m), Base.OneTo(n)) end -# Kronecker factor decomposition: SectorMatrix = sector ⊗ data +# ---- accessors ---- + +# Primitive accessors: sector(sm) and data(sm). # sector() returns the structural delta factor (SectorIdentity), not the stored SectorRange. # Access the stored SectorRange via sm.sector or sectoraxes(sm)[1]. sector(sm::SectorMatrix) = SectorIdentity{eltype(sm)}(sm.sector) dataaxes(sm::SectorMatrix) = axes(data(sm)) +# Derived accessors: sectoraxes and axes are written in terms of sector and data. +sectoraxes(sm::SectorMatrix) = axes(sector(sm)) + sectortype(::Type{<:SectorMatrix{T, D, S}}) where {T, D, S} = S datatype(::Type{SectorMatrix{T, D, S}}) where {T, D, S} = D function Base.axes(sm::SectorMatrix) - return ( - SectorOneTo(sm.sector, size(data(sm), 1)), - SectorOneTo(dual(sm.sector), size(data(sm), 2)), - ) + return map(SectorOneTo, sectoraxes(sm), dataaxes(sm)) end Base.copy(sm::SectorMatrix) = SectorMatrix(sm.sector, copy(data(sm))) @@ -50,11 +69,13 @@ function Base.convert( x::SectorMatrix{T₂, E, S} )::SectorMatrix{T₁, D, S} where {T₁, T₂, D, E, S} D === E && return x - return SectorMatrix(x.sector, convert(D, data(x))) + return SectorMatrix{T₁, D, S}(x.sector, convert(D, data(x))) end -function Base.similar(sm::SectorMatrix, ::Type{T}) where {T} - return SectorMatrix(sm.sector, similar(data(sm), T)) +function Base.similar(sm::SectorMatrix{<:Any, <:Any, S}, ::Type{T}) where {T, S} + new_data = similar(data(sm), T) + D = typeof(new_data) + return SectorMatrix{T, D, S}(sm.sector, new_data) end function KroneckerArrays.:(⊗)(s::SectorIdentity, data::AbstractMatrix) diff --git a/src/sectoroneto.jl b/src/sectoroneto.jl index 3126b2b8..22fcfd7f 100644 --- a/src/sectoroneto.jl +++ b/src/sectoroneto.jl @@ -14,6 +14,7 @@ end # Convenience: SectorRange with default data length SectorOneTo(s::SectorRange) = SectorOneTo(s, 1) +SectorOneTo(s::SectorRange, r::Base.OneTo) = SectorOneTo(s, last(r)) # Primitive accessors sector(r::SectorOneTo) = r.sector diff --git a/test/test_abelianarray.jl b/test/test_abelianarray.jl index 8d00ddfb..260a1b96 100644 --- a/test/test_abelianarray.jl +++ b/test/test_abelianarray.jl @@ -1,4 +1,4 @@ -using BlockArrays: BlockArrays, Block, blocklength +using BlockArrays: BlockArrays, Block, blockedrange, blocklength using BlockSparseArrays: eachblockstoredindex using GradedArrays: GradedArrays, AbelianGradedArray, AbelianSectorArray, AbstractGradedArray, FusedGradedMatrix, GradedOneTo, SU2, SectorRange, U1, data, @@ -213,3 +213,53 @@ using Test: @test, @test_throws, @testset @test all(iszero, a.blockdata[(1, 1)]) end end + +@testset "FusedGradedMatrix undef constructor" begin + sectors = [U1(0), U1(1)] + cod_bls = [2, 3] + dom_bls = [1, 2] + + @testset "Convenience constructor (defaults D = Matrix{T})" begin + m = FusedGradedMatrix{Float64}(undef, sectors, cod_bls, dom_bls) + @test m isa FusedGradedMatrix{Float64, Matrix{Float64}, U1} + @test length(m.sectors) == 2 + @test m.sectors == sectors + @test size(m.blocks[1]) == (2, 1) + @test size(m.blocks[2]) == (3, 2) + end + + @testset "Fully parameterized constructor" begin + m = FusedGradedMatrix{Float64, Matrix{Float64}, U1}( + undef, sectors, (blockedrange(cod_bls), blockedrange(dom_bls)) + ) + @test m isa FusedGradedMatrix{Float64, Matrix{Float64}, U1} + @test size(m.blocks[1]) == (2, 1) + end + + @testset "Tuple BlockedOneTo form" begin + m = FusedGradedMatrix{Float64}( + undef, sectors, (blockedrange([2, 3]), blockedrange([1, 2])) + ) + @test m isa FusedGradedMatrix{Float64, Matrix{Float64}, U1} + @test size(m.blocks[1]) == (2, 1) + @test size(m.blocks[2]) == (3, 2) + end + + @testset "Rejects mismatched lengths" begin + @test_throws ArgumentError FusedGradedMatrix{Float64}( + undef, sectors, cod_bls, [1] + ) + end + + @testset "Rejects unsorted sectors" begin + @test_throws ArgumentError FusedGradedMatrix{Float64}( + undef, [U1(1), U1(0)], cod_bls, dom_bls + ) + end + + @testset "Rejects non-unique sectors" begin + @test_throws ArgumentError FusedGradedMatrix{Float64}( + undef, [U1(0), U1(0)], [2, 3], [1, 2] + ) + end +end diff --git a/test/test_data.jl b/test/test_data.jl new file mode 100644 index 00000000..d742914e --- /dev/null +++ b/test/test_data.jl @@ -0,0 +1,108 @@ +using BlockArrays: Block, blockedrange +using BlockSparseArrays: eachblockstoredindex +using GradedArrays: GradedArrays, AbelianGradedArray, AbelianSectorArray, Data, + FusedGradedMatrix, GradedOneTo, SectorMatrix, U1, data, dual, gradedrange, sectoraxes, + sectors +using Test: @test, @test_throws, @testset + +@testset "Data indexing" begin + @testset "FusedGradedMatrix" begin + sectors = [U1(0), U1(1)] + m = FusedGradedMatrix(sectors, [ones(2, 3), 2 * ones(4, 5)]) + + @testset "Data getindex returns copy of raw block" begin + d = m[Data(1, 1)] + @test d isa Matrix{Float64} + @test d == ones(2, 3) + # Verify it's a copy, not a view + d[1, 1] = 999.0 + @test m.blocks[1][1, 1] == 1.0 + end + + @testset "Data getindex second block" begin + d = m[Data(2, 2)] + @test d == 2 * ones(4, 5) + end + + @testset "Data setindex! copies into raw block" begin + sectors2 = [U1(0), U1(1)] + m2 = FusedGradedMatrix(sectors2, [zeros(2, 3), zeros(4, 5)]) + new_data = 7 * ones(2, 3) + m2[Data(1, 1)] = new_data + @test m2.blocks[1] == new_data + # Verify it's a copy, not aliased + new_data[1, 1] = 0.0 + @test m2.blocks[1][1, 1] == 7.0 + end + + @testset "Data setindex! size mismatch errors" begin + m3 = FusedGradedMatrix([U1(0)], [zeros(2, 3)]) + @test_throws DimensionMismatch (m3[Data(1, 1)] = ones(5, 5)) + end + + @testset "Data off-diagonal errors" begin + @test_throws ErrorException m[Data(1, 2)] + @test_throws ErrorException (m[Data(1, 2)] = ones(2, 5)) + end + + @testset "Block setindex! with SectorMatrix" begin + sectors2 = [U1(0), U1(1)] + m4 = FusedGradedMatrix(sectors2, [zeros(2, 3), zeros(4, 5)]) + sm = SectorMatrix(U1(0), 7 * ones(2, 3)) + m4[Block(1, 1)] = sm + @test m4.blocks[1] == 7 * ones(2, 3) + end + + @testset "Block setindex! verifies sector" begin + sectors2 = [U1(0), U1(1)] + m5 = FusedGradedMatrix(sectors2, [zeros(2, 3), zeros(4, 5)]) + sm_wrong = SectorMatrix(U1(2), ones(2, 3)) + @test_throws ErrorException (m5[Block(1, 1)] = sm_wrong) + end + + @testset "Block setindex! off-diagonal errors" begin + sectors2 = [U1(0), U1(1)] + m6 = FusedGradedMatrix(sectors2, [zeros(2, 3), zeros(4, 5)]) + sm = SectorMatrix(U1(0), ones(2, 3)) + @test_throws ErrorException (m6[Block(1, 2)] = sm) + end + end + + @testset "AbelianGradedArray" begin + g1 = gradedrange([U1(0) => 2, U1(1) => 3]) + g2 = gradedrange([U1(0) => 1, U1(-1) => 2]) + a = AbelianGradedArray{Float64}(undef, g1, g2) + a[Block(1, 1)] = ones(2, 1) + a[Block(2, 2)] = 2 * ones(3, 2) + + @testset "Data getindex returns copy of raw block" begin + d = a[Data(1, 1)] + @test d isa Matrix{Float64} + @test d == ones(2, 1) + # Verify it's a copy + d[1, 1] = 999.0 + @test data(a[Block(1, 1)]) == ones(2, 1) + end + + @testset "Data getindex second block" begin + d = a[Data(2, 2)] + @test d == 2 * ones(3, 2) + end + + @testset "Data setindex! copies into raw block" begin + a2 = AbelianGradedArray{Float64}(undef, g1, g2) + a2[Block(1, 1)] = zeros(2, 1) + new_data = 5 * ones(2, 1) + a2[Data(1, 1)] = new_data + @test data(a2[Block(1, 1)]) == 5 * ones(2, 1) + # Verify it's a copy + new_data[1, 1] = 0.0 + @test data(a2[Block(1, 1)])[1, 1] == 5.0 + end + + @testset "Data access on unstored block errors" begin + @test_throws Exception a[Data(1, 2)] + @test_throws Exception (a[Data(1, 2)] = ones(2, 2)) + end + end +end diff --git a/test/test_exports.jl b/test/test_exports.jl index 4217e081..7ae2f09c 100644 --- a/test/test_exports.jl +++ b/test/test_exports.jl @@ -13,6 +13,7 @@ using Test: @test, @testset :AbelianSectorDelta, :AbelianSectorMatrix, :AbelianSectorVector, + :Data, :FusedGradedMatrix, :Fib, :GradedArrays, diff --git a/test/test_sectormatrix.jl b/test/test_sectormatrix.jl index b06113ae..0020e7bd 100644 --- a/test/test_sectormatrix.jl +++ b/test/test_sectormatrix.jl @@ -2,7 +2,7 @@ using GradedArrays: GradedArrays, AbelianSectorArray, SU2, SectorIdentity, Secto SectorOneTo, SectorRange, U1, data, dataaxes, dual, isdual, sector, sectoraxes, sectortype, ⊗ using TensorKitSectors: TensorKitSectors as TKS -using Test: @test, @testset +using Test: @test, @test_throws, @testset @testset "SectorMatrix" begin @testset "Construction from SectorRange + data" begin @@ -136,4 +136,26 @@ using Test: @test, @testset sa .= sm @test data(sa) ≈ data(sm) end + + @testset "Undef constructor (Int dims)" begin + sm = SectorMatrix{Float64}(undef, U1(0), 3, 4) + @test sm isa SectorMatrix{Float64, Matrix{Float64}, U1} + @test size(data(sm)) == (3, 4) + @test sectoraxes(sm, 1) == U1(0) + end + + @testset "Undef constructor (AbstractUnitRange dims)" begin + sm = SectorMatrix{Float64}(undef, U1(1), Base.OneTo(2), Base.OneTo(5)) + @test sm isa SectorMatrix{Float64, Matrix{Float64}, U1} + @test size(data(sm)) == (2, 5) + @test sectoraxes(sm, 1) == U1(1) + end + + @testset "Undef constructor (fully parameterized)" begin + sm = SectorMatrix{Float64, Matrix{Float64}, U1}( + undef, U1(0), Base.OneTo(3), Base.OneTo(4) + ) + @test sm isa SectorMatrix{Float64, Matrix{Float64}, U1} + @test size(data(sm)) == (3, 4) + end end From 6d3a0e21f6897639021f4bed221a852be7478ec0 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Sat, 11 Apr 2026 13:09:37 -0400 Subject: [PATCH 02/12] Clean up AbelianGradedArray undef constructor with eachblockaxis Co-Authored-By: Claude Opus 4.6 (1M context) --- src/abeliangradedarray.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/abeliangradedarray.jl b/src/abeliangradedarray.jl index f7835b6d..123276fa 100644 --- a/src/abeliangradedarray.jl +++ b/src/abeliangradedarray.jl @@ -27,14 +27,14 @@ end function AbelianGradedArray{T, N, D, S}( ::UndefInitializer, axs::NTuple{N, GradedOneTo{S}} ) where {T, N, D <: AbstractArray{T, N}, S <: SectorRange} + block_axes = map(ax -> data.(eachblockaxis(ax)), axs) + function allocate_block(bk) + bk_inds = Int.(Tuple(bk)) + return similar(D, ntuple(d -> block_axes[d][bk_inds[d]], Val(N))) + end bks = allowedblocks(axs) blockdata = Dict{NTuple{N, Int}, D}( - Int.(Tuple(bk)) => - similar( - D, - ntuple(d -> Base.OneTo(blocklengths(axs[d])[Int(Tuple(bk)[d])]), Val(N)) - ) - for bk in bks + Int.(Tuple(bk)) => allocate_block(bk) for bk in bks ) return AbelianGradedArray{T, N, D, S}(blockdata, axs) end From 3ea55952e0fd22fd1359ee6bb395cfb0fdd7c17e Mon Sep 17 00:00:00 2001 From: mtfishman Date: Sat, 11 Apr 2026 13:11:43 -0400 Subject: [PATCH 03/12] Add parameterized AbelianSectorArray delta constructor, Tuple(::Data), use accessor in convert Co-Authored-By: Claude Opus 4.6 (1M context) --- src/abeliansectorarray.jl | 8 +++++++- src/data.jl | 1 + 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/src/abeliansectorarray.jl b/src/abeliansectorarray.jl index e48726b4..d672e7dd 100644 --- a/src/abeliansectorarray.jl +++ b/src/abeliansectorarray.jl @@ -36,6 +36,12 @@ function AbelianSectorArray( ) where {N} return AbelianSectorArray(delta.sectors, data) end +function AbelianSectorArray{T, N, A, S}( + delta::AbelianSectorDelta{<:Any, N, S}, + data::A + ) where {T, N, A <: AbstractArray{T, N}, S <: SectorRange} + return AbelianSectorArray{T, N, A, S}(delta.sectors, data) +end const AbelianSectorVector{T, A <: AbstractVector{T}, S <: SectorRange} = AbelianSectorArray{T, 1, A, S} @@ -84,7 +90,7 @@ function Base.convert( x::AbelianSectorArray{T₂, N, B, S} )::AbelianSectorArray{T₁, N, A, S} where {T₁, T₂, N, A, B, S} A === B && return x - return AbelianSectorArray{T₁, N, A, S}(x.sectors, convert(A, data(x))) + return AbelianSectorArray{T₁, N, A, S}(sector(x), convert(A, data(x))) end # ======================== permutedims ======================== diff --git a/src/data.jl b/src/data.jl index 103f9fd2..3795349e 100644 --- a/src/data.jl +++ b/src/data.jl @@ -11,3 +11,4 @@ end Data(n::Vararg{Int, N}) where {N} = Data{N}(n) BlockArrays.Block(I::Data) = Block(I.n) Data(I::Block) = Data(Int.(Tuple(I))) +Base.Tuple(I::Data) = I.n From 1e3ea4578fd2259ef1637498eb7a34d081ca9f1f Mon Sep 17 00:00:00 2001 From: mtfishman Date: Sat, 11 Apr 2026 13:34:32 -0400 Subject: [PATCH 04/12] Refactor: generic Block/Data interface, FusedGradedMatrix inner constructor, cleanup - Move Block/Data view/getindex/setindex! to AbstractGradedArray generic interface - Concrete types only implement view(a, ::Block{N}) as primitive - Add parameterized FusedGradedMatrix{T,D,S} inner constructor - Use AbstractGradedMatrix{T} as supertype for FusedGradedMatrix - Simplify FusedGradedMatrix undef constructor with eachblockaxis - Add AbelianSectorArray{T,N,A,S}(::AbelianSectorDelta, data) constructor - Use accessor in AbelianSectorArray convert - Add Tuple(::Data), SectorOneTo(::SectorRange, ::Base.OneTo) - Use view(m, I) .= value for Block setindex! on FusedGradedMatrix - Clean up AbelianGradedArray undef constructor with eachblockaxis Co-Authored-By: Claude Opus 4.6 (1M context) --- src/abeliangradedarray.jl | 53 +---------------------- src/abstractgradedarray.jl | 51 ++++++++++++++++++++++ src/fusedgradedmatrix.jl | 86 ++++++++------------------------------ test/test_abelianarray.jl | 2 +- test/test_data.jl | 4 +- 5 files changed, 73 insertions(+), 123 deletions(-) diff --git a/src/abeliangradedarray.jl b/src/abeliangradedarray.jl index 123276fa..4ed696dc 100644 --- a/src/abeliangradedarray.jl +++ b/src/abeliangradedarray.jl @@ -69,16 +69,12 @@ BlockSparseArrays.blocktype(a::AbelianGradedArray) = BlockSparseArrays.blocktype # view (primitive): returns AbelianSectorArray sharing data with blockdata # --------------------------------------------------------------------------- -# view: returns a AbelianSectorArray sharing data (errors for unstored blocks) -function Base.view(a::AbelianGradedArray{T, N}, I::Vararg{Block{1}, N}) where {T, N} - bk = ntuple(d -> Int(I[d]), Val(N)) +function Base.view(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 -function Base.view(a::AbelianGradedArray{T, N}, I::Block{N}) where {T, N} - return view(a, Tuple(I)...) -end # --------------------------------------------------------------------------- # blocks — lazy view delegating to view (following BlockArrays convention) @@ -111,51 +107,6 @@ function Base.setindex!( return b end -# --------------------------------------------------------------------------- -# getindex / setindex! on AbelianGradedArray with Block — convenience wrappers -# --------------------------------------------------------------------------- - -# getindex: returns a copy (errors for unstored blocks) -function Base.getindex(a::AbelianGradedArray{T, N}, I::Vararg{Block{1}, N}) where {T, N} - return copy(view(a, I...)) -end -function Base.getindex(a::AbelianGradedArray{T, N}, I::Block{N}) where {T, N} - return a[Tuple(I)...] -end - -# setindex!: Block{N} unpacks to Vararg{Block{1}, N} (following BlockArrays convention) -function Base.setindex!(a::AbelianGradedArray{<:Any, N}, value, I::Block{N}) where {N} - return setindex!(a, value, Tuple(I)...) -end - -# Primitive: view existing block, then broadcast in. -# Handles both AbelianSectorArray and raw data values. -function Base.setindex!( - a::AbelianGradedArray{<:Any, N}, value::AbstractArray{<:Any, N}, I::Vararg{Block{1}, N} - ) where {N} - view(a, I...) .= value - return a -end - -# --------------------------------------------------------------------------- -# Data indexing — raw block data without sector wrappers -# --------------------------------------------------------------------------- - -function Base.view(a::AbelianGradedArray{T, N}, I::Data{N}) where {T, N} - return data(view(a, Block(I))) -end - -function Base.getindex(a::AbelianGradedArray{T, N}, I::Data{N}) where {T, N} - return copy(view(a, I)) -end - -function Base.setindex!( - a::AbelianGradedArray{<:Any, N}, value::AbstractArray{<:Any, N}, I::Data{N} - ) where {N} - view(a, I) .= value - return a -end - # --------------------------------------------------------------------------- # Splitting getindex: each I[d][k] = Block(b)[r] means dest block k comes # from source block b at subrange r. Inverse of the merging getindex. diff --git a/src/abstractgradedarray.jl b/src/abstractgradedarray.jl index 32a52158..f89d16f4 100644 --- a/src/abstractgradedarray.jl +++ b/src/abstractgradedarray.jl @@ -26,3 +26,54 @@ function Base.setindex!(::AbstractGradedArray, _, ::Vararg{Int}) "Scalar indexing is not supported for AbstractGradedArray. Use block indexing." ) end + +# --------------------------------------------------------------------------- +# Block indexing interface +# +# Concrete subtypes must implement: +# view(a::ConcreteType, ::Block{N}) → sector-wrapped view (e.g. SectorMatrix) +# +# Everything else is derived here. +# --------------------------------------------------------------------------- + +function Base.view(a::AbstractGradedArray{T, N}, I::Vararg{Block{1}, N}) where {T, N} + return view(a, Block(Int.(I))) +end + +function Base.getindex(a::AbstractGradedArray{T, N}, I::Block{N}) where {T, N} + return copy(view(a, I)) +end +function Base.getindex(a::AbstractGradedArray{T, N}, I::Vararg{Block{1}, N}) where {T, N} + return a[Block(Int.(I))] +end + +function Base.setindex!(a::AbstractGradedArray{<:Any, N}, value, I::Block{N}) where {N} + return setindex!(a, value, Tuple(I)...) +end +function Base.setindex!( + a::AbstractGradedArray{<:Any, N}, value, I::Vararg{Block{1}, N} + ) where {N} + view(a, I...) .= value + return a +end + +# --------------------------------------------------------------------------- +# Data indexing — raw block data without sector wrappers +# +# Built on top of Block view: view(a, Data(I)) = data(view(a, Block(I))) +# --------------------------------------------------------------------------- + +function Base.view(a::AbstractGradedArray{T, N}, I::Data{N}) where {T, N} + return data(view(a, Block(I))) +end + +function Base.getindex(a::AbstractGradedArray{T, N}, I::Data{N}) where {T, N} + return copy(view(a, I)) +end + +function Base.setindex!( + a::AbstractGradedArray{<:Any, N}, value::AbstractArray{<:Any, N}, I::Data{N} + ) where {N} + view(a, I) .= value + return a +end diff --git a/src/fusedgradedmatrix.jl b/src/fusedgradedmatrix.jl index 7f95ae97..f048b9d8 100644 --- a/src/fusedgradedmatrix.jl +++ b/src/fusedgradedmatrix.jl @@ -18,13 +18,13 @@ derived from `size(blocks[i], 1)`. The domain (column) axis is dual with sectors `dual(sectors[i])` and sector lengths from `size(blocks[i], 2)`. """ struct FusedGradedMatrix{T, D <: AbstractMatrix{T}, S <: SectorRange} <: - AbstractGradedArray{T, 2} + AbstractGradedMatrix{T} sectors::Vector{S} blocks::Vector{D} - function FusedGradedMatrix( + function FusedGradedMatrix{T, D, S}( sectors::Vector{S}, blocks::Vector{D} - ) where {T, S <: SectorRange, D <: AbstractMatrix{T}} + ) where {T, D <: AbstractMatrix{T}, S <: SectorRange} length(sectors) == length(blocks) || throw(ArgumentError("sectors and blocks must have the same length")) issorted(sectors) || @@ -34,34 +34,26 @@ struct FusedGradedMatrix{T, D <: AbstractMatrix{T}, S <: SectorRange} <: return new{T, D, S}(sectors, blocks) end end +function FusedGradedMatrix( + sectors::Vector{S}, + blocks::Vector{D} + ) where {T, S <: SectorRange, D <: AbstractMatrix{T}} + return FusedGradedMatrix{T, D, S}(sectors, blocks) +end # ======================== undef constructors ======================== -function _check_undef_args(sectors, axes::Tuple{BlockedOneTo, BlockedOneTo}) - n = length(sectors) - length(blocklengths(axes[1])) == n && length(blocklengths(axes[2])) == n || - throw( - ArgumentError( - "sectors length ($n) must match block count in axes" - ) - ) - return nothing -end - -# Innermost: fully parameterized, takes (BlockedOneTo, BlockedOneTo) axes. function FusedGradedMatrix{T, D, S}( ::UndefInitializer, sectors::Vector{S}, axes::Tuple{BlockedOneTo, BlockedOneTo} ) where {T, D <: AbstractMatrix{T}, S <: SectorRange} - _check_undef_args(sectors, axes) - cod_bls = blocklengths(axes[1]) - dom_bls = blocklengths(axes[2]) - blks = [ - similar(D, (Base.OneTo(cod_bls[i]), Base.OneTo(dom_bls[i]))) - for i in eachindex(sectors) - ] - return FusedGradedMatrix(sectors, blks) + cod_axes = eachblockaxis(axes[1]) + dom_axes = eachblockaxis(axes[2]) + length(cod_axes) == length(dom_axes) == length(sectors) || + throw(ArgumentError("axes block counts must match sectors length")) + blks = [similar(D, (cod_axes[i], dom_axes[i])) for i in eachindex(sectors)] + return FusedGradedMatrix{T, D, S}(sectors, blks) end # Convenience: default D = Matrix{T}. @@ -82,7 +74,7 @@ function FusedGradedMatrix{T}( ) where {T} return FusedGradedMatrix{T}( undef, sectors, - (blockedrange(codomain_blocklengths), blockedrange(domain_blocklengths)) + blockedrange.((codomain_blocklengths, domain_blocklengths)) ) end @@ -105,7 +97,7 @@ Base.size(m::FusedGradedMatrix) = map(length, axes(m)) Base.eltype(::Type{FusedGradedMatrix{T}}) where {T} = T Base.eltype(::Type{<:FusedGradedMatrix{T}}) where {T} = T -# ======================== Block indexing ======================== +# ======================== Block indexing (primitive) ======================== function Base.view(m::FusedGradedMatrix, I::Block{2}) i, j = Int.(Tuple(I)) @@ -113,50 +105,6 @@ function Base.view(m::FusedGradedMatrix, I::Block{2}) error("Off-diagonal access not supported for block-diagonal FusedGradedMatrix") return SectorMatrix(m.sectors[i], m.blocks[i]) end -function Base.view(m::FusedGradedMatrix, i::Block{1}, j::Block{1}) - return view(m, Block(Int(i), Int(j))) -end - -function Base.getindex(m::FusedGradedMatrix, I::Block{2}) - return copy(view(m, I)) -end -function Base.getindex(m::FusedGradedMatrix, i::Block{1}, j::Block{1}) - return m[Block(Int(i), Int(j))] -end - -function Base.setindex!(m::FusedGradedMatrix, value::SectorMatrix, I::Block{2}) - i, j = Int.(Tuple(I)) - i == j || - error("Off-diagonal access not supported for block-diagonal FusedGradedMatrix") - value.sector == m.sectors[i] || - error("Sector mismatch: block has $(m.sectors[i]) but value has $(value.sector)") - size(data(value)) == size(m.blocks[i]) || - throw( - DimensionMismatch( - "data size $(size(data(value))) does not match block size $(size(m.blocks[i]))" - ) - ) - copyto!(m.blocks[i], data(value)) - return m -end -function Base.setindex!(m::FusedGradedMatrix, value::SectorMatrix, i::Block{1}, j::Block{1}) - return setindex!(m, value, Block(Int(i), Int(j))) -end - -# ======================== Data indexing ======================== - -function Base.view(m::FusedGradedMatrix, I::Data{2}) - return data(view(m, Block(I))) -end - -function Base.getindex(m::FusedGradedMatrix, I::Data{2}) - return copy(view(m, I)) -end - -function Base.setindex!(m::FusedGradedMatrix, value::AbstractMatrix, I::Data{2}) - view(m, I) .= value - return m -end # ======================== eachblockstoredindex ======================== diff --git a/test/test_abelianarray.jl b/test/test_abelianarray.jl index 260a1b96..e55e1642 100644 --- a/test/test_abelianarray.jl +++ b/test/test_abelianarray.jl @@ -246,7 +246,7 @@ end end @testset "Rejects mismatched lengths" begin - @test_throws ArgumentError FusedGradedMatrix{Float64}( + @test_throws Exception FusedGradedMatrix{Float64}( undef, sectors, cod_bls, [1] ) end diff --git a/test/test_data.jl b/test/test_data.jl index d742914e..c07b5a3f 100644 --- a/test/test_data.jl +++ b/test/test_data.jl @@ -57,14 +57,14 @@ using Test: @test, @test_throws, @testset sectors2 = [U1(0), U1(1)] m5 = FusedGradedMatrix(sectors2, [zeros(2, 3), zeros(4, 5)]) sm_wrong = SectorMatrix(U1(2), ones(2, 3)) - @test_throws ErrorException (m5[Block(1, 1)] = sm_wrong) + @test_throws DimensionMismatch (m5[Block(1, 1)] = sm_wrong) end @testset "Block setindex! off-diagonal errors" begin sectors2 = [U1(0), U1(1)] m6 = FusedGradedMatrix(sectors2, [zeros(2, 3), zeros(4, 5)]) sm = SectorMatrix(U1(0), ones(2, 3)) - @test_throws ErrorException (m6[Block(1, 2)] = sm) + @test_throws Exception (m6[Block(1, 2)] = sm) end end From 15d9e208d2684fc9774b1d699dc1f2bef66398a5 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Sat, 11 Apr 2026 13:46:58 -0400 Subject: [PATCH 05/12] Refactor FusedGradedMatrix mul!, *, and conversions to use Data/Block interface MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Add check_mul_axes for FusedGradedMatrix - Use blockdiagindices + view(m, Data(I)) in mul! and * - Rewrite FusedGradedMatrix(::AbelianGradedMatrix) with undef constructor + Data indexing - Rewrite AbelianGradedArray(::FusedGradedMatrix) with Data indexing - Fix Array(a[bI]) → copy(view(a, Data(bI))) - Use eachblockstoredindex in blocks (blockdiagindices depends on blocks) Co-Authored-By: Claude Opus 4.6 (1M context) --- src/fusedgradedmatrix.jl | 45 +++++++++++++++++++++------------------- 1 file changed, 24 insertions(+), 21 deletions(-) diff --git a/src/fusedgradedmatrix.jl b/src/fusedgradedmatrix.jl index f048b9d8..6452e915 100644 --- a/src/fusedgradedmatrix.jl +++ b/src/fusedgradedmatrix.jl @@ -116,7 +116,7 @@ end using LinearAlgebra: Diagonal function BlockArrays.blocks(m::FusedGradedMatrix) - diagblocks = map(b -> view(m, b), blockdiagindices(m)) + diagblocks = map(I -> view(m, I), eachblockstoredindex(m)) return Diagonal(collect(diagblocks)) end @@ -138,22 +138,31 @@ end # ======================== mul! ======================== +function check_mul_axes( + C::FusedGradedMatrix, A::FusedGradedMatrix, B::FusedGradedMatrix + ) + axes(A, 2) == dual(axes(B, 1)) || + throw(DimensionMismatch("sector mismatch in contracted dimension")) + axes(C, 1) == axes(A, 1) || throw(DimensionMismatch()) + axes(C, 2) == axes(B, 2) || throw(DimensionMismatch()) + return nothing +end + function LinearAlgebra.mul!( C::FusedGradedMatrix, A::FusedGradedMatrix, B::FusedGradedMatrix, α::Number, β::Number ) - C.sectors == A.sectors == B.sectors || - throw(DimensionMismatch("FusedGradedMatrix sectors must match")) - for i in eachindex(C.blocks) - mul!(C.blocks[i], A.blocks[i], B.blocks[i], α, β) + check_mul_axes(C, A, B) + for I in blockdiagindices(C) + mul!(view(C, Data(I)), view(A, Data(I)), view(B, Data(I)), α, β) end return C end function Base.:(*)(A::FusedGradedMatrix{T₁}, B::FusedGradedMatrix{T₂}) where {T₁, T₂} - A.sectors == B.sectors || throw(DimensionMismatch("sectors must match")) + axes(A, 2) == dual(axes(B, 1)) || throw(DimensionMismatch("sectors must match")) T = Base.promote_op(LinearAlgebra.matprod, T₁, T₂) - result_blocks = [A.blocks[i] * B.blocks[i] for i in eachindex(A.blocks)] + result_blocks = [view(A, Data(I)) * view(B, Data(I)) for I in blockdiagindices(A)] return FusedGradedMatrix(copy(A.sectors), result_blocks) end @@ -214,15 +223,12 @@ function FusedGradedMatrix(a::AbelianGradedMatrix{T}) where {T} ) fused_sectors = collect(row_sectors) - # Default: zero blocks with the right dimensions from row/col axes - row_bls = blocklengths(row_ax) - col_bls = blocklengths(col_ax) - diag_blocks = [zeros(T, row_bls[i], col_bls[i]) for i in 1:n] + fused_axes = blockedrange.((datalengths(row_ax), datalengths(col_ax))) + m = FI.zero!(FusedGradedMatrix{T}(undef, fused_sectors, fused_axes)) for bI in eachblockstoredindex(a) - i = Int(Tuple(bI)[1]) - diag_blocks[i] = Array(a[bI]) + m[Data(bI)] = view(a, Data(bI)) end - return FusedGradedMatrix(fused_sectors, diag_blocks) + return m end """ @@ -232,13 +238,10 @@ Convert a `FusedGradedMatrix` to a 2D `AbelianGradedArray`. Inverse of `FusedGradedMatrix(::AbelianGradedArray)`. """ function AbelianGradedArray(m::FusedGradedMatrix{T}) where {T} - codomain, domain = axes(m) - a = similar(m, (codomain, domain)) - for (i, block) in enumerate(m.blocks) - iszero(block) && continue - row_sector = sectors(codomain)[i] - col_sector = sectors(domain)[i] - a[Block(i, i)] = AbelianSectorArray((row_sector, col_sector), block) + a = FI.zero!(similar(m, axes(m))) + for I in blockdiagindices(m) + iszero(view(m, Data(I))) && continue + a[Data(I)] = view(m, Data(I)) end return a end From 29a826c324b3aa883441bc2a55df015b746e7b87 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Sat, 11 Apr 2026 13:48:01 -0400 Subject: [PATCH 06/12] Simplify FusedGradedMatrix/AbelianGradedArray conversions Co-Authored-By: Claude Opus 4.6 (1M context) --- src/fusedgradedmatrix.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/fusedgradedmatrix.jl b/src/fusedgradedmatrix.jl index 6452e915..9989a42e 100644 --- a/src/fusedgradedmatrix.jl +++ b/src/fusedgradedmatrix.jl @@ -224,9 +224,9 @@ function FusedGradedMatrix(a::AbelianGradedMatrix{T}) where {T} fused_sectors = collect(row_sectors) fused_axes = blockedrange.((datalengths(row_ax), datalengths(col_ax))) - m = FI.zero!(FusedGradedMatrix{T}(undef, fused_sectors, fused_axes)) - for bI in eachblockstoredindex(a) - m[Data(bI)] = view(a, Data(bI)) + m = FusedGradedMatrix{T}(undef, fused_sectors, fused_axes) + for I in eachblockstoredindex(a) + m[Data(I)] = view(a, Data(I)) end return m end @@ -238,9 +238,8 @@ Convert a `FusedGradedMatrix` to a 2D `AbelianGradedArray`. Inverse of `FusedGradedMatrix(::AbelianGradedArray)`. """ function AbelianGradedArray(m::FusedGradedMatrix{T}) where {T} - a = FI.zero!(similar(m, axes(m))) + a = similar(m, axes(m)) for I in blockdiagindices(m) - iszero(view(m, Data(I))) && continue a[Data(I)] = view(m, Data(I)) end return a From 6f76f0a9ca0feb9aa067b1299b4e4f927f4c2e7c Mon Sep 17 00:00:00 2001 From: mtfishman Date: Sat, 11 Apr 2026 13:48:43 -0400 Subject: [PATCH 07/12] Use blockdiagindices in FusedGradedMatrix(::AbelianGradedMatrix) conversion Co-Authored-By: Claude Opus 4.6 (1M context) --- src/fusedgradedmatrix.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fusedgradedmatrix.jl b/src/fusedgradedmatrix.jl index 9989a42e..09727cc3 100644 --- a/src/fusedgradedmatrix.jl +++ b/src/fusedgradedmatrix.jl @@ -225,7 +225,7 @@ function FusedGradedMatrix(a::AbelianGradedMatrix{T}) where {T} fused_sectors = collect(row_sectors) fused_axes = blockedrange.((datalengths(row_ax), datalengths(col_ax))) m = FusedGradedMatrix{T}(undef, fused_sectors, fused_axes) - for I in eachblockstoredindex(a) + for I in blockdiagindices(m) m[Data(I)] = view(a, Data(I)) end return m From b35ffcad672a86c3f20e64b03dfc28e8f0c9a09e Mon Sep 17 00:00:00 2001 From: mtfishman Date: Sat, 11 Apr 2026 13:50:25 -0400 Subject: [PATCH 08/12] Simplify FusedGradedMatrix(::AbelianGradedMatrix) validation Co-Authored-By: Claude Opus 4.6 (1M context) --- src/fusedgradedmatrix.jl | 18 ++---------------- 1 file changed, 2 insertions(+), 16 deletions(-) diff --git a/src/fusedgradedmatrix.jl b/src/fusedgradedmatrix.jl index 09727cc3..70a1d0aa 100644 --- a/src/fusedgradedmatrix.jl +++ b/src/fusedgradedmatrix.jl @@ -203,26 +203,12 @@ Convert a 2D block-diagonal `AbelianGradedArray` (as produced by `matricize`) in """ function FusedGradedMatrix(a::AbelianGradedMatrix{T}) where {T} row_ax, col_ax = axes(a) - n = blocklength(row_ax) - blocklength(col_ax) == n || - throw( - ArgumentError("AbelianGradedMatrix must have matching row/column block counts") - ) - - row_sectors = sectors(row_ax) - col_sectors = sectors(col_ax) - row_sectors == dual.(col_sectors) || throw( + sectors(row_ax) == dual.(sectors(col_ax)) || throw( ArgumentError( "AbelianGradedMatrix axes must be canonical duals to convert to FusedGradedMatrix" ) ) - BlockSparseArrays.isblockdiagonal(a) || throw( - ArgumentError( - "AbelianGradedMatrix must be block-diagonal to convert to FusedGradedMatrix" - ) - ) - - fused_sectors = collect(row_sectors) + fused_sectors = collect(sectors(row_ax)) fused_axes = blockedrange.((datalengths(row_ax), datalengths(col_ax))) m = FusedGradedMatrix{T}(undef, fused_sectors, fused_axes) for I in blockdiagindices(m) From 7b0a93c8e85e34aaea8901f2ef61768c4d14c322 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Sat, 11 Apr 2026 13:58:38 -0400 Subject: [PATCH 09/12] Refactor FusedGradedMatrix * to delegate to mul! via allocate_output Co-Authored-By: Claude Opus 4.6 (1M context) --- src/fusedgradedmatrix.jl | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/src/fusedgradedmatrix.jl b/src/fusedgradedmatrix.jl index 70a1d0aa..f85bf0ca 100644 --- a/src/fusedgradedmatrix.jl +++ b/src/fusedgradedmatrix.jl @@ -159,13 +159,23 @@ function LinearAlgebra.mul!( return C end -function Base.:(*)(A::FusedGradedMatrix{T₁}, B::FusedGradedMatrix{T₂}) where {T₁, T₂} - axes(A, 2) == dual(axes(B, 1)) || throw(DimensionMismatch("sectors must match")) - T = Base.promote_op(LinearAlgebra.matprod, T₁, T₂) - result_blocks = [view(A, Data(I)) * view(B, Data(I)) for I in blockdiagindices(A)] +function allocate_output(::typeof(*), A::FusedGradedMatrix, B::FusedGradedMatrix) + cod_axes = data.(eachblockaxis(axes(A, 1))) + dom_axes = data.(eachblockaxis(axes(B, 2))) + result_blocks = [ + similar( + Base.promote_op(*, typeof(view(A, Data(I))), typeof(view(B, Data(I)))), + (cod_axes[Int(Tuple(I)[1])], dom_axes[Int(Tuple(I)[2])]) + ) for I in blockdiagindices(A) + ] return FusedGradedMatrix(copy(A.sectors), result_blocks) end +function Base.:(*)(A::FusedGradedMatrix, B::FusedGradedMatrix) + C = allocate_output(*, A, B) + return mul!(C, A, B, true, false) +end + # ======================== similar ======================== function Base.similar(m::FusedGradedMatrix, ::Type{T}) where {T} From cb3d1a554ed201b1693601ded26e219f61a2bfab Mon Sep 17 00:00:00 2001 From: mtfishman Date: Sat, 11 Apr 2026 13:59:37 -0400 Subject: [PATCH 10/12] Use mul!(C, A, B) in *, fix flaky broadcast test with undef data Co-Authored-By: Claude Opus 4.6 (1M context) --- src/fusedgradedmatrix.jl | 2 +- test/test_tensoralgebra.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/fusedgradedmatrix.jl b/src/fusedgradedmatrix.jl index f85bf0ca..ef246e58 100644 --- a/src/fusedgradedmatrix.jl +++ b/src/fusedgradedmatrix.jl @@ -173,7 +173,7 @@ end function Base.:(*)(A::FusedGradedMatrix, B::FusedGradedMatrix) C = allocate_output(*, A, B) - return mul!(C, A, B, true, false) + return mul!(C, A, B) end # ======================== similar ======================== diff --git a/test/test_tensoralgebra.jl b/test/test_tensoralgebra.jl index 45d6324f..e0b9fb56 100644 --- a/test/test_tensoralgebra.jl +++ b/test/test_tensoralgebra.jl @@ -103,8 +103,8 @@ end @testset "AbelianGradedArray linear broadcasting" begin g1 = gradedrange([U1(0) => 2, U1(1) => 3]) g2 = gradedrange([U1(0) => 1, U1(-1) => 2]) - a = AbelianGradedArray{Float64}(undef, g1, g2) - b = AbelianGradedArray{Float64}(undef, g1, g2) + a = zeros(Float64, g1, g2) + b = zeros(Float64, g1, g2) # Use allowed block (2,2): U1(1) × U1(-1) = 0 block_a = randn!(Matrix{Float64}(undef, 3, 2)) From 64e164e91eca7b682b4391c3522321e980dedba7 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Sat, 11 Apr 2026 14:05:22 -0400 Subject: [PATCH 11/12] Add eachdataaxis and eachsectoraxis accessors Co-Authored-By: Claude Opus 4.6 (1M context) --- src/GradedArrays.jl | 1 + src/abeliangradedarray.jl | 2 +- src/fusedgradedmatrix.jl | 4 ++-- src/gradedoneto.jl | 2 ++ src/sectoroneto.jl | 2 ++ test/test_exports.jl | 2 ++ 6 files changed, 10 insertions(+), 3 deletions(-) diff --git a/src/GradedArrays.jl b/src/GradedArrays.jl index 2cbd5489..0b879100 100644 --- a/src/GradedArrays.jl +++ b/src/GradedArrays.jl @@ -14,6 +14,7 @@ export FusedGradedMatrix export dual, flip, gradedrange, isdual, data, dataaxes, dataaxes1, datalength, datalengths, + eachdataaxis, eachsectoraxis, sector, sectoraxes, sectoraxes1, sectorlength, sectorlengths, sectors, sectortype, Data diff --git a/src/abeliangradedarray.jl b/src/abeliangradedarray.jl index 4ed696dc..693844bf 100644 --- a/src/abeliangradedarray.jl +++ b/src/abeliangradedarray.jl @@ -27,7 +27,7 @@ end function AbelianGradedArray{T, N, D, S}( ::UndefInitializer, axs::NTuple{N, GradedOneTo{S}} ) where {T, N, D <: AbstractArray{T, N}, S <: SectorRange} - block_axes = map(ax -> data.(eachblockaxis(ax)), axs) + block_axes = map(eachdataaxis, axs) function allocate_block(bk) bk_inds = Int.(Tuple(bk)) return similar(D, ntuple(d -> block_axes[d][bk_inds[d]], Val(N))) diff --git a/src/fusedgradedmatrix.jl b/src/fusedgradedmatrix.jl index ef246e58..a6cc27cf 100644 --- a/src/fusedgradedmatrix.jl +++ b/src/fusedgradedmatrix.jl @@ -160,8 +160,8 @@ function LinearAlgebra.mul!( end function allocate_output(::typeof(*), A::FusedGradedMatrix, B::FusedGradedMatrix) - cod_axes = data.(eachblockaxis(axes(A, 1))) - dom_axes = data.(eachblockaxis(axes(B, 2))) + cod_axes = eachdataaxis(axes(A, 1)) + dom_axes = eachdataaxis(axes(B, 2)) result_blocks = [ similar( Base.promote_op(*, typeof(view(A, Data(I))), typeof(view(B, Data(I)))), diff --git a/src/gradedoneto.jl b/src/gradedoneto.jl index c76c1a88..f70c5fa6 100644 --- a/src/gradedoneto.jl +++ b/src/gradedoneto.jl @@ -77,6 +77,8 @@ function BlockSparseArrays.eachblockaxis(g::GradedOneTo) for (s, m) in zip(sectors(g), datalengths(g)) ] end +eachdataaxis(g::GradedOneTo) = data.(eachblockaxis(g)) +eachsectoraxis(g::GradedOneTo) = sector.(eachblockaxis(g)) function BlockSparseArrays.mortar_axis(axs::AbstractVector{<:SectorOneTo}) isempty(axs) && throw( diff --git a/src/sectoroneto.jl b/src/sectoroneto.jl index 22fcfd7f..70b3ce59 100644 --- a/src/sectoroneto.jl +++ b/src/sectoroneto.jl @@ -69,6 +69,8 @@ to_gradedrange(r::SectorOneTo) = gradedrange([sector(r) => datalength(r)]) # ======================== BlockSparseArrays interface ======================== BlockSparseArrays.eachblockaxis(r::SectorOneTo) = [r] +eachdataaxis(r::SectorOneTo) = [data(r)] +eachsectoraxis(r::SectorOneTo) = [sector(r)] # ======================== tensor_product ======================== diff --git a/test/test_exports.jl b/test/test_exports.jl index 7ae2f09c..09fc4a57 100644 --- a/test/test_exports.jl +++ b/test/test_exports.jl @@ -24,6 +24,8 @@ using Test: @test, @testset :dataaxes1, :datalength, :datalengths, + :eachdataaxis, + :eachsectoraxis, :O2, :SectorIdentity, :SectorMatrix, From e9770975f779a26792ba228654c690b60275cdf2 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Sat, 11 Apr 2026 14:10:01 -0400 Subject: [PATCH 12/12] Simplify FusedGradedMatrix(::AbelianGradedMatrix) with datalengths.(axes(a)) Co-Authored-By: Claude Opus 4.6 (1M context) --- src/fusedgradedmatrix.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/fusedgradedmatrix.jl b/src/fusedgradedmatrix.jl index a6cc27cf..35269b0f 100644 --- a/src/fusedgradedmatrix.jl +++ b/src/fusedgradedmatrix.jl @@ -212,14 +212,13 @@ Convert a 2D block-diagonal `AbelianGradedArray` (as produced by `matricize`) in `FusedGradedMatrix`. Extracts diagonal blocks from the stored entries. """ function FusedGradedMatrix(a::AbelianGradedMatrix{T}) where {T} - row_ax, col_ax = axes(a) - sectors(row_ax) == dual.(sectors(col_ax)) || throw( + sectors(axes(a, 1)) == dual.(sectors(axes(a, 2))) || throw( ArgumentError( "AbelianGradedMatrix axes must be canonical duals to convert to FusedGradedMatrix" ) ) - fused_sectors = collect(sectors(row_ax)) - fused_axes = blockedrange.((datalengths(row_ax), datalengths(col_ax))) + fused_sectors = collect(sectors(axes(a, 1))) + fused_axes = blockedrange.(datalengths.(axes(a))) m = FusedGradedMatrix{T}(undef, fused_sectors, fused_axes) for I in blockdiagindices(m) m[Data(I)] = view(a, Data(I))