From d4365f0904be3481c29e7e60e871e97acae7a76f Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Tue, 29 Jul 2025 09:11:33 -0400
Subject: [PATCH 1/6] Generalize matricize

---
 Project.toml                                  |  2 +-
 .../BlockSparseArraysTensorAlgebraExt.jl      | 30 ++++++++++++++++---
 .../BlockArraysExtensions.jl                  | 12 ++++----
 .../abstractblocksparsearray.jl               |  3 +-
 src/blocksparsearray/blocksparsearray.jl      | 10 +++++--
 .../getunstoredblock.jl                       |  2 +-
 src/factorizations/svd.jl                     | 12 ++++----
 7 files changed, 52 insertions(+), 19 deletions(-)

diff --git a/Project.toml b/Project.toml
index 7f0208a3..db5a78bd 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.8.1"
+version = "0.8.2"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
diff --git a/ext/BlockSparseArraysTensorAlgebraExt/BlockSparseArraysTensorAlgebraExt.jl b/ext/BlockSparseArraysTensorAlgebraExt/BlockSparseArraysTensorAlgebraExt.jl
index d5551f43..fc4d5940 100644
--- a/ext/BlockSparseArraysTensorAlgebraExt/BlockSparseArraysTensorAlgebraExt.jl
+++ b/ext/BlockSparseArraysTensorAlgebraExt/BlockSparseArraysTensorAlgebraExt.jl
@@ -15,19 +15,41 @@ function TensorAlgebra.FusionStyle(::Type{<:AbstractBlockSparseArray})
   return BlockReshapeFusion()
 end
 
+using BlockArrays: Block, blocklength, blocks
+using BlockSparseArrays: blocksparse
+using SparseArraysBase: eachstoredindex
+using TensorAlgebra: TensorAlgebra, matricize, unmatricize
 function TensorAlgebra.matricize(
   ::BlockReshapeFusion, a::AbstractArray, biperm::BlockedTrivialPermutation{2}
 )
-  new_axes = fuseaxes(axes(a), biperm)
-  return blockreshape(a, new_axes)
+  ax = fuseaxes(axes(a), biperm)
+  reshaped_blocks_a = reshape(blocks(a), map(blocklength, ax))
+  bs = Dict(
+    Block(Tuple(I)) => matricize(reshaped_blocks_a[I], biperm) for
+    I in eachstoredindex(reshaped_blocks_a)
+  )
+  return blocksparse(bs, ax)
 end
 
+using BlockArrays: blocklengths
 function TensorAlgebra.unmatricize(
   ::BlockReshapeFusion,
   m::AbstractMatrix,
-  blocked_axes::BlockedTuple{2,<:Any,<:Tuple{Vararg{AbstractUnitRange}}},
+  blocked_ax::BlockedTuple{2,<:Any,<:Tuple{Vararg{AbstractUnitRange}}},
 )
-  return blockreshape(m, Tuple(blocked_axes)...)
+  ax = Tuple(blocked_ax)
+  reshaped_blocks_m = reshape(blocks(m), map(blocklength, ax))
+  function f(I)
+    block_axes_I = BlockedTuple(
+      map(ntuple(identity, length(ax))) do i
+        return Base.axes1(ax[i][Block(I[i])])
+      end,
+      blocklengths(blocked_ax),
+    )
+    return unmatricize(reshaped_blocks_m[I], block_axes_I)
+  end
+  bs = Dict(Block(Tuple(I)) => f(I) for I in eachstoredindex(reshaped_blocks_m))
+  return blocksparse(bs, ax)
 end
 
 end
diff --git a/src/BlockArraysExtensions/BlockArraysExtensions.jl b/src/BlockArraysExtensions/BlockArraysExtensions.jl
index 58155403..8b2a1eff 100644
--- a/src/BlockArraysExtensions/BlockArraysExtensions.jl
+++ b/src/BlockArraysExtensions/BlockArraysExtensions.jl
@@ -267,13 +267,15 @@ tuple_oneto(n) = ntuple(identity, n)
 
 function _blockreshape(a::AbstractArray, axes::Tuple{Vararg{AbstractUnitRange}})
   reshaped_blocks_a = reshape(blocks(a), blocklength.(axes))
-  reshaped_a = similar(a, axes)
-  for I in eachstoredindex(reshaped_blocks_a)
-    block_size_I = map(i -> length(axes[i][Block(I[i])]), tuple_oneto(length(axes)))
+  function f(I)
+    block_axes_I = map(ntuple(identity, length(axes))) do i
+      return Base.axes1(axes[i][Block(I[i])])
+    end
     # TODO: Better converter here.
-    reshaped_a[Block(Tuple(I))] = reshape(reshaped_blocks_a[I], block_size_I)
+    return reshape(reshaped_blocks_a[I], block_axes_I)
   end
-  return reshaped_a
+  bs = Dict(Block(Tuple(I)) => f(I) for I in eachstoredindex(reshaped_blocks_a))
+  return blocksparse(bs, axes)
 end
 
 function blockreshape(
diff --git a/src/abstractblocksparsearray/abstractblocksparsearray.jl b/src/abstractblocksparsearray/abstractblocksparsearray.jl
index 3d688e45..1824e830 100644
--- a/src/abstractblocksparsearray/abstractblocksparsearray.jl
+++ b/src/abstractblocksparsearray/abstractblocksparsearray.jl
@@ -92,7 +92,8 @@ function Base.setindex!(
   # in Julia v1.10 (https://github.com/JuliaLang/julia/pull/48895,
   # https://github.com/JuliaLang/julia/pull/52487).
   # TODO: Delete once we drop support for Julia v1.10.
-  blocks(a)[Int.(I)...] = _convert(blocktype(a), value)
+  aI = @view! a[I...]
+  copyto!(aI, value)
   return a
 end
 
diff --git a/src/blocksparsearray/blocksparsearray.jl b/src/blocksparsearray/blocksparsearray.jl
index 15af2f77..030178ab 100644
--- a/src/blocksparsearray/blocksparsearray.jl
+++ b/src/blocksparsearray/blocksparsearray.jl
@@ -262,13 +262,19 @@ function blocksparsezeros(::BlockType{A}, axes...) where {A<:AbstractArray}
   # to make a bit more generic.
   return BlockSparseArray{eltype(A),ndims(A),A}(undef, axes...)
 end
-function blocksparse(d::Dict{<:Block,<:AbstractArray}, axes...)
-  a = blocksparsezeros(BlockType(valtype(d)), axes...)
+function blocksparse(d::Dict{<:Block,<:AbstractArray}, ax::Tuple)
+  blockaxtype = Tuple{map(eltype ∘ eachblockaxis, ax)...}
+  # TODO: Catch if inference fails and use `valtype(d)` instead.
+  blockt = Base.promote_op(similar, Type{valtype(d)}, blockaxtype)
+  a = blocksparsezeros(BlockType(blockt), ax)
   for I in eachindex(d)
     a[I] = d[I]
   end
   return a
 end
+function blocksparse(d::Dict{<:Block,<:AbstractArray}, ax::AbstractUnitRange...)
+  return blocksparse(d, ax)
+end
 
 # Base `AbstractArray` interface
 Base.axes(a::BlockSparseArray) = a.axes
diff --git a/src/blocksparsearrayinterface/getunstoredblock.jl b/src/blocksparsearrayinterface/getunstoredblock.jl
index d704f05a..c5ae1dca 100644
--- a/src/blocksparsearrayinterface/getunstoredblock.jl
+++ b/src/blocksparsearrayinterface/getunstoredblock.jl
@@ -20,7 +20,7 @@ end
 @inline function Base.getindex(a::ZeroBlocks{N,A}, I::Vararg{Int,N}) where {N,A}
   # TODO: Use `BlockArrays.eachblockaxes`.
   ax = ntuple(N) do d
-    return only(axes(a.parentaxes[d][Block(I[d])]))
+    return eachblockaxis(a.parentaxes[d])[I[d]]
   end
   !isconcretetype(A) && return zero!(similar(Array{eltype(A),N}, ax))
   return zero!(similar(A, ax))
diff --git a/src/factorizations/svd.jl b/src/factorizations/svd.jl
index 8c616d49..d03d16b9 100644
--- a/src/factorizations/svd.jl
+++ b/src/factorizations/svd.jl
@@ -49,12 +49,12 @@ function MatrixAlgebraKit.initialize_output(
   ::typeof(svd_compact!), A::AbstractBlockSparseMatrix, alg::BlockPermutedDiagonalAlgorithm
 )
   bm, bn = blocksize(A)
-  bmn = min(bm, bn)
+  (bmn, mindim) = findmin((bm, bn))
 
   brows = eachblockaxis(axes(A, 1))
   bcols = eachblockaxis(axes(A, 2))
   u_axes = similar(brows, bmn)
-  v_axes = similar(brows, bmn)
+  v_axes = similar(bcols, bmn)
 
   # fill in values for blocks that are present
   bIs = collect(eachblockstoredindex(A))
@@ -62,8 +62,9 @@ function MatrixAlgebraKit.initialize_output(
   bcolIs = Int.(last.(Tuple.(bIs)))
   for bI in eachblockstoredindex(A)
     row, col = Int.(Tuple(bI))
-    u_axes[col] = infimum(brows[row], bcols[col])
-    v_axes[col] = infimum(bcols[col], brows[row])
+    dim = (row, col)[mindim]
+    u_axes[dim] = infimum(brows[row], bcols[col])
+    v_axes[dim] = infimum(bcols[col], brows[row])
   end
 
   # fill in values for blocks that aren't present, pairing them in order of occurence
@@ -83,9 +84,10 @@ function MatrixAlgebraKit.initialize_output(
   # allocate output
   for bI in eachblockstoredindex(A)
     brow, bcol = Tuple(bI)
+    bdim = (brow, bcol)[mindim]
     block = @view!(A[bI])
     block_alg = block_algorithm(alg, block)
-    U[brow, bcol], S[bcol, bcol], Vt[bcol, bcol] = MatrixAlgebraKit.initialize_output(
+    U[brow, bdim], S[bdim, bdim], Vt[bdim, bcol] = MatrixAlgebraKit.initialize_output(
       svd_compact!, block, block_alg
     )
   end

From ad8c6b2422b1381a6bb8bce07e0f8d2c9f974e73 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Wed, 30 Jul 2025 18:18:31 -0400
Subject: [PATCH 2/6] Revert changes to SVD

---
 Project.toml                                       |  3 +++
 .../BlockSparseArraysTensorProductsExt.jl          | 14 ++++++++++++++
 .../abstractblocksparsearray.jl                    |  2 +-
 src/factorizations/svd.jl                          | 12 +++++-------
 4 files changed, 23 insertions(+), 8 deletions(-)
 create mode 100644 ext/BlockSparseArraysTensorProductsExt/BlockSparseArraysTensorProductsExt.jl

diff --git a/Project.toml b/Project.toml
index db5a78bd..0ce971cc 100644
--- a/Project.toml
+++ b/Project.toml
@@ -22,9 +22,11 @@ TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138"
 
 [weakdeps]
 TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a"
+TensorProducts = "decf83d6-1968-43f4-96dc-fdb3fe15fc6d"
 
 [extensions]
 BlockSparseArraysTensorAlgebraExt = "TensorAlgebra"
+BlockSparseArraysTensorProductsExt = "TensorProducts"
 
 [compat]
 Adapt = "4.1.1"
@@ -43,6 +45,7 @@ MatrixAlgebraKit = "0.2.2"
 SparseArraysBase = "0.7.1"
 SplitApplyCombine = "1.2.3"
 TensorAlgebra = "0.3.2"
+TensorProducts = "0.1.7"
 Test = "1.10"
 TypeParameterAccessors = "0.4.1"
 julia = "1.10"
diff --git a/ext/BlockSparseArraysTensorProductsExt/BlockSparseArraysTensorProductsExt.jl b/ext/BlockSparseArraysTensorProductsExt/BlockSparseArraysTensorProductsExt.jl
new file mode 100644
index 00000000..d3733d89
--- /dev/null
+++ b/ext/BlockSparseArraysTensorProductsExt/BlockSparseArraysTensorProductsExt.jl
@@ -0,0 +1,14 @@
+module BlockSparseArraysTensorProductsExt
+
+using BlockSparseArrays: BlockUnitRange, blockrange, eachblockaxis
+using TensorProducts: TensorProducts, tensor_product
+# TODO: Dispatch on `FusionStyle` to allow different kinds of products,
+# for example to allow merging common symmetry sectors.
+function TensorProducts.tensor_product(a1::BlockUnitRange, a2::BlockUnitRange)
+  new_blockaxes = vec(
+    map(splat(tensor_product), Iterators.product(eachblockaxis(a1), eachblockaxis(a2)))
+  )
+  return blockrange(new_blockaxes)
+end
+
+end
diff --git a/src/abstractblocksparsearray/abstractblocksparsearray.jl b/src/abstractblocksparsearray/abstractblocksparsearray.jl
index 1824e830..5bc25c4b 100644
--- a/src/abstractblocksparsearray/abstractblocksparsearray.jl
+++ b/src/abstractblocksparsearray/abstractblocksparsearray.jl
@@ -93,7 +93,7 @@ function Base.setindex!(
   # https://github.com/JuliaLang/julia/pull/52487).
   # TODO: Delete once we drop support for Julia v1.10.
   aI = @view! a[I...]
-  copyto!(aI, value)
+  aI .= value
   return a
 end
 
diff --git a/src/factorizations/svd.jl b/src/factorizations/svd.jl
index d03d16b9..8c616d49 100644
--- a/src/factorizations/svd.jl
+++ b/src/factorizations/svd.jl
@@ -49,12 +49,12 @@ function MatrixAlgebraKit.initialize_output(
   ::typeof(svd_compact!), A::AbstractBlockSparseMatrix, alg::BlockPermutedDiagonalAlgorithm
 )
   bm, bn = blocksize(A)
-  (bmn, mindim) = findmin((bm, bn))
+  bmn = min(bm, bn)
 
   brows = eachblockaxis(axes(A, 1))
   bcols = eachblockaxis(axes(A, 2))
   u_axes = similar(brows, bmn)
-  v_axes = similar(bcols, bmn)
+  v_axes = similar(brows, bmn)
 
   # fill in values for blocks that are present
   bIs = collect(eachblockstoredindex(A))
@@ -62,9 +62,8 @@ function MatrixAlgebraKit.initialize_output(
   bcolIs = Int.(last.(Tuple.(bIs)))
   for bI in eachblockstoredindex(A)
     row, col = Int.(Tuple(bI))
-    dim = (row, col)[mindim]
-    u_axes[dim] = infimum(brows[row], bcols[col])
-    v_axes[dim] = infimum(bcols[col], brows[row])
+    u_axes[col] = infimum(brows[row], bcols[col])
+    v_axes[col] = infimum(bcols[col], brows[row])
   end
 
   # fill in values for blocks that aren't present, pairing them in order of occurence
@@ -84,10 +83,9 @@ function MatrixAlgebraKit.initialize_output(
   # allocate output
   for bI in eachblockstoredindex(A)
     brow, bcol = Tuple(bI)
-    bdim = (brow, bcol)[mindim]
     block = @view!(A[bI])
     block_alg = block_algorithm(alg, block)
-    U[brow, bdim], S[bdim, bdim], Vt[bdim, bcol] = MatrixAlgebraKit.initialize_output(
+    U[brow, bcol], S[bcol, bcol], Vt[bcol, bcol] = MatrixAlgebraKit.initialize_output(
       svd_compact!, block, block_alg
     )
   end

From 5944fb66abb6a38a80a7c1011c78e50ac8657728 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Thu, 31 Jul 2025 12:13:23 -0400
Subject: [PATCH 3/6] Fix tests

---
 .../BlockArraysExtensions.jl                  |  5 ++++-
 .../abstractblocksparsearray.jl               | 20 ++++++++++++-------
 src/blocksparsearray/blocksparsearray.jl      |  6 ++++--
 .../blocksparsearrayinterface.jl              |  4 ++++
 4 files changed, 25 insertions(+), 10 deletions(-)

diff --git a/src/BlockArraysExtensions/BlockArraysExtensions.jl b/src/BlockArraysExtensions/BlockArraysExtensions.jl
index 8b2a1eff..4a86a980 100644
--- a/src/BlockArraysExtensions/BlockArraysExtensions.jl
+++ b/src/BlockArraysExtensions/BlockArraysExtensions.jl
@@ -180,7 +180,10 @@ end
 # ```
 # but includes `BlockIndices`, where the blocks aren't contiguous.
 const BlockSliceCollection = Union{
-  Base.Slice,BlockSlice{<:BlockRange{1}},BlockIndices{<:Vector{<:Block{1}}}
+  Base.Slice,
+  BlockSlice{<:Block{1}},
+  BlockSlice{<:BlockRange{1}},
+  BlockIndices{<:Vector{<:Block{1}}},
 }
 const BlockIndexRangeSlice = BlockSlice{
   <:BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexRange{1}}}
diff --git a/src/abstractblocksparsearray/abstractblocksparsearray.jl b/src/abstractblocksparsearray/abstractblocksparsearray.jl
index 5bc25c4b..b58d2805 100644
--- a/src/abstractblocksparsearray/abstractblocksparsearray.jl
+++ b/src/abstractblocksparsearray/abstractblocksparsearray.jl
@@ -87,13 +87,19 @@ function Base.setindex!(
       ),
     )
   end
-  # Custom `_convert` works around the issue that
-  # `convert(::Type{<:Diagonal}, ::AbstractMatrix)` isnt' defined
-  # in Julia v1.10 (https://github.com/JuliaLang/julia/pull/48895,
-  # https://github.com/JuliaLang/julia/pull/52487).
-  # TODO: Delete once we drop support for Julia v1.10.
-  aI = @view! a[I...]
-  aI .= value
+  if !isconcretetype(blocktype(a))
+    # Custom `_convert` works around the issue that
+    # `convert(::Type{<:Diagonal}, ::AbstractMatrix)` isnt' defined
+    # in Julia v1.10 (https://github.com/JuliaLang/julia/pull/48895,
+    # https://github.com/JuliaLang/julia/pull/52487).
+    # TODO: Delete once we drop support for Julia v1.10.
+    blocks(a)[Int.(I)...] = _convert(blocktype(a), value)
+  else
+    # This writes into existing blocks, or constructs blocks
+    # using the axes.
+    aI = @view! a[I...]
+    aI .= value
+  end
   return a
 end
 
diff --git a/src/blocksparsearray/blocksparsearray.jl b/src/blocksparsearray/blocksparsearray.jl
index 030178ab..1ede2964 100644
--- a/src/blocksparsearray/blocksparsearray.jl
+++ b/src/blocksparsearray/blocksparsearray.jl
@@ -272,8 +272,10 @@ function blocksparse(d::Dict{<:Block,<:AbstractArray}, ax::Tuple)
   end
   return a
 end
-function blocksparse(d::Dict{<:Block,<:AbstractArray}, ax::AbstractUnitRange...)
-  return blocksparse(d, ax)
+function blocksparse(
+  d::Dict{<:Block,<:AbstractArray}, blocklens::AbstractVector{<:Integer}...
+)
+  return blocksparse(d, map(blockedrange, blocklens))
 end
 
 # Base `AbstractArray` interface
diff --git a/src/blocksparsearrayinterface/blocksparsearrayinterface.jl b/src/blocksparsearrayinterface/blocksparsearrayinterface.jl
index 30d75baf..e017e286 100644
--- a/src/blocksparsearrayinterface/blocksparsearrayinterface.jl
+++ b/src/blocksparsearrayinterface/blocksparsearrayinterface.jl
@@ -537,7 +537,11 @@ function SparseArraysBase.getunstoredindex(
   return error("Not implemented.")
 end
 
+# Convert a blockwise slice on a block sparse array
+# to an elementwise slice on the corresponding sparse array
+# of blocks.
 to_blocks_indices(I::BlockSlice{<:BlockRange{1}}) = Int.(I.block)
+to_blocks_indices(I::BlockSlice{<:Block{1}}) = Int(I.block):Int(I.block)
 to_blocks_indices(I::BlockIndices{<:Vector{<:Block{1}}}) = Int.(I.blocks)
 to_blocks_indices(I::Base.Slice) = Base.OneTo(blocklength(I.indices))
 

From bfbc49a0936582cc345f58333c4cc8c22e7ecf6e Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Fri, 8 Aug 2025 17:09:44 -0400
Subject: [PATCH 4/6] Add missing BlockUnitRange functionality

---
 src/BlockArraysExtensions/blockrange.jl            |  8 ++++++++
 .../abstractblocksparsearray.jl                    | 14 +++++++-------
 src/blocksparsearray/blocksparsearray.jl           | 13 ++++++++-----
 src/blocksparsearrayinterface/getunstoredblock.jl  |  7 +++++++
 4 files changed, 30 insertions(+), 12 deletions(-)

diff --git a/src/BlockArraysExtensions/blockrange.jl b/src/BlockArraysExtensions/blockrange.jl
index a1f291eb..ea991e13 100644
--- a/src/BlockArraysExtensions/blockrange.jl
+++ b/src/BlockArraysExtensions/blockrange.jl
@@ -8,6 +8,9 @@ end
 function blockrange(eachblockaxis)
   return BlockUnitRange(blockedrange(length.(eachblockaxis)), eachblockaxis)
 end
+function blockrange(first::Integer, eachblockaxis)
+  return BlockUnitRange(blockedrange(first, length.(eachblockaxis)), eachblockaxis)
+end
 Base.first(r::BlockUnitRange) = first(r.r)
 Base.last(r::BlockUnitRange) = last(r.r)
 BlockArrays.blocklasts(r::BlockUnitRange) = blocklasts(r.r)
@@ -15,6 +18,11 @@ eachblockaxis(r::BlockUnitRange) = r.eachblockaxis
 function Base.getindex(r::BlockUnitRange, I::Block{1})
   return eachblockaxis(r)[Int(I)] .+ (first(r.r[I]) - 1)
 end
+function Base.getindex(r::BlockUnitRange, I::BlockRange{1})
+  return blockrange(first(r), eachblockaxis(r)[Int.(I)])
+end
+Base.axes(r::BlockUnitRange) = (blockrange(eachblockaxis(r)),)
+Base.axes1(r::BlockUnitRange) = blockrange(eachblockaxis(r))
 
 using BlockArrays: BlockedOneTo
 const BlockOneTo{T<:Integer,B,CS,R<:BlockedOneTo{T,CS}} = BlockUnitRange{T,B,CS,R}
diff --git a/src/abstractblocksparsearray/abstractblocksparsearray.jl b/src/abstractblocksparsearray/abstractblocksparsearray.jl
index b58d2805..5d5805cd 100644
--- a/src/abstractblocksparsearray/abstractblocksparsearray.jl
+++ b/src/abstractblocksparsearray/abstractblocksparsearray.jl
@@ -87,18 +87,18 @@ function Base.setindex!(
       ),
     )
   end
-  if !isconcretetype(blocktype(a))
+  if isstored(a, I...)
+    # This writes into existing blocks, or constructs blocks
+    # using the axes.
+    aI = @view! a[I...]
+    aI .= value
+  else
     # Custom `_convert` works around the issue that
     # `convert(::Type{<:Diagonal}, ::AbstractMatrix)` isnt' defined
     # in Julia v1.10 (https://github.com/JuliaLang/julia/pull/48895,
     # https://github.com/JuliaLang/julia/pull/52487).
-    # TODO: Delete once we drop support for Julia v1.10.
+    # TODO: Delete `_convert` once we drop support for Julia v1.10.
     blocks(a)[Int.(I)...] = _convert(blocktype(a), value)
-  else
-    # This writes into existing blocks, or constructs blocks
-    # using the axes.
-    aI = @view! a[I...]
-    aI .= value
   end
   return a
 end
diff --git a/src/blocksparsearray/blocksparsearray.jl b/src/blocksparsearray/blocksparsearray.jl
index 1ede2964..97426861 100644
--- a/src/blocksparsearray/blocksparsearray.jl
+++ b/src/blocksparsearray/blocksparsearray.jl
@@ -262,11 +262,14 @@ function blocksparsezeros(::BlockType{A}, axes...) where {A<:AbstractArray}
   # to make a bit more generic.
   return BlockSparseArray{eltype(A),ndims(A),A}(undef, axes...)
 end
-function blocksparse(d::Dict{<:Block,<:AbstractArray}, ax::Tuple)
-  blockaxtype = Tuple{map(eltype ∘ eachblockaxis, ax)...}
-  # TODO: Catch if inference fails and use `valtype(d)` instead.
-  blockt = Base.promote_op(similar, Type{valtype(d)}, blockaxtype)
-  a = blocksparsezeros(BlockType(blockt), ax)
+function blocksparse(d::Dict, ax::Tuple)
+  a = blocksparsezeros(BlockType(valtype(d)), ax...)
+
+  ## blockaxtype = Tuple{map(eltype ∘ eachblockaxis, ax)...}
+  ## # TODO: Catch if inference fails and use `valtype(d)` instead.
+  ## blockt = Base.promote_op(similar, Type{valtype(d)}, blockaxtype)
+  ## a = blocksparsezeros(BlockType(blockt), ax)
+
   for I in eachindex(d)
     a[I] = d[I]
   end
diff --git a/src/blocksparsearrayinterface/getunstoredblock.jl b/src/blocksparsearrayinterface/getunstoredblock.jl
index c5ae1dca..0669096f 100644
--- a/src/blocksparsearrayinterface/getunstoredblock.jl
+++ b/src/blocksparsearrayinterface/getunstoredblock.jl
@@ -18,10 +18,17 @@ function Base.AbstractArray{A}(a::ZeroBlocks{N}) where {N,A}
 end
 
 @inline function Base.getindex(a::ZeroBlocks{N,A}, I::Vararg{Int,N}) where {N,A}
+
+  #foreach(display, a.parentaxes)
+  #foreach(display ∘ eachblockaxis, a.parentaxes)
+
   # TODO: Use `BlockArrays.eachblockaxes`.
   ax = ntuple(N) do d
     return eachblockaxis(a.parentaxes[d])[I[d]]
   end
+
+  #@show ax
+
   !isconcretetype(A) && return zero!(similar(Array{eltype(A),N}, ax))
   return zero!(similar(A, ax))
 end

From da70733cfba95fb173d0d37f2a5c9e9cc09e1b00 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Mon, 11 Aug 2025 21:05:35 -0400
Subject: [PATCH 5/6] Bump to v0.9

---
 Project.toml                                   |  2 +-
 docs/Project.toml                              |  2 +-
 examples/Project.toml                          |  2 +-
 .../BlockSparseArraysTensorAlgebraExt.jl       | 18 ++++++++++++++----
 src/blocksparsearray/blocksparsearray.jl       | 14 ++++++++------
 test/Project.toml                              |  2 +-
 6 files changed, 26 insertions(+), 14 deletions(-)

diff --git a/Project.toml b/Project.toml
index 0ce971cc..1e1903ec 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.8.2"
+version = "0.9.0"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
diff --git a/docs/Project.toml b/docs/Project.toml
index 7885be60..edbaa0bc 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -6,6 +6,6 @@ Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 
 [compat]
 BlockArrays = "1"
-BlockSparseArrays = "0.8"
+BlockSparseArrays = "0.9"
 Documenter = "1"
 Literate = "2"
diff --git a/examples/Project.toml b/examples/Project.toml
index 6c17bd70..7e02e0de 100644
--- a/examples/Project.toml
+++ b/examples/Project.toml
@@ -5,5 +5,5 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 BlockArrays = "1"
-BlockSparseArrays = "0.8"
+BlockSparseArrays = "0.9"
 Test = "1"
diff --git a/ext/BlockSparseArraysTensorAlgebraExt/BlockSparseArraysTensorAlgebraExt.jl b/ext/BlockSparseArraysTensorAlgebraExt/BlockSparseArraysTensorAlgebraExt.jl
index fc4d5940..7912ef51 100644
--- a/ext/BlockSparseArraysTensorAlgebraExt/BlockSparseArraysTensorAlgebraExt.jl
+++ b/ext/BlockSparseArraysTensorAlgebraExt/BlockSparseArraysTensorAlgebraExt.jl
@@ -24,10 +24,20 @@ function TensorAlgebra.matricize(
 )
   ax = fuseaxes(axes(a), biperm)
   reshaped_blocks_a = reshape(blocks(a), map(blocklength, ax))
-  bs = Dict(
-    Block(Tuple(I)) => matricize(reshaped_blocks_a[I], biperm) for
-    I in eachstoredindex(reshaped_blocks_a)
-  )
+  key(I) = Block(Tuple(I))
+  value(I) = matricize(reshaped_blocks_a[I], biperm)
+  Is = eachstoredindex(reshaped_blocks_a)
+  bs = if isempty(Is)
+    # Catch empty case and make sure the type is constrained properly.
+    # This seems to only be necessary in Julia versions below v1.11,
+    # try removing it when we drop support for those versions.
+    keytype = Base.promote_op(key, eltype(Is))
+    valtype = Base.promote_op(value, eltype(Is))
+    valtype′ = !isconcretetype(valtype) ? AbstractMatrix{eltype(a)} : valtype
+    Dict{keytype,valtype′}()
+  else
+    Dict(key(I) => value(I) for I in Is)
+  end
   return blocksparse(bs, ax)
 end
 
diff --git a/src/blocksparsearray/blocksparsearray.jl b/src/blocksparsearray/blocksparsearray.jl
index 97426861..c1b98e95 100644
--- a/src/blocksparsearray/blocksparsearray.jl
+++ b/src/blocksparsearray/blocksparsearray.jl
@@ -263,13 +263,15 @@ function blocksparsezeros(::BlockType{A}, axes...) where {A<:AbstractArray}
   return BlockSparseArray{eltype(A),ndims(A),A}(undef, axes...)
 end
 function blocksparse(d::Dict, ax::Tuple)
+  d_narrow = if isempty(d)
+    Dict{Block{length(ax),Int},AbstractArray{Any,length(ax)}}()
+  else
+    Dict(Pair(kv) for kv in d)
+  end
+  return blocksparse(d_narrow, ax)
+end
+function blocksparse(d::Dict{<:Block{N},<:AbstractArray{<:Any,N}}, ax::Tuple) where {N}
   a = blocksparsezeros(BlockType(valtype(d)), ax...)
-
-  ## blockaxtype = Tuple{map(eltype ∘ eachblockaxis, ax)...}
-  ## # TODO: Catch if inference fails and use `valtype(d)` instead.
-  ## blockt = Base.promote_op(similar, Type{valtype(d)}, blockaxtype)
-  ## a = blocksparsezeros(BlockType(blockt), ax)
-
   for I in eachindex(d)
     a[I] = d[I]
   end
diff --git a/test/Project.toml b/test/Project.toml
index 52ef8620..12fb6f03 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -25,7 +25,7 @@ Adapt = "4"
 Aqua = "0.8"
 ArrayLayouts = "1"
 BlockArrays = "1"
-BlockSparseArrays = "0.8"
+BlockSparseArrays = "0.9"
 DerivableInterfaces = "0.5"
 DiagonalArrays = "0.3"
 GPUArraysCore = "0.2"

From a6efc60e6e768a2767cfc6b7acfc38254b683c4a Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Mon, 11 Aug 2025 22:41:21 -0400
Subject: [PATCH 6/6] Stricter blocksparse constructor

---
 src/blocksparsearray/blocksparsearray.jl | 10 +---------
 1 file changed, 1 insertion(+), 9 deletions(-)

diff --git a/src/blocksparsearray/blocksparsearray.jl b/src/blocksparsearray/blocksparsearray.jl
index c1b98e95..f708f345 100644
--- a/src/blocksparsearray/blocksparsearray.jl
+++ b/src/blocksparsearray/blocksparsearray.jl
@@ -262,15 +262,7 @@ function blocksparsezeros(::BlockType{A}, axes...) where {A<:AbstractArray}
   # to make a bit more generic.
   return BlockSparseArray{eltype(A),ndims(A),A}(undef, axes...)
 end
-function blocksparse(d::Dict, ax::Tuple)
-  d_narrow = if isempty(d)
-    Dict{Block{length(ax),Int},AbstractArray{Any,length(ax)}}()
-  else
-    Dict(Pair(kv) for kv in d)
-  end
-  return blocksparse(d_narrow, ax)
-end
-function blocksparse(d::Dict{<:Block{N},<:AbstractArray{<:Any,N}}, ax::Tuple) where {N}
+function blocksparse(d::Dict{<:Block,<:AbstractArray}, ax::Tuple)
   a = blocksparsezeros(BlockType(valtype(d)), ax...)
   for I in eachindex(d)
     a[I] = d[I]