diff --git a/NDTensors/Project.toml b/NDTensors/Project.toml
index dead86e738..bd56eda80c 100644
--- a/NDTensors/Project.toml
+++ b/NDTensors/Project.toml
@@ -1,7 +1,7 @@
 name = "NDTensors"
 uuid = "23ae76d9-e61a-49c4-8f12-3f1a16adf9cf"
 authors = ["Matthew Fishman <mfishman@flatironinstitute.org>"]
-version = "0.3.35"
+version = "0.3.36"
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
diff --git a/NDTensors/src/lib/BlockSparseArrays/ext/BlockSparseArraysGradedAxesExt/src/BlockSparseArraysGradedAxesExt.jl b/NDTensors/src/lib/BlockSparseArrays/ext/BlockSparseArraysGradedAxesExt/src/BlockSparseArraysGradedAxesExt.jl
index 9a92e8b08d..adf7b35de0 100644
--- a/NDTensors/src/lib/BlockSparseArrays/ext/BlockSparseArraysGradedAxesExt/src/BlockSparseArraysGradedAxesExt.jl
+++ b/NDTensors/src/lib/BlockSparseArrays/ext/BlockSparseArraysGradedAxesExt/src/BlockSparseArraysGradedAxesExt.jl
@@ -85,29 +85,6 @@ function Base.axes(a::Adjoint{<:Any,<:AbstractBlockSparseMatrix})
   return dual.(reverse(axes(a')))
 end
 
-# TODO: Delete this definition in favor of the one in
-# GradedAxes once https://github.com/JuliaArrays/BlockArrays.jl/pull/405 is merged.
-# TODO: Make a special definition for `BlockedVector{<:Block{1}}` in order
-# to merge blocks.
-function GradedAxes.blockedunitrange_getindices(
-  a::AbstractBlockedUnitRange, indices::AbstractVector{<:Union{Block{1},BlockIndexRange{1}}}
-)
-  # Without converting `indices` to `Vector`,
-  # mapping `indices` outputs a `BlockVector`
-  # which is harder to reason about.
-  blocks = map(index -> a[index], Vector(indices))
-  # We pass `length.(blocks)` to `mortar` in order
-  # to pass block labels to the axes of the output,
-  # if they exist. This makes it so that
-  # `only(axes(a[indices])) isa `GradedUnitRange`
-  # if `a isa `GradedUnitRange`, for example.
-  # TODO: Remove `unlabel` once `BlockArray` axes
-  # type is generalized in BlockArrays.jl.
-  # TODO: Support using `BlockSparseVector`, need
-  # to make more `BlockSparseArray` constructors.
-  return BlockSparseArray(blocks, (blockedrange(length.(blocks)),))
-end
-
 # This definition is only needed since calls like
 # `a[[Block(1), Block(2)]]` where `a isa AbstractGradedUnitRange`
 # returns a `BlockSparseVector` instead of a `BlockVector`
diff --git a/NDTensors/src/lib/BlockSparseArrays/ext/BlockSparseArraysGradedAxesExt/test/runtests.jl b/NDTensors/src/lib/BlockSparseArrays/ext/BlockSparseArraysGradedAxesExt/test/runtests.jl
index 638dae82db..a7a54fb9cb 100644
--- a/NDTensors/src/lib/BlockSparseArrays/ext/BlockSparseArraysGradedAxesExt/test/runtests.jl
+++ b/NDTensors/src/lib/BlockSparseArrays/ext/BlockSparseArraysGradedAxesExt/test/runtests.jl
@@ -56,8 +56,7 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
       @test size(b) == (4, 4, 4, 4)
       @test blocksize(b) == (2, 2, 2, 2)
       @test blocklengths.(axes(b)) == ([2, 2], [2, 2], [2, 2], [2, 2])
-      # TODO: Fix this for `BlockedArray`.
-      @test_broken nstored(b) == 256
+      @test nstored(b) == 256
       # TODO: Fix this for `BlockedArray`.
       @test_broken block_nstored(b) == 16
       for i in 1:ndims(a)
diff --git a/NDTensors/src/lib/BlockSparseArrays/src/BlockArraysExtensions/BlockArraysExtensions.jl b/NDTensors/src/lib/BlockSparseArrays/src/BlockArraysExtensions/BlockArraysExtensions.jl
index bb2634ce14..fb015e7ee4 100644
--- a/NDTensors/src/lib/BlockSparseArrays/src/BlockArraysExtensions/BlockArraysExtensions.jl
+++ b/NDTensors/src/lib/BlockSparseArrays/src/BlockArraysExtensions/BlockArraysExtensions.jl
@@ -4,6 +4,7 @@ using BlockArrays:
   AbstractBlockVector,
   Block,
   BlockRange,
+  BlockedOneTo,
   BlockedUnitRange,
   BlockVector,
   BlockSlice,
@@ -19,19 +20,6 @@ using Dictionaries: Dictionary, Indices
 using ..GradedAxes: blockedunitrange_getindices
 using ..SparseArrayInterface: stored_indices
 
-# GenericBlockSlice works around an issue that the indices of BlockSlice
-# are restricted to Int element type.
-# TODO: Raise an issue/make a pull request in BlockArrays.jl.
-struct GenericBlockSlice{B,T<:Integer,I<:AbstractUnitRange{T}} <: AbstractUnitRange{T}
-  block::B
-  indices::I
-end
-BlockArrays.Block(bs::GenericBlockSlice{<:Block}) = bs.block
-for f in (:axes, :unsafe_indices, :axes1, :first, :last, :size, :length, :unsafe_length)
-  @eval Base.$f(S::GenericBlockSlice) = Base.$f(S.indices)
-end
-Base.getindex(S::GenericBlockSlice, i::Integer) = getindex(S.indices, i)
-
 # BlockIndices works around an issue that the indices of BlockSlice
 # are restricted to AbstractUnitRange{Int}.
 struct BlockIndices{B,T<:Integer,I<:AbstractVector{T}} <: AbstractVector{T}
@@ -42,6 +30,63 @@ for f in (:axes, :unsafe_indices, :axes1, :first, :last, :size, :length, :unsafe
   @eval Base.$f(S::BlockIndices) = Base.$f(S.indices)
 end
 Base.getindex(S::BlockIndices, i::Integer) = getindex(S.indices, i)
+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
+  # subslice of a block correctly (i.e. it ignores `i.indices`).
+  @assert length(S.indices[Block(i)]) == length(i.indices)
+  return BlockSlice(S.blocks[Int(Block(i))], S.indices[Block(i)])
+end
+function Base.getindex(S::BlockIndices, i::BlockSlice{<:BlockRange{1}})
+  # TODO: Check that `i.indices` is consistent with `S.indices`.
+  # TODO: Turn this into a `blockedunitrange_getindices` definition.
+  subblocks = S.blocks[Int.(i.block)]
+  subindices = mortar(
+    map(1:length(i.block)) do I
+      r = blocks(i.indices)[I]
+      return S.indices[first(r)]:S.indices[last(r)]
+    end,
+  )
+  return BlockIndices(subblocks, subindices)
+end
+
+# TODO: This is type piracy. This is used in `reindex` when making
+# views of blocks of sliced block arrays, for example:
+# ```julia
+# a = BlockSparseArray{elt}(undef, ([2, 3], [2, 3]))
+# b = @view a[[Block(1)[1:1], Block(2)[1:2]], [Block(1)[1:1], Block(2)[1:2]]]
+# b[Block(1, 1)]
+# ```
+# Without this change, BlockArrays has the slicing behavior:
+# ```julia
+# julia> mortar([Block(1)[1:1], Block(2)[1:2]])[BlockSlice(Block(2), 2:3)]
+# 2-element Vector{BlockIndex{1, Tuple{Int64}, Tuple{Int64}}}:
+#  Block(2)[1]
+#  Block(2)[2]
+# ```
+# while with this change it has the slicing behavior:
+# ```julia
+# julia> mortar([Block(1)[1:1], Block(2)[1:2]])[BlockSlice(Block(2), 2:3)]
+# Block(2)[1:2]
+# ```
+# i.e. it preserves the types of the blocks better. Upstream this fix to
+# BlockArrays.jl. Also consider overloading `reindex` so that it calls
+# a custom `getindex` function to avoid type piracy in the meantime.
+# Also fix this in BlockArrays:
+# ```julia
+# julia> mortar([Block(1)[1:1], Block(2)[1:2]])[Block(2)]
+# 2-element Vector{BlockIndex{1, Tuple{Int64}, Tuple{Int64}}}:
+#  Block(2)[1]
+#  Block(2)[2]
+# ```
+function Base.getindex(
+  a::BlockVector{<:BlockIndex{1},<:AbstractVector{<:BlockIndexRange{1}}},
+  I::BlockSlice{<:Block{1}},
+)
+  # Check that the block slice corresponds to the correct block.
+  @assert I.indices == only(axes(a))[Block(I)]
+  return blocks(a)[Int(Block(I))]
+end
 
 # Outputs a `BlockUnitRange`.
 function sub_axis(a::AbstractUnitRange, indices)
@@ -185,15 +230,12 @@ function blockrange(axis::AbstractUnitRange, r::AbstractVector{<:Block{1}})
   return r
 end
 
-using BlockArrays: BlockSlice
-function blockrange(axis::AbstractUnitRange, r::BlockSlice)
-  return blockrange(axis, r.block)
+function blockrange(axis::BlockedOneTo{<:Integer}, r::BlockVector{<:Integer})
+  return error("Slicing not implemented for range of type `$(typeof(r))`.")
 end
 
-# GenericBlockSlice works around an issue that the indices of BlockSlice
-# are restricted to Int element type.
-# TODO: Raise an issue/make a pull request in BlockArrays.jl.
-function blockrange(axis::AbstractUnitRange, r::GenericBlockSlice)
+using BlockArrays: BlockSlice
+function blockrange(axis::AbstractUnitRange, r::BlockSlice)
   return blockrange(axis, r.block)
 end
 
diff --git a/NDTensors/src/lib/BlockSparseArrays/src/BlockSparseArrays.jl b/NDTensors/src/lib/BlockSparseArrays/src/BlockSparseArrays.jl
index 09064999e4..a31d0269d5 100644
--- a/NDTensors/src/lib/BlockSparseArrays/src/BlockSparseArrays.jl
+++ b/NDTensors/src/lib/BlockSparseArrays/src/BlockSparseArrays.jl
@@ -8,7 +8,7 @@ include("abstractblocksparsearray/abstractblocksparsearray.jl")
 include("abstractblocksparsearray/wrappedabstractblocksparsearray.jl")
 include("abstractblocksparsearray/abstractblocksparsematrix.jl")
 include("abstractblocksparsearray/abstractblocksparsevector.jl")
-include("abstractblocksparsearray/view.jl")
+include("abstractblocksparsearray/views.jl")
 include("abstractblocksparsearray/arraylayouts.jl")
 include("abstractblocksparsearray/sparsearrayinterface.jl")
 include("abstractblocksparsearray/broadcast.jl")
diff --git a/NDTensors/src/lib/BlockSparseArrays/src/abstractblocksparsearray/broadcast.jl b/NDTensors/src/lib/BlockSparseArrays/src/abstractblocksparsearray/broadcast.jl
index c00f3f9e1c..5887673036 100644
--- a/NDTensors/src/lib/BlockSparseArrays/src/abstractblocksparsearray/broadcast.jl
+++ b/NDTensors/src/lib/BlockSparseArrays/src/abstractblocksparsearray/broadcast.jl
@@ -12,7 +12,7 @@ function Broadcast.BroadcastStyle(
       <:Any,
       <:Any,
       <:AbstractBlockSparseArray,
-      <:Tuple{BlockSlice{<:Any,<:AbstractBlockedUnitRange},Vararg{Any}},
+      <:Tuple{BlockSlice{<:Any,<:Any,<:AbstractBlockedUnitRange},Vararg{Any}},
     },
   },
 )
@@ -25,8 +25,8 @@ function Broadcast.BroadcastStyle(
       <:Any,
       <:AbstractBlockSparseArray,
       <:Tuple{
-        BlockSlice{<:Any,<:AbstractBlockedUnitRange},
-        BlockSlice{<:Any,<:AbstractBlockedUnitRange},
+        BlockSlice{<:Any,<:Any,<:AbstractBlockedUnitRange},
+        BlockSlice{<:Any,<:Any,<:AbstractBlockedUnitRange},
         Vararg{Any},
       },
     },
@@ -40,7 +40,7 @@ function Broadcast.BroadcastStyle(
       <:Any,
       <:Any,
       <:AbstractBlockSparseArray,
-      <:Tuple{Any,BlockSlice{<:Any,<:AbstractBlockedUnitRange},Vararg{Any}},
+      <:Tuple{Any,BlockSlice{<:Any,<:Any,<:AbstractBlockedUnitRange},Vararg{Any}},
     },
   },
 )
