[BlockSparseArrays] BlockSparseArray functionality #1336
Comments
ogauthe
added
enhancement
|
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 |
|
In terms of I have a protype of a QR decomposition of a |
mtfishman
changed the title
|
Also note that slicing like this should work right now: a[Block(1, 1)[1:2, 1:2]]i.e. you can take slices within a specified block. See BlockArrays.jl for a reference on that slicing notation. |
|
new feature request: I updated the first comment. |
|
new issue: |
|
new issue:
|
|
@ogauthe a number of these issues were fixed by #1332, I've updated the list in the first post accordingly. I added regression tests in #1360 for ones that still need to be fixed, and additionally added placeholder tests that I've marked as broken in the BlockSparseArrays tests. Please continue to update this post with new issues you find, and/or make PRs with broken behavior marked with |
|
Feature request: |
|
I think ideally Alternatively, Good question about whether or not the axes should get dualed if |
|
The solution to accept any Axes<:Tuple{Vararg{<:AbstractUnitRange,N}},Then one can construct a g1 = gradedrange([U1(0) => 1])
m1 = BlockSparseArray{Float64}(dual(g1), g1,)outputs |
|
Thanks for investigating. That seems like the right move to generalize the axes in that way. Hopefully that error is easy enough to circumvent. |
|
I continue in exploring the effect of |
mtfishman
changed the title
|
issue: I cannot write a slice of a block a[BlockArrays.Block(1,1)][1:2,1:2] = ones((2,2))does not write |
mtfishman
changed the title
|
issue: a[BlockArrays.Block(1,1)] = ones((2,2))
println(LinearAlgebra.norm(a)) # 2.0
a[BlockArrays.Block(1,1)][1, 1] = NaN
println(LinearAlgebra.norm(a[BlockArrays.Block(1,1)])) # NaN
println(LinearAlgebra.norm(a)) # AssertionErroroutputs I just checked that replacing |
|
issue: a block can be written with an invalid shape. An error should be raised. a = BlockSparseArray{Float64}([2, 3], [2, 3])
println(size(a)) # (5,5)
b = BlockArrays.Block(1,1)
println(size(a[b])) # (2,2)
a[b] = ones((3,3))
println(size(a)) # (5,5)
println(size(a[b])) # (3,3) |
This issue lists functionalities and feature requests for
BlockSparseArray.Feature requests/issues
BlockSparseArrayto anyAbstractUnitRange(e.g. fixr = gradedrange([U1(0) => 1]); a = BlockSparseArray{Float64}(dual(r), r,)).BlockSparseMatrix, i.e.a'(in progress in1).Base.adjoint(::BlockSparseMatrix).display(::Adjoint{T, BlockSparseArray)anddisplay(::Transpose{T, BlockSparseArray)are not implementeda[Block(2), Block(2)] = randn(3, 3)for initializing a block (currently onlya[Block(2, 2)] = randn(3, 3)works).a[Block(2, 2)] .= 1.0as well as from a slicea[Block(1,1)][1:1,1:1] = ones((1,1))LinearAlgebra.norm(a)crashes whenacontainsNaNsvd,qr, etc. 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 calledblock_svd,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 2 for a prototype of a block-wise QR.block_stored_indicestostored_blocks, for outputting a listVector{<:Block}representing which blocks are nonzero (or zero but stored). We can keepblock_stored_indicesaround but have it outputVector{<:CartesianIndex}, which is the same list of blocks but in the formatCartesianIndexinstead ofBlock. That can be helpful when interpreting theBlockSparseArrayas aSparseArrayDOKof blocks, i.e. the data structureblocks(a::BlockSparseArray). However, that can easily be accessed withstored_indices(blocks(a::BlockSparseArray))so I don't know if it is needed (it follows the general pattern inBlockArraysof definingblock*(a)functions that act as*(blocks(a)), i.e.blocksize(a) == size(blocks(a)), etc.).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 in1).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 in1).nstored(blocks(a))to get the number of initialized/stored blocks. Ideally there would be shorthands for this likeblock_nstored(a)(in progress in1)..*=and./=, such asa .*= 2, are broken (in progress in1).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.Footnotes
[BlockSparseArrays] More general broadcasting and slicing #1332 ↩ ↩2 ↩3 ↩4 ↩5
https://github.com/ITensor/ITensors.jl/blob/v0.3.57/NDTensors/src/lib/BlockSparseArrays/src/backup/LinearAlgebraExt/qr.jl ↩
The text was updated successfully, but these errors were encountered: