From ec8c96ef54868f31bfa613937ddf7541d1d13ced Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Mon, 3 Mar 2025 12:26:15 -0500
Subject: [PATCH 01/14] Require undef in BlockSparseArray constructors

---
 Project.toml                                 |   2 +-
 src/abstractblocksparsearray/arraylayouts.jl |   2 +-
 src/blocksparsearray/blocksparsearray.jl     | 102 ++++++++--------
 test/test_basics.jl                          | 116 +++++++++----------
 4 files changed, 104 insertions(+), 118 deletions(-)

diff --git a/Project.toml b/Project.toml
index f36743c..0a2cd55 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "BlockSparseArrays"
 uuid = "2c9a651f-6452-4ace-a6ac-809f4280fbb4"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
-version = "0.2.27"
+version = "0.3.0"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
diff --git a/src/abstractblocksparsearray/arraylayouts.jl b/src/abstractblocksparsearray/arraylayouts.jl
index e02eb97..8b77120 100644
--- a/src/abstractblocksparsearray/arraylayouts.jl
+++ b/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
diff --git a/src/blocksparsearray/blocksparsearray.jl b/src/blocksparsearray/blocksparsearray.jl
index 44b8d77..f87d943 100644
--- a/src/blocksparsearray/blocksparsearray.jl
+++ b/src/blocksparsearray/blocksparsearray.jl
@@ -3,9 +3,6 @@ using DerivableInterfaces: @interface
 using Dictionaries: Dictionary
 using SparseArraysBase: SparseArrayDOK
 
-# TODO: Delete this.
-## using BlockArrays: blocks
-
 struct BlockSparseArray{
   T,
   N,
@@ -27,40 +24,48 @@ const BlockSparseVector{T,A<:AbstractVector{T},Blocks<:AbstractVector{A},Axes<:T
   T,1,A,Blocks,Axes
 }
 
+# TODO: Rename to `sparsemortar`.
 function BlockSparseArray(
   block_data::Dictionary{<:Block{N},<:AbstractArray{<:Any,N}},
   axes::Tuple{Vararg{AbstractUnitRange,N}},
 ) where {N}
   blocks = default_blocks(block_data, axes)
+  # TODO: Rename to `sparsemortar`.
   return BlockSparseArray(blocks, axes)
 end
 
+# TODO: Rename to `sparsemortar`.
 function BlockSparseArray(
   block_indices::Vector{<:Block{N}},
   block_data::Vector{<:AbstractArray{<:Any,N}},
   axes::Tuple{Vararg{AbstractUnitRange,N}},
 ) where {N}
+  # TODO: Rename to `sparsemortar`.
   return BlockSparseArray(Dictionary(block_indices, block_data), axes)
 end
 
+# TODO: Rename to `sparsemortar`.
 function BlockSparseArray{T,N,A,Blocks}(
   blocks::AbstractArray{<:AbstractArray{T,N},N}, axes::Tuple{Vararg{AbstractUnitRange,N}}
 ) where {T,N,A<:AbstractArray{T,N},Blocks<:AbstractArray{A,N}}
   return BlockSparseArray{T,N,A,Blocks,typeof(axes)}(blocks, axes)
 end
 
+# TODO: Rename to `sparsemortar`.
 function BlockSparseArray{T,N,A}(
   blocks::AbstractArray{<:AbstractArray{T,N},N}, axes::Tuple{Vararg{AbstractUnitRange,N}}
 ) where {T,N,A<:AbstractArray{T,N}}
   return BlockSparseArray{T,N,A,typeof(blocks)}(blocks, axes)
 end
 
+# TODO: Rename to `sparsemortar`.
 function BlockSparseArray{T,N}(
   blocks::AbstractArray{<:AbstractArray{T,N},N}, axes::Tuple{Vararg{AbstractUnitRange,N}}
 ) where {T,N}
   return BlockSparseArray{T,N,eltype(blocks),typeof(blocks),typeof(axes)}(blocks, axes)
 end
 
+# TODO: Rename to `sparsemortar`.
 function BlockSparseArray{T,N}(
   block_data::Dictionary{Block{N,Int},<:AbstractArray{T,N}},
   axes::Tuple{Vararg{AbstractUnitRange,N}},
@@ -70,99 +75,88 @@ function BlockSparseArray{T,N}(
 end
 
 function BlockSparseArray{T,N,A}(
-  axes::Tuple{Vararg{AbstractUnitRange,N}}
+  ::UndefInitializer, axes::Tuple{Vararg{AbstractUnitRange,N}}
 ) where {T,N,A<:AbstractArray{T,N}}
   blocks = default_blocks(A, axes)
   return BlockSparseArray{T,N,A}(blocks, axes)
 end
 
 function BlockSparseArray{T,N,A}(
-  axes::Vararg{AbstractUnitRange,N}
+  ::UndefInitializer, axes::Vararg{AbstractUnitRange,N}
 ) where {T,N,A<:AbstractArray{T,N}}
-  return BlockSparseArray{T,N,A}(axes)
+  return BlockSparseArray{T,N,A}(undef, axes)
 end
 
 function BlockSparseArray{T,N,A}(
-  dims::Tuple{Vararg{Vector{Int},N}}
+  ::UndefInitializer, dims::Tuple{Vararg{Vector{Int},N}}
 ) where {T,N,A<:AbstractArray{T,N}}
-  return BlockSparseArray{T,N,A}(blockedrange.(dims))
+  return BlockSparseArray{T,N,A}(undef, blockedrange.(dims))
 end
 
 # Fix ambiguity error.
-function BlockSparseArray{T,0,A}(axes::Tuple{}) where {T,A<:AbstractArray{T,0}}
+function BlockSparseArray{T,0,A}(
+  ::UndefInitializer, axes::Tuple{}
+) where {T,A<:AbstractArray{T,0}}
   blocks = default_blocks(A, axes)
   return BlockSparseArray{T,0,A}(blocks, axes)
 end
 
 function BlockSparseArray{T,N,A}(
-  dims::Vararg{Vector{Int},N}
+  ::UndefInitializer, dims::Vararg{Vector{Int},N}
 ) where {T,N,A<:AbstractArray{T,N}}
-  return BlockSparseArray{T,N,A}(dims)
-end
-
-function BlockSparseArray{T,N}(axes::Tuple{Vararg{AbstractUnitRange,N}}) where {T,N}
-  return BlockSparseArray{T,N,default_arraytype(T, axes)}(axes)
-end
-
-function BlockSparseArray{T,N}(axes::Vararg{AbstractUnitRange,N}) where {T,N}
-  return BlockSparseArray{T,N}(axes)
+  return BlockSparseArray{T,N,A}(undef, dims)
 end
 
-function BlockSparseArray{T,0}(axes::Tuple{}) where {T}
-  return BlockSparseArray{T,0,default_arraytype(T, axes)}(axes)
-end
-
-function BlockSparseArray{T,N}(dims::Tuple{Vararg{Vector{Int},N}}) where {T,N}
-  return BlockSparseArray{T,N}(blockedrange.(dims))
-end
-
-function BlockSparseArray{T,N}(dims::Vararg{Vector{Int},N}) where {T,N}
-  return BlockSparseArray{T,N}(dims)
+function BlockSparseArray{T,N}(
+  ::UndefInitializer, axes::Tuple{Vararg{AbstractUnitRange,N}}
+) where {T,N}
+  return BlockSparseArray{T,N,default_arraytype(T, axes)}(undef, axes)
 end
 
-function BlockSparseArray{T}(dims::Tuple{Vararg{Vector{Int}}}) where {T}
-  return BlockSparseArray{T,length(dims)}(dims)
+function BlockSparseArray{T,N}(
+  ::UndefInitializer, axes::Vararg{AbstractUnitRange,N}
+) where {T,N}
+  return BlockSparseArray{T,N}(undef, axes)
 end
 
-function BlockSparseArray{T}(axes::Tuple{Vararg{AbstractUnitRange}}) where {T}
-  return BlockSparseArray{T,length(axes)}(axes)
+function BlockSparseArray{T,0}(::UndefInitializer, axes::Tuple{}) where {T}
+  return BlockSparseArray{T,0,default_arraytype(T, axes)}(undef, axes)
 end
 
-function BlockSparseArray{T}(axes::Tuple{}) where {T}
-  return BlockSparseArray{T,length(axes)}(axes)
+function BlockSparseArray{T,N}(
+  ::UndefInitializer, dims::Tuple{Vararg{Vector{Int},N}}
+) where {T,N}
+  return BlockSparseArray{T,N}(undef, blockedrange.(dims))
 end
 
-function BlockSparseArray{T}(dims::Vararg{Vector{Int}}) where {T}
-  return BlockSparseArray{T}(dims)
+function BlockSparseArray{T,N}(::UndefInitializer, dims::Vararg{Vector{Int},N}) where {T,N}
+  return BlockSparseArray{T,N}(undef, dims)
 end
 
-function BlockSparseArray{T}(axes::Vararg{AbstractUnitRange}) where {T}
-  return BlockSparseArray{T}(axes)
+function BlockSparseArray{T}(::UndefInitializer, dims::Tuple{Vararg{Vector{Int}}}) where {T}
+  return BlockSparseArray{T,length(dims)}(undef, dims)
 end
 
-function BlockSparseArray{T}() where {T}
-  return BlockSparseArray{T}(())
+function BlockSparseArray{T}(
+  ::UndefInitializer, axes::Tuple{Vararg{AbstractUnitRange}}
+) where {T}
+  return BlockSparseArray{T,length(axes)}(undef, axes)
 end
 
-# undef
-function BlockSparseArray{T,N,A,Blocks}(
-  ::UndefInitializer, args...
-) where {T,N,A<:AbstractArray{T,N},Blocks<:AbstractArray{A,N}}
-  return BlockSparseArray{T,N,A,Blocks}(args...)
+function BlockSparseArray{T}(::UndefInitializer, axes::Tuple{}) where {T}
+  return BlockSparseArray{T,length(axes)}(undef, axes)
 end
 
-function BlockSparseArray{T,N,A}(
-  ::UndefInitializer, args...
-) where {T,N,A<:AbstractArray{T,N}}
-  return BlockSparseArray{T,N,A}(args...)
+function BlockSparseArray{T}(::UndefInitializer, dims::Vararg{Vector{Int}}) where {T}
+  return BlockSparseArray{T}(undef, dims)
 end
 
-function BlockSparseArray{T,N}(::UndefInitializer, args...) where {T,N}
-  return BlockSparseArray{T,N}(args...)
+function BlockSparseArray{T}(::UndefInitializer, axes::Vararg{AbstractUnitRange}) where {T}
+  return BlockSparseArray{T}(undef, axes)
 end
 
-function BlockSparseArray{T}(::UndefInitializer, args...) where {T}
-  return BlockSparseArray{T}(args...)
+function BlockSparseArray{T}(::UndefInitializer) where {T}
+  return BlockSparseArray{T}(undef, ())
 end
 
 # Base `AbstractArray` interface
diff --git a/test/test_basics.jl b/test/test_basics.jl
index da93d22..2ae7e16 100644
--- a/test/test_basics.jl
+++ b/test/test_basics.jl
@@ -48,20 +48,20 @@ arrayts = (Array, JLArray)
   dev(a) = adapt(arrayt, a)
   @testset "Broken" begin
     # TODO: Fix this and turn it into a proper test.
-    a = dev(BlockSparseArray{elt}([2, 3], [2, 3]))
+    a = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))
     a[Block(1, 1)] = dev(randn(elt, 2, 2))
     a[Block(2, 2)] = dev(randn(elt, 3, 3))
     @test_broken a[:, 4]
 
     # TODO: Fix this and turn it into a proper test.
-    a = dev(BlockSparseArray{elt}([2, 3], [2, 3]))
+    a = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))
     a[Block(1, 1)] = dev(randn(elt, 2, 2))
     a[Block(2, 2)] = dev(randn(elt, 3, 3))
     @test_broken a[:, [2, 4]]
     @test_broken a[[3, 5], [2, 4]]
 
     # TODO: Fix this and turn it into a proper test.
-    a = dev(BlockSparseArray{elt}([2, 3], [2, 3]))
+    a = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))
     a[Block(1, 1)] = dev(randn(elt, 2, 2))
     a[Block(2, 2)] = dev(randn(elt, 3, 3))
     @allowscalar @test a[2:4, 4] == Array(a)[2:4, 4]