diff --git a/NDTensors/src/lib/BlockSparseArrays/src/abstractblocksparsearray/map.jl b/NDTensors/src/lib/BlockSparseArrays/src/abstractblocksparsearray/map.jl
index 582d4193a3..c376ff631a 100644
--- a/NDTensors/src/lib/BlockSparseArrays/src/abstractblocksparsearray/map.jl
+++ b/NDTensors/src/lib/BlockSparseArrays/src/abstractblocksparsearray/map.jl
@@ -104,10 +104,10 @@ end
 
 # TODO: Why isn't this calling `mapreduce` already?
 function Base.iszero(a::BlockSparseArrayLike)
-  return sparse_iszero(a)
+  return sparse_iszero(blocks(a))
 end
 
 # TODO: Why isn't this calling `mapreduce` already?
 function Base.isreal(a::BlockSparseArrayLike)
-  return sparse_isreal(a)
+  return sparse_isreal(blocks(a))
 end
diff --git a/NDTensors/src/lib/BlockSparseArrays/src/abstractblocksparsearray/sparsearrayinterface.jl b/NDTensors/src/lib/BlockSparseArrays/src/abstractblocksparsearray/sparsearrayinterface.jl
index b92c82ae81..ac866afc9a 100644
--- a/NDTensors/src/lib/BlockSparseArrays/src/abstractblocksparsearray/sparsearrayinterface.jl
+++ b/NDTensors/src/lib/BlockSparseArrays/src/abstractblocksparsearray/sparsearrayinterface.jl
@@ -32,3 +32,7 @@ end
 function SparseArrayInterface.sparse_storage(a::AbstractBlockSparseArray)
   return BlockSparseStorage(a)
 end
+
+function SparseArrayInterface.nstored(a::BlockSparseArrayLike)
+  return sum(nstored, sparse_storage(blocks(a)); init=zero(Int))
+end
diff --git a/NDTensors/src/lib/BlockSparseArrays/src/abstractblocksparsearray/view.jl b/NDTensors/src/lib/BlockSparseArrays/src/abstractblocksparsearray/view.jl
deleted file mode 100644
index 61f303384f..0000000000
--- a/NDTensors/src/lib/BlockSparseArrays/src/abstractblocksparsearray/view.jl
+++ /dev/null
@@ -1,29 +0,0 @@
-using BlockArrays: BlockIndexRange, BlockRange, BlockSlice, block
-
-# TODO: Define `AnyBlockSparseVector`.
-function Base.view(a::BlockSparseArrayLike{<:Any,N}, index::Block{N}) where {N}
-  return blocksparse_view(a, index)
-end
-
-# Fix ambiguity error with `BlockArrays`.
-function Base.view(
-  a::SubArray{
-    <:Any,
-    N,
-    <:AbstractBlockSparseArray,
-    <:Tuple{
-      Vararg{
-        Union{Base.Slice,BlockSlice{<:BlockRange{1,<:Tuple{AbstractUnitRange{Int}}}}},N
-      },
-    },
-  },
-  index::Block{N},
-) where {N}
-  return blocksparse_view(a, index)
-end
-
-# Fix ambiguity error with `BlockArrays`.
-# TODO: Define `AnyBlockSparseVector`.
-function Base.view(a::BlockSparseArrayLike{<:Any,1}, index::Block{1})
-  return blocksparse_view(a, index)
-end
diff --git a/NDTensors/src/lib/BlockSparseArrays/src/abstractblocksparsearray/views.jl b/NDTensors/src/lib/BlockSparseArrays/src/abstractblocksparsearray/views.jl
new file mode 100644
index 0000000000..c09ac27890
--- /dev/null
+++ b/NDTensors/src/lib/BlockSparseArrays/src/abstractblocksparsearray/views.jl
@@ -0,0 +1,24 @@
+using BlockArrays: Block, BlockSlices
+
+function blocksparse_view(a, I...)
+  return Base.invoke(view, Tuple{AbstractArray,Vararg{Any}}, a, I...)
+end
+
+# These definitions circumvent some generic definitions in BlockArrays.jl:
+# https://github.com/JuliaArrays/BlockArrays.jl/blob/master/src/views.jl
+# which don't handle subslices of blocks properly.
+function Base.view(
+  a::SubArray{<:Any,N,<:BlockSparseArrayLike,<:NTuple{N,BlockSlices}}, I::Block{N}
+) where {N}
+  return blocksparse_view(a, I)
+end
+function Base.view(
+  a::SubArray{<:Any,N,<:BlockSparseArrayLike,<:NTuple{N,BlockSlices}}, I::Vararg{Block{1},N}
+) where {N}
+  return blocksparse_view(a, I...)
+end
+function Base.view(
+  V::SubArray{<:Any,1,<:BlockSparseArrayLike,<:Tuple{BlockSlices}}, I::Block{1}
+)
+  return blocksparse_view(a, I)
+end
diff --git a/NDTensors/src/lib/BlockSparseArrays/src/abstractblocksparsearray/wrappedabstractblocksparsearray.jl b/NDTensors/src/lib/BlockSparseArrays/src/abstractblocksparsearray/wrappedabstractblocksparsearray.jl
index 0a81605262..aec6783308 100644
--- a/NDTensors/src/lib/BlockSparseArrays/src/abstractblocksparsearray/wrappedabstractblocksparsearray.jl
+++ b/NDTensors/src/lib/BlockSparseArrays/src/abstractblocksparsearray/wrappedabstractblocksparsearray.jl
@@ -18,86 +18,32 @@ const BlockSparseArrayLike{T,N} = Union{
   <:AbstractBlockSparseArray{T,N},<:WrappedAbstractBlockSparseArray{T,N}
 }
 
-# Used when making views.
-# TODO: Move to blocksparsearrayinterface.
-function blocksparse_to_indices(a, inds, I)
-  return (unblock(a, inds, I), to_indices(a, BlockArrays._maybetail(inds), Base.tail(I))...)
-end
-
-# TODO: Move to blocksparsearrayinterface.
-function blocksparse_to_indices(a, I)
-  return to_indices(a, axes(a), I)
-end
-
-# Used when making views.
+# a[1:2, 1:2]
 function Base.to_indices(
-  a::BlockSparseArrayLike, inds, I::Tuple{AbstractVector{<:Block{1}},Vararg{Any}}
+  a::BlockSparseArrayLike, inds, I::Tuple{UnitRange{<:Integer},Vararg{Any}}
 )
   return blocksparse_to_indices(a, inds, I)
 end
 
+# a[[Block(2), Block(1)], [Block(2), Block(1)]]
 function Base.to_indices(
-  a::BlockSparseArrayLike, inds, I::Tuple{AbstractUnitRange{<:Integer},Vararg{Any}}
+  a::BlockSparseArrayLike, inds, I::Tuple{Vector{<:Block{1}},Vararg{Any}}
 )
   return blocksparse_to_indices(a, inds, I)
 end
 
-# Fixes ambiguity error with BlockArrays.
-function Base.to_indices(a::BlockSparseArrayLike, inds, I::Tuple{BlockRange{1},Vararg{Any}})
-  return blocksparse_to_indices(a, inds, I)
-end
-
+# a[[Block(1)[1:2], Block(2)[1:2]], [Block(1)[1:2], Block(2)[1:2]]]
 function Base.to_indices(
-  a::BlockSparseArrayLike, I::Tuple{AbstractVector{<:Block{1}},Vararg{Any}}
-)
-  return blocksparse_to_indices(a, I)
-end
-
-# Handle case of indexing with `[Block(1)[1:2], Block(2)[1:2]]`
-# by converting it to a `BlockVector` with
-# `mortar([Block(1)[1:2], Block(2)[1:2]])`.
-function Base.to_indices(
-  a::BlockSparseArrayLike, inds, I::Tuple{AbstractVector{<:BlockIndexRange{1}},Vararg{Any}}
+  a::BlockSparseArrayLike, inds, I::Tuple{Vector{<:BlockIndexRange{1}},Vararg{Any}}
 )
   return to_indices(a, inds, (mortar(I[1]), Base.tail(I)...))
 end
 
-# Fixes ambiguity error with BlockArrays.
-function Base.to_indices(a::BlockSparseArrayLike, I::Tuple{BlockRange{1},Vararg{Any}})
-  return blocksparse_to_indices(a, I)
-end
-
+# a[BlockVector([Block(2), Block(1)], [2]), BlockVector([Block(2), Block(1)], [2])]
 function Base.to_indices(
-  a::BlockSparseArrayLike, I::Tuple{AbstractUnitRange{<:Integer},Vararg{Any}}
-)
-  return blocksparse_to_indices(a, I)
-end
-
-# Used inside `Base.to_indices` when making views.
-# TODO: Move to blocksparsearrayinterface.
-# TODO: Make a special definition for `BlockedVector{<:Block{1}}` in order
-# to merge blocks.
-function blocksparse_unblock(a, inds, I::Tuple{AbstractVector{<:Block{1}},Vararg{Any}})
-  return BlockIndices(I[1], blockedunitrange_getindices(inds[1], I[1]))
-end
-
-# TODO: Move to blocksparsearrayinterface.
-function blocksparse_unblock(a, inds, I::Tuple{AbstractUnitRange{<:Integer},Vararg{Any}})
-  bs = blockrange(inds[1], I[1])
-  # GenericBlockSlice works around an issue that the indices of BlockSlice
-  # are restricted to Int element type.
-  # TODO: Raise an issue/make a pull request in BlockArrays.jl.
-  return GenericBlockSlice(bs, blockedunitrange_getindices(inds[1], I[1]))
-end
-
-function BlockArrays.unblock(a, inds, I::Tuple{AbstractVector{<:Block{1}},Vararg{Any}})
-  return blocksparse_unblock(a, inds, I)
-end
-
-function BlockArrays.unblock(
-  a::BlockSparseArrayLike, inds, I::Tuple{AbstractUnitRange{<:Integer},Vararg{Any}}
+  a::BlockSparseArrayLike, inds, I::Tuple{BlockVector{<:Block{1}},Vararg{Any}}
 )
-  return blocksparse_unblock(a, inds, I)
+  return blocksparse_to_indices(a, inds, I)
 end
 
 # BlockArrays `AbstractBlockArray` interface
@@ -135,15 +81,6 @@ function Base.getindex(
   return ArrayLayouts.layout_getindex(a, I...)
 end
 
-function Base.getindex(a::BlockSparseArrayLike{<:Any,N}, block::Block{N}) where {N}
-  return blocksparse_getindex(a, block)
-end
-function Base.getindex(
-  a::BlockSparseArrayLike{<:Any,N}, block::Vararg{Block{1},N}
-) where {N}
-  return blocksparse_getindex(a, block...)
-end
-
 # TODO: Define `blocksparse_isassigned`.
 function Base.isassigned(
   a::BlockSparseArrayLike{<:Any,N}, index::Vararg{Block{1},N}
@@ -167,20 +104,8 @@ function Base.setindex!(a::BlockSparseArrayLike{<:Any,N}, value, I::BlockIndex{N
   blocksparse_setindex!(a, value, I)
   return a
 end
-
-function Base.setindex!(
-  a::BlockSparseArrayLike{<:Any,N}, value, I::Vararg{Block{1},N}
-) where {N}
-  a[Block(Int.(I))] = value
-  return a
-end
-function Base.setindex!(a::BlockSparseArrayLike{<:Any,N}, value, I::Block{N}) where {N}
-  blocksparse_setindex!(a, value, I)
-  return a
-end
-
-# Fix ambiguity error
-function Base.setindex!(a::BlockSparseArrayLike{<:Any,1}, value, I::Block{1})
+# Fixes ambiguity error with BlockArrays.jl
+function Base.setindex!(a::BlockSparseArrayLike{<:Any,1}, value, I::BlockIndex{1})
   blocksparse_setindex!(a, value, I)
   return a
 end
@@ -205,12 +130,6 @@ function Base.fill!(a::BlockSparseArrayLike, value)
   return a
 end
 
-# `BlockArrays` interface
-# TODO: Is this needed if `blocks` is defined?
-function BlockArrays.viewblock(a::BlockSparseArrayLike{<:Any,N}, I::Block{N,Int}) where {N}
-  return blocksparse_viewblock(a, I)
-end
-
 # Needed by `BlockArrays` matrix multiplication interface
 function Base.similar(
   arraytype::Type{<:BlockSparseArrayLike}, axes::Tuple{Vararg{AbstractUnitRange{<:Integer}}}
diff --git a/NDTensors/src/lib/BlockSparseArrays/src/blocksparsearrayinterface/blocksparsearrayinterface.jl b/NDTensors/src/lib/BlockSparseArrays/src/blocksparsearrayinterface/blocksparsearrayinterface.jl
index 29b0ed2f3f..39e684a55f 100644
--- a/NDTensors/src/lib/BlockSparseArrays/src/blocksparsearrayinterface/blocksparsearrayinterface.jl
+++ b/NDTensors/src/lib/BlockSparseArrays/src/blocksparsearrayinterface/blocksparsearrayinterface.jl
@@ -1,12 +1,13 @@
 using BlockArrays:
   AbstractBlockVector,
   Block,
-  BlockedUnitRange,
   BlockIndex,
+  BlockVector,
+  BlockedUnitRange,
   block,
   blockcheckbounds,
-  blocks,
   blocklengths,
+  blocks,
   findblockindex
 using LinearAlgebra: Adjoint, Transpose
 using ..SparseArrayInterface: perm, iperm, nstored, sparse_zero!
@@ -19,42 +20,30 @@ function blocksparse_getindex(a::AbstractArray{<:Any,N}, I::Vararg{Int,N}) where
   return a[findblockindex.(axes(a), I)...]
 end
 
-function blocksparse_getindex(a::AbstractArray{<:Any,N}, I::Block{N}) where {N}
-  return blocksparse_getindex(a, Tuple(I)...)
-end
-function blocksparse_getindex(a::AbstractArray{<:Any,N}, I::Vararg{Block{1},N}) where {N}
-  # TODO: Avoid copy if the block isn't stored.
-  return copy(blocks(a)[Int.(I)...])
+# a[1:2, 1:2]
+# TODO: This definition means that the result of slicing a block sparse array
+# with a non-blocked unit range is blocked. We may want to change that behavior,
+# and make that explicit with `@blocked a[1:2, 1:2]`. See the discussion in
+# https://github.com/JuliaArrays/BlockArrays.jl/issues/347 and also
+# https://github.com/ITensor/ITensors.jl/issues/1336.
+function blocksparse_to_indices(a, inds, I::Tuple{UnitRange{<:Integer},Vararg{Any}})
+  bs1 = blockrange(inds[1], I[1])
+  I1 = BlockSlice(bs1, blockedunitrange_getindices(inds[1], I[1]))
+  return (I1, to_indices(a, Base.tail(inds), Base.tail(I))...)
 end
 
-# TODO: Implement as `copy(@view a[I...])`, which is then implemented
-# through `ArrayLayouts.sub_materialize`.
-using ..SparseArrayInterface: set_getindex_zero_function
-function blocksparse_getindex(
-  a::AbstractArray{<:Any,N}, I::Vararg{AbstractVector{<:Block{1}},N}
-) where {N}
-  blocks_a = blocks(a)
-  # Convert to cartesian indices of the underlying sparse array of blocks.
-  CI = map(i -> Int.(i), I)
-  subblocks_a = blocks_a[CI...]
-  subaxes = ntuple(ndims(a)) do i
-    return only(axes(axes(a, i)[I[i]]))
-  end
-  subblocks_a = set_getindex_zero_function(subblocks_a, BlockZero(subaxes))
-  return typeof(a)(subblocks_a, subaxes)
+# a[[Block(2), Block(1)], [Block(2), Block(1)]]
+function blocksparse_to_indices(a, inds, I::Tuple{Vector{<:Block{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
 
-# Slice by block and merge according to the blocking structure of the indices.
-function blocksparse_getindex(
-  a::AbstractArray{<:Any,N}, I::Vararg{AbstractBlockVector{<:Block{1}},N}
-) where {N}
-  a_sub = a[Vector.(I)...]
-  # TODO: Define `blocklengths(::AbstractBlockVector)`? Maybe make a PR
-  # to `BlockArrays`.
-  blockmergers = blockedrange.(blocklengths.(only.(axes.(I))))
-  # TODO: Need to implement this!
-  a_merged = block_merge(a_sub, blockmergers...)
-  return a_merged
+# a[BlockVector([Block(2), Block(1)], [2]), BlockVector([Block(2), Block(1)], [2])]
+# Permute and merge blocks.
+# TODO: This isn't merging blocks yet, that needs to be implemented that.
+function blocksparse_to_indices(a, inds, I::Tuple{BlockVector{<:Block{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
 
 # TODO: Need to implement this!
@@ -75,26 +64,6 @@ function blocksparse_setindex!(a::AbstractArray{<:Any,N}, value, I::BlockIndex{N
   return a
 end
 
-function blocksparse_setindex!(a::AbstractArray{<:Any,N}, value, I::Block{N}) where {N}
-  blocksparse_setindex!(a, value, Tuple(I)...)
-  return a
-end
-function blocksparse_setindex!(
-  a::AbstractArray{<:Any,N}, value, I::Vararg{Block{1},N}
-) where {N}
-  i = Int.(I)
-  @boundscheck blockcheckbounds(a, i...)
-  # TODO: Use `blocksizes(a)[i...]` when we upgrade to
-  # BlockArrays.jl v1.
-  if size(value) ≠ size(view(a, I...))
-    return throw(
-      DimensionMismatch("Trying to set a block with an array of the wrong size.")
-    )
-  end
-  blocks(a)[i...] = value
-  return a
-end
-
 function blocksparse_fill!(a::AbstractArray, value)
   for b in BlockRange(a)
     # We can't use:
@@ -119,20 +88,6 @@ function blocksparse_fill!(a::AbstractArray, value)
   return a
 end
 
-function blocksparse_view(a::AbstractArray{<:Any,N}, I::Block{N}) where {N}
-  return blocksparse_view(a, Tuple(I)...)
-end
-function blocksparse_view(a::AbstractArray{<:Any,N}, I::Vararg{Block{1},N}) where {N}
-  return SubArray(a, to_indices(a, I))
-end
-
-function blocksparse_viewblock(a::AbstractArray{<:Any,N}, I::Block{N}) where {N}
-  # TODO: Create a conversion function, say `CartesianIndex(Int.(Tuple(I)))`.
-  i = I.n
-  @boundscheck blockcheckbounds(a, i...)
-  return blocks(a)[i...]
-end
-
 function block_nstored(a::AbstractArray)
   return nstored(blocks(a))
 end
@@ -303,8 +258,12 @@ end
 # `SparseArrayInterface.sparse_storage`, which is used
 # to defined this.
 SparseArrayInterface.nstored(a::SparseSubArrayBlocks) = length(stored_indices(a))
+
+## struct SparseSubArrayBlocksStorage{Array<:SparseSubArrayBlocks}
+##   array::Array
+## end
 function SparseArrayInterface.sparse_storage(a::SparseSubArrayBlocks)
-  return error("Not implemented")
+  return map(I -> a[I], stored_indices(a))
 end
 
 function blocksparse_blocks(a::SubArray)
diff --git a/NDTensors/src/lib/BlockSparseArrays/test/test_basics.jl b/NDTensors/src/lib/BlockSparseArrays/test/test_basics.jl
index 08246c363a..2f3b8d6b92 100644
--- a/NDTensors/src/lib/BlockSparseArrays/test/test_basics.jl
+++ b/NDTensors/src/lib/BlockSparseArrays/test/test_basics.jl
@@ -23,23 +23,6 @@ using Test: @test, @test_broken, @test_throws, @testset
 include("TestBlockSparseArraysUtils.jl")
 @testset "BlockSparseArrays (eltype=$elt)" for elt in
                                                (Float32, Float64, ComplexF32, ComplexF64)
-  @testset "Broken" begin
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
-    b = @view a[[Block(2), Block(1)], [Block(2), Block(1)]]
-    @test b isa SubArray{<:Any,<:Any,<:BlockSparseArray}
-    @test_broken b[2:4, 2:4]
-
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
-    b = @views a[[Block(2), Block(1)], [Block(2), Block(1)]][Block(1, 1)]
-    @test_broken b isa SubArray{<:Any,<:Any,<:BlockSparseArray}
-
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
-    b = @views a[Block(1, 1)][1:2, 1:1]
-    @test b isa SubArray{<:Any,<:Any,<:BlockSparseArray}
-    for i in parentindices(b)
-      @test_broken i isa BlockSlice{<:BlockIndexRange{1}}
-    end
-  end
   @testset "Basics" begin
     a = BlockSparseArray{elt}([2, 3], [2, 3])
     @test a == BlockSparseArray{elt}(blockedrange([2, 3]), blockedrange([2, 3]))
@@ -114,6 +97,26 @@ include("TestBlockSparseArraysUtils.jl")
     @test size(b) == size(a)
     @test blocksize(b) == blocksize(a)
 
+    a = BlockSparseArray{elt}([2, 3], [3, 4])
+    b = @view a[[Block(2), Block(1)], [Block(2), Block(1)]]
+    c = @view b[Block(1, 1)]
+    @test iszero(a)
+    @test iszero(nstored(a))
+    @test iszero(b)
+    @test iszero(nstored(b))
+    @test iszero(c)
+    @test iszero(nstored(c))
+    a[5, 7] = 1
+    @test !iszero(a)
+    @test nstored(a) == 3 * 4
+    @test !iszero(b)
+    @test nstored(b) == 3 * 4
+    @test !iszero(c)
+    @test nstored(c) == 3 * 4
+    d = @view a[1:4, 1:6]
+    @test iszero(d)
+    @test nstored(d) == 2 * 3
+
     a = BlockSparseArray{elt}(undef, ([2, 3], [3, 4]))
     @views for b in [Block(1, 2), Block(2, 1)]
       a[b] = randn(elt, size(a[b]))
@@ -247,12 +250,23 @@ include("TestBlockSparseArraysUtils.jl")
     @views for b in [Block(1, 2), Block(2, 1)]
       a[b] = randn(elt, size(a[b]))
     end
-    b = a[2:4, 2:4]
-    @test b == Array(a)[2:4, 2:4]
-    @test size(b) == (3, 3)
-    @test blocksize(b) == (2, 2)
-    @test nstored(b) == 1 * 1 + 2 * 2
-    @test block_nstored(b) == 2
+    for b in (a[2:4, 2:4], @view(a[2:4, 2:4]))
+      @test b == Array(a)[2:4, 2:4]
+      @test size(b) == (3, 3)
+      @test blocksize(b) == (2, 2)
+      @test nstored(b) == 1 * 1 + 2 * 2
+      @test block_nstored(b) == 2
+      for f in (getindex, view)
+        @test size(f(b, Block(1, 1))) == (1, 2)
+        @test size(f(b, Block(2, 1))) == (2, 2)
+        @test size(f(b, Block(1, 2))) == (1, 1)
+        @test size(f(b, Block(2, 2))) == (2, 1)
+        @test f(b, Block(1, 1)) == a[Block(1, 1)[2:2, 2:3]]
+        @test f(b, Block(2, 1)) == a[Block(2, 1)[1:2, 2:3]]
+        @test f(b, Block(1, 2)) == a[Block(1, 2)[2:2, 1:1]]
+        @test f(b, Block(2, 2)) == a[Block(2, 2)[1:2, 1:1]]
+      end
+    end
 
     a = BlockSparseArray{elt}(undef, ([2, 3], [3, 4]))
     @views for b in [Block(1, 2), Block(2, 1)]
@@ -357,6 +371,18 @@ include("TestBlockSparseArraysUtils.jl")
     @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])
+    b = @views a[Block(2, 2)][1:2, 2:3]
+    @test b isa SubArray{<:Any,<:Any,<:BlockSparseArray}
+    for i in parentindices(b)
+      @test i isa BlockSlice{<:BlockIndexRange{1}}
+    end
+    x = randn(elt, 2, 2)
+    b .= x
+    @test a[Block(2, 2)[1:2, 2:3]] == x
+    @test a[Block(2, 2)[1:2, 2:3]] == b
+    @test block_nstored(a) == 1
+
     a = BlockSparseArray{elt}([2, 3], [2, 3])
     @views for b in [Block(1, 1), Block(2, 2)]
       a[b] = randn(elt, size(a[b]))
@@ -529,6 +555,7 @@ include("TestBlockSparseArraysUtils.jl")
       @test b isa SubArray{<:Any,<:Any,<:BlockSparseArray}
       @test block_nstored(b) == 1
       @test b[Block(1, 1)] == x
+      @test @view(b[Block(1, 1)]) isa SubArray{<:Any,<:Any,<:BlockSparseArray}
       for blck in [Block(2, 1), Block(1, 2), Block(2, 2)]
         @test iszero(b[blck])
       end
@@ -550,6 +577,32 @@ include("TestBlockSparseArraysUtils.jl")
     @test a[Block(1, 2)] == x
     @test block_nstored(a) == 1
 
+    a = BlockSparseArray{elt}([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
+    b = @view a[[Block(3), Block(2), Block(1)], [Block(3), Block(2), Block(1)]]
+    @test b isa SubArray{<:Any,<:Any,<:BlockSparseArray}
+    c = @view b[4:8, 4:8]
+    @test c isa SubArray{<:Any,<:Any,<:BlockSparseArray}
+    @test size(c) == (5, 5)
+    @test block_nstored(c) == 2
+    @test blocksize(c) == (2, 2)
+    @test blocklengths.(axes(c)) == ([2, 3], [2, 3])
+    @test size(c[Block(1, 1)]) == (2, 2)
+    @test c[Block(1, 1)] == a[Block(2, 2)[2:3, 2:3]]
+    @test size(c[Block(2, 2)]) == (3, 3)
+    @test c[Block(2, 2)] == a[Block(1, 1)[1:3, 1:3]]
+    @test size(c[Block(2, 1)]) == (3, 2)
+    @test iszero(c[Block(2, 1)])
+    @test size(c[Block(1, 2)]) == (2, 3)
+    @test iszero(c[Block(1, 2)])
+
+    x = randn(elt, 3, 3)
+    c[Block(2, 2)] = x
+    @test c[Block(2, 2)] == x
+    @test a[Block(1, 1)[1:3, 1:3]] == x
+
     a = BlockSparseArray{elt}([2, 3], [3, 4])
     b = @view a[[Block(2), Block(1)], [Block(2), Block(1)]]
     for index in parentindices(@view(b[Block(1, 1)]))
diff --git a/NDTensors/src/lib/GradedAxes/src/blockedunitrange.jl b/NDTensors/src/lib/GradedAxes/src/blockedunitrange.jl
index ac7d489d60..f420fae614 100644
--- a/NDTensors/src/lib/GradedAxes/src/blockedunitrange.jl
+++ b/NDTensors/src/lib/GradedAxes/src/blockedunitrange.jl
@@ -88,26 +88,23 @@ function blockedunitrange_getindices(
   return map(index -> a[index], indices)
 end
 
-## # TODO: Bring back this definition and remove the one in
-## # BlockSparseArraysGradedAxesExt once
-## # https://github.com/JuliaArrays/BlockArrays.jl/pull/405 is merged.
-## # TODO: Move this to a `BlockArraysExtensions` library.
-## # TODO: Make a special definition for `BlockedVector{<:Block{1}}` in order
-## # to merge blocks.
-## function blockedunitrange_getindices(
-##   a::AbstractBlockedUnitRange, indices::AbstractVector{<:Union{Block{1},BlockIndexRange{1}}}
-## )
-##   # Without converting `indices` to `Vector`,
-##   # mapping `indices` outputs a `BlockVector`
-##   # which is harder to reason about.
-##   blocks = map(index -> a[index], Vector(indices))
-##   # We pass `length.(blocks)` to `mortar` in order
-##   # to pass block labels to the axes of the output,
-##   # if they exist. This makes it so that
-##   # `only(axes(a[indices])) isa `GradedUnitRange`
-##   # if `a isa `GradedUnitRange`, for example.
-##   return mortar(blocks, length.(blocks))
-## end
+# TODO: Move this to a `BlockArraysExtensions` library.
+# TODO: Make a special definition for `BlockedVector{<:Block{1}}` in order
+# to merge blocks.
+function blockedunitrange_getindices(
+  a::AbstractBlockedUnitRange, indices::AbstractVector{<:Union{Block{1},BlockIndexRange{1}}}
+)
+  # Without converting `indices` to `Vector`,
+  # mapping `indices` outputs a `BlockVector`
+  # which is harder to reason about.
+  blocks = map(index -> a[index], Vector(indices))
+  # We pass `length.(blocks)` to `mortar` in order
+  # to pass block labels to the axes of the output,
+  # if they exist. This makes it so that
+  # `only(axes(a[indices])) isa `GradedUnitRange`
+  # if `a isa `GradedUnitRange`, for example.
+  return mortar(blocks, length.(blocks))
+end
 
 # TODO: Move this to a `BlockArraysExtensions` library.
 function blockedunitrange_getindices(a::AbstractBlockedUnitRange, indices::Block{1})
diff --git a/NDTensors/src/lib/GradedAxes/test/test_basics.jl b/NDTensors/src/lib/GradedAxes/test/test_basics.jl
index 87d6910550..e1dcc67174 100644
--- a/NDTensors/src/lib/GradedAxes/test/test_basics.jl
+++ b/NDTensors/src/lib/GradedAxes/test/test_basics.jl
@@ -9,7 +9,6 @@ using BlockArrays:
   blocklength,
   blocklengths,
   blocks
-using NDTensors.BlockSparseArrays: BlockSparseVector
 using NDTensors.GradedAxes: GradedOneTo, GradedUnitRange, blocklabels, gradedrange
 using NDTensors.LabelledNumbers: LabelledUnitRange, islabelled, label, labelled, unlabel
 using Test: @test, @test_broken, @testset
@@ -159,14 +158,7 @@ using Test: @test, @test_broken, @testset
 
   x = gradedrange(["x" => 2, "y" => 3, "z" => 4])
   a = x[[Block(3), Block(2)]]
-  # This is a BlockSparseArray since BlockArray
-  # doesn't support axes with general integer
-  # element types. That is being fixed in:
-  # https://github.com/JuliaArrays/BlockArrays.jl/pull/405
-  # TODO: Change to `BlockVector` once we update
-  # `GradedAxes` once the axes of `BlockArray`
-  # are generalized.
-  @test a isa BlockSparseVector
+  @test a isa BlockVector
   @test length(a) == 7
   @test blocklength(a) == 2
   # TODO: `BlockArrays` doesn't define `blocklengths`
@@ -187,14 +179,7 @@ using Test: @test, @test_broken, @testset
 
   x = gradedrange(["x" => 2, "y" => 3, "z" => 4])
   a = x[[Block(3)[2:3], Block(2)[2:3]]]
-  # This is a BlockSparseArray since BlockArray
-  # doesn't support axes with general integer
-  # element types. That is being fixed in:
-  # https://github.com/JuliaArrays/BlockArrays.jl/pull/405
-  # TODO: Change to `BlockVector` once we update
-  # `GradedAxes` once the axes of `BlockArray`
-  # are generalized.
-  @test a isa BlockSparseVector
+  @test a isa BlockVector
   @test length(a) == 4
   @test blocklength(a) == 2
   # TODO: `BlockArrays` doesn't define `blocklengths`
diff --git a/NDTensors/src/lib/SparseArrayInterface/src/abstractsparsearray/sparsearrayinterface.jl b/NDTensors/src/lib/SparseArrayInterface/src/abstractsparsearray/sparsearrayinterface.jl
index 897a77b960..98095225f7 100644
--- a/NDTensors/src/lib/SparseArrayInterface/src/abstractsparsearray/sparsearrayinterface.jl
+++ b/NDTensors/src/lib/SparseArrayInterface/src/abstractsparsearray/sparsearrayinterface.jl
@@ -16,3 +16,32 @@ function SparseArrayInterface.setindex_notstored!(
   iszero(value) && return a
   return error("Setting the specified unstored index is not supported.")
 end
+
+# TODO: Make this into a generic definition of all `AbstractArray`?
+# TODO: Check if this is efficient, or determine if this mapping should
+# be performed in `storage_index_to_index` and/or `index_to_storage_index`.
+function SparseArrayInterface.sparse_storage(a::SubArray{<:Any,<:Any,<:AbstractSparseArray})
+  parent_storage = sparse_storage(parent(a))
+  all_sliced_storage_indices = map(keys(parent_storage)) do I
+    return map_index(a.indices, I)
+  end
+  sliced_storage_indices = filter(!isnothing, all_sliced_storage_indices)
+  sliced_parent_storage = map(I -> parent_storage[I], keys(sliced_storage_indices))
+  return typeof(parent_storage)(sliced_storage_indices, sliced_parent_storage)
+end
+
+# TODO: Make this into a generic definition of all `AbstractArray`?
+function SparseArrayInterface.stored_indices(
+  a::PermutedDimsArray{<:Any,<:Any,<:Any,<:Any,<:AbstractSparseArray}
+)
+  return Iterators.map(
+    I -> CartesianIndex(map(i -> I[i], perm(a))), stored_indices(parent(a))
+  )
+end
+
+# TODO: Make this into a generic definition of all `AbstractArray`?
+function SparseArrayInterface.sparse_storage(
+  a::PermutedDimsArray{<:Any,<:Any,<:Any,<:Any,<:AbstractSparseArray}
+)
+  return sparse_storage(parent(a))
+end
diff --git a/NDTensors/src/lib/SparseArrayInterface/src/sparsearrayinterface/interface.jl b/NDTensors/src/lib/SparseArrayInterface/src/sparsearrayinterface/interface.jl
index ef2bc93ad2..ed5f7cb1d8 100644
--- a/NDTensors/src/lib/SparseArrayInterface/src/sparsearrayinterface/interface.jl
+++ b/NDTensors/src/lib/SparseArrayInterface/src/sparsearrayinterface/interface.jl
@@ -1,13 +1,12 @@
 # Also look into:
 # https://juliaarrays.github.io/ArrayInterface.jl/stable/sparsearrays/
 
-# Minimal interface
-# Data structure of the stored (generally nonzero) values
-sparse_storage(a::AbstractArray) = error("Not implemented")
-
-sparse_storage(a::Array) = a
+# Minimal sparse array interface.
+# Data structure of the stored (generally nonzero) values.
+# By default assume it is dense, so all values are stored.
+sparse_storage(a::AbstractArray) = a
 
-# Minimal interface
+# Minimal sparse array interface.
 # Map an index into the stored data to a CartesianIndex of the
 # outer array.
 storage_index_to_index(a::AbstractArray, I) = I
diff --git a/NDTensors/src/lib/SparseArrayInterface/src/sparsearrayinterface/wrappers.jl b/NDTensors/src/lib/SparseArrayInterface/src/sparsearrayinterface/wrappers.jl
index 351cb8f134..b4364b80ed 100644
--- a/NDTensors/src/lib/SparseArrayInterface/src/sparsearrayinterface/wrappers.jl
+++ b/NDTensors/src/lib/SparseArrayInterface/src/sparsearrayinterface/wrappers.jl
@@ -8,9 +8,6 @@ iperm(::PermutedDimsArray{<:Any,<:Any,<:Any,IP}) where {IP} = IP
 genperm(v, perm) = map(j -> v[j], perm)
 genperm(v::CartesianIndex, perm) = CartesianIndex(map(j -> Tuple(v)[j], perm))
 
-# TODO: Should the keys get permuted?
-sparse_storage(a::PermutedDimsArray) = sparse_storage(parent(a))
-
 function storage_index_to_index(a::PermutedDimsArray, I)
   return genperm(storage_index_to_index(parent(a), I), perm(a))
 end
@@ -36,18 +33,6 @@ function map_index(
   return CartesianIndex(new_index)
 end
 
-# TODO: Check if this is efficient, or determine if this mapping should
-# be performed in `storage_index_to_index` and/or `index_to_storage_index`.
-function sparse_storage(a::SubArray)
-  parent_storage = sparse_storage(parent(a))
-  all_sliced_storage_indices = map(keys(parent_storage)) do I
-    return map_index(a.indices, I)
-  end
-  sliced_storage_indices = filter(!isnothing, all_sliced_storage_indices)
-  sliced_parent_storage = map(I -> parent_storage[I], keys(sliced_storage_indices))
-  return typeof(parent_storage)(sliced_storage_indices, sliced_parent_storage)
-end
-
 function storage_index_to_index(a::SubArray, I)
   return storage_index_to_index(parent(a), I)
 end
diff --git a/NDTensors/src/lib/SparseArrayInterface/test/SparseArrayInterfaceTestUtils/Project.toml b/NDTensors/src/lib/SparseArrayInterface/test/SparseArrayInterfaceTestUtils/Project.toml
new file mode 100644
index 0000000000..9b1d5ccd25
--- /dev/null
+++ b/NDTensors/src/lib/SparseArrayInterface/test/SparseArrayInterfaceTestUtils/Project.toml
@@ -0,0 +1,2 @@
+[deps]
+NDTensors = "23ae76d9-e61a-49c4-8f12-3f1a16adf9cf"
diff --git a/NDTensors/src/lib/SparseArrayInterface/test/SparseArrayInterfaceTestUtils/SparseArrays.jl b/NDTensors/src/lib/SparseArrayInterface/test/SparseArrayInterfaceTestUtils/SparseArrays.jl
index abb531303c..05b2cb3072 100644
--- a/NDTensors/src/lib/SparseArrayInterface/test/SparseArrayInterfaceTestUtils/SparseArrays.jl
+++ b/NDTensors/src/lib/SparseArrayInterface/test/SparseArrayInterfaceTestUtils/SparseArrays.jl
@@ -65,6 +65,16 @@ function SparseArrayInterface.setindex_notstored!(
   return a
 end
 
+# TODO: Make this into a generic definition of all `AbstractArray`?
+using NDTensors.SparseArrayInterface: perm, stored_indices
+function SparseArrayInterface.stored_indices(
+  a::PermutedDimsArray{<:Any,<:Any,<:Any,<:Any,<:SparseArray}
+)
+  return Iterators.map(
+    I -> CartesianIndex(map(i -> I[i], perm(a))), stored_indices(parent(a))
+  )
+end
+
 # Empty the storage, helps with efficiency in `map!` to drop
 # zeros.
 function SparseArrayInterface.dropall!(a::SparseArray)