Require undef in BlockSparseArray constructors

Project.toml

@@ -1,7 +1,7 @@
 name = "BlockSparseArrays"
 uuid = "2c9a651f-6452-4ace-a6ac-809f4280fbb4"
 authors = ["ITensor developers <[email protected]> and contributors"]
-version = "0.2.28"
+version = "0.3.0"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

README.md

@@ -34,129 +34,17 @@ julia> Pkg.add("BlockSparseArrays")
 ## Examples
 
 ````julia
-using BlockArrays: BlockArrays, BlockedVector, Block, blockedrange
+using BlockArrays: Block
 using BlockSparseArrays: BlockSparseArray, blockstoredlength
-using Test: @test, @test_broken
-
-function main()
-  # Block dimensions
-  i1 = [2, 3]
-  i2 = [2, 3]
-
-  i_axes = (blockedrange(i1), blockedrange(i2))
-
-  function block_size(axes, block)
-    return length.(getindex.(axes, Block.(block.n)))
-  end
-
-  # Data
-  nz_blocks = Block.([(1, 1), (2, 2)])
-  nz_block_sizes = [block_size(i_axes, nz_block) for nz_block in nz_blocks]
-  nz_block_lengths = prod.(nz_block_sizes)
-
-  # Blocks with contiguous underlying data
-  d_data = BlockedVector(randn(sum(nz_block_lengths)), nz_block_lengths)
-  d_blocks = [
-    reshape(@view(d_data[Block(i)]), block_size(i_axes, nz_blocks[i])) for
-    i in 1:length(nz_blocks)
-  ]
-  b = BlockSparseArray(nz_blocks, d_blocks, i_axes)
-
-  @test blockstoredlength(b) == 2
-
-  # Blocks with discontiguous underlying data
-  d_blocks = randn.(nz_block_sizes)
-  b = BlockSparseArray(nz_blocks, d_blocks, i_axes)
-
-  @test blockstoredlength(b) == 2
-
-  # Access a block
-  @test b[Block(1, 1)] == d_blocks[1]
-
-  # Access a zero block, returns a zero matrix
-  @test b[Block(1, 2)] == zeros(2, 3)
-
-  # Set a zero block
-  a₁₂ = randn(2, 3)
-  b[Block(1, 2)] = a₁₂
-  @test b[Block(1, 2)] == a₁₂
-
-  # Matrix multiplication
-  @test b * b ≈ Array(b) * Array(b)
-
-  permuted_b = permutedims(b, (2, 1))
-  @test permuted_b isa BlockSparseArray
-  @test permuted_b == permutedims(Array(b), (2, 1))
-
-  @test b + b ≈ Array(b) + Array(b)
-  @test b + b isa BlockSparseArray
-  # TODO: Fix this, broken.
-  @test_broken blockstoredlength(b + b) == 2
-
-  scaled_b = 2b
-  @test scaled_b ≈ 2Array(b)
-  @test scaled_b isa BlockSparseArray
-
-  # TODO: Fix this, broken.
-  @test_broken reshape(b, ([4, 6, 6, 9],)) isa BlockSparseArray{<:Any,1}
-
-  return nothing
-end
-
-main()
-````
-
-# BlockSparseArrays.jl and BlockArrays.jl interface
-
-````julia
-using BlockArrays: BlockArrays, Block
-using BlockSparseArrays: BlockSparseArray
-
-i1 = [2, 3]
-i2 = [2, 3]
-B = BlockSparseArray{Float64}(i1, i2)
-B[Block(1, 1)] = randn(2, 2)
-B[Block(2, 2)] = randn(3, 3)
-
-# Minimal interface
-
-# Specifies the block structure
-@show collect.(BlockArrays.blockaxes(axes(B, 1)))
-
-# Index range of a block
-@show axes(B, 1)[Block(1)]
-
-# Last index of each block
-@show BlockArrays.blocklasts(axes(B, 1))
-
-# Find the block containing the index
-@show BlockArrays.findblock(axes(B, 1), 3)
-
-# Retrieve a block
-@show B[Block(1, 1)]
-@show BlockArrays.viewblock(B, Block(1, 1))
-
-# Check block bounds
-@show BlockArrays.blockcheckbounds(B, 2, 2)
-@show BlockArrays.blockcheckbounds(B, Block(2, 2))
-
-# Derived interface
-
-# Specifies the block structure
-@show collect(Iterators.product(BlockArrays.blockaxes(B)...))
-
-# Iterate over block views
-@show sum.(BlockArrays.eachblock(B))
-
-# Reshape into 1-d
-# TODO: Fix this, broken.
-# @show BlockArrays.blockvec(B)[Block(1)]
-
-# Array-of-array view
-@show BlockArrays.blocks(B)[1, 1] == B[Block(1, 1)]
-
-# Access an index within a block
-@show B[Block(1, 1)[1, 1]] == B[1, 1]
+using Test: @test
+
+a = BlockSparseArray{Float64}(undef, [2, 3], [2, 3])
+a[Block(1, 2)] = randn(2, 3)
+a[Block(2, 1)] = randn(3, 2)
+@test blockstoredlength(a) == 2
+b = a .+ 2 .* a'
+@test Array(b) ≈ Array(a) + 2 * Array(a')
+@test blockstoredlength(b) == 2
 ````
 
 ---

docs/Project.toml

@@ -6,6 +6,6 @@ Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 
 [compat]
 BlockArrays = "1"
-BlockSparseArrays = "0.2"
+BlockSparseArrays = "0.3"
 Documenter = "1"
 Literate = "2"

examples/Project.toml

@@ -5,5 +5,5 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 BlockArrays = "1"
-BlockSparseArrays = "0.2"
+BlockSparseArrays = "0.3"
 Test = "1"

examples/README.jl

@@ -39,124 +39,14 @@ julia> Pkg.add("BlockSparseArrays")
 
 # ## Examples
 
-using BlockArrays: BlockArrays, BlockedVector, Block, blockedrange
+using BlockArrays: Block
 using BlockSparseArrays: BlockSparseArray, blockstoredlength
-using Test: @test, @test_broken
-
-function main()
-  ## Block dimensions
-  i1 = [2, 3]
-  i2 = [2, 3]
-
-  i_axes = (blockedrange(i1), blockedrange(i2))
-
-  function block_size(axes, block)
-    return length.(getindex.(axes, Block.(block.n)))
-  end
-
-  ## Data
-  nz_blocks = Block.([(1, 1), (2, 2)])
-  nz_block_sizes = [block_size(i_axes, nz_block) for nz_block in nz_blocks]
-  nz_block_lengths = prod.(nz_block_sizes)
-
-  ## Blocks with contiguous underlying data
-  d_data = BlockedVector(randn(sum(nz_block_lengths)), nz_block_lengths)
-  d_blocks = [
-    reshape(@view(d_data[Block(i)]), block_size(i_axes, nz_blocks[i])) for
-    i in 1:length(nz_blocks)
-  ]
-  b = BlockSparseArray(nz_blocks, d_blocks, i_axes)
-
-  @test blockstoredlength(b) == 2
-
-  ## Blocks with discontiguous underlying data
-  d_blocks = randn.(nz_block_sizes)
-  b = BlockSparseArray(nz_blocks, d_blocks, i_axes)
-
-  @test blockstoredlength(b) == 2
-
-  ## Access a block
-  @test b[Block(1, 1)] == d_blocks[1]
-
-  ## Access a zero block, returns a zero matrix
-  @test b[Block(1, 2)] == zeros(2, 3)
-
-  ## Set a zero block
-  a₁₂ = randn(2, 3)
-  b[Block(1, 2)] = a₁₂
-  @test b[Block(1, 2)] == a₁₂
-
-  ## Matrix multiplication
-  @test b * b ≈ Array(b) * Array(b)
-
-  permuted_b = permutedims(b, (2, 1))
-  @test permuted_b isa BlockSparseArray
-  @test permuted_b == permutedims(Array(b), (2, 1))
-
-  @test b + b ≈ Array(b) + Array(b)
-  @test b + b isa BlockSparseArray
-  ## TODO: Fix this, broken.
-  @test_broken blockstoredlength(b + b) == 2
-
-  scaled_b = 2b
-  @test scaled_b ≈ 2Array(b)
-  @test scaled_b isa BlockSparseArray
-
-  ## TODO: Fix this, broken.
-  @test_broken reshape(b, ([4, 6, 6, 9],)) isa BlockSparseArray{<:Any,1}
-
-  return nothing
-end
-
-main()
-
-# # BlockSparseArrays.jl and BlockArrays.jl interface
-
-using BlockArrays: BlockArrays, Block
-using BlockSparseArrays: BlockSparseArray
-
-i1 = [2, 3]
-i2 = [2, 3]
-B = BlockSparseArray{Float64}(i1, i2)
-B[Block(1, 1)] = randn(2, 2)
-B[Block(2, 2)] = randn(3, 3)
-
-## Minimal interface
-
-## Specifies the block structure
-@show collect.(BlockArrays.blockaxes(axes(B, 1)))
-
-## Index range of a block
-@show axes(B, 1)[Block(1)]
-
-## Last index of each block
-@show BlockArrays.blocklasts(axes(B, 1))
-
-## Find the block containing the index
-@show BlockArrays.findblock(axes(B, 1), 3)
-
-## Retrieve a block
-@show B[Block(1, 1)]
-@show BlockArrays.viewblock(B, Block(1, 1))
-
-## Check block bounds
-@show BlockArrays.blockcheckbounds(B, 2, 2)
-@show BlockArrays.blockcheckbounds(B, Block(2, 2))
-
-## Derived interface
-
-## Specifies the block structure
-@show collect(Iterators.product(BlockArrays.blockaxes(B)...))
-
-## Iterate over block views
-@show sum.(BlockArrays.eachblock(B))
-
-## Reshape into 1-d
-## TODO: Fix this, broken.
-## @show BlockArrays.blockvec(B)[Block(1)]
-
-## Array-of-array view
-@show BlockArrays.blocks(B)[1, 1] == B[Block(1, 1)]
-
-## Access an index within a block
-@show B[Block(1, 1)[1, 1]] == B[1, 1]
+using Test: @test
+
+a = BlockSparseArray{Float64}(undef, [2, 3], [2, 3])
+a[Block(1, 2)] = randn(2, 3)
+a[Block(2, 1)] = randn(3, 2)
+@test blockstoredlength(a) == 2
+b = a .+ 2 .* a'
+@test Array(b) ≈ Array(a) + 2 * Array(a')
+@test blockstoredlength(b) == 2

ext/BlockSparseArraysGradedUnitRangesExt/BlockSparseArraysGradedUnitRangesExt.jl

@@ -17,7 +17,7 @@ function similar_blocksparse(
   return BlockSparseArray{
     elt,length(axes),similartype(unwrap_array_type(blocktype(a)), elt, axes)
   }(
-    axes
+    undef, axes
   )
 end
 

src/BlockSparseArrays.jl

@@ -5,7 +5,8 @@ export BlockSparseArray,
   BlockSparseVector,
   blockstoredlength,
   eachblockstoredindex,
-  eachstoredblock
+  eachstoredblock,
+  sparsemortar
 
 # factorizations
 include("factorizations/svd.jl")
@@ -38,7 +39,6 @@ include("abstractblocksparsearray/cat.jl")
 include("abstractblocksparsearray/adapt.jl")
 
 # functions specifically for BlockSparseArray
-include("blocksparsearray/defaults.jl")
 include("blocksparsearray/blocksparsearray.jl")
 include("blocksparsearray/blockdiagonalarray.jl")
 

src/abstractblocksparsearray/arraylayouts.jl

@@ -38,7 +38,7 @@ function ArrayLayouts.sub_materialize(layout::BlockLayout{<:SparseLayout}, a, ax
   # TODO: Define `blocktype`/`blockstype` for `SubArray` wrapping `BlockSparseArray`.
   # TODO: Use `similar`?
   blocktype_a = blocktype(parent(a))
-  a_dest = BlockSparseArray{eltype(a),length(axes),blocktype_a}(axes)
+  a_dest = BlockSparseArray{eltype(a),length(axes),blocktype_a}(undef, axes)
   a_dest .= a
   return a_dest
 end

src/blocksparsearray/blockdiagonalarray.jl

@@ -1,7 +1,8 @@
-# type alias for block-diagonal
-using LinearAlgebra: Diagonal
+using BlockArrays: blockedrange
 using DiagonalArrays: DiagonalArrays, diagonal
+using LinearAlgebra: Diagonal
 
+# type alias for block-diagonal
 const BlockDiagonal{T,A,Axes,V<:AbstractVector{A}} = BlockSparseMatrix{
   T,A,Diagonal{A,V},Axes
 }
@@ -12,8 +13,8 @@ const BlockSparseDiagonal{T,A<:AbstractBlockSparseVector{T}} = Diagonal{T,A}
 end
 
 function BlockDiagonal(blocks::AbstractVector{<:AbstractMatrix})
-  return BlockSparseArray(
-    Diagonal(blocks), (blockedrange(size.(blocks, 1)), blockedrange(size.(blocks, 2)))
+  return _BlockSparseArray(
+    Diagonal(blocks), blockedrange.((size.(blocks, 1), size.(blocks, 2)))
   )
 end