@@ -85,10 +85,6 @@ arrayts = (Array, JLArray)
       ## BlockSparseMatrix{elt,Matrix{elt},SparseMatrixDOK{Matrix{elt}}}, # TODO
     )
       for args in (
-        bs,
-        (bs,),
-        blockedrange.(bs),
-        (blockedrange.(bs),),
         (undef, bs),
         (undef, bs...),
         (undef, blockedrange.(bs)),
@@ -117,10 +113,6 @@ arrayts = (Array, JLArray)
       ## BlockSparseVector{elt,Vector{elt},SparseVectorDOK{Vector{elt}}}, # TODO
     )
       for args in (
-        bs,
-        (bs,),
-        blockedrange.(bs),
-        (blockedrange.(bs),),
         (undef, bs),
         (undef, bs...),
         (undef, blockedrange.(bs)),
@@ -146,7 +138,7 @@ arrayts = (Array, JLArray)
     # TODO: This is difficult to determine just from type information.
     @test_broken blocktype(typeof(a)) <: SubArray{elt,2,arrayt{elt,2}}
 
-    a = BlockSparseMatrix{elt,arrayt{elt,2}}([1, 1], [1, 1])
+    a = BlockSparseMatrix{elt,arrayt{elt,2}}(undef, [1, 1], [1, 1])
     @test (@constinferred blockstype(a)) <: SparseMatrixDOK{arrayt{elt,2}}
     @test (@constinferred blockstype(typeof(a))) <: SparseMatrixDOK{arrayt{elt,2}}
     @test (@constinferred blocktype(a)) <: arrayt{elt,2}
@@ -167,9 +159,9 @@ arrayts = (Array, JLArray)
     @test_broken blocktype(typeof(a)) <: SubArray{elt,2,arrayt{elt,2}}
   end
   @testset "Basics" begin
-    a = dev(BlockSparseArray{elt}([2, 3], [2, 3]))
+    a = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))
     @allowscalar @test a == dev(
-      BlockSparseArray{elt}(blockedrange([2, 3]), blockedrange([2, 3]))
+      BlockSparseArray{elt}(undef, blockedrange([2, 3]), blockedrange([2, 3]))
     )
     @test eltype(a) === elt
     @test axes(a) == (1:5, 1:5)
@@ -182,7 +174,7 @@ arrayts = (Array, JLArray)
     @allowscalar @test all(I -> iszero(a[I]), eachindex(a))
     @test_throws DimensionMismatch a[Block(1, 1)] = randn(elt, 2, 3)
 
-    a = BlockSparseArray{elt}([2, 3], [2, 3])
+    a = BlockSparseArray{elt}(undef, [2, 3], [2, 3])
     a[3, 3] = 33
     @test eltype(a) === elt
     @test axes(a) == (1:5, 1:5)
@@ -201,7 +193,7 @@ arrayts = (Array, JLArray)
       end
     end
 
-    a = BlockSparseArray{elt}([2, 3], [2, 3])
+    a = BlockSparseArray{elt}(undef, [2, 3], [2, 3])
     a[Block(2, 1)] = randn(elt, 3, 2)
     a[Block(1, 2)] = randn(elt, 2, 3)
     @test issetequal(eachstoredblock(a), [a[Block(2, 1)], a[Block(1, 2)]])
@@ -211,7 +203,7 @@ arrayts = (Array, JLArray)
     @test isnan(norm(a))
 
     # Empty constructor
-    for a in (dev(BlockSparseArray{elt}()), dev(BlockSparseArray{elt}(undef)))
+    for a in (dev(BlockSparseArray{elt}(undef)),)
       @test size(a) == ()
       @test isone(length(a))
       @test blocksize(a) == ()
@@ -247,7 +239,7 @@ arrayts = (Array, JLArray)
     end
 
     @testset "Transpose" begin
-      a = dev(BlockSparseArray{elt}([2, 2], [3, 3, 1]))
+      a = dev(BlockSparseArray{elt}(undef, [2, 2], [3, 3, 1]))
       a[Block(1, 1)] = dev(randn(elt, 2, 3))
       a[Block(2, 3)] = dev(randn(elt, 2, 1))
 
@@ -271,7 +263,7 @@ arrayts = (Array, JLArray)
     end
 
     @testset "Adjoint" begin
-      a = dev(BlockSparseArray{elt}([2, 2], [3, 3, 1]))
+      a = dev(BlockSparseArray{elt}(undef, [2, 2], [3, 3, 1]))
       a[Block(1, 1)] = dev(randn(elt, 2, 3))
       a[Block(2, 3)] = dev(randn(elt, 2, 1))
 
@@ -305,12 +297,12 @@ arrayts = (Array, JLArray)
 
     # TODO: Broken on GPU.
     if arrayt ≠ Array
-      a = dev(BlockSparseArray{elt}([2, 3], [3, 4]))
+      a = dev(BlockSparseArray{elt}(undef, [2, 3], [3, 4]))
       @test_broken a[Block(1, 2)] .= 2
     end
 
     # TODO: Broken on GPU.
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
+    a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
     a[Block(1, 2)] .= 2
     @test eltype(a) == elt
     @test all(==(2), a[Block(1, 2)])
@@ -322,12 +314,12 @@ arrayts = (Array, JLArray)
 
     # TODO: Broken on GPU.
     if arrayt ≠ Array
-      a = dev(BlockSparseArray{elt}([2, 3], [3, 4]))
+      a = dev(BlockSparseArray{elt}(undef, [2, 3], [3, 4]))
       @test_broken a[Block(1, 2)] .= 0
     end
 
     # TODO: Broken on GPU.
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
+    a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
     a[Block(1, 2)] .= 0
     @test eltype(a) == elt
     @test iszero(a[Block(1, 1)])
@@ -367,7 +359,7 @@ arrayts = (Array, JLArray)
     @test size(b) == size(a)
     @test blocksize(b) == blocksize(a)
 
-    a = dev(BlockSparseArray{elt}([2, 3], [3, 4]))
+    a = dev(BlockSparseArray{elt}(undef, [2, 3], [3, 4]))
     b = @view a[[Block(2), Block(1)], [Block(2), Block(1)]]
     c = @view b[Block(1, 1)]
     @test iszero(a)
@@ -468,7 +460,7 @@ arrayts = (Array, JLArray)
     @test blockstoredlength(b) == 2
     @test storedlength(b) == 2 * 4 + 3 * 3
 
-    a = dev(BlockSparseArray{elt}([1, 1, 1], [1, 2, 3], [2, 2, 1], [1, 2, 1]))
+    a = dev(BlockSparseArray{elt}(undef, [1, 1, 1], [1, 2, 3], [2, 2, 1], [1, 2, 1]))
     a[Block(3, 2, 2, 3)] = dev(randn(elt, 1, 2, 2, 1))
     perm = (2, 3, 4, 1)
     for b in (PermutedDimsArray(a, perm), permutedims(a, perm))
@@ -623,7 +615,7 @@ arrayts = (Array, JLArray)
     b[Block(1, 1)] .= 1
     @test b[Block(1, 1)] == trues(blocksizes(b)[1, 1])
 
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
+    a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
     b = @view a[Block(2, 2)]
     @test size(b) == (3, 4)
     for i in parentindices(b)
@@ -632,7 +624,7 @@ arrayts = (Array, JLArray)
     @test parentindices(b)[1] == 1:3
     @test parentindices(b)[2] == 1:4
 
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
+    a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
     b = @view a[Block(2, 2)[1:2, 2:2]]
     @test size(b) == (2, 1)
     for i in parentindices(b)
@@ -655,7 +647,7 @@ arrayts = (Array, JLArray)
     @test @view(a[Block(2, 2)])[1:1, 1:2] == x
     @test a[3:3, 4:5] == x
 
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
+    a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
     b = @views a[Block(2, 2)][1:2, 2:3]
     @test b isa SubArray{<:Any,<:Any,<:BlockView}
     for i in parentindices(b)
@@ -667,7 +659,7 @@ arrayts = (Array, JLArray)
     @test a[Block(2, 2)[1:2, 2:3]] == b
     @test blockstoredlength(a) == 1
 
-    a = BlockSparseArray{elt}([2, 3], [2, 3])
+    a = BlockSparseArray{elt}(undef, [2, 3], [2, 3])
     @views for b in [Block(1, 1), Block(2, 2)]
       a[b] = randn(elt, size(a[b]))
     end
@@ -681,7 +673,7 @@ arrayts = (Array, JLArray)
       end
     end
 
-    a = dev(BlockSparseArray{elt}([2, 3], [2, 3]))
+    a = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))
     @views for b in [Block(1, 1), Block(2, 2)]
       # TODO: Use `blocksizes(a)[Int.(Tuple(b))...]` once available.
       a[b] = dev(randn(elt, size(a[b])))
@@ -768,7 +760,7 @@ arrayts = (Array, JLArray)
     @test !isassigned(a, Block(1)[1], Block(0)[1])
     @test !isassigned(a, Block(3)[3], Block(2)[4])
 
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
+    a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
     @test iszero(a)
     @test iszero(blockstoredlength(a))
     fill!(a, 0)
@@ -782,7 +774,7 @@ arrayts = (Array, JLArray)
     @test iszero(a)
     @test iszero(blockstoredlength(a))
 
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
+    a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
     @test iszero(a)
     @test iszero(blockstoredlength(a))
     zero!(a)
@@ -796,7 +788,7 @@ arrayts = (Array, JLArray)
     @test iszero(a)
     @test iszero(blockstoredlength(a))
 
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
+    a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
     @test iszero(a)
     @test iszero(blockstoredlength(a))
     a .= 0
@@ -811,7 +803,7 @@ arrayts = (Array, JLArray)
     @test iszero(blockstoredlength(a))
 
     # TODO: Broken on GPU.
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
+    a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
     for I in (Block.(1:2), [Block(1), Block(2)])
       b = @view a[I, I]
       x = randn(elt, 3, 4)
@@ -831,14 +823,14 @@ arrayts = (Array, JLArray)
     end
 
     function f1()
-      a = BlockSparseArray{elt}([2, 3], [3, 4])
+      a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
       b = @view a[[Block(2), Block(1)], [Block(2), Block(1)]]
       x = randn(elt, 3, 4)
       b[Block(1, 1)] .= x
       return (; a, b, x)
     end
     function f2()
-      a = BlockSparseArray{elt}([2, 3], [3, 4])
+      a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
       b = @view a[[Block(2), Block(1)], [Block(2), Block(1)]]
       x = randn(elt, 3, 4)
       b[Block(1, 1)] = x
@@ -861,7 +853,7 @@ arrayts = (Array, JLArray)
       @test_throws DimensionMismatch b[Block(1, 1)] .= randn(2, 3)
     end
 
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
+    a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
     b = @views a[[Block(2), Block(1)], [Block(2), Block(1)]][Block(2, 1)]
     @test iszero(b)
     @test size(b) == (2, 4)
@@ -871,7 +863,7 @@ arrayts = (Array, JLArray)
     @test a[Block(1, 2)] == x
     @test blockstoredlength(a) == 1
 
-    a = BlockSparseArray{elt}([4, 3, 2], [4, 3, 2])
+    a = BlockSparseArray{elt}(undef, [4, 3, 2], [4, 3, 2])
     @views for B in [Block(1, 1), Block(2, 2), Block(3, 3)]
       a[B] = randn(elt, size(a[B]))
     end
@@ -897,13 +889,13 @@ arrayts = (Array, JLArray)
     @test c[Block(2, 2)] == x
     @test a[Block(1, 1)[1:3, 1:3]] == x
 
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
+    a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
     b = @view a[[Block(2), Block(1)], [Block(2), Block(1)]]
     for index in parentindices(@view(b[Block(1, 1)]))
       @test index isa Base.OneTo{Int}
     end
 
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
+    a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
     a[Block(1, 1)] = randn(elt, 2, 3)
     b = @view a[Block(1, 1)[1:2, 1:1]]
     @test b isa SubArray{elt,2,Matrix{elt}}
@@ -911,7 +903,7 @@ arrayts = (Array, JLArray)
       @test i isa UnitRange{Int}
     end
 
-    a = BlockSparseArray{elt}([2, 2, 2, 2], [2, 2, 2, 2])
+    a = BlockSparseArray{elt}(undef, [2, 2, 2, 2], [2, 2, 2, 2])
     @views for I in [Block(1, 1), Block(2, 2), Block(3, 3), Block(4, 4)]
       a[I] = randn(elt, size(a[I]))
     end
@@ -932,7 +924,7 @@ arrayts = (Array, JLArray)
       @test c == a[Block.(3:4), Block.(3:4)]
     end
 
-    a = BlockSparseArray{elt}([2, 3], [2, 3])
+    a = BlockSparseArray{elt}(undef, [2, 3], [2, 3])
     a[Block(1, 1)] = randn(elt, 2, 2)
     a[Block(2, 2)] = randn(elt, 3, 3)
     for I in (mortar([Block(1)[2:2], Block(2)[2:3]]), [Block(1)[2:2], Block(2)[2:3]])
@@ -941,7 +933,7 @@ arrayts = (Array, JLArray)
     end
 
     # Merge and permute blocks.
-    a = BlockSparseArray{elt}([2, 2, 2, 2], [2, 2, 2, 2])
+    a = BlockSparseArray{elt}(undef, [2, 2, 2, 2], [2, 2, 2, 2])
     @views for I in [Block(1, 1), Block(2, 2), Block(3, 3), Block(4, 4)]
       a[I] = randn(elt, size(a[I]))
     end
@@ -966,7 +958,7 @@ arrayts = (Array, JLArray)
 
     # TODO: Add more tests of this, it may
     # only be working accidentally.
-    a = BlockSparseArray{elt}([2, 3], [2, 3])
+    a = BlockSparseArray{elt}(undef, [2, 3], [2, 3])
     a[Block(1, 1)] = randn(elt, 2, 2)
     a[Block(2, 2)] = randn(elt, 3, 3)
     @test a[2:4, 4] == Array(a)[2:4, 4]
@@ -975,7 +967,7 @@ arrayts = (Array, JLArray)
   end
   @testset "view!" begin
     for blk in ((Block(2, 2),), (Block(2), Block(2)))
-      a = BlockSparseArray{elt}([2, 3], [2, 3])
+      a = BlockSparseArray{elt}(undef, [2, 3], [2, 3])
       b = view!(a, blk...)
       x = randn(elt, 3, 3)
       b .= x
@@ -986,7 +978,7 @@ arrayts = (Array, JLArray)
       @test @view!(a[blk...]) == x
     end
     for blk in ((Block(2, 2),), (Block(2), Block(2)))
-      a = BlockSparseArray{elt}([2, 3], [2, 3])
+      a = BlockSparseArray{elt}(undef, [2, 3], [2, 3])
       b = @view! a[blk...]
       x = randn(elt, 3, 3)
       b .= x
@@ -997,7 +989,7 @@ arrayts = (Array, JLArray)
       @test @view!(a[blk...]) == x
     end
     for blk in ((Block(2, 2)[2:3, 1:2],), (Block(2)[2:3], Block(2)[1:2]))
-      a = BlockSparseArray{elt}([2, 3], [2, 3])
+      a = BlockSparseArray{elt}(undef, [2, 3], [2, 3])
       b = view!(a, blk...)
       x = randn(elt, 2, 2)
       b .= x
@@ -1008,7 +1000,7 @@ arrayts = (Array, JLArray)
       @test @view!(a[blk...]) == x
     end
     for blk in ((Block(2, 2)[2:3, 1:2],), (Block(2)[2:3], Block(2)[1:2]))
-      a = BlockSparseArray{elt}([2, 3], [2, 3])
+      a = BlockSparseArray{elt}(undef, [2, 3], [2, 3])
       b = @view! a[blk...]
       x = randn(elt, 2, 2)
       b .= x
@@ -1020,9 +1012,9 @@ arrayts = (Array, JLArray)
     end
   end
   @testset "LinearAlgebra" begin
-    a1 = dev(BlockSparseArray{elt}([2, 3], [2, 3]))
+    a1 = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))
     a1[Block(1, 1)] = dev(randn(elt, size(@view(a1[Block(1, 1)]))))
-    a2 = dev(BlockSparseArray{elt}([2, 3], [2, 3]))
+    a2 = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))
     a2[Block(1, 1)] = dev(randn(elt, size(@view(a1[Block(1, 1)]))))
     a_dest = a1 * a2
     @allowscalar @test Array(a_dest) ≈ Array(a1) * Array(a2)
@@ -1030,10 +1022,10 @@ arrayts = (Array, JLArray)
     @test blockstoredlength(a_dest) == 1
   end
   @testset "Matrix multiplication" begin
-    a1 = dev(BlockSparseArray{elt}([2, 3], [2, 3]))
+    a1 = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))
     a1[Block(1, 2)] = dev(randn(elt, size(@view(a1[Block(1, 2)]))))
     a1[Block(2, 1)] = dev(randn(elt, size(@view(a1[Block(2, 1)]))))
-    a2 = dev(BlockSparseArray{elt}([2, 3], [2, 3]))
+    a2 = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))
     a2[Block(1, 2)] = dev(randn(elt, size(@view(a2[Block(1, 2)]))))
     a2[Block(2, 1)] = dev(randn(elt, size(@view(a2[Block(2, 1)]))))
     for (a1′, a2′) in ((a1, a2), (a1', a2), (a1, a2'), (a1', a2'))
@@ -1042,19 +1034,19 @@ arrayts = (Array, JLArray)
     end
   end
   @testset "Dot product" begin
-    a1 = dev(BlockSparseArray{elt}([2, 3, 4]))
+    a1 = dev(BlockSparseArray{elt}(undef, [2, 3, 4]))
     a1[Block(1)] = dev(randn(elt, size(@view(a1[Block(1)]))))
     a1[Block(3)] = dev(randn(elt, size(@view(a1[Block(3)]))))
-    a2 = dev(BlockSparseArray{elt}([2, 3, 4]))
+    a2 = dev(BlockSparseArray{elt}(undef, [2, 3, 4]))
     a2[Block(2)] = dev(randn(elt, size(@view(a1[Block(2)]))))
     a2[Block(3)] = dev(randn(elt, size(@view(a1[Block(3)]))))
     @test a1' * a2 ≈ Array(a1)' * Array(a2)
     @test dot(a1, a2) ≈ a1' * a2
   end
   @testset "cat" begin
-    a1 = dev(BlockSparseArray{elt}([2, 3], [2, 3]))
+    a1 = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))
     a1[Block(2, 1)] = dev(randn(elt, size(@view(a1[Block(2, 1)]))))
-    a2 = dev(BlockSparseArray{elt}([2, 3], [2, 3]))
+    a2 = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))
     a2[Block(1, 2)] = dev(randn(elt, size(@view(a2[Block(1, 2)]))))
 
     a_dest = cat(a1, a2; dims=1)
@@ -1079,9 +1071,9 @@ arrayts = (Array, JLArray)
     @test a_dest[Block(3, 4)] == a2[Block(1, 2)]
   end
   @testset "TensorAlgebra" begin
-    a1 = dev(BlockSparseArray{elt}([2, 3], [2, 3]))
+    a1 = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))
     a1[Block(1, 1)] = dev(randn(elt, size(@view(a1[Block(1, 1)]))))
-    a2 = dev(BlockSparseArray{elt}([2, 3], [2, 3]))
+    a2 = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))
     a2[Block(1, 1)] = dev(randn(elt, size(@view(a1[Block(1, 1)]))))
     # TODO: Make this work, requires customization of `TensorAlgebra.fusedims` and
     # `TensorAlgebra.splitdims` in terms of `BlockSparseArrays.blockreshape`,
@@ -1107,7 +1099,7 @@ arrayts = (Array, JLArray)
     matrixt_elt = arrayt{elt,2}
     arrayt_elt = arrayt{elt,3}
 
-    a = BlockSparseVector{elt,arrayt{elt,1}}([2, 2])
+    a = BlockSparseVector{elt,arrayt{elt,1}}(undef, [2, 2])
     res = sprint(summary, a)
     function ref_vec(elt, arrayt, prefix="")
       return "2-blocked 4-element $(prefix)BlockSparseVector{$(elt), $(arrayt), …, …}"
@@ -1116,7 +1108,7 @@ arrayts = (Array, JLArray)
     @test (res == ref_vec(elt, vectort_elt)) ||
       (res == ref_vec(elt, vectort_elt, "BlockSparseArrays."))
 
-    a = BlockSparseMatrix{elt,arrayt{elt,2}}([2, 2], [2, 2])
+    a = BlockSparseMatrix{elt,arrayt{elt,2}}(undef, [2, 2], [2, 2])
     res = sprint(summary, a)
     function ref_mat(elt, arrayt, prefix="")
       return "2×2-blocked 4×4 $(prefix)BlockSparseMatrix{$(elt), $(arrayt), …, …}"
@@ -1125,7 +1117,7 @@ arrayts = (Array, JLArray)
     @test (res == ref_mat(elt, matrixt_elt)) ||
       (res == ref_mat(elt, matrixt_elt, "BlockSparseArrays."))
 
-    a = BlockSparseArray{elt,3,arrayt{elt,3}}([2, 2], [2, 2], [2, 2])
+    a = BlockSparseArray{elt,3,arrayt{elt,3}}(undef, [2, 2], [2, 2], [2, 2])
     res = sprint(summary, a)
     function ref_arr(elt, arrayt, prefix="")
       return "2×2×2-blocked 4×4×4 $(prefix)BlockSparseArray{$(elt), 3, $(arrayt), …, …}"
@@ -1136,7 +1128,7 @@ arrayts = (Array, JLArray)
     if elt === Float64
       # Not testing other element types since they change the
       # spacing so it isn't easy to make the test general.
-      a = BlockSparseMatrix{elt,arrayt{elt,2}}([2, 2], [2, 2])
+      a = BlockSparseMatrix{elt,arrayt{elt,2}}(undef, [2, 2], [2, 2])
       @allowscalar a[1, 2] = 12
       @test sprint(show, "text/plain", a) ==
         "$(summary(a)):\n $(zero(eltype(a)))  $(eltype(a)(12))  │   ⋅    ⋅ \n $(zero(eltype(a)))   $(zero(eltype(a)))  │   ⋅    ⋅ \n ───────────┼──────────\n  ⋅     ⋅   │   ⋅    ⋅ \n  ⋅     ⋅   │   ⋅    ⋅ "

From acef25f86f9ed7f34b750c22aec86ff16e675c68 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Mon, 3 Mar 2025 13:01:19 -0500
Subject: [PATCH 02/14] Fix more tests

---
 .../BlockSparseArraysGradedUnitRangesExt.jl   |  2 +-
 src/blocksparsearray/blocksparsearray.jl      | 77 ++++++++++---------
 src/blocksparsearrayinterface/map.jl          |  2 +-
 test/test_basics.jl                           |  1 +
 test/test_gradedunitrangesext.jl              | 22 +++---
 test/test_svd.jl                              |  2 +-
 test/test_tensoralgebraext.jl                 |  2 +-
 7 files changed, 55 insertions(+), 53 deletions(-)

diff --git a/ext/BlockSparseArraysGradedUnitRangesExt/BlockSparseArraysGradedUnitRangesExt.jl b/ext/BlockSparseArraysGradedUnitRangesExt/BlockSparseArraysGradedUnitRangesExt.jl
index 376f53e..49762ed 100644
--- a/ext/BlockSparseArraysGradedUnitRangesExt/BlockSparseArraysGradedUnitRangesExt.jl
+++ b/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
 
diff --git a/src/blocksparsearray/blocksparsearray.jl b/src/blocksparsearray/blocksparsearray.jl
index f87d943..1de153a 100644
--- a/src/blocksparsearray/blocksparsearray.jl
+++ b/src/blocksparsearray/blocksparsearray.jl
@@ -24,61 +24,62 @@ const BlockSparseVector{T,A<:AbstractVector{T},Blocks<:AbstractVector{A},Axes<:T
   T,1,A,Blocks,Axes
 }
 
-# TODO: Rename to `sparsemortar`.
-function BlockSparseArray(
+"""
+    sparsemortar(blocks::AbstractArray{<:AbstractArray{T,N},N}, axes) -> ::BlockSparseArray{T,N}
+
+Construct a block sparse array from a sparse array of arrays and specified blocked axes.
+The block sizes must be commensurate with the blocks of the axes.
+"""
+function sparsemortar(
+  blocks::AbstractArray{<:AbstractArray{T,N},N}, axes::Tuple{Vararg{AbstractUnitRange,N}}
+) where {T,N}
+  return BlockSparseArray{T,N,eltype(blocks),typeof(blocks),typeof(axes)}(blocks, axes)
+end
+
+"""
+    sparsemortar(blocks::Dictionary{<:Block{N},<:AbstractArray{T,N}}, axes) -> ::BlockSparseArray{T,N}
+
+Construct a block sparse array from a dictionary mapping the locations of the stored/nonzero
+blocks to the data of those blocks, along with a set of blocked axes.
+The block sizes must be commensurate with the blocks of the specified axes.
+"""
+function sparsemortar(
   block_data::Dictionary{<:Block{N},<:AbstractArray{<:Any,N}},
   axes::Tuple{Vararg{AbstractUnitRange,N}},
 ) where {N}
   blocks = default_blocks(block_data, axes)
-  # TODO: Rename to `sparsemortar`.
-  return BlockSparseArray(blocks, axes)
+  return sparsemortar(blocks, axes)
 end
 
-# TODO: Rename to `sparsemortar`.
-function BlockSparseArray(
+"""
+    sparsemortar(block_indices::Vector{<:Block{N}}, block_data::Vector{<:AbstractArray{T,N}}, axes) -> ::BlockSparseArray{T,N}
+
+Construct a block sparse array from a list of locations of the stored/nonzero blocks,
+a corresponding list of the data of those blocks, along with a set of blocked axes.
+The block sizes must be commensurate with the blocks of the specified axes.
+"""
+function sparsemortar(
   block_indices::Vector{<:Block{N}},
   block_data::Vector{<:AbstractArray{<:Any,N}},
   axes::Tuple{Vararg{AbstractUnitRange,N}},
 ) where {N}
-  # TODO: Rename to `sparsemortar`.
-  return BlockSparseArray(Dictionary(block_indices, block_data), axes)
+  return sparsemortar(Dictionary(block_indices, block_data), axes)
 end
 
-# TODO: Rename to `sparsemortar`.
-function BlockSparseArray{T,N,A,Blocks}(
-  blocks::AbstractArray{<:AbstractArray{T,N},N}, axes::Tuple{Vararg{AbstractUnitRange,N}}
-) where {T,N,A<:AbstractArray{T,N},Blocks<:AbstractArray{A,N}}
-  return BlockSparseArray{T,N,A,Blocks,typeof(axes)}(blocks, axes)
-end
+@doc """
+    BlockSparseArray{T}(undef, dims)
+    BlockSparseArray{T,N}(undef, dims)
+    BlockSparseArray{T,N,A}(undef, dims)
 
-# TODO: Rename to `sparsemortar`.
-function BlockSparseArray{T,N,A}(
-  blocks::AbstractArray{<:AbstractArray{T,N},N}, axes::Tuple{Vararg{AbstractUnitRange,N}}
-) where {T,N,A<:AbstractArray{T,N}}
-  return BlockSparseArray{T,N,A,typeof(blocks)}(blocks, axes)
-end
-
-# TODO: Rename to `sparsemortar`.
-function BlockSparseArray{T,N}(
-  blocks::AbstractArray{<:AbstractArray{T,N},N}, axes::Tuple{Vararg{AbstractUnitRange,N}}
-) where {T,N}
-  return BlockSparseArray{T,N,eltype(blocks),typeof(blocks),typeof(axes)}(blocks, axes)
-end
-
-# TODO: Rename to `sparsemortar`.
-function BlockSparseArray{T,N}(
-  block_data::Dictionary{Block{N,Int},<:AbstractArray{T,N}},
-  axes::Tuple{Vararg{AbstractUnitRange,N}},
-) where {T,N}
-  blocks = default_blocks(block_data, axes)
-  return BlockSparseArray{T,N}(blocks, axes)
-end
+Construct an uninitialized N-dimensional BlockSparseArray containing elements of type T. `dims` should be a list
+of block lengths in each dimension or a list of blocked ranges representing the axes.
+""" BlockSparseArray
 
 function BlockSparseArray{T,N,A}(
   ::UndefInitializer, axes::Tuple{Vararg{AbstractUnitRange,N}}
 ) where {T,N,A<:AbstractArray{T,N}}
   blocks = default_blocks(A, axes)
-  return BlockSparseArray{T,N,A}(blocks, axes)
+  return sparsemortar(blocks, axes)
 end
 
 function BlockSparseArray{T,N,A}(
@@ -98,7 +99,7 @@ function BlockSparseArray{T,0,A}(
   ::UndefInitializer, axes::Tuple{}
 ) where {T,A<:AbstractArray{T,0}}
   blocks = default_blocks(A, axes)
-  return BlockSparseArray{T,0,A}(blocks, axes)
+  return sparsemortar(blocks, axes)
 end
 
 function BlockSparseArray{T,N,A}(
diff --git a/src/blocksparsearrayinterface/map.jl b/src/blocksparsearrayinterface/map.jl
index b97bff8..6b5a446 100644
--- a/src/blocksparsearrayinterface/map.jl
+++ b/src/blocksparsearrayinterface/map.jl
@@ -112,5 +112,5 @@ function map_stored_blocks(f, a::AbstractArray)
   # since they can't necessarily by `Diagonal` if there are rectangular blocks.
   mapped_blocks = Dictionary{eltype(bs),eltype(ds)}(bs, ds)
   # TODO: Use `similartype(typeof(a), eltype(eltype(mapped_blocks)))(...)`.
-  return BlockSparseArray(mapped_blocks, axes(a))
+  return sparsemortar(mapped_blocks, axes(a))
 end
diff --git a/test/test_basics.jl b/test/test_basics.jl
index 2ae7e16..a52f1dc 100644
--- a/test/test_basics.jl
+++ b/test/test_basics.jl
@@ -30,6 +30,7 @@ using BlockSparseArrays:
   eachstoredblock,
   blockstype,
   blocktype,
+  sparsemortar,
   view!
 using GPUArraysCore: @allowscalar
 using JLArrays: JLArray
diff --git a/test/test_gradedunitrangesext.jl b/test/test_gradedunitrangesext.jl
index 4988028..fbbd029 100644
--- a/test/test_gradedunitrangesext.jl
+++ b/test/test_gradedunitrangesext.jl
@@ -20,7 +20,7 @@ using TensorAlgebra: fusedims, splitdims
 using LinearAlgebra: adjoint
 using Random: randn!
 function randn_blockdiagonal(elt::Type, axes::Tuple)
-  a = BlockSparseArray{elt}(axes)
+  a = BlockSparseArray{elt}(undef, axes)
   blockdiaglength = minimum(blocksize(a))
   for i in 1:blockdiaglength
     b = Block(ntuple(Returns(i), ndims(a)))
@@ -135,7 +135,7 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
   @testset "dual axes" begin
     r = gradedrange([U1(0) => 2, U1(1) => 2])
     for ax in ((r, r), (dual(r), r), (r, dual(r)), (dual(r), dual(r)))
-      a = BlockSparseArray{elt}(ax...)
+      a = BlockSparseArray{elt}(undef, ax...)
       @views for b in [Block(1, 1), Block(2, 2)]
         a[b] = randn(elt, size(a[b]))
       end
@@ -178,7 +178,7 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
 
     @testset "GradedOneTo" begin
       r = gradedrange([U1(0) => 2, U1(1) => 2])
-      a = BlockSparseArray{elt}(r, r)
+      a = BlockSparseArray{elt}(undef, r, r)
       @views for i in [Block(1, 1), Block(2, 2)]
         a[i] = randn(elt, size(a[i]))
       end
@@ -199,7 +199,7 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
 
     @testset "GradedUnitRange" begin
       r = gradedrange([U1(0) => 2, U1(1) => 2])[1:3]
-      a = BlockSparseArray{elt}(r, r)
+      a = BlockSparseArray{elt}(undef, r, r)
       @views for i in [Block(1, 1), Block(2, 2)]
         a[i] = randn(elt, size(a[i]))
       end
@@ -226,7 +226,7 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
     # Test case when all axes are dual.
     @testset "dual GradedOneTo" begin
       r = gradedrange([U1(-1) => 2, U1(1) => 2])
-      a = BlockSparseArray{elt}(dual(r), dual(r))
+      a = BlockSparseArray{elt}(undef, dual(r), dual(r))
       @views for i in [Block(1, 1), Block(2, 2)]
         a[i] = randn(elt, size(a[i]))
       end
@@ -251,7 +251,7 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
 
     @testset "dual GradedUnitRange" begin
       r = gradedrange([U1(0) => 2, U1(1) => 2])[1:3]
-      a = BlockSparseArray{elt}(dual(r), dual(r))
+      a = BlockSparseArray{elt}(undef, dual(r), dual(r))
       @views for i in [Block(1, 1), Block(2, 2)]
         a[i] = randn(elt, size(a[i]))
       end
@@ -277,7 +277,7 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
 
     @testset "dual BlockedUnitRange" begin    # self dual
       r = blockedrange([2, 2])
-      a = BlockSparseArray{elt}(dual(r), dual(r))
+      a = BlockSparseArray{elt}(undef, dual(r), dual(r))
       @views for i in [Block(1, 1), Block(2, 2)]
         a[i] = randn(elt, size(a[i]))
       end
@@ -302,7 +302,7 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
       gradedrange([U1(0) => 2, U1(1) => 2]),
       gradedrange([U1(0) => 2, U1(1) => 2])[begin:end],
     )
-      a = BlockSparseArray{elt}(r, r)
+      a = BlockSparseArray{elt}(undef, r, r)
       @views for i in [Block(1, 1), Block(2, 2)]
         a[i] = randn(elt, size(a[i]))
       end
@@ -330,16 +330,16 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
   end
   @testset "Matrix multiplication" begin
     r = gradedrange([U1(0) => 2, U1(1) => 3])
-    a1 = BlockSparseArray{elt}(dual(r), r)
+    a1 = BlockSparseArray{elt}(undef, dual(r), r)
     a1[Block(1, 2)] = randn(elt, size(@view(a1[Block(1, 2)])))
     a1[Block(2, 1)] = randn(elt, size(@view(a1[Block(2, 1)])))
-    a2 = BlockSparseArray{elt}(dual(r), r)
+    a2 = BlockSparseArray{elt}(undef, dual(r), r)
     a2[Block(1, 2)] = randn(elt, size(@view(a2[Block(1, 2)])))
     a2[Block(2, 1)] = randn(elt, size(@view(a2[Block(2, 1)])))
     @test Array(a1 * a2) ≈ Array(a1) * Array(a2)
     @test Array(a1' * a2') ≈ Array(a1') * Array(a2')
 
-    a2 = BlockSparseArray{elt}(r, dual(r))
+    a2 = BlockSparseArray{elt}(undef, r, dual(r))
     a2[Block(1, 2)] = randn(elt, size(@view(a2[Block(1, 2)])))
     a2[Block(2, 1)] = randn(elt, size(@view(a2[Block(2, 1)])))
     @test Array(a1' * a2) ≈ Array(a1') * Array(a2)
diff --git a/test/test_svd.jl b/test/test_svd.jl
index 8d8dfe7..73efa16 100644
--- a/test/test_svd.jl
+++ b/test/test_svd.jl
@@ -47,7 +47,7 @@ end
 # -----------
 @testset "($m, $n) BlockSparseMatrix{$T}" for ((m, n), T) in
                                               Iterators.product(blockszs, eltypes)
-  a = BlockSparseArray{T}(m, n)
+  a = BlockSparseArray{T}(undef, m, n)
 
   # test empty matrix
   usv_empty = svd(a)
diff --git a/test/test_tensoralgebraext.jl b/test/test_tensoralgebraext.jl
index 470efaa..3fe0b86 100644
--- a/test/test_tensoralgebraext.jl
+++ b/test/test_tensoralgebraext.jl
@@ -9,7 +9,7 @@ using SymmetrySectors: U1
 using TensorAlgebra: contract
 
 function randn_blockdiagonal(elt::Type, axes::Tuple)
-  a = BlockSparseArray{elt}(axes)
+  a = BlockSparseArray{elt}(undef, axes)
   blockdiaglength = minimum(blocksize(a))
   for i in 1:blockdiaglength
     b = Block(ntuple(Returns(i), ndims(a)))

From 16671a0b00a13353ab95f1a9b96d0716880512f2 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Mon, 3 Mar 2025 13:05:56 -0500
Subject: [PATCH 03/14] Update README

---
 README.md          | 8 ++++----
 examples/README.jl | 8 ++++----
 2 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/README.md b/README.md
index e2f415e..edfee83 100644
--- a/README.md
+++ b/README.md
@@ -35,7 +35,7 @@ julia> Pkg.add("BlockSparseArrays")
 
 ````julia
 using BlockArrays: BlockArrays, BlockedVector, Block, blockedrange
-using BlockSparseArrays: BlockSparseArray, blockstoredlength
+using BlockSparseArrays: BlockSparseArray, blockstoredlength, sparsemortar
 using Test: @test, @test_broken
 
 function main()
@@ -60,13 +60,13 @@ function main()
     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)
+  b = sparsemortar(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)
+  b = sparsemortar(nz_blocks, d_blocks, i_axes)
 
   @test blockstoredlength(b) == 2
 
@@ -114,7 +114,7 @@ using BlockSparseArrays: BlockSparseArray
 
 i1 = [2, 3]
 i2 = [2, 3]
-B = BlockSparseArray{Float64}(i1, i2)
+B = BlockSparseArray{Float64}(undef, i1, i2)
 B[Block(1, 1)] = randn(2, 2)
 B[Block(2, 2)] = randn(3, 3)
 
diff --git a/examples/README.jl b/examples/README.jl
index 2d00e5b..a1be77f 100644
--- a/examples/README.jl
+++ b/examples/README.jl
@@ -40,7 +40,7 @@ julia> Pkg.add("BlockSparseArrays")
 # ## Examples
 
 using BlockArrays: BlockArrays, BlockedVector, Block, blockedrange
-using BlockSparseArrays: BlockSparseArray, blockstoredlength
+using BlockSparseArrays: BlockSparseArray, blockstoredlength, sparsemortar
 using Test: @test, @test_broken
 
 function main()
@@ -65,13 +65,13 @@ function main()
     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)
+  b = sparsemortar(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)
+  b = sparsemortar(nz_blocks, d_blocks, i_axes)
 
   @test blockstoredlength(b) == 2
 
@@ -117,7 +117,7 @@ using BlockSparseArrays: BlockSparseArray
 
 i1 = [2, 3]
 i2 = [2, 3]
-B = BlockSparseArray{Float64}(i1, i2)
+B = BlockSparseArray{Float64}(undef, i1, i2)
 B[Block(1, 1)] = randn(2, 2)
 B[Block(2, 2)] = randn(3, 3)
 

From 81e2328f601fb1fc767695a78a9b385634a9f1d8 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Mon, 3 Mar 2025 13:12:58 -0500
Subject: [PATCH 04/14] Update exports

---
 src/BlockSparseArrays.jl | 3 ++-
 test/test_exports.jl     | 1 +
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/BlockSparseArrays.jl b/src/BlockSparseArrays.jl
index d58a9b7..b9ac7b0 100644
--- a/src/BlockSparseArrays.jl
+++ b/src/BlockSparseArrays.jl
@@ -5,7 +5,8 @@ export BlockSparseArray,
   BlockSparseVector,
   blockstoredlength,
   eachblockstoredindex,
-  eachstoredblock
+  eachstoredblock,
+  sparsemortar
 
 # factorizations
 include("factorizations/svd.jl")
diff --git a/test/test_exports.jl b/test/test_exports.jl
index 0e0d20e..7121c79 100644
--- a/test/test_exports.jl
+++ b/test/test_exports.jl
@@ -9,6 +9,7 @@ using Test: @test, @testset
     :blockstoredlength,
     :eachblockstoredindex,
     :eachstoredblock,
+    :sparsemortar,
   ]
   @test issetequal(names(BlockSparseArrays), exports)
 end

From 452d6749f6bb32c76f6e6fb8e2d91e20c26f7700 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Mon, 3 Mar 2025 21:21:21 -0500
Subject: [PATCH 05/14] compat versions

---
 docs/Project.toml     | 2 +-
 examples/Project.toml | 2 +-
 test/Project.toml     | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index 6dca2c9..1c7114e 100644
--- a/docs/Project.toml
+++ b/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"
diff --git a/examples/Project.toml b/examples/Project.toml
index 48948bf..cf876f1 100644
--- a/examples/Project.toml
+++ b/examples/Project.toml
@@ -5,5 +5,5 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 BlockArrays = "1"
-BlockSparseArrays = "0.2"
+BlockSparseArrays = "0.3"
 Test = "1"
diff --git a/test/Project.toml b/test/Project.toml
index 4c93abd..817fd10 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -26,7 +26,7 @@ Adapt = "4"
 Aqua = "0.8"
 ArrayLayouts = "1"
 BlockArrays = "1"
-BlockSparseArrays = "0.2"
+BlockSparseArrays = "0.3"
 DiagonalArrays = "0.2"
 GPUArraysCore = "0.2"
 GradedUnitRanges = "0.1"

From a63764fd937d7709eb99ffcd3a1d31aca4d69250 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Tue, 4 Mar 2025 12:29:29 -0500
Subject: [PATCH 06/14] Better constructor code logic

---
 README.md                                  | 134 ++-------------------
 examples/README.jl                         | 132 ++------------------
 src/blocksparsearray/blockdiagonalarray.jl |   2 +-
 src/blocksparsearray/blocksparsearray.jl   |  58 ++++-----
 src/blocksparsearrayinterface/map.jl       |  21 +---
 5 files changed, 52 insertions(+), 295 deletions(-)

diff --git a/README.md b/README.md
index edfee83..72a1fb9 100644
--- a/README.md
+++ b/README.md
@@ -34,129 +34,17 @@ julia> Pkg.add("BlockSparseArrays")
 ## Examples
 
 ````julia
-using BlockArrays: BlockArrays, BlockedVector, Block, blockedrange
-using BlockSparseArrays: BlockSparseArray, blockstoredlength, sparsemortar
-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 = sparsemortar(nz_blocks, d_blocks, i_axes)
-
-  @test blockstoredlength(b) == 2
-
-  # Blocks with discontiguous underlying data
-  d_blocks = randn.(nz_block_sizes)
-  b = sparsemortar(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}(undef, 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 BlockArrays: Block
+using BlockSparseArrays: BlockSparseArray, blockstoredlength
+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
 ````
 
 ---
diff --git a/examples/README.jl b/examples/README.jl
index a1be77f..ba0467b 100644
--- a/examples/README.jl
+++ b/examples/README.jl
@@ -39,124 +39,14 @@ julia> Pkg.add("BlockSparseArrays")
 
 # ## Examples
 
-using BlockArrays: BlockArrays, BlockedVector, Block, blockedrange
-using BlockSparseArrays: BlockSparseArray, blockstoredlength, sparsemortar
-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 = sparsemortar(nz_blocks, d_blocks, i_axes)
-
-  @test blockstoredlength(b) == 2
-
-  ## Blocks with discontiguous underlying data
-  d_blocks = randn.(nz_block_sizes)
-  b = sparsemortar(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}(undef, 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 BlockArrays: Block
+using BlockSparseArrays: BlockSparseArray, blockstoredlength
+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
diff --git a/src/blocksparsearray/blockdiagonalarray.jl b/src/blocksparsearray/blockdiagonalarray.jl
index 7d0f281..0926431 100644
--- a/src/blocksparsearray/blockdiagonalarray.jl
+++ b/src/blocksparsearray/blockdiagonalarray.jl
@@ -12,7 +12,7 @@ const BlockSparseDiagonal{T,A<:AbstractBlockSparseVector{T}} = Diagonal{T,A}
 end
 
 function BlockDiagonal(blocks::AbstractVector{<:AbstractMatrix})
-  return BlockSparseArray(
+  return sparsemortar(
     Diagonal(blocks), (blockedrange(size.(blocks, 1)), blockedrange(size.(blocks, 2)))
   )
 end
diff --git a/src/blocksparsearray/blocksparsearray.jl b/src/blocksparsearray/blocksparsearray.jl
index 1de153a..ab7fb22 100644
--- a/src/blocksparsearray/blocksparsearray.jl
+++ b/src/blocksparsearray/blocksparsearray.jl
@@ -3,24 +3,34 @@ using DerivableInterfaces: @interface
 using Dictionaries: Dictionary
 using SparseArraysBase: SparseArrayDOK
 
+function _BlockSparseArray end
+
 struct BlockSparseArray{
   T,
   N,
   A<:AbstractArray{T,N},
   Blocks<:AbstractArray{A,N},
-  Axes<:Tuple{Vararg{AbstractUnitRange,N}},
+  Axes<:Tuple{Vararg{AbstractUnitRange{<:Integer},N}},
 } <: AbstractBlockSparseArray{T,N}
   blocks::Blocks
   axes::Axes
+  global @inline function _BlockSparseArray(
+    blocks::AbstractArray{<:AbstractArray{T,N},N},
+    axes::Tuple{Vararg{AbstractUnitRange{<:Integer},N}},
+  ) where {T,N}
+    Base.require_one_based_indexing(axes...)
+    Base.require_one_based_indexing(blocks)
+    return new{T,N,eltype(blocks),typeof(blocks),typeof(axes)}(blocks, axes)
+  end
 end
 
 # TODO: Can this definition be shortened?
-const BlockSparseMatrix{T,A<:AbstractMatrix{T},Blocks<:AbstractMatrix{A},Axes<:Tuple{AbstractUnitRange,AbstractUnitRange}} = BlockSparseArray{
+const BlockSparseMatrix{T,A<:AbstractMatrix{T},Blocks<:AbstractMatrix{A},Axes<:Tuple{AbstractUnitRange{<:Integer},AbstractUnitRange{<:Integer}}} = BlockSparseArray{
   T,2,A,Blocks,Axes
 }
 
 # TODO: Can this definition be shortened?
-const BlockSparseVector{T,A<:AbstractVector{T},Blocks<:AbstractVector{A},Axes<:Tuple{AbstractUnitRange}} = BlockSparseArray{
+const BlockSparseVector{T,A<:AbstractVector{T},Blocks<:AbstractVector{A},Axes<:Tuple{AbstractUnitRange{<:Integer}}} = BlockSparseArray{
   T,1,A,Blocks,Axes
 }
 
@@ -31,39 +41,17 @@ Construct a block sparse array from a sparse array of arrays and specified block
 The block sizes must be commensurate with the blocks of the axes.
 """
 function sparsemortar(
-  blocks::AbstractArray{<:AbstractArray{T,N},N}, axes::Tuple{Vararg{AbstractUnitRange,N}}
+  blocks::AbstractArray{<:AbstractArray{T,N},N},
+  axes::Tuple{Vararg{AbstractUnitRange{<:Integer},N}},
 ) where {T,N}
-  return BlockSparseArray{T,N,eltype(blocks),typeof(blocks),typeof(axes)}(blocks, axes)
-end
-
-"""
-    sparsemortar(blocks::Dictionary{<:Block{N},<:AbstractArray{T,N}}, axes) -> ::BlockSparseArray{T,N}
-
-Construct a block sparse array from a dictionary mapping the locations of the stored/nonzero
-blocks to the data of those blocks, along with a set of blocked axes.
-The block sizes must be commensurate with the blocks of the specified axes.
-"""
-function sparsemortar(
-  block_data::Dictionary{<:Block{N},<:AbstractArray{<:Any,N}},
-  axes::Tuple{Vararg{AbstractUnitRange,N}},
-) where {N}
-  blocks = default_blocks(block_data, axes)
-  return sparsemortar(blocks, axes)
+  return _BlockSparseArray(blocks, axes)
 end
 
-"""
-    sparsemortar(block_indices::Vector{<:Block{N}}, block_data::Vector{<:AbstractArray{T,N}}, axes) -> ::BlockSparseArray{T,N}
-
-Construct a block sparse array from a list of locations of the stored/nonzero blocks,
-a corresponding list of the data of those blocks, along with a set of blocked axes.
-The block sizes must be commensurate with the blocks of the specified axes.
-"""
-function sparsemortar(
-  block_indices::Vector{<:Block{N}},
-  block_data::Vector{<:AbstractArray{<:Any,N}},
-  axes::Tuple{Vararg{AbstractUnitRange,N}},
-) where {N}
-  return sparsemortar(Dictionary(block_indices, block_data), axes)
+function BlockArrays.mortar(
+  blocks::SparseArrayDOK{<:AbstractArray{T,N},N},
+  axes::Tuple{Vararg{AbstractUnitRange{<:Integer},N}},
+) where {T,N}
+  return _BlockSparseArray(blocks, axes)
 end
 
 @doc """
@@ -79,7 +67,7 @@ function BlockSparseArray{T,N,A}(
   ::UndefInitializer, axes::Tuple{Vararg{AbstractUnitRange,N}}
 ) where {T,N,A<:AbstractArray{T,N}}
   blocks = default_blocks(A, axes)
-  return sparsemortar(blocks, axes)
+  return _BlockSparseArray(blocks, axes)
 end
 
 function BlockSparseArray{T,N,A}(
@@ -99,7 +87,7 @@ function BlockSparseArray{T,0,A}(
   ::UndefInitializer, axes::Tuple{}
 ) where {T,A<:AbstractArray{T,0}}
   blocks = default_blocks(A, axes)
-  return sparsemortar(blocks, axes)
+  return _BlockSparseArray(blocks, axes)
 end
 
 function BlockSparseArray{T,N,A}(
diff --git a/src/blocksparsearrayinterface/map.jl b/src/blocksparsearrayinterface/map.jl
index 6b5a446..6d892d3 100644
--- a/src/blocksparsearrayinterface/map.jl
+++ b/src/blocksparsearrayinterface/map.jl
@@ -96,21 +96,12 @@ end
 
 # TODO: Decide what to do with these.
 function map_stored_blocks(f, a::AbstractArray)
-  # TODO: Implement this as:
-  # ```julia
-  # mapped_blocks = SparseArraysInterface.map_stored(f, blocks(a))
-  # BlockSparseArray(mapped_blocks, axes(a))
-  # ```
-  # TODO: `block_stored_indices` should output `Indices` storing
-  # the stored Blocks, not a `Dictionary` from cartesian indices
-  # to Blocks.
   bs = collect(eachblockstoredindex(a))
   ds = map(b -> f(@view(a[b])), bs)
-  # We manually specify the block type using `Base.promote_op`
-  # since `a[b]` may not be inferrable. For example, if `blocktype(a)`
-  # is `Diagonal{Float64,Vector{Float64}}`, the non-stored blocks are `Matrix{Float64}`
-  # since they can't necessarily by `Diagonal` if there are rectangular blocks.
-  mapped_blocks = Dictionary{eltype(bs),eltype(ds)}(bs, ds)
-  # TODO: Use `similartype(typeof(a), eltype(eltype(mapped_blocks)))(...)`.
-  return sparsemortar(mapped_blocks, axes(a))
+  # TODO: Use `similartype` instead?
+  a′ = BlockSparseArray{eltype(eltype(ds)),ndims(a),eltype(ds)}(undef, axes(a))
+  for (b, d) in zip(bs, ds)
+    a′[b] = d
+  end
+  return a′
 end

From 307437dc39c9092ba5cc5dd509d61e9b6b2f6329 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Tue, 4 Mar 2025 15:33:19 -0500
Subject: [PATCH 07/14] undef_blocks constructor for SparseArrayDOK, more
 sparsemortar constructors

---
 src/BlockSparseArrays.jl                   |   1 -
 src/blocksparsearray/blockdiagonalarray.jl |   4 +-
 src/blocksparsearray/blocksparsearray.jl   | 120 ++++++++++++++++-----
 src/blocksparsearray/defaults.jl           |  42 --------
 test/test_basics.jl                        |  21 +++-
 5 files changed, 117 insertions(+), 71 deletions(-)
 delete mode 100644 src/blocksparsearray/defaults.jl

diff --git a/src/BlockSparseArrays.jl b/src/BlockSparseArrays.jl
index b9ac7b0..a350dff 100644
--- a/src/BlockSparseArrays.jl
+++ b/src/BlockSparseArrays.jl
@@ -39,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")
 
diff --git a/src/blocksparsearray/blockdiagonalarray.jl b/src/blocksparsearray/blockdiagonalarray.jl
index 0926431..ce43cee 100644
--- a/src/blocksparsearray/blockdiagonalarray.jl
+++ b/src/blocksparsearray/blockdiagonalarray.jl
@@ -12,9 +12,7 @@ const BlockSparseDiagonal{T,A<:AbstractBlockSparseVector{T}} = Diagonal{T,A}
 end
 
 function BlockDiagonal(blocks::AbstractVector{<:AbstractMatrix})
-  return sparsemortar(
-    Diagonal(blocks), (blockedrange(size.(blocks, 1)), blockedrange(size.(blocks, 2)))
-  )
+  return sparsemortar(Diagonal(blocks), size.(blocks, 1), size.(blocks, 2))
 end
 
 function DiagonalArrays.diagonal(S::BlockSparseVector)
diff --git a/src/blocksparsearray/blocksparsearray.jl b/src/blocksparsearray/blocksparsearray.jl
index ab7fb22..e3cc65a 100644
--- a/src/blocksparsearray/blocksparsearray.jl
+++ b/src/blocksparsearray/blocksparsearray.jl
@@ -1,8 +1,64 @@
-using BlockArrays: BlockArrays, Block, BlockedUnitRange, blockedrange, blocklength
+using BlockArrays:
+  BlockArrays,
+  Block,
+  BlockedUnitRange,
+  UndefBlocksInitializer,
+  blockedrange,
+  blocklength,
+  undef_blocks
 using DerivableInterfaces: @interface
 using Dictionaries: Dictionary
 using SparseArraysBase: SparseArrayDOK
 
+"""
+    SparseArrayDOK{T}(::UndefBlocksInitializer, axes)
+    SparseArrayDOK{T,N}(::UndefBlocksInitializer, axes)
+
+Construct the block structure of an undefined BlockSparseArray that will have
+blocked axes `axes`.
+"""
+function SparseArraysBase.SparseArrayDOK{T,N}(
+  ::UndefBlocksInitializer, ax::Tuple{Vararg{AbstractUnitRange{<:Integer},N}}
+) where {T,N}
+  return SparseArrayDOK{T,N}(undef, blocklength.(ax), GetUnstoredBlock(ax))
+end
+function SparseArraysBase.SparseArrayDOK{T,N}(
+  ::UndefBlocksInitializer, ax::Vararg{AbstractUnitRange{<:Integer},N}
+) where {T,N}
+  return SparseArrayDOK{T,N}(undef_blocks, ax)
+end
+function SparseArraysBase.SparseArrayDOK{T,N}(
+  ::UndefBlocksInitializer, dims::Tuple{Vararg{AbstractVector{<:Integer},N}}
+) where {T,N}
+  return SparseArrayDOK{T,N}(undef_blocks, blockedrange.(dims))
+end
+function SparseArraysBase.SparseArrayDOK{T,N}(
+  ::UndefBlocksInitializer, dims::Vararg{AbstractVector{<:Integer},N}
+) where {T,N}
+  return SparseArrayDOK{T,N}(undef_blocks, blockedrange.(dims))
+end
+
+function SparseArraysBase.SparseArrayDOK{T}(
+  ::UndefBlocksInitializer, ax::Tuple{Vararg{AbstractUnitRange{<:Integer},N}}
+) where {T,N}
+  return SparseArrayDOK{T,N}(undef_blocks, ax)
+end
+function SparseArraysBase.SparseArrayDOK{T}(
+  ::UndefBlocksInitializer, ax::Vararg{AbstractUnitRange{<:Integer},N}
+) where {T,N}
+  return SparseArrayDOK{T,N}(undef_blocks, ax)
+end
+function SparseArraysBase.SparseArrayDOK{T}(
+  ::UndefBlocksInitializer, dims::Tuple{Vararg{AbstractVector{<:Integer},N}}
+) where {T,N}
+  return SparseArrayDOK{T,N}(undef_blocks, blockedrange.(dims))
+end
+function SparseArraysBase.SparseArrayDOK{T}(
+  ::UndefBlocksInitializer, dims::Vararg{AbstractVector{<:Integer},N}
+) where {T,N}
+  return SparseArrayDOK{T,N}(undef_blocks, blockedrange.(dims))
+end
+
 function _BlockSparseArray end
 
 struct BlockSparseArray{
@@ -46,12 +102,22 @@ function sparsemortar(
 ) where {T,N}
   return _BlockSparseArray(blocks, axes)
 end
-
-function BlockArrays.mortar(
-  blocks::SparseArrayDOK{<:AbstractArray{T,N},N},
-  axes::Tuple{Vararg{AbstractUnitRange{<:Integer},N}},
+function sparsemortar(
+  blocks::AbstractArray{<:AbstractArray{T,N},N},
+  axes::Vararg{AbstractUnitRange{<:Integer},N},
 ) where {T,N}
-  return _BlockSparseArray(blocks, axes)
+  return sparsemortar(blocks, axes)
+end
+function sparsemortar(
+  blocks::AbstractArray{<:AbstractArray{T,N},N},
+  dims::Tuple{Vararg{AbstractVector{<:Integer},N}},
+) where {T,N}
+  return sparsemortar(blocks, blockedrange.(dims))
+end
+function sparsemortar(
+  blocks::AbstractArray{<:AbstractArray{T,N},N}, dims::Vararg{AbstractVector{<:Integer},N}
+) where {T,N}
+  return sparsemortar(blocks, dims)
 end
 
 @doc """
@@ -64,20 +130,19 @@ of block lengths in each dimension or a list of blocked ranges representing the
 """ BlockSparseArray
 
 function BlockSparseArray{T,N,A}(
-  ::UndefInitializer, axes::Tuple{Vararg{AbstractUnitRange,N}}
+  ::UndefInitializer, axes::Tuple{Vararg{AbstractUnitRange{<:Integer},N}}
 ) where {T,N,A<:AbstractArray{T,N}}
-  blocks = default_blocks(A, axes)
-  return _BlockSparseArray(blocks, axes)
+  return _BlockSparseArray(SparseArrayDOK{A}(undef_blocks, axes), axes)
 end
 
 function BlockSparseArray{T,N,A}(
-  ::UndefInitializer, axes::Vararg{AbstractUnitRange,N}
+  ::UndefInitializer, axes::Vararg{AbstractUnitRange{<:Integer},N}
 ) where {T,N,A<:AbstractArray{T,N}}
   return BlockSparseArray{T,N,A}(undef, axes)
 end
 
 function BlockSparseArray{T,N,A}(
-  ::UndefInitializer, dims::Tuple{Vararg{Vector{Int},N}}
+  ::UndefInitializer, dims::Tuple{Vararg{AbstractVector{<:Integer},N}}
 ) where {T,N,A<:AbstractArray{T,N}}
   return BlockSparseArray{T,N,A}(undef, blockedrange.(dims))
 end
@@ -86,48 +151,51 @@ end
 function BlockSparseArray{T,0,A}(
   ::UndefInitializer, axes::Tuple{}
 ) where {T,A<:AbstractArray{T,0}}
-  blocks = default_blocks(A, axes)
-  return _BlockSparseArray(blocks, axes)
+  return _BlockSparseArray(SparseArrayDOK{A}(undef_blocks, axes), axes)
 end
 
 function BlockSparseArray{T,N,A}(
-  ::UndefInitializer, dims::Vararg{Vector{Int},N}
+  ::UndefInitializer, dims::Vararg{AbstractVector{<:Integer},N}
 ) where {T,N,A<:AbstractArray{T,N}}
   return BlockSparseArray{T,N,A}(undef, dims)
 end
 
 function BlockSparseArray{T,N}(
-  ::UndefInitializer, axes::Tuple{Vararg{AbstractUnitRange,N}}
+  ::UndefInitializer, axes::Tuple{Vararg{AbstractUnitRange{<:Integer},N}}
 ) where {T,N}
-  return BlockSparseArray{T,N,default_arraytype(T, axes)}(undef, axes)
+  return BlockSparseArray{T,N,Array{T,N}}(undef, axes)
 end
 
 function BlockSparseArray{T,N}(
-  ::UndefInitializer, axes::Vararg{AbstractUnitRange,N}
+  ::UndefInitializer, axes::Vararg{AbstractUnitRange{<:Integer},N}
 ) where {T,N}
   return BlockSparseArray{T,N}(undef, axes)
 end
 
 function BlockSparseArray{T,0}(::UndefInitializer, axes::Tuple{}) where {T}
-  return BlockSparseArray{T,0,default_arraytype(T, axes)}(undef, axes)
+  return BlockSparseArray{T,0,Array{T,0}}(undef, axes)
 end
 
 function BlockSparseArray{T,N}(
-  ::UndefInitializer, dims::Tuple{Vararg{Vector{Int},N}}
+  ::UndefInitializer, dims::Tuple{Vararg{AbstractVector{<:Integer},N}}
 ) where {T,N}
   return BlockSparseArray{T,N}(undef, blockedrange.(dims))
 end
 
-function BlockSparseArray{T,N}(::UndefInitializer, dims::Vararg{Vector{Int},N}) where {T,N}
+function BlockSparseArray{T,N}(
+  ::UndefInitializer, dims::Vararg{AbstractVector{<:Integer},N}
+) where {T,N}
   return BlockSparseArray{T,N}(undef, dims)
 end
 
-function BlockSparseArray{T}(::UndefInitializer, dims::Tuple{Vararg{Vector{Int}}}) where {T}
+function BlockSparseArray{T}(
+  ::UndefInitializer, dims::Tuple{Vararg{AbstractVector{<:Integer}}}
+) where {T}
   return BlockSparseArray{T,length(dims)}(undef, dims)
 end
 
 function BlockSparseArray{T}(
-  ::UndefInitializer, axes::Tuple{Vararg{AbstractUnitRange}}
+  ::UndefInitializer, axes::Tuple{Vararg{AbstractUnitRange{<:Integer}}}
 ) where {T}
   return BlockSparseArray{T,length(axes)}(undef, axes)
 end
@@ -136,11 +204,15 @@ function BlockSparseArray{T}(::UndefInitializer, axes::Tuple{}) where {T}
   return BlockSparseArray{T,length(axes)}(undef, axes)
 end
 
-function BlockSparseArray{T}(::UndefInitializer, dims::Vararg{Vector{Int}}) where {T}
+function BlockSparseArray{T}(
+  ::UndefInitializer, dims::Vararg{AbstractVector{<:Integer}}
+) where {T}
   return BlockSparseArray{T}(undef, dims)
 end
 
-function BlockSparseArray{T}(::UndefInitializer, axes::Vararg{AbstractUnitRange}) where {T}
+function BlockSparseArray{T}(
+  ::UndefInitializer, axes::Vararg{AbstractUnitRange{<:Integer}}
+) where {T}
   return BlockSparseArray{T}(undef, axes)
 end
 
diff --git a/src/blocksparsearray/defaults.jl b/src/blocksparsearray/defaults.jl
deleted file mode 100644
index 788a975..0000000
--- a/src/blocksparsearray/defaults.jl
+++ /dev/null
@@ -1,42 +0,0 @@
-using BlockArrays: Block
-using Dictionaries: Dictionary
-using SparseArraysBase: SparseArrayDOK
-
-# Construct the sparse structure storing the blocks
-function default_blockdata(
-  block_data::Dictionary{<:CartesianIndex{N},<:AbstractArray{<:Any,N}},
-  axes::Tuple{Vararg{AbstractUnitRange,N}},
-) where {N}
-  return error()
-end
-
-function default_blocks(
-  block_indices::Vector{<:Block{N}},
-  block_data::Vector{<:AbstractArray{<:Any,N}},
-  axes::Tuple{Vararg{AbstractUnitRange,N}},
-) where {N}
-  return default_blocks(Dictionary(block_indices, block_data), axes)
-end
-
-function default_blocks(
-  block_data::Dictionary{<:Block{N},<:AbstractArray{<:Any,N}},
-  axes::Tuple{Vararg{AbstractUnitRange,N}},
-) where {N}
-  return default_blocks(blocks_to_cartesianindices(block_data), axes)
-end
-
-function default_arraytype(elt::Type, axes::Tuple{Vararg{AbstractUnitRange}})
-  return Array{elt,length(axes)}
-end
-
-function default_blocks(blocktype::Type, axes::Tuple{Vararg{AbstractUnitRange}})
-  block_data = Dictionary{Block{length(axes),Int},blocktype}()
-  return default_blocks(block_data, axes)
-end
-
-function default_blocks(
-  block_data::Dictionary{<:CartesianIndex{N},<:AbstractArray{<:Any,N}},
-  axes::Tuple{Vararg{AbstractUnitRange,N}},
-) where {N}
-  return SparseArrayDOK(block_data, blocklength.(axes), GetUnstoredBlock(axes))
-end
diff --git a/test/test_basics.jl b/test/test_basics.jl
index a52f1dc..e39e266 100644
--- a/test/test_basics.jl
+++ b/test/test_basics.jl
@@ -17,7 +17,8 @@ using BlockArrays:
   blocklengths,
   blocksize,
   blocksizes,
-  mortar
+  mortar,
+  undef_blocks
 using BlockSparseArrays:
   @view!,
   BlockSparseArray,
@@ -158,6 +159,24 @@ arrayts = (Array, JLArray)
     @test (@constinferred blocktype(a)) <: SubArray{elt,2,arrayt{elt,2}}
     # TODO: This is difficult to determine just from type information.
     @test_broken blocktype(typeof(a)) <: SubArray{elt,2,arrayt{elt,2}}
+
+    # sparsemortar
+    for ax in (
+      ([2, 3], [2, 3]),
+      (([2, 3], [2, 3]),),
+      blockedrange.(([2, 3], [2, 3])),
+      (blockedrange.(([2, 3], [2, 3])),),
+    )
+      blocks = SparseArrayDOK{arrayt{elt,2}}(undef_blocks, ax...)
+      blocks[2, 1] = arrayt(randn(elt, 3, 2))
+      blocks[1, 2] = arrayt(randn(elt, 2, 3))
+      a = sparsemortar(blocks, ax...)
+      @test a isa BlockSparseArray{elt,2,arrayt{elt,2}}
+      @test iszero(a[Block(1, 1)])
+      @test a[Block(2, 1)] == blocks[2, 1]
+      @test a[Block(1, 2)] == blocks[1, 2]
+      @test iszero(a[Block(2, 2)])
+    end
   end
   @testset "Basics" begin
     a = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))

From 96d3ac27485d9f1083e3deca7338ac6ce37c33e3 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Tue, 4 Mar 2025 16:08:10 -0500
Subject: [PATCH 08/14] Rewrite map_stored_blocks

---
 src/blocksparsearrayinterface/map.jl | 21 +++++++++++++--------
 1 file changed, 13 insertions(+), 8 deletions(-)

diff --git a/src/blocksparsearrayinterface/map.jl b/src/blocksparsearrayinterface/map.jl
index 6d892d3..8664279 100644
--- a/src/blocksparsearrayinterface/map.jl
+++ b/src/blocksparsearrayinterface/map.jl
@@ -1,5 +1,7 @@
+using BlockArrays: blocks, eachstoredindex, undef_blocks
 using DerivableInterfaces: @interface, AbstractArrayInterface, interface
 using GPUArraysCore: @allowscalar
+using SparseArraysBase: SparseArrayDOK
 
 # TODO: Rewrite this so that it takes the blocking structure
 # made by combining the blocking of the axes (i.e. the blocking that
@@ -94,14 +96,17 @@ function map_zero_dim! end
   return a_dest
 end
 
-# TODO: Decide what to do with these.
+# TODO: Do we need this function or can we just use `map`?
+# Probably it should be a special version of `map` where we
+# specify the function preserves zeros, i.e.
+# `map(f, a; preserves_zeros=true)` or `@preserves_zeros map(f, a)`.
 function map_stored_blocks(f, a::AbstractArray)
-  bs = collect(eachblockstoredindex(a))
-  ds = map(b -> f(@view(a[b])), bs)
-  # TODO: Use `similartype` instead?
-  a′ = BlockSparseArray{eltype(eltype(ds)),ndims(a),eltype(ds)}(undef, axes(a))
-  for (b, d) in zip(bs, ds)
-    a′[b] = d
+  blocks_a = blocks(a)
+  stored_indices = collect(eachstoredindex(a))
+  stored_blocks = map(I -> f(blocks_a[I]), stored_indices)
+  blocks_a′ = SparseArrayDOK{eltype(stored_blocks)}(undef_blocks, axes(a))
+  for (I, b) in zip(stored_indices, stored_blocks)
+    blocks_a′[I] = b
   end
-  return a′
+  return sparsemortar(blocks_a′, axes(a))
 end

From 47a5eb5bf5f259adf8e8f4e6554125f81b1fdc86 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Tue, 4 Mar 2025 16:19:45 -0500
Subject: [PATCH 09/14] Fix tests

---
 src/blocksparsearrayinterface/map.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/blocksparsearrayinterface/map.jl b/src/blocksparsearrayinterface/map.jl
index 8664279..de67adc 100644
--- a/src/blocksparsearrayinterface/map.jl
+++ b/src/blocksparsearrayinterface/map.jl
@@ -102,7 +102,7 @@ end
 # `map(f, a; preserves_zeros=true)` or `@preserves_zeros map(f, a)`.
 function map_stored_blocks(f, a::AbstractArray)
   blocks_a = blocks(a)
-  stored_indices = collect(eachstoredindex(a))
+  stored_indices = collect(eachstoredindex(blocks_a))
   stored_blocks = map(I -> f(blocks_a[I]), stored_indices)
   blocks_a′ = SparseArrayDOK{eltype(stored_blocks)}(undef_blocks, axes(a))
   for (I, b) in zip(stored_indices, stored_blocks)

From 760cf9cc53f773fbfa9c1f2401ae35466fbe82b7 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Tue, 4 Mar 2025 16:21:13 -0500
Subject: [PATCH 10/14] Fix tests

---
 src/blocksparsearrayinterface/map.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/blocksparsearrayinterface/map.jl b/src/blocksparsearrayinterface/map.jl
index de67adc..715689a 100644
--- a/src/blocksparsearrayinterface/map.jl
+++ b/src/blocksparsearrayinterface/map.jl
@@ -1,7 +1,7 @@
-using BlockArrays: blocks, eachstoredindex, undef_blocks
+using BlockArrays: blocks, undef_blocks
 using DerivableInterfaces: @interface, AbstractArrayInterface, interface
 using GPUArraysCore: @allowscalar
-using SparseArraysBase: SparseArrayDOK
+using SparseArraysBase: SparseArrayDOK, eachstoredindex
 
 # TODO: Rewrite this so that it takes the blocking structure
 # made by combining the blocking of the axes (i.e. the blocking that

From b8aa0cdc531c3a66c44a805c6f3124651f0f9805 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Tue, 4 Mar 2025 22:44:44 -0500
Subject: [PATCH 11/14] Various improvements

---
 src/blocksparsearray/blockdiagonalarray.jl |  2 +-
 src/blocksparsearray/blocksparsearray.jl   | 84 +++++++++++-----------
 src/blocksparsearrayinterface/map.jl       | 20 +++---
 3 files changed, 54 insertions(+), 52 deletions(-)

diff --git a/src/blocksparsearray/blockdiagonalarray.jl b/src/blocksparsearray/blockdiagonalarray.jl
index ce43cee..ebec0c1 100644
--- a/src/blocksparsearray/blockdiagonalarray.jl
+++ b/src/blocksparsearray/blockdiagonalarray.jl
@@ -12,7 +12,7 @@ const BlockSparseDiagonal{T,A<:AbstractBlockSparseVector{T}} = Diagonal{T,A}
 end
 
 function BlockDiagonal(blocks::AbstractVector{<:AbstractMatrix})
-  return sparsemortar(Diagonal(blocks), size.(blocks, 1), size.(blocks, 2))
+  return _BlockSparseArray(Diagonal(blocks), size.(blocks, 1), size.(blocks, 2))
 end
 
 function DiagonalArrays.diagonal(S::BlockSparseVector)
diff --git a/src/blocksparsearray/blocksparsearray.jl b/src/blocksparsearray/blocksparsearray.jl
index e3cc65a..aa280c0 100644
--- a/src/blocksparsearray/blocksparsearray.jl
+++ b/src/blocksparsearray/blocksparsearray.jl
@@ -11,11 +11,13 @@ using Dictionaries: Dictionary
 using SparseArraysBase: SparseArrayDOK
 
 """
-    SparseArrayDOK{T}(::UndefBlocksInitializer, axes)
-    SparseArrayDOK{T,N}(::UndefBlocksInitializer, axes)
+    SparseArrayDOK{T}(undef_blocks, axes)
+    SparseArrayDOK{T,N}(undef_blocks, axes)
 
 Construct the block structure of an undefined BlockSparseArray that will have
 blocked axes `axes`.
+
+See also: [`undef_blocks`](@ref BlockArrays.undef_blocks)`
 """
 function SparseArraysBase.SparseArrayDOK{T,N}(
   ::UndefBlocksInitializer, ax::Tuple{Vararg{AbstractUnitRange{<:Integer},N}}
@@ -28,14 +30,17 @@ function SparseArraysBase.SparseArrayDOK{T,N}(
   return SparseArrayDOK{T,N}(undef_blocks, ax)
 end
 function SparseArraysBase.SparseArrayDOK{T,N}(
-  ::UndefBlocksInitializer, dims::Tuple{Vararg{AbstractVector{<:Integer},N}}
+  ::UndefBlocksInitializer,
+  dims::Tuple{AbstractVector{<:Integer},Vararg{AbstractVector{<:Integer}}},
 ) where {T,N}
   return SparseArrayDOK{T,N}(undef_blocks, blockedrange.(dims))
 end
 function SparseArraysBase.SparseArrayDOK{T,N}(
-  ::UndefBlocksInitializer, dims::Vararg{AbstractVector{<:Integer},N}
+  ::UndefBlocksInitializer,
+  dim1::AbstractVector{<:Integer},
+  dim_rest::AbstractVector{<:Integer}...,
 ) where {T,N}
-  return SparseArrayDOK{T,N}(undef_blocks, blockedrange.(dims))
+  return SparseArrayDOK{T,N}(undef_blocks, (dim1, dim_rest...))
 end
 
 function SparseArraysBase.SparseArrayDOK{T}(
@@ -49,14 +54,17 @@ function SparseArraysBase.SparseArrayDOK{T}(
   return SparseArrayDOK{T,N}(undef_blocks, ax)
 end
 function SparseArraysBase.SparseArrayDOK{T}(
-  ::UndefBlocksInitializer, dims::Tuple{Vararg{AbstractVector{<:Integer},N}}
-) where {T,N}
-  return SparseArrayDOK{T,N}(undef_blocks, blockedrange.(dims))
+  ::UndefBlocksInitializer,
+  dims::Tuple{AbstractVector{<:Integer},Vararg{AbstractVector{<:Integer}}},
+) where {T}
+  return SparseArrayDOK{T}(undef_blocks, blockedrange.(dims))
 end
 function SparseArraysBase.SparseArrayDOK{T}(
-  ::UndefBlocksInitializer, dims::Vararg{AbstractVector{<:Integer},N}
-) where {T,N}
-  return SparseArrayDOK{T,N}(undef_blocks, blockedrange.(dims))
+  ::UndefBlocksInitializer,
+  dim1::AbstractVector{<:Integer},
+  dim_rest::AbstractVector{<:Integer}...,
+) where {T}
+  return SparseArrayDOK{T}(undef_blocks, (dim1, dim_rest...))
 end
 
 function _BlockSparseArray end
@@ -110,14 +118,16 @@ function sparsemortar(
 end
 function sparsemortar(
   blocks::AbstractArray{<:AbstractArray{T,N},N},
-  dims::Tuple{Vararg{AbstractVector{<:Integer},N}},
+  dims::Tuple{AbstractVector{<:Integer},Vararg{AbstractVector{<:Integer}}},
 ) where {T,N}
   return sparsemortar(blocks, blockedrange.(dims))
 end
 function sparsemortar(
-  blocks::AbstractArray{<:AbstractArray{T,N},N}, dims::Vararg{AbstractVector{<:Integer},N}
+  blocks::AbstractArray{<:AbstractArray{T,N},N},
+  dim1::AbstractVector{<:Integer},
+  dim_rest::AbstractVector{<:Integer}...,
 ) where {T,N}
-  return sparsemortar(blocks, dims)
+  return sparsemortar(blocks, (dim1, dim_rest...))
 end
 
 @doc """
@@ -142,22 +152,18 @@ function BlockSparseArray{T,N,A}(
 end
 
 function BlockSparseArray{T,N,A}(
-  ::UndefInitializer, dims::Tuple{Vararg{AbstractVector{<:Integer},N}}
+  ::UndefInitializer,
+  dims::Tuple{AbstractVector{<:Integer},Vararg{AbstractVector{<:Integer}}},
 ) where {T,N,A<:AbstractArray{T,N}}
   return BlockSparseArray{T,N,A}(undef, blockedrange.(dims))
 end
 
-# Fix ambiguity error.
-function BlockSparseArray{T,0,A}(
-  ::UndefInitializer, axes::Tuple{}
-) where {T,A<:AbstractArray{T,0}}
-  return _BlockSparseArray(SparseArrayDOK{A}(undef_blocks, axes), axes)
-end
-
 function BlockSparseArray{T,N,A}(
-  ::UndefInitializer, dims::Vararg{AbstractVector{<:Integer},N}
+  ::UndefInitializer,
+  dim1::AbstractVector{<:Integer},
+  dim_rest::AbstractVector{<:Integer}...,
 ) where {T,N,A<:AbstractArray{T,N}}
-  return BlockSparseArray{T,N,A}(undef, dims)
+  return BlockSparseArray{T,N,A}(undef, (dim1, dim_rest...))
 end
 
 function BlockSparseArray{T,N}(
@@ -172,24 +178,24 @@ function BlockSparseArray{T,N}(
   return BlockSparseArray{T,N}(undef, axes)
 end
 
-function BlockSparseArray{T,0}(::UndefInitializer, axes::Tuple{}) where {T}
-  return BlockSparseArray{T,0,Array{T,0}}(undef, axes)
-end
-
 function BlockSparseArray{T,N}(
-  ::UndefInitializer, dims::Tuple{Vararg{AbstractVector{<:Integer},N}}
+  ::UndefInitializer,
+  dims::Tuple{AbstractVector{<:Integer},Vararg{AbstractVector{<:Integer}}},
 ) where {T,N}
   return BlockSparseArray{T,N}(undef, blockedrange.(dims))
 end
 
 function BlockSparseArray{T,N}(
-  ::UndefInitializer, dims::Vararg{AbstractVector{<:Integer},N}
+  ::UndefInitializer,
+  dim1::AbstractVector{<:Integer},
+  dim_rest::AbstractVector{<:Integer}...,
 ) where {T,N}
-  return BlockSparseArray{T,N}(undef, dims)
+  return BlockSparseArray{T,N}(undef, (dim1, dim_rest...))
 end
 
 function BlockSparseArray{T}(
-  ::UndefInitializer, dims::Tuple{Vararg{AbstractVector{<:Integer}}}
+  ::UndefInitializer,
+  dims::Tuple{AbstractVector{<:Integer},Vararg{AbstractVector{<:Integer}}},
 ) where {T}
   return BlockSparseArray{T,length(dims)}(undef, dims)
 end
@@ -200,14 +206,12 @@ function BlockSparseArray{T}(
   return BlockSparseArray{T,length(axes)}(undef, axes)
 end
 
-function BlockSparseArray{T}(::UndefInitializer, axes::Tuple{}) where {T}
-  return BlockSparseArray{T,length(axes)}(undef, axes)
-end
-
 function BlockSparseArray{T}(
-  ::UndefInitializer, dims::Vararg{AbstractVector{<:Integer}}
+  ::UndefInitializer,
+  dim1::AbstractVector{<:Integer},
+  dim_rest::AbstractVector{<:Integer}...,
 ) where {T}
-  return BlockSparseArray{T}(undef, dims)
+  return BlockSparseArray{T}(undef, (dim1, dim_rest...))
 end
 
 function BlockSparseArray{T}(
@@ -216,10 +220,6 @@ function BlockSparseArray{T}(
   return BlockSparseArray{T}(undef, axes)
 end
 
-function BlockSparseArray{T}(::UndefInitializer) where {T}
-  return BlockSparseArray{T}(undef, ())
-end
-
 # Base `AbstractArray` interface
 Base.axes(a::BlockSparseArray) = a.axes
 
diff --git a/src/blocksparsearrayinterface/map.jl b/src/blocksparsearrayinterface/map.jl
index 715689a..2df7133 100644
--- a/src/blocksparsearrayinterface/map.jl
+++ b/src/blocksparsearrayinterface/map.jl
@@ -1,7 +1,5 @@
-using BlockArrays: blocks, undef_blocks
 using DerivableInterfaces: @interface, AbstractArrayInterface, interface
 using GPUArraysCore: @allowscalar
-using SparseArraysBase: SparseArrayDOK, eachstoredindex
 
 # TODO: Rewrite this so that it takes the blocking structure
 # made by combining the blocking of the axes (i.e. the blocking that
@@ -101,12 +99,16 @@ end
 # specify the function preserves zeros, i.e.
 # `map(f, a; preserves_zeros=true)` or `@preserves_zeros map(f, a)`.
 function map_stored_blocks(f, a::AbstractArray)
-  blocks_a = blocks(a)
-  stored_indices = collect(eachstoredindex(blocks_a))
-  stored_blocks = map(I -> f(blocks_a[I]), stored_indices)
-  blocks_a′ = SparseArrayDOK{eltype(stored_blocks)}(undef_blocks, axes(a))
-  for (I, b) in zip(stored_indices, stored_blocks)
-    blocks_a′[I] = b
+  block_stored_indices = collect(eachblockstoredindex(a))
+  if isempty(block_stored_indices)
+    blocktype_a′ = Base.promote_op(f, blocktype(a))
+    return BlockSparseArray{eltype(blocktype_a′),ndims(a),blocktype_a′}(undef, axes(a))
   end
-  return sparsemortar(blocks_a′, axes(a))
+  stored_blocks = map(B -> f(@view!(a[B])), block_stored_indices)
+  blocktype_a′ = eltype(stored_blocks)
+  a′ = BlockSparseArray{eltype(blocktype_a′),ndims(a),blocktype_a′}(undef, axes(a))
+  for (B, b) in zip(block_stored_indices, stored_blocks)
+    a′[B] = b
+  end
+  return a′
 end

From 5affd8da9ea26ed681259c97b2223475fc903067 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Tue, 4 Mar 2025 22:52:43 -0500
Subject: [PATCH 12/14] Update docs

---
 src/blocksparsearray/blocksparsearray.jl | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/blocksparsearray/blocksparsearray.jl b/src/blocksparsearray/blocksparsearray.jl
index aa280c0..9d7e570 100644
--- a/src/blocksparsearray/blocksparsearray.jl
+++ b/src/blocksparsearray/blocksparsearray.jl
@@ -17,7 +17,9 @@ using SparseArraysBase: SparseArrayDOK
 Construct the block structure of an undefined BlockSparseArray that will have
 blocked axes `axes`.
 
-See also: [`undef_blocks`](@ref BlockArrays.undef_blocks)`
+Note that `undef_blocks` is defined in
+[BlockArrays.jl](https://juliaarrays.github.io/BlockArrays.jl/stable/lib/public/#BlockArrays.undef_blocks)
+and should be imported from that package to use it as an input to this constructor.
 """
 function SparseArraysBase.SparseArrayDOK{T,N}(
   ::UndefBlocksInitializer, ax::Tuple{Vararg{AbstractUnitRange{<:Integer},N}}

From dd31ac38ba4f10f07a8e3e270d2563124280c85a Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Tue, 4 Mar 2025 23:30:00 -0500
Subject: [PATCH 13/14] Fix tests

---
 src/blocksparsearray/blockdiagonalarray.jl | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/blocksparsearray/blockdiagonalarray.jl b/src/blocksparsearray/blockdiagonalarray.jl
index ebec0c1..e8668fb 100644
--- a/src/blocksparsearray/blockdiagonalarray.jl
+++ b/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,7 +13,7 @@ const BlockSparseDiagonal{T,A<:AbstractBlockSparseVector{T}} = Diagonal{T,A}
 end
 
 function BlockDiagonal(blocks::AbstractVector{<:AbstractMatrix})
-  return _BlockSparseArray(Diagonal(blocks), size.(blocks, 1), size.(blocks, 2))
+  return _BlockSparseArray(Diagonal(blocks), blockedrange.((size.(blocks, 1), size.(blocks, 2))))
 end
 
 function DiagonalArrays.diagonal(S::BlockSparseVector)

From 74d02d070795c195f004f553ed725c0d8d6cb6c2 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Tue, 4 Mar 2025 23:32:31 -0500
Subject: [PATCH 14/14] Format

---
 src/blocksparsearray/blockdiagonalarray.jl | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/blocksparsearray/blockdiagonalarray.jl b/src/blocksparsearray/blockdiagonalarray.jl
index e8668fb..05345ac 100644
--- a/src/blocksparsearray/blockdiagonalarray.jl
+++ b/src/blocksparsearray/blockdiagonalarray.jl
@@ -13,7 +13,9 @@ const BlockSparseDiagonal{T,A<:AbstractBlockSparseVector{T}} = Diagonal{T,A}
 end
 
 function BlockDiagonal(blocks::AbstractVector{<:AbstractMatrix})
-  return _BlockSparseArray(Diagonal(blocks), blockedrange.((size.(blocks, 1), size.(blocks, 2))))
+  return _BlockSparseArray(
+    Diagonal(blocks), blockedrange.((size.(blocks, 1), size.(blocks, 2)))
+  )
 end
 
 function DiagonalArrays.diagonal(S::BlockSparseVector)