From 4e4429260413d3d6583f983f65af279e140312f3 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Tue, 1 Jul 2025 18:15:05 -0400
Subject: [PATCH 01/11] Towards more general truncation and slicing

---
 Project.toml                                  |  2 +-
 src/BlockArraysExtensions/blockedunitrange.jl |  8 +++-
 src/abstractblocksparsearray/linearalgebra.jl | 20 ++++++++-
 src/factorizations/truncation.jl              | 45 +++++++++++--------
 4 files changed, 53 insertions(+), 22 deletions(-)

diff --git a/Project.toml b/Project.toml
index 17a00d07..318969d5 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.7.21"
+version = "0.7.22"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
diff --git a/src/BlockArraysExtensions/blockedunitrange.jl b/src/BlockArraysExtensions/blockedunitrange.jl
index 2adf198f..cdd90161 100644
--- a/src/BlockArraysExtensions/blockedunitrange.jl
+++ b/src/BlockArraysExtensions/blockedunitrange.jl
@@ -314,6 +314,8 @@ end
 # `Base.getindex(a::Block, b...)`.
 _getindex(a::Block{N}, b::Vararg{Any,N}) where {N} = GenericBlockIndex(a, b)
 _getindex(a::Block{N}, b::Vararg{Integer,N}) where {N} = a[b...]
+_getindex(a::Block{N}, b::Vararg{AbstractUnitRange{<:Integer},N}) where {N} = a[b...]
+_getindex(a::Block{N}, b::Vararg{AbstractVector,N}) where {N} = BlockIndexVector(a, b)
 # Fix ambiguity.
 _getindex(a::Block{0}) = a[]
 
@@ -372,7 +374,11 @@ function blockedunitrange_getindices(
   a::AbstractBlockedUnitRange,
   indices::BlockVector{<:GenericBlockIndex{1},<:Vector{<:BlockIndexVector{1}}},
 )
-  return mortar(map(b -> a[b], blocks(indices)))
+  blks = map(b -> a[b], blocks(indices))
+  # Preserve any extra structure in the axes, like a
+  # Kronecker structure, symmetry sectors, etc.
+  ax = mortar_axis(map(b -> axis(a[b]), blocks(indices)))
+  return mortar(blks, (ax,))
 end
 
 # This is a specialization of `BlockArrays.unblock`:
diff --git a/src/abstractblocksparsearray/linearalgebra.jl b/src/abstractblocksparsearray/linearalgebra.jl
index 0c9483fb..2b03de9e 100644
--- a/src/abstractblocksparsearray/linearalgebra.jl
+++ b/src/abstractblocksparsearray/linearalgebra.jl
@@ -1,4 +1,4 @@
-using LinearAlgebra: LinearAlgebra, Adjoint, Transpose, norm, tr
+using LinearAlgebra: LinearAlgebra, Adjoint, Transpose, diag, norm, tr
 
 # Like: https://github.com/JuliaLang/julia/blob/v1.11.1/stdlib/LinearAlgebra/src/transpose.jl#L184
 # but also takes the dual of the axes.
@@ -33,6 +33,24 @@ function LinearAlgebra.tr(a::AnyAbstractBlockSparseMatrix)
   return tr_a
 end
 
+# TODO: Define in DiagonalArrays.jl.
+function diagaxis(a::AbstractArray)
+  LinearAlgebra.checksquare(a)
+  return axes(a, 1)
+end
+function LinearAlgebra.diag(a::AnyAbstractBlockSparseMatrix)
+  # TODO: Add `checkblocksquare` to also check it is square blockwise.
+  LinearAlgebra.checksquare(a)
+  diagaxes = map(blockdiagindices(a)) do I
+    return diagaxis(@view(a[I]))
+  end
+  r = blockrange(diagaxes)
+  stored_blocks = Dict((
+    Tuple(I)[1] => diag(@view!(a[I])) for I in eachstoredblockdiagindex(a)
+  ))
+  return blocksparse(stored_blocks, (r,))
+end
+
 # TODO: Define `SparseArraysBase.isdiag`, define as
 # `isdiag(blocks(a))`.
 function blockisdiag(a::AbstractArray)
diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl
index 1ad4e72e..07b6924e 100644
--- a/src/factorizations/truncation.jl
+++ b/src/factorizations/truncation.jl
@@ -1,14 +1,11 @@
-using MatrixAlgebraKit: TruncationStrategy, diagview, eig_trunc!, eigh_trunc!, svd_trunc!
-
-function MatrixAlgebraKit.diagview(A::BlockSparseMatrix{T,Diagonal{T,Vector{T}}}) where {T}
-  D = BlockSparseVector{T}(undef, axes(A, 1))
-  for I in eachblockstoredindex(A)
-    if ==(Int.(Tuple(I))...)
-      D[Tuple(I)[1]] = diagview(A[I])
-    end
-  end
-  return D
-end
+using MatrixAlgebraKit:
+  TruncationStrategy,
+  diagview,
+  eig_trunc!,
+  eigh_trunc!,
+  findtruncated,
+  svd_trunc!,
+  truncate!
 
 """
     BlockPermutedDiagonalTruncationStrategy(strategy::TruncationStrategy)
@@ -27,7 +24,7 @@ function MatrixAlgebraKit.truncate!(
   strategy::TruncationStrategy,
 )
   # TODO assert blockdiagonal
-  return MatrixAlgebraKit.truncate!(
+  return truncate!(
     svd_trunc!, (U, S, Vᴴ), BlockPermutedDiagonalTruncationStrategy(strategy)
   )
 end
@@ -38,9 +35,7 @@ for f in [:eig_trunc!, :eigh_trunc!]
       (D, V)::NTuple{2,AbstractBlockSparseMatrix},
       strategy::TruncationStrategy,
     )
-      return MatrixAlgebraKit.truncate!(
-        $f, (D, V), BlockPermutedDiagonalTruncationStrategy(strategy)
-      )
+      return truncate!($f, (D, V), BlockPermutedDiagonalTruncationStrategy(strategy))
     end
   end
 end
@@ -50,10 +45,22 @@ end
 function MatrixAlgebraKit.findtruncated(
   values::AbstractVector, strategy::BlockPermutedDiagonalTruncationStrategy
 )
-  ind = MatrixAlgebraKit.findtruncated(values, strategy.strategy)
+  ind = findtruncated(Vector(values), strategy.strategy)
   indexmask = falses(length(values))
   indexmask[ind] .= true
-  return indexmask
+  return to_truncated_indices(values, indexmask)
+end
+
+# Allow customizing the indices output by `findtruncated`
+# based on the type of `values`, for example to preserve
+# a block or Kronecker structure.
+to_truncated_indices(values, I) = I
+function to_truncated_indices(values::AbstractBlockVector, I::AbstractVector{Bool})
+  I′ = BlockedVector(I, blocklengths(axis(values)))
+  blocks = map(BlockRange(values)) do b
+    return _getindex(b, to_truncated_indices(values[b], I′[b]))
+  end
+  return blocks
 end
 
 function MatrixAlgebraKit.truncate!(
@@ -61,7 +68,7 @@ function MatrixAlgebraKit.truncate!(
   (U, S, Vᴴ)::NTuple{3,AbstractBlockSparseMatrix},
   strategy::BlockPermutedDiagonalTruncationStrategy,
 )
-  I = MatrixAlgebraKit.findtruncated(diagview(S), strategy)
+  I = MatrixAlgebraKit.findtruncated(diag(S), strategy)
   return (U[:, I], S[I, I], Vᴴ[I, :])
 end
 for f in [:eig_trunc!, :eigh_trunc!]
@@ -71,7 +78,7 @@ for f in [:eig_trunc!, :eigh_trunc!]
       (D, V)::NTuple{2,AbstractBlockSparseMatrix},
       strategy::BlockPermutedDiagonalTruncationStrategy,
     )
-      I = MatrixAlgebraKit.findtruncated(diagview(D), strategy)
+      I = MatrixAlgebraKit.findtruncated(diag(D), strategy)
       return (D[I, I], V[:, I])
     end
   end

From 5f4563d7d04d442582c982232ab9e19febb962e2 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Wed, 2 Jul 2025 10:37:14 -0400
Subject: [PATCH 02/11] Preserve sector information better in slicing

---
 .../BlockArraysExtensions.jl                  |  6 ++
 src/BlockArraysExtensions/blockedunitrange.jl |  6 +-
 test/test_factorizations.jl                   | 64 ++++++++++++-------
 3 files changed, 51 insertions(+), 25 deletions(-)

diff --git a/src/BlockArraysExtensions/BlockArraysExtensions.jl b/src/BlockArraysExtensions/BlockArraysExtensions.jl
index c4589782..4a194d24 100644
--- a/src/BlockArraysExtensions/BlockArraysExtensions.jl
+++ b/src/BlockArraysExtensions/BlockArraysExtensions.jl
@@ -68,6 +68,12 @@ for f in (:axes, :unsafe_indices, :axes1, :first, :last, :size, :length, :unsafe
 end
 Base.getindex(S::BlockIndices, i::Integer) = getindex(S.indices, i)
 
+# TODO: Move this to a `BlockArraysExtensions` library.
+function blockedunitrange_getindices(a::AbstractBlockedUnitRange, indices::BlockIndices)
+  # TODO: Is this a good definition? It ignores `indices.indices`.
+  return a[indices.blocks]
+end
+
 # Generalization of to `BlockArrays._blockslice`:
 # https://github.com/JuliaArrays/BlockArrays.jl/blob/v1.6.3/src/views.jl#L13-L14
 # Used by `BlockArrays.unblock`, which is used in `to_indices`
diff --git a/src/BlockArraysExtensions/blockedunitrange.jl b/src/BlockArraysExtensions/blockedunitrange.jl
index cdd90161..604bef09 100644
--- a/src/BlockArraysExtensions/blockedunitrange.jl
+++ b/src/BlockArraysExtensions/blockedunitrange.jl
@@ -368,7 +368,11 @@ function blockedunitrange_getindices(
   a::AbstractBlockedUnitRange,
   indices::BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexVector{1}}},
 )
-  return mortar(map(b -> a[b], blocks(indices)))
+  blks = map(b -> a[b], blocks(indices))
+  # Preserve any extra structure in the axes, like a
+  # Kronecker structure, symmetry sectors, etc.
+  ax = mortar_axis(map(b -> axis(a[b]), blocks(indices)))
+  return mortar(blks, (ax,))
 end
 function blockedunitrange_getindices(
   a::AbstractBlockedUnitRange,
diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl
index 5113497b..651a8f65 100644
--- a/test/test_factorizations.jl
+++ b/test/test_factorizations.jl
@@ -146,20 +146,23 @@ test_params = Iterators.product(blockszs, eltypes)
   @test test_svd(a, usv_empty)
 
   # test blockdiagonal
+  rng = StableRNG(123)
   for i in LinearAlgebra.diagind(blocks(a))
     I = CartesianIndices(blocks(a))[i]
-    a[Block(I.I...)] = rand(T, size(blocks(a)[i]))
+    a[Block(I.I...)] = rand(rng, T, size(blocks(a)[i]))
   end
   usv = svd_compact(a)
   @test test_svd(a, usv)
 
-  perm = Random.randperm(length(m))
+  rng = StableRNG(123)
+  perm = Random.randperm(rng, length(m))
   b = a[Block.(perm), Block.(1:length(n))]
   usv = svd_compact(b)
   @test test_svd(b, usv)
 
   # test permuted blockdiagonal with missing row/col
-  I_removed = rand(eachblockstoredindex(b))
+  rng = StableRNG(123)
+  I_removed = rand(rng, eachblockstoredindex(b))
   c = copy(b)
   delete!(blocks(c).storage, CartesianIndex(Int.(Tuple(I_removed))))
   usv = svd_compact(c)
@@ -176,20 +179,23 @@ end
   @test test_svd(a, usv_empty; full=true)
 
   # test blockdiagonal
+  rng = StableRNG(123)
   for i in LinearAlgebra.diagind(blocks(a))
     I = CartesianIndices(blocks(a))[i]
-    a[Block(I.I...)] = rand(T, size(blocks(a)[i]))
+    a[Block(I.I...)] = rand(rng, T, size(blocks(a)[i]))
   end
   usv = svd_full(a)
   @test test_svd(a, usv; full=true)
 
-  perm = Random.randperm(length(m))
+  rng = StableRNG(123)
+  perm = Random.randperm(rng, length(m))
   b = a[Block.(perm), Block.(1:length(n))]
   usv = svd_full(b)
   @test test_svd(b, usv; full=true)
 
   # test permuted blockdiagonal with missing row/col
-  I_removed = rand(eachblockstoredindex(b))
+  rng = StableRNG(123)
+  I_removed = rand(rng, eachblockstoredindex(b))
   c = copy(b)
   delete!(blocks(c).storage, CartesianIndex(Int.(Tuple(I_removed))))
   usv = svd_full(c)
@@ -203,9 +209,10 @@ end
   a = BlockSparseArray{T}(undef, m, n)
 
   # test blockdiagonal
+  rng = StableRNG(123)
   for i in LinearAlgebra.diagind(blocks(a))
     I = CartesianIndices(blocks(a))[i]
-    a[Block(I.I...)] = rand(T, size(blocks(a)[i]))
+    a[Block(I.I...)] = rand(rng, T, size(blocks(a)[i]))
   end
 
   minmn = min(size(a)...)
@@ -236,7 +243,8 @@ end
   @test (V1ᴴ * V1ᴴ' ≈ LinearAlgebra.I)
 
   # test permuted blockdiagonal
-  perm = Random.randperm(length(m))
+  rng = StableRNG(123)
+  perm = Random.randperm(rng, length(m))
   b = a[Block.(perm), Block.(1:length(n))]
   for trunc in (truncrank(r), trunctol(atol))
     U1, S1, V1ᴴ = svd_trunc(b; trunc)
@@ -270,8 +278,9 @@ end
 @testset "qr_compact (T=$T)" for T in (Float32, Float64, ComplexF32, ComplexF64)
   for i in [2, 3], j in [2, 3], k in [2, 3], l in [2, 3]
     A = BlockSparseArray{T}(undef, ([i, j], [k, l]))
-    A[Block(1, 1)] = randn(T, i, k)
-    A[Block(2, 2)] = randn(T, j, l)
+    rng = StableRNG(123)
+    A[Block(1, 1)] = randn(rng, T, i, k)
+    A[Block(2, 2)] = randn(rng, T, j, l)
     Q, R = qr_compact(A)
     @test Matrix(Q'Q) ≈ LinearAlgebra.I
     @test A ≈ Q * R
@@ -281,8 +290,9 @@ end
 @testset "qr_full (T=$T)" for T in (Float32, Float64, ComplexF32, ComplexF64)
   for i in [2, 3], j in [2, 3], k in [2, 3], l in [2, 3]
     A = BlockSparseArray{T}(undef, ([i, j], [k, l]))
-    A[Block(1, 1)] = randn(T, i, k)
-    A[Block(2, 2)] = randn(T, j, l)
+    rng = StableRNG(123)
+    A[Block(1, 1)] = randn(rng, T, i, k)
+    A[Block(2, 2)] = randn(rng, T, j, l)
     Q, R = qr_full(A)
     Q′, R′ = qr_full(Matrix(A))
     @test size(Q) == size(Q′)
@@ -296,8 +306,9 @@ end
 @testset "lq_compact" for T in (Float32, Float64, ComplexF32, ComplexF64)
   for i in [2, 3], j in [2, 3], k in [2, 3], l in [2, 3]
     A = BlockSparseArray{T}(undef, ([i, j], [k, l]))
-    A[Block(1, 1)] = randn(T, i, k)
-    A[Block(2, 2)] = randn(T, j, l)
+    rng = StableRNG(123)
+    A[Block(1, 1)] = randn(rng, T, i, k)
+    A[Block(2, 2)] = randn(rng, T, j, l)
     L, Q = lq_compact(A)
     @test Matrix(Q * Q') ≈ LinearAlgebra.I
     @test A ≈ L * Q
@@ -307,8 +318,9 @@ end
 @testset "lq_full" for T in (Float32, Float64, ComplexF32, ComplexF64)
   for i in [2, 3], j in [2, 3], k in [2, 3], l in [2, 3]
     A = BlockSparseArray{T}(undef, ([i, j], [k, l]))
-    A[Block(1, 1)] = randn(T, i, k)
-    A[Block(2, 2)] = randn(T, j, l)
+    rng = StableRNG(123)
+    A[Block(1, 1)] = randn(rng, T, i, k)
+    A[Block(2, 2)] = randn(rng, T, j, l)
     L, Q = lq_full(A)
     L′, Q′ = lq_full(Matrix(A))
     @test size(L) == size(L′)
@@ -321,8 +333,9 @@ end
 
 @testset "left_polar (T=$T)" for T in (Float32, Float64, ComplexF32, ComplexF64)
   A = BlockSparseArray{T}(undef, ([3, 4], [2, 3]))
-  A[Block(1, 1)] = randn(T, 3, 2)
-  A[Block(2, 2)] = randn(T, 4, 3)
+  rng = StableRNG(123)
+  A[Block(1, 1)] = randn(rng, T, 3, 2)
+  A[Block(2, 2)] = randn(rng, T, 4, 3)
 
   U, C = left_polar(A)
   @test U * C ≈ A
@@ -331,8 +344,9 @@ end
 
 @testset "right_polar (T=$T)" for T in (Float32, Float64, ComplexF32, ComplexF64)
   A = BlockSparseArray{T}(undef, ([2, 3], [3, 4]))
-  A[Block(1, 1)] = randn(T, 2, 3)
-  A[Block(2, 2)] = randn(T, 3, 4)
+  rng = StableRNG(123)
+  A[Block(1, 1)] = randn(rng, T, 2, 3)
+  A[Block(2, 2)] = randn(rng, T, 3, 4)
 
   C, U = right_polar(A)
   @test C * U ≈ A
@@ -341,8 +355,9 @@ end
 
 @testset "left_orth (T=$T)" for T in (Float32, Float64, ComplexF32, ComplexF64)
   A = BlockSparseArray{T}(undef, ([3, 4], [2, 3]))
-  A[Block(1, 1)] = randn(T, 3, 2)
-  A[Block(2, 2)] = randn(T, 4, 3)
+  rng = StableRNG(123)
+  A[Block(1, 1)] = randn(rng, T, 3, 2)
+  A[Block(2, 2)] = randn(rng, T, 4, 3)
 
   for kind in (:polar, :qr, :svd)
     U, C = left_orth(A; kind)
@@ -358,8 +373,9 @@ end
 
 @testset "right_orth (T=$T)" for T in (Float32, Float64, ComplexF32, ComplexF64)
   A = BlockSparseArray{T}(undef, ([2, 3], [3, 4]))
-  A[Block(1, 1)] = randn(T, 2, 3)
-  A[Block(2, 2)] = randn(T, 3, 4)
+  rng = StableRNG(123)
+  A[Block(1, 1)] = randn(rng, T, 2, 3)
+  A[Block(2, 2)] = randn(rng, T, 3, 4)
 
   for kind in (:lq, :polar, :svd)
     C, U = right_orth(A; kind)

From 9d4edc1b5f900ee9583064d28b01f3b8243e8b68 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Wed, 2 Jul 2025 10:43:59 -0400
Subject: [PATCH 03/11] Namespace

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

diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl
index 07b6924e..d173f88d 100644
--- a/src/factorizations/truncation.jl
+++ b/src/factorizations/truncation.jl
@@ -68,7 +68,7 @@ function MatrixAlgebraKit.truncate!(
   (U, S, Vᴴ)::NTuple{3,AbstractBlockSparseMatrix},
   strategy::BlockPermutedDiagonalTruncationStrategy,
 )
-  I = MatrixAlgebraKit.findtruncated(diag(S), strategy)
+  I = findtruncated(diag(S), strategy)
   return (U[:, I], S[I, I], Vᴴ[I, :])
 end
 for f in [:eig_trunc!, :eigh_trunc!]
@@ -78,7 +78,7 @@ for f in [:eig_trunc!, :eigh_trunc!]
       (D, V)::NTuple{2,AbstractBlockSparseMatrix},
       strategy::BlockPermutedDiagonalTruncationStrategy,
     )
-      I = MatrixAlgebraKit.findtruncated(diag(D), strategy)
+      I = findtruncated(diag(D), strategy)
       return (D[I, I], V[:, I])
     end
   end

From ac9b71cf46879c9722f5b702e1b278d4f918cf37 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Thu, 17 Jul 2025 18:10:37 -0400
Subject: [PATCH 04/11] Better support for slicing with Colon

---
 Project.toml                                       | 2 +-
 src/BlockArraysExtensions/BlockArraysExtensions.jl | 1 +
 src/BlockArraysExtensions/blockedunitrange.jl      | 9 +++++++++
 src/BlockArraysExtensions/blockrange.jl            | 5 +++++
 src/abstractblocksparsearray/views.jl              | 8 ++++++++
 5 files changed, 24 insertions(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 318969d5..1bb211fd 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.7.22"
+version = "0.7.23"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
diff --git a/src/BlockArraysExtensions/BlockArraysExtensions.jl b/src/BlockArraysExtensions/BlockArraysExtensions.jl
index 4a194d24..3b0ee09f 100644
--- a/src/BlockArraysExtensions/BlockArraysExtensions.jl
+++ b/src/BlockArraysExtensions/BlockArraysExtensions.jl
@@ -185,6 +185,7 @@ const GenericBlockIndexVectorSlices = BlockIndices{
   <:BlockVector{<:GenericBlockIndex{1},<:Vector{<:BlockIndexVector}}
 }
 const SubBlockSliceCollection = Union{
+  Base.Slice,
   BlockIndexRangeSlice,
   BlockIndexRangeSlices,
   BlockIndexVectorSlices,
diff --git a/src/BlockArraysExtensions/blockedunitrange.jl b/src/BlockArraysExtensions/blockedunitrange.jl
index 604bef09..b111c83f 100644
--- a/src/BlockArraysExtensions/blockedunitrange.jl
+++ b/src/BlockArraysExtensions/blockedunitrange.jl
@@ -349,6 +349,15 @@ BlockArrays.Block(b::BlockIndexVector) = b.block
 
 Base.copy(a::BlockIndexVector) = BlockIndexVector(a.block, copy.(a.indices))
 
+# Copied from BlockArrays.BlockIndexRange.
+function Base.show(io::IO, B::BlockIndexVector)
+  show(io, Block(B))
+  print(io, "[")
+  print_tuple_elements(io, B.indices)
+  print(io, "]")
+end
+Base.show(io::IO, ::MIME"text/plain", B::BlockIndexVector) = show(io, B)
+
 function Base.getindex(b::AbstractBlockedUnitRange, Kkr::BlockIndexVector{1})
   return b[block(Kkr)][Kkr.indices...]
 end
diff --git a/src/BlockArraysExtensions/blockrange.jl b/src/BlockArraysExtensions/blockrange.jl
index 7edd0133..106bd10b 100644
--- a/src/BlockArraysExtensions/blockrange.jl
+++ b/src/BlockArraysExtensions/blockrange.jl
@@ -16,6 +16,11 @@ function Base.getindex(r::BlockUnitRange, I::Block{1})
   return eachblockaxis(r)[Int(I)] .+ (first(r.r[I]) - 1)
 end
 
+const BlockOneTo{T<:Integer,B,CS,R<:BlockedOneTo{T,CS}} = BlockUnitRange{T,B,CS,R}
+Base.axes(S::Base.Slice{<:BlockOneTo}) = (S.indices,)
+Base.axes1(S::Base.Slice{<:BlockOneTo}) = S.indices
+Base.unsafe_indices(S::Base.Slice{<:BlockOneTo}) = (S.indices,)
+
 function BlockArrays.combine_blockaxes(r1::BlockUnitRange, r2::BlockUnitRange)
   if eachblockaxis(r1) ≠ eachblockaxis(r2)
     return throw(ArgumentError("BlockUnitRanges must have the same block axes"))
diff --git a/src/abstractblocksparsearray/views.jl b/src/abstractblocksparsearray/views.jl
index 084cb2cf..9db3a8b9 100644
--- a/src/abstractblocksparsearray/views.jl
+++ b/src/abstractblocksparsearray/views.jl
@@ -302,6 +302,14 @@ end
 
 blockedslice_blocks(x::BlockSlice) = x.block
 blockedslice_blocks(x::BlockIndices) = x.blocks
+# Reinterpret the slice blockwise.
+function blockedslice_blocks(x::Base.Slice)
+  return mortar(
+    map(BlockRange(x.indices)) do b
+      return BlockIndexRange(b, Base.Slice(Base.axes1(x.indices[b])))
+    end,
+  )
+end
 
 # TODO: Define `@interface interface(a) viewblock`.
 function BlockArrays.viewblock(

From 4f70c04354f23533ade3c080d72407141eb42598 Mon Sep 17 00:00:00 2001
From: Matt Fishman <mtfishman@users.noreply.github.com>
Date: Thu, 17 Jul 2025 18:20:18 -0400
Subject: [PATCH 05/11] Update Project.toml

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 1bb211fd..318969d5 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.7.23"
+version = "0.7.22"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

From f5c47e09dbf4650879fab3f7cd5e5560a592d0ea Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Thu, 17 Jul 2025 18:27:45 -0400
Subject: [PATCH 06/11] Namespace

---
 src/BlockArraysExtensions/blockrange.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/BlockArraysExtensions/blockrange.jl b/src/BlockArraysExtensions/blockrange.jl
index 106bd10b..a1f291eb 100644
--- a/src/BlockArraysExtensions/blockrange.jl
+++ b/src/BlockArraysExtensions/blockrange.jl
@@ -16,6 +16,7 @@ function Base.getindex(r::BlockUnitRange, I::Block{1})
   return eachblockaxis(r)[Int(I)] .+ (first(r.r[I]) - 1)
 end
 
+using BlockArrays: BlockedOneTo
 const BlockOneTo{T<:Integer,B,CS,R<:BlockedOneTo{T,CS}} = BlockUnitRange{T,B,CS,R}
 Base.axes(S::Base.Slice{<:BlockOneTo}) = (S.indices,)
 Base.axes1(S::Base.Slice{<:BlockOneTo}) = S.indices

From 41947d1163e9ff771656d35dc9fc439e9fcdd5a6 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Thu, 17 Jul 2025 18:29:14 -0400
Subject: [PATCH 07/11] Bump version

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 318969d5..1bb211fd 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.7.22"
+version = "0.7.23"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

From 0b4e871b14133d1e426c7bb3a154e2444f805a9c Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Fri, 18 Jul 2025 14:34:00 -0400
Subject: [PATCH 08/11] Start fixing ambiguities

---
 src/abstractblocksparsearray/views.jl | 17 +++++++++++++++++
 1 file changed, 17 insertions(+)

diff --git a/src/abstractblocksparsearray/views.jl b/src/abstractblocksparsearray/views.jl
index 9db3a8b9..77788cd4 100644
--- a/src/abstractblocksparsearray/views.jl
+++ b/src/abstractblocksparsearray/views.jl
@@ -189,6 +189,18 @@ function Base.view(
 ) where {T,N}
   return viewblock(a, block...)
 end
+# Fix ambiguity error with BlockArrays.jl.
+function Base.view(
+  a::SubArray{
+    T,
+    N,
+    <:AbstractBlockSparseArray{T,N},
+    <:Tuple{Vararg{Union{Base.Slice,BlockSlice{Union{},<:Integer}},N}},
+  },
+  block::Block{N},
+) where {T,N}
+  return viewblock(a, block)
+end
 
 # XXX: TODO: Distinguish if a sub-view of the block needs to be taken!
 # Define a new `SubBlockSlice` which is used in:
@@ -327,6 +339,11 @@ function BlockArrays.viewblock(
   end
   return @view parent(a)[brs...]
 end
+
+## function BlockArrays.viewblock
+##   a::SubArray{
+##     T,N,<:AbstractBlockSparseArray{T,N},<:Tuple{Vararg{Union{Base.Slice, BlockArrays.BlockSlice{Union{}, T} where T<:Integer, BlockSparseArrays.BlockIndices{Union{}, T} where T<:Integer},N}}, ::Vararg{Block{1}, N}) where {T, N}
+
 # TODO: Define `@interface interface(a) viewblock`.
 function BlockArrays.viewblock(
   a::SubArray{

From f3ff0d7ea13c3d6e0f6b1f98f4854e6034e3bac3 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Tue, 22 Jul 2025 10:55:53 -0400
Subject: [PATCH 09/11] Fix ambiguities

---
 src/abstractblocksparsearray/views.jl | 39 ++++++++++++++++++---------
 1 file changed, 26 insertions(+), 13 deletions(-)

diff --git a/src/abstractblocksparsearray/views.jl b/src/abstractblocksparsearray/views.jl
index 77788cd4..e585937d 100644
--- a/src/abstractblocksparsearray/views.jl
+++ b/src/abstractblocksparsearray/views.jl
@@ -189,19 +189,27 @@ function Base.view(
 ) where {T,N}
   return viewblock(a, block...)
 end
-# Fix ambiguity error with BlockArrays.jl.
+
+# Disambiguate between block reindexing of blockwise views
+# (`BlockSliceCollection`) and subblockwise views (`SubBlockSliceCollection`),
+# which both include `Base.Slice`.
 function Base.view(
-  a::SubArray{
-    T,
-    N,
-    <:AbstractBlockSparseArray{T,N},
-    <:Tuple{Vararg{Union{Base.Slice,BlockSlice{Union{},<:Integer}},N}},
-  },
+  a::SubArray{T,N,<:AbstractBlockSparseArray{T,N},<:Tuple{Vararg{Base.Slice,N}}},
   block::Block{N},
 ) where {T,N}
   return viewblock(a, block)
 end
 
+# Block reindexing of blockwise views (`BlockSliceCollection`).
+function viewblock_blockslice(a::SubArray{<:Any,N}, block::Vararg{Block{1},N}) where {N}
+  I = CartesianIndex(Int.(block))
+  # TODO: Use `eachblockstoredindex`.
+  if I ∈ eachstoredindex(blocks(a))
+    return blocks(a)[I]
+  end
+  return BlockView(parent(a), Block.(Base.reindex(parentindices(blocks(a)), Tuple(I))))
+end
+
 # XXX: TODO: Distinguish if a sub-view of the block needs to be taken!
 # Define a new `SubBlockSlice` which is used in:
 # `@interface interface(a) to_indices(a, inds, I::Tuple{UnitRange{<:Integer},Vararg{Any}})`
@@ -211,12 +219,17 @@ function BlockArrays.viewblock(
   a::SubArray{T,N,<:AbstractBlockSparseArray{T,N},<:Tuple{Vararg{BlockSliceCollection,N}}},
   block::Vararg{Block{1},N},
 ) where {T,N}
-  I = CartesianIndex(Int.(block))
-  # TODO: Use `eachblockstoredindex`.
-  if I ∈ eachstoredindex(blocks(a))
-    return blocks(a)[I]
-  end
-  return BlockView(parent(a), Block.(Base.reindex(parentindices(blocks(a)), Tuple(I))))
+  return viewblock_blockslice(a, block...)
+end
+
+# Disambiguate between block reindexing of blockwise views
+# (`BlockSliceCollection`) and subblockwise views (`SubBlockSliceCollection`),
+# which both include `Base.Slice`.
+function BlockArrays.viewblock(
+  a::SubArray{T,N,<:AbstractBlockSparseArray{T,N},<:Tuple{Vararg{Base.Slice,N}}},
+  block::Vararg{Block{1},N},
+) where {T,N}
+  return viewblock_blockslice(a, block...)
 end
 
 function to_blockindexrange(

From b4991a9667d489a74b8e7bebff2a35c7c54f50e8 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Tue, 22 Jul 2025 15:36:05 -0400
Subject: [PATCH 10/11] Fix tests

---
 .../BlockArraysExtensions.jl                  | 10 ++++++++
 src/abstractblocksparsearray/views.jl         |  9 +++++++
 .../blocksparsearrayinterface.jl              | 24 ++++++++++++-------
 3 files changed, 34 insertions(+), 9 deletions(-)

diff --git a/src/BlockArraysExtensions/BlockArraysExtensions.jl b/src/BlockArraysExtensions/BlockArraysExtensions.jl
index 3b0ee09f..58155403 100644
--- a/src/BlockArraysExtensions/BlockArraysExtensions.jl
+++ b/src/BlockArraysExtensions/BlockArraysExtensions.jl
@@ -88,6 +88,9 @@ function _blockslice(x, y::AbstractVector)
   return BlockIndices(x, y)
 end
 
+# TODO: Constrain the type of `BlockIndices` more, this seems
+# to assume that `S.blocks` is a list of blocks as opposed to
+# a flat list of block indices like the definition below.
 function Base.getindex(S::BlockIndices, i::BlockSlice{<:Block{1}})
   # TODO: Check that `i.indices` is consistent with `S.indices`.
   # It seems like this isn't handling the case where `i` is a
@@ -96,6 +99,13 @@ function Base.getindex(S::BlockIndices, i::BlockSlice{<:Block{1}})
   return _blockslice(S.blocks[Int(Block(i))], S.indices[Block(i)])
 end
 
+function Base.getindex(
+  S::BlockIndices{<:AbstractBlockVector{<:BlockIndex{1}}}, i::BlockSlice{<:Block{1}}
+)
+  @assert length(S.indices[Block(i)]) == length(i.indices)
+  return _blockslice(S.blocks[Block(i)], S.indices[Block(i)])
+end
+
 # This is used in slicing like:
 # a = BlockSparseArray{Float64}([2, 2, 2, 2], [2, 2, 2, 2])
 # I = BlockedVector([Block(4), Block(3), Block(2), Block(1)], [2, 2])
diff --git a/src/abstractblocksparsearray/views.jl b/src/abstractblocksparsearray/views.jl
index e585937d..ad9de691 100644
--- a/src/abstractblocksparsearray/views.jl
+++ b/src/abstractblocksparsearray/views.jl
@@ -316,6 +316,15 @@ function Base.view(
 ) where {T,N}
   return viewblock(a, block...)
 end
+# Fix ambiguity error.
+function Base.view(
+  a::SubArray{
+    T,N,<:AbstractBlockSparseArray{T,N},<:Tuple{Vararg{SubBlockSliceCollection,N}}
+  },
+  block::Vararg{Block{1},N},
+) where {T,N}
+  return viewblock(a, block...)
+end
 function BlockArrays.viewblock(
   a::SubArray{
     T,N,<:AbstractBlockSparseArray{T,N},<:Tuple{Vararg{SubBlockSliceCollection,N}}
diff --git a/src/blocksparsearrayinterface/blocksparsearrayinterface.jl b/src/blocksparsearrayinterface/blocksparsearrayinterface.jl
index 9ca44e5a..30d75baf 100644
--- a/src/blocksparsearrayinterface/blocksparsearrayinterface.jl
+++ b/src/blocksparsearrayinterface/blocksparsearrayinterface.jl
@@ -222,19 +222,25 @@ end
 
 # a[mortar([Block(1)[1:2], Block(2)[1:3]]), mortar([Block(1)[1:2], Block(2)[1:3]])]
 # a[[Block(1)[1:2], Block(2)[1:3]], [Block(1)[1:2], Block(2)[1:3]]]
-@interface ::AbstractBlockSparseArrayInterface function Base.to_indices(
-  a, inds, I::Tuple{BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexRange{1}}},Vararg{Any}}
-)
-  I1 = BlockIndices(I[1], blockedunitrange_getindices(inds[1], I[1]))
-  return (I1, to_indices(a, Base.tail(inds), Base.tail(I))...)
-end
-
 @interface ::AbstractBlockSparseArrayInterface function Base.to_indices(
   a,
   inds,
-  I::Tuple{BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexVector{1}}},Vararg{Any}},
+  I::Tuple{
+    BlockVector{<:BlockIndex{1},<:Vector{<:Union{BlockIndexRange{1},BlockIndexVector{1}}}},
+    Vararg{Any},
+  },
 )
-  I1 = BlockIndices(I[1], blockedunitrange_getindices(inds[1], I[1]))
+  # Index the `inds` by the `BlockIndexRange`/`BlockIndexVector`s on each block
+  # in order to canonicalize the indices and preserve metadata,
+  # such as sector data for symmetric tensors.
+  bs = mortar(
+    map(blocks(I[1])) do bi
+      b = Block(bi)
+      binds = only(bi.indices)
+      return BlockIndexVector(b, Base.axes1(inds[1][b])[binds])
+    end,
+  )
+  I1 = BlockIndices(bs, blockedunitrange_getindices(inds[1], I[1]))
   return (I1, to_indices(a, Base.tail(inds), Base.tail(I))...)
 end
 @interface ::AbstractBlockSparseArrayInterface function Base.to_indices(

From 2b006570fe99e53e41cb55e2bac8d0f378a284cf Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Tue, 22 Jul 2025 15:37:52 -0400
Subject: [PATCH 11/11] Cleanup

---
 src/abstractblocksparsearray/views.jl | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/src/abstractblocksparsearray/views.jl b/src/abstractblocksparsearray/views.jl
index ad9de691..19de8a31 100644
--- a/src/abstractblocksparsearray/views.jl
+++ b/src/abstractblocksparsearray/views.jl
@@ -362,10 +362,6 @@ function BlockArrays.viewblock(
   return @view parent(a)[brs...]
 end
 
-## function BlockArrays.viewblock
-##   a::SubArray{
-##     T,N,<:AbstractBlockSparseArray{T,N},<:Tuple{Vararg{Union{Base.Slice, BlockArrays.BlockSlice{Union{}, T} where T<:Integer, BlockSparseArrays.BlockIndices{Union{}, T} where T<:Integer},N}}, ::Vararg{Block{1}, N}) where {T, N}
-
 # TODO: Define `@interface interface(a) viewblock`.
 function BlockArrays.viewblock(
   a::SubArray{