mtfishman
commented
6 months ago Thanks. This currently works:
julia> using NDTensors.BlockSparseArrays: Block, BlockSparseArray, blocks
julia> using LinearAlgebra: I
julia> a = BlockSparseArray{Float64}([2, 2], [2, 2])
2×2-blocked 4×4 BlockSparseArray{Float64, 2, Matrix{Float64}, NDTensors.SparseArrayDOKs.SparseArrayDOK{Matrix{Float64}, 2, NDTensors.BlockSparseArrays.BlockZero{Tuple{BlockArrays.BlockedUnitRange{Vector{Int64}}, BlockArrays.BlockedUnitRange{Vector{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{Vector{Int64}}, BlockArrays.BlockedUnitRange{Vector{Int64}}}}:
0.0 0.0 │ 0.0 0.0
0.0 0.0 │ 0.0 0.0
──────────┼──────────
0.0 0.0 │ 0.0 0.0
0.0 0.0 │ 0.0 0.0
julia> a[Block(2, 2)] = I(3)
3×3 Diagonal{Bool, Vector{Bool}}:
1 ⋅ ⋅
⋅ 1 ⋅
⋅ ⋅ 1
julia> a
2×2-blocked 4×4 BlockSparseArray{Float64, 2, Matrix{Float64}, NDTensors.SparseArrayDOKs.SparseArrayDOK{Matrix{Float64}, 2, NDTensors.BlockSparseArrays.BlockZero{Tuple{BlockArrays.BlockedUnitRange{Vector{Int64}}, BlockArrays.BlockedUnitRange{Vector{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{Vector{Int64}}, BlockArrays.BlockedUnitRange{Vector{Int64}}}}:
0.0 0.0 │ 0.0 0.0
0.0 0.0 │ 0.0 0.0
──────────┼──────────
0.0 0.0 │ 1.0 0.0
0.0 0.0 │ 0.0 1.0
julia> using NDTensors.SparseArrayInterface: stored_indices
julia> stored_indices(blocks(a))
1-element Dictionaries.MappedDictionary{CartesianIndex{2}, CartesianIndex{2}, NDTensors.SparseArrayInterface.var"#1#2"{NDTensors.SparseArrayDOKs.SparseArrayDOK{Matrix{Float64}, 2, NDTensors.BlockSparseArrays.BlockZero{Tuple{BlockArrays.BlockedUnitRange{Vector{Int64}}, BlockArrays.BlockedUnitRange{Vector{Int64}}}}}}, Tuple{Dictionaries.Indices{CartesianIndex{2}}}}
CartesianIndex(2, 2) │ CartesianIndex(2, 2)though using this alternative syntax is currently broken:
julia> a = BlockSparseArray{Float64}([2, 2], [2, 2])
2×2-blocked 4×4 BlockSparseArray{Float64, 2, Matrix{Float64}, NDTensors.SparseArrayDOKs.SparseArrayDOK{Matrix{Float64}, 2, NDTensors.BlockSparseArrays.BlockZero{Tuple{BlockArrays.BlockedUnitRange{Vector{Int64}}, BlockArrays.BlockedUnitRange{Vector{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{Vector{Int64}}, BlockArrays.BlockedUnitRange{Vector{Int64}}}}:
0.0 0.0 │ 0.0 0.0
0.0 0.0 │ 0.0 0.0
──────────┼──────────
0.0 0.0 │ 0.0 0.0
0.0 0.0 │ 0.0 0.0
julia> a[Block(2), Block(2)] = I(3)
ERROR: DimensionMismatch: tried to assign (3, 3) array to (2, 2) block
Stacktrace:
[1] setindex!(::BlockSparseArray{…}, ::Diagonal{…}, ::Block{…}, ::Block{…})
@ BlockArrays ~/.julia/packages/BlockArrays/L5yjb/src/abstractblockarray.jl:165
[2] top-level scope
@ REPL[30]:1
Some type information was truncated. Use `show(err)` to see complete types.I would have to think about if it makes sense to support a[Block(2, 2)] .= 1.0 or a[Block(2), Block(2)] .= 1.0 if Block(2, 2) isn't initialized yet (probably we should support that).
ogauthe
emstoudenmire
This issue lists functionalities and feature requests for
BlockSparseArray.Feature requests/issues
block_stored_indices(::LinearAlgebra.Adjoint{T, BlockSparseArray})does not transpose its indicesLinearAlgebra.Adjoint{T, NDTensors.BlockSparseArrays.BlockSparseArray}returns aBlockedArraywhen sliced withBlock.Base.catand related functions.a[1:2, 1:2]) to output non-blocked arrays, and define@blocked a[1:2, 1:2]to explicitly preserve blocking. See the discussion in https://github.com/JuliaArrays/BlockArrays.jl/issues/347.LinearAlgebra.norm(a)crashes whenacontainsNaN.svd,qr, etc. See #1515. These are well defined if the block sparse matrix has a block structure (i.e. the sparsity pattern of the sparse array of arraysblocks(a)) corresponding to a generalized permutation matrix. Probably they should be called something likeblock_svd,block_eigen,block_qr, etc. to distinguish that they are meant to be used on block sparse matrices with those structures (and error if they don't have that structure). See ^linalg for a prototype of a blockwise QR. See also BlockDiagonals.jl for an example in Julia of blockwise factorizations, they use a naming schemesvd_blockwise,eigen_blockwise, etc. The slicing operation introduced in #1489 will be useful for performing block-wise truncated factorizations.Strided.@stridedfails when called withview(::BlockSparseArray, ::Block). As a workaround it will work if you useview!/@view!instroduced in #1498.BlockSparseArrayLiketoAnyBlockSparseArray, which is the naming convention used in other Julia packages for a similar concept^anycuarray.block_nstoredtoblock_stored_lengthandnstoredtostored_length.Array, for exampleDiagonalArrays.DiagonalArray,SparseArrayDOKs.SparseArrayDOK,LinearAlgebra.Diagonal, etc.BlockSparseArraycan have blocks that are AbstractArray subtype, however some operations don't preserve those types properly (i.e. implicitly convert toArrayblocks) or don't work.block_stored_indicestostored_blocks, for outputting a listVector{<:Block}representing which blocks are stored.a = BlockSparseArray{Float64}([2, 3], [2, 3]); @view a[Block(1, 1)]returns aSubArraywhere the last type parameter which marks whether or not the slice supports faster linear indexing isfalse, while it should betrueif that is the case for that block ofa(this is addressed by #1513,@view a[Block(1, 1)]no longer outputs aSubArray, but rather either the block data directly or aBlockViewobject if the block doesn't exist yet).TensorAlgebra.contractfails when called withview(::BlockSparseArray, ::Block)orreshape(view(::BlockSparseArray, ::Block), ...). As a workaround it will work if you useview!/@view!instroduced in #1498.a = BlockSparseArray{Float64}([2, 3], [2, 3]); b = @view a[Block.(1:2), Block.(1:2)]; b[Block(1, 1)] = randn(2, 2)doesn't set the blockBlock(1, 1)(it remains uninitialized, i.e. structurally zero). I think the issue is that@view b[Block(1, 1)]makes two layers ofSubArraywrappers instead of flattening down to a single layer, and those two layers are not being dispatched on properly (in general we only catch if something is aBlockSparseArrayor aBlockSparseArraywrapped in a single wrapper layer).BlockArraysv1.1, see CI for #1363. Fixed by #1503.r = gradedrange([U1(0) => 1]); a = BlockSparseArray{Float64}(r, r); size(view(a, Block(1,1))[1:1,1:1])returns a tuple ofLabelledIntegerinstead ofInt(see discussion, keep it that way at least for now).r = gradedrange([U1(0) => 1]); a = BlockSparseArray{Float64}(dual(r), r); @view(a[Block(1, 1)])[1:1, 1:1]and other combinations ofduallead to method ambiguity errors.dualis not preserved when adding/subtractingBlockSparseArrays, i.e.g = gradedrange([U1(0) => 1]); m = BlockSparseArray{Float64}(dual(g), g); isdual(axes(m + m, 1))should betruebut isfalse.r = gradedrange([U1(0) => 1]); a = BlockSparseArray{Float64}(r, r); @view a[Block(1, 1)].a[2:4, 2:4], by usingBlockArrays.BlockSlice.a[Block(2), Block(2)] = randn(3, 3).a[Block(2, 2)] .= 1.@view(a[Block(1, 1)])[1:1, 1:1] = 1.a[Block(1, 1)] = bifsize(a[Block(1, 1)]) != size(b).BlockSparseMatrixinvolving dual axes.BlockSparseMatrix, i.e.a' * aanda * a', with and without dual axes.adjoint(::BlockSparseMatrix). Can be implemented by overloadingaxes(::Adjoint{<:Any,<:AbstractBlockSparseMatrix}).show(::Adjoint{<:Any,<:BlockSparseMatrix})andshow(::Transpose{<:Any,<:BlockSparseMatrix))are broken.eachindex(::BlockSparseArray)involving dual axes.BlockSparseMatrix, i.e.a'(in progress in^adjoint).Base.similar(a::BlockSparseArray, eltype::type)andBase.similar(a::BlockSparseArray, eltype::type, size::NTuple{N,AbstractUnitRange})do not seteltypecopy(::BlockSparseArray)copy the blocks.a[1:2, 1:2]is not implemented yet and needs to be implemented (in progress in^slicing).stored_indices(blocks(a)))to get a list ofBlockcorresponding to initialized/stored blocks. Ideally there would be shorthands for this likeblock_stored_indices(a)(in progress in^slicing).nstored(blocks(a))to get the number of initialized/stored blocks. Ideally there would be shorthands for this likeblock_nstored(a)(in progress in^slicing)..*=and./=, such asa .*= 2, are broken (in progress in[^1]).Base.:*(::BlockSparseArray, x::Number)andBase.:/(::BlockSparseArray, x::Number)are not definedBase.:*(::ComplexF64, ::BlockSparseArray{Float64})does not change data type for empty array and crashes ifacontains data.[^blockarrays_v1]: https://github.com/ITensor/ITensors.jl/pull/1452, https://github.com/JuliaArrays/BlockArrays.jl/pull/255