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..0b879100 100644 --- a/src/GradedArrays.jl +++ b/src/GradedArrays.jl @@ -11,19 +11,22 @@ export AbstractSectorArray, export AbstractGradedArray, AbstractGradedMatrix export AbelianGradedArray, AbelianGradedVector, AbelianGradedMatrix export FusedGradedMatrix -export gradedrange export dual, flip, gradedrange, isdual, data, dataaxes, dataaxes1, datalength, datalengths, + eachdataaxis, eachsectoraxis, 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 +36,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..693844bf 100644 --- a/src/abeliangradedarray.jl +++ b/src/abeliangradedarray.jl @@ -27,11 +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(eachdataaxis, 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)) => - D(undef, ntuple(d -> 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 @@ -66,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) @@ -108,32 +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 - # --------------------------------------------------------------------------- # 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..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} @@ -55,9 +61,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 +69,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 +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(sector(x), convert(A, data(x))) + return AbelianSectorArray{T₁, N, A, S}(sector(x), convert(A, data(x))) end # ======================== permutedims ======================== 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/data.jl b/src/data.jl new file mode 100644 index 00000000..3795349e --- /dev/null +++ b/src/data.jl @@ -0,0 +1,14 @@ +""" + 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))) +Base.Tuple(I::Data) = I.n diff --git a/src/fusedgradedmatrix.jl b/src/fusedgradedmatrix.jl index a402ff56..35269b0f 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,6 +34,49 @@ 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 FusedGradedMatrix{T, D, S}( + ::UndefInitializer, + sectors::Vector{S}, + axes::Tuple{BlockedOneTo, BlockedOneTo} + ) where {T, D <: AbstractMatrix{T}, S <: SectorRange} + 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}. +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, domain_blocklengths)) + ) +end # ======================== Accessors ======================== @@ -54,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)) @@ -62,16 +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 # ======================== eachblockstoredindex ======================== @@ -83,8 +116,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(I -> view(m, I), eachblockstoredindex(m)) + return Diagonal(collect(diagblocks)) end # ======================== fill! / zero! ======================== @@ -105,25 +138,44 @@ 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")) - T = Base.promote_op(LinearAlgebra.matprod, T₁, T₂) - result_blocks = [A.blocks[i] * B.blocks[i] for i in eachindex(A.blocks)] +function allocate_output(::typeof(*), A::FusedGradedMatrix, B::FusedGradedMatrix) + 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)))), + (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) +end + # ======================== similar ======================== function Base.similar(m::FusedGradedMatrix, ::Type{T}) where {T} @@ -160,36 +212,18 @@ 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) - 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(axes(a, 1)) == dual.(sectors(axes(a, 2))) || 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) - # 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] - for bI in eachblockstoredindex(a) - i = Int(Tuple(bI)[1]) - diag_blocks[i] = Array(a[bI]) + 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)) end - return FusedGradedMatrix(fused_sectors, diag_blocks) + return m end """ @@ -199,13 +233,9 @@ 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 = similar(m, axes(m)) + for I in blockdiagindices(m) + a[Data(I)] = view(m, Data(I)) end return a end 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/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..70b3ce59 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 @@ -68,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_abelianarray.jl b/test/test_abelianarray.jl index 8d00ddfb..e55e1642 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 Exception 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..c07b5a3f --- /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 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 Exception (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..09fc4a57 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, @@ -23,6 +24,8 @@ using Test: @test, @testset :dataaxes1, :datalength, :datalengths, + :eachdataaxis, + :eachsectoraxis, :O2, :SectorIdentity, :SectorMatrix, 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 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))