77[ ![ Code Style: Blue] ( https://img.shields.io/badge/code%20style-blue-4495d1.svg )] ( https://github.com/invenia/BlueStyle )
88[ ![ Aqua] ( https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg )] ( https://github.com/JuliaTesting/Aqua.jl )
99
10+ A block sparse array type in Julia based on the [ ` BlockArrays.jl ` ] ( https://github.com/JuliaArrays/BlockArrays.jl ) interface.
11+
1012## Installation instructions
1113
1214This package resides in the ` ITensor/ITensorRegistry ` local registry.
@@ -32,10 +34,131 @@ julia> Pkg.add("BlockSparseArrays")
3234## Examples
3335
3436```` julia
35- using BlockSparseArrays: BlockSparseArrays
37+ using BlockArrays: BlockArrays, BlockedVector, Block, blockedrange
38+ using BlockSparseArrays: BlockSparseArray, block_stored_length
39+ using Test: @test , @test_broken
40+
41+ function main ()
42+ # Block dimensions
43+ i1 = [2 , 3 ]
44+ i2 = [2 , 3 ]
45+
46+ i_axes = (blockedrange (i1), blockedrange (i2))
47+
48+ function block_size (axes, block)
49+ return length .(getindex .(axes, Block .(block. n)))
50+ end
51+
52+ # Data
53+ nz_blocks = Block .([(1 , 1 ), (2 , 2 )])
54+ nz_block_sizes = [block_size (i_axes, nz_block) for nz_block in nz_blocks]
55+ nz_block_lengths = prod .(nz_block_sizes)
56+
57+ # Blocks with contiguous underlying data
58+ d_data = BlockedVector (randn (sum (nz_block_lengths)), nz_block_lengths)
59+ d_blocks = [
60+ reshape (@view (d_data[Block (i)]), block_size (i_axes, nz_blocks[i])) for
61+ i in 1 : length (nz_blocks)
62+ ]
63+ b = BlockSparseArray (nz_blocks, d_blocks, i_axes)
64+
65+ @test block_stored_length (b) == 2
66+
67+ # Blocks with discontiguous underlying data
68+ d_blocks = randn .(nz_block_sizes)
69+ b = BlockSparseArray (nz_blocks, d_blocks, i_axes)
70+
71+ @test block_stored_length (b) == 2
72+
73+ # Access a block
74+ @test b[Block (1 , 1 )] == d_blocks[1 ]
75+
76+ # Access a zero block, returns a zero matrix
77+ @test b[Block (1 , 2 )] == zeros (2 , 3 )
78+
79+ # Set a zero block
80+ a₁₂ = randn (2 , 3 )
81+ b[Block (1 , 2 )] = a₁₂
82+ @test b[Block (1 , 2 )] == a₁₂
83+
84+ # Matrix multiplication
85+ # TODO : Fix this, broken.
86+ @test_broken b * b ≈ Array (b) * Array (b)
87+
88+ permuted_b = permutedims (b, (2 , 1 ))
89+ @test permuted_b isa BlockSparseArray
90+ @test permuted_b == permutedims (Array (b), (2 , 1 ))
91+
92+ @test b + b ≈ Array (b) + Array (b)
93+ @test b + b isa BlockSparseArray
94+ # TODO : Fix this, broken.
95+ @test_broken block_stored_length (b + b) == 2
96+
97+ scaled_b = 2 b
98+ @test scaled_b ≈ 2 Array (b)
99+ @test scaled_b isa BlockSparseArray
100+
101+ # TODO : Fix this, broken.
102+ @test_broken reshape (b, ([4 , 6 , 6 , 9 ],)) isa BlockSparseArray{<: Any ,1 }
103+
104+ return nothing
105+ end
106+
107+ main ()
36108````
37109
38- Examples go here.
110+ # BlockSparseArrays.jl and BlockArrays.jl interface
111+
112+ ```` julia
113+ using BlockArrays: BlockArrays, Block
114+ using BlockSparseArrays: BlockSparseArray
115+
116+ i1 = [2 , 3 ]
117+ i2 = [2 , 3 ]
118+ B = BlockSparseArray {Float64} (i1, i2)
119+ B[Block (1 , 1 )] = randn (2 , 2 )
120+ B[Block (2 , 2 )] = randn (3 , 3 )
121+
122+ # Minimal interface
123+
124+ # Specifies the block structure
125+ @show collect .(BlockArrays. blockaxes (axes (B, 1 )))
126+
127+ # Index range of a block
128+ @show axes (B, 1 )[Block (1 )]
129+
130+ # Last index of each block
131+ @show BlockArrays. blocklasts (axes (B, 1 ))
132+
133+ # Find the block containing the index
134+ @show BlockArrays. findblock (axes (B, 1 ), 3 )
135+
136+ # Retrieve a block
137+ @show B[Block (1 , 1 )]
138+ @show BlockArrays. viewblock (B, Block (1 , 1 ))
139+
140+ # Check block bounds
141+ @show BlockArrays. blockcheckbounds (B, 2 , 2 )
142+ @show BlockArrays. blockcheckbounds (B, Block (2 , 2 ))
143+
144+ # Derived interface
145+
146+ # Specifies the block structure
147+ @show collect (Iterators. product (BlockArrays. blockaxes (B)... ))
148+
149+ # Iterate over block views
150+ @show sum .(BlockArrays. eachblock (B))
151+
152+ # Reshape into 1-d
153+ # TODO : Fix this, broken.
154+ # @show BlockArrays.blockvec(B)[Block(1)]
155+
156+ # Array-of-array view
157+ @show BlockArrays. blocks (B)[1 , 1 ] == B[Block (1 , 1 )]
158+
159+ # Access an index within a block
160+ @show B[Block (1 , 1 )[1 , 1 ]] == B[1 , 1 ]
161+ ````
39162
40163---
41164
0 commit comments