diff --git a/NDTensors/src/NDTensors.jl b/NDTensors/src/NDTensors.jl
index 43139bbd56..f9b6dcd154 100644
--- a/NDTensors/src/NDTensors.jl
+++ b/NDTensors/src/NDTensors.jl
@@ -31,6 +31,7 @@ for lib in [
   :SparseArrayDOKs,
   :DiagonalArrays,
   :BlockSparseArrays,
+  :GradedAxes,
   :NamedDimsArrays,
   :SmallVectors,
   :SortedSets,
diff --git a/NDTensors/src/arraystorage/diagonalarray/storage/contract.jl b/NDTensors/src/arraystorage/diagonalarray/storage/contract.jl
index ac372eee03..2665141cc5 100644
--- a/NDTensors/src/arraystorage/diagonalarray/storage/contract.jl
+++ b/NDTensors/src/arraystorage/diagonalarray/storage/contract.jl
@@ -1,3 +1,6 @@
+using .SparseArrayInterface: densearray
+using .DiagonalArrays: DiagIndex, diaglength
+
 # TODO: Move to a different file.
 Unwrap.parenttype(::Type{<:DiagonalArray{<:Any,<:Any,P}}) where {P} = P
 
@@ -99,7 +102,7 @@ function contract!(
       coffset += ii * custride[i]
     end
     c = zero(eltype(C))
-    for j in 1:DiagonalArrays.diaglength(A)
+    for j in 1:diaglength(A)
       # With α == 0 && β == 1
       C[cstart + j * c_cstride + coffset] +=
         A[DiagIndex(j)] * B[bstart + j * b_cstride + boffset]
diff --git a/NDTensors/src/dims.jl b/NDTensors/src/dims.jl
index 14f8187394..e3f57ff520 100644
--- a/NDTensors/src/dims.jl
+++ b/NDTensors/src/dims.jl
@@ -1,3 +1,5 @@
+using .DiagonalArrays: DiagonalArrays
+
 export dense, dims, dim, mindim, diaglength
 
 # dim and dims are used in the Tensor interface, overload 
@@ -26,7 +28,7 @@ mindim(inds::Tuple) = minimum(dims(inds))
 
 mindim(::Tuple{}) = 1
 
-diaglength(inds::Tuple) = mindim(inds)
+DiagonalArrays.diaglength(inds::Tuple) = mindim(inds)
 
 """
     dim_to_strides(ds)
@@ -94,4 +96,3 @@ dim(T::Tensor) = dim(inds(T))
 dim(T::Tensor, i::Int) = dim(inds(T), i)
 maxdim(T::Tensor) = maxdim(inds(T))
 mindim(T::Tensor) = mindim(inds(T))
-diaglength(T::Tensor) = mindim(T)
diff --git a/NDTensors/src/lib/DiagonalArrays/README.md b/NDTensors/src/lib/DiagonalArrays/README.md
index 07cc9e40c8..1043d095f7 100644
--- a/NDTensors/src/lib/DiagonalArrays/README.md
+++ b/NDTensors/src/lib/DiagonalArrays/README.md
@@ -3,12 +3,14 @@
 A Julia `DiagonalArray` type.
 
 ````julia
-using NDTensors.DiagonalArrays: DiagonalArray, DiagonalMatrix, DiagIndex, DiagIndices, isdiagindex
+using NDTensors.DiagonalArrays:
+  DiagonalArray, DiagonalMatrix, DiagIndex, DiagIndices, diaglength, isdiagindex
 using Test
 
 function main()
   d = DiagonalMatrix([1.0, 2.0, 3.0])
   @test eltype(d) == Float64
+  @test diaglength(d) == 3
   @test size(d) == (3, 3)
   @test d[1, 1] == 1
   @test d[2, 2] == 2
@@ -17,6 +19,7 @@ function main()
 
   d = DiagonalArray([1.0, 2.0, 3.0], 3, 4, 5)
   @test eltype(d) == Float64
+  @test diaglength(d) == 3
   @test d[1, 1, 1] == 1
   @test d[2, 2, 2] == 2
   @test d[3, 3, 3] == 3
diff --git a/NDTensors/src/lib/DiagonalArrays/examples/README.jl b/NDTensors/src/lib/DiagonalArrays/examples/README.jl
index 204ac937e5..074adec5c9 100644
--- a/NDTensors/src/lib/DiagonalArrays/examples/README.jl
+++ b/NDTensors/src/lib/DiagonalArrays/examples/README.jl
@@ -3,12 +3,13 @@
 # A Julia `DiagonalArray` type.
 
 using NDTensors.DiagonalArrays:
-  DiagonalArray, DiagonalMatrix, DiagIndex, DiagIndices, isdiagindex
+  DiagonalArray, DiagonalMatrix, DiagIndex, DiagIndices, diaglength, isdiagindex
 using Test
 
 function main()
   d = DiagonalMatrix([1.0, 2.0, 3.0])
   @test eltype(d) == Float64
+  @test diaglength(d) == 3
   @test size(d) == (3, 3)
   @test d[1, 1] == 1
   @test d[2, 2] == 2
@@ -17,6 +18,7 @@ function main()
 
   d = DiagonalArray([1.0, 2.0, 3.0], 3, 4, 5)
   @test eltype(d) == Float64
+  @test diaglength(d) == 3
   @test d[1, 1, 1] == 1
   @test d[2, 2, 2] == 2
   @test d[3, 3, 3] == 3
diff --git a/NDTensors/src/lib/DiagonalArrays/src/diaginterface.jl b/NDTensors/src/lib/DiagonalArrays/src/diaginterface.jl
index de7c8ad51c..19092010a4 100644
--- a/NDTensors/src/lib/DiagonalArrays/src/diaginterface.jl
+++ b/NDTensors/src/lib/DiagonalArrays/src/diaginterface.jl
@@ -1,5 +1,11 @@
 using Compat: allequal
 
+diaglength(a::AbstractArray{<:Any,0}) = 1
+
+function diaglength(a::AbstractArray)
+  return minimum(size(a))
+end
+
 function isdiagindex(a::AbstractArray{<:Any,N}, I::CartesianIndex{N}) where {N}
   @boundscheck checkbounds(a, I)
   return allequal(Tuple(I))
@@ -16,8 +22,7 @@ function diagstride(a::AbstractArray)
 end
 
 function diagindices(a::AbstractArray)
-  diaglength = minimum(size(a))
-  maxdiag = LinearIndices(a)[CartesianIndex(ntuple(Returns(diaglength), ndims(a)))]
+  maxdiag = LinearIndices(a)[CartesianIndex(ntuple(Returns(diaglength(a)), ndims(a)))]
   return 1:diagstride(a):maxdiag
 end
 
diff --git a/NDTensors/src/lib/DiagonalArrays/test/Project.toml b/NDTensors/src/lib/DiagonalArrays/test/Project.toml
new file mode 100644
index 0000000000..9b1d5ccd25
--- /dev/null
+++ b/NDTensors/src/lib/DiagonalArrays/test/Project.toml
@@ -0,0 +1,2 @@
+[deps]
+NDTensors = "23ae76d9-e61a-49c4-8f12-3f1a16adf9cf"
diff --git a/NDTensors/src/lib/DiagonalArrays/test/runtests.jl b/NDTensors/src/lib/DiagonalArrays/test/runtests.jl
index a0064b5deb..1baef0ec29 100644
--- a/NDTensors/src/lib/DiagonalArrays/test/runtests.jl
+++ b/NDTensors/src/lib/DiagonalArrays/test/runtests.jl
@@ -1,6 +1,6 @@
-using Test
-using NDTensors.DiagonalArrays
-
+@eval module $(gensym())
+using Test: @test, @testset
+using NDTensors.DiagonalArrays: DiagonalArrays
 @testset "Test NDTensors.DiagonalArrays" begin
   @testset "README" begin
     @test include(
@@ -9,4 +9,12 @@ using NDTensors.DiagonalArrays
       ),
     ) isa Any
   end
+  @testset "Basics" begin
+    using NDTensors.DiagonalArrays: diaglength
+    a = fill(1.0, 2, 3)
+    @test diaglength(a) == 2
+    a = fill(1.0)
+    @test diaglength(a) == 1
+  end
+end
 end
diff --git a/NDTensors/src/lib/GradedAxes/src/GradedAxes.jl b/NDTensors/src/lib/GradedAxes/src/GradedAxes.jl
new file mode 100644
index 0000000000..4b0163f885
--- /dev/null
+++ b/NDTensors/src/lib/GradedAxes/src/GradedAxes.jl
@@ -0,0 +1,6 @@
+module GradedAxes
+include("groupsortperm.jl")
+include("tensor_product.jl")
+include("abstractgradedunitrange.jl")
+include("gradedunitrange.jl")
+end
diff --git a/NDTensors/src/lib/GradedAxes/src/abstractgradedunitrange.jl b/NDTensors/src/lib/GradedAxes/src/abstractgradedunitrange.jl
new file mode 100644
index 0000000000..3a1ea50b6e
--- /dev/null
+++ b/NDTensors/src/lib/GradedAxes/src/abstractgradedunitrange.jl
@@ -0,0 +1,104 @@
+using BlockArrays:
+  BlockArrays,
+  Block,
+  BlockRange,
+  BlockedUnitRange,
+  blockaxes,
+  blockedrange,
+  blockfirsts,
+  blocklasts,
+  blocklengths,
+  findblock
+using Dictionaries: Dictionary
+
+# Fuse two symmetry labels
+fuse(l1, l2) = error("Not implemented")
+
+abstract type AbstractGradedUnitRange{T,G} <: AbstractUnitRange{Int} end
+
+BlockArrays.blockedrange(a::AbstractGradedUnitRange) = error("Not implemented")
+sectors(a::AbstractGradedUnitRange) = error("Not implemented")
+scale_factor(a::AbstractGradedUnitRange) = error("Not implemented")
+
+# BlockArrays block axis interface
+BlockArrays.blockaxes(a::AbstractGradedUnitRange) = blockaxes(blockedrange(a))
+Base.getindex(a::AbstractGradedUnitRange, b::Block{1}) = blockedrange(a)[b]
+BlockArrays.blockfirsts(a::AbstractGradedUnitRange) = blockfirsts(blockedrange(a))
+BlockArrays.blocklasts(a::AbstractGradedUnitRange) = blocklasts(blockedrange(a))
+function BlockArrays.findblock(a::AbstractGradedUnitRange, k::Integer)
+  return findblock(blockedrange(a), k)
+end
+
+# Base axis interface
+Base.getindex(a::AbstractGradedUnitRange, I::Integer) = blockedrange(a)[I]
+Base.first(a::AbstractGradedUnitRange) = first(blockedrange(a))
+Base.last(a::AbstractGradedUnitRange) = last(blockedrange(a))
+Base.length(a::AbstractGradedUnitRange) = length(blockedrange(a))
+Base.step(a::AbstractGradedUnitRange) = step(blockedrange(a))
+Base.unitrange(b::AbstractGradedUnitRange) = first(b):last(b)
+
+sector(a::AbstractGradedUnitRange, b::Block{1}) = sectors(a)[only(b.n)]
+sector(a::AbstractGradedUnitRange, I::Integer) = sector(a, findblock(a, I))
+
+# Tensor product, no sorting
+function tensor_product(a1::AbstractGradedUnitRange, a2::AbstractGradedUnitRange)
+  a = tensor_product(blockedrange(a1), blockedrange(a2))
+  sectors_a = map(Iterators.product(sectors(a1), sectors(a2))) do (l1, l2)
+    return fuse(scale_factor(a1) * l1, scale_factor(a2) * l2)
+  end
+  return gradedrange(a, vec(sectors_a))
+end
+
+function Base.show(io::IO, mimetype::MIME"text/plain", a::AbstractGradedUnitRange)
+  show(io, mimetype, sectors(a))
+  println(io)
+  println(io, "Scale factor = ", scale_factor(a))
+  return show(io, mimetype, blockedrange(a))
+end
+
+function blockmerge(a::AbstractGradedUnitRange, grouped_perm::Vector{Vector{Int}})
+  merged_sectors = map(group -> sector(a, Block(first(group))), grouped_perm)
+  lengths = blocklengths(a)
+  merged_lengths = map(group -> sum(@view(lengths[group])), grouped_perm)
+  return gradedrange(merged_sectors, merged_lengths)
+end
+
+# Sort and merge by the grade of the blocks.
+function blockmergesort(a::AbstractGradedUnitRange)
+  grouped_perm = blockmergesortperm(a)
+  return blockmerge(a, grouped_perm)
+end
+
+# Get the permutation for sorting, then group by common elements.
+# groupsortperm([2, 1, 2, 3]) == [[2], [1, 3], [4]]
+function blockmergesortperm(a::AbstractGradedUnitRange)
+  return groupsortperm(sectors(a))
+end
+
+function sub_axis(a::AbstractGradedUnitRange, blocks)
+  a_sub = sub_axis(blockedrange(a), blocks)
+  sectors_sub = map(b -> sector(a, b), Indices(blocks))
+  return AbstractGradedUnitRange(a_sub, sectors_sub)
+end
+
+function fuse(a1::AbstractGradedUnitRange, a2::AbstractGradedUnitRange)
+  a = tensor_product(a1, a2)
+  return blockmergesort(a)
+end
+
+## TODO: Add this back.
+## # Slicing
+## ## using BlockArrays: BlockRange, _BlockedUnitRange
+## Base.@propagate_inbounds function Base.getindex(
+##   b::AbstractGradedUnitRange, KR::BlockRange{1}
+## )
+##   cs = blocklasts(b)
+##   isempty(KR) && return _BlockedUnitRange(1, cs[1:0])
+##   K, J = first(KR), last(KR)
+##   k, j = Integer(K), Integer(J)
+##   bax = blockaxes(b, 1)
+##   @boundscheck K in bax || throw(BlockBoundsError(b, K))
+##   @boundscheck J in bax || throw(BlockBoundsError(b, J))
+##   K == first(bax) && return _BlockedUnitRange(first(b), cs[k:j])
+##   return _BlockedUnitRange(cs[k - 1] + 1, cs[k:j])
+## end
diff --git a/NDTensors/src/lib/GradedAxes/src/gradedunitrange.jl b/NDTensors/src/lib/GradedAxes/src/gradedunitrange.jl
new file mode 100644
index 0000000000..ce5578964c
--- /dev/null
+++ b/NDTensors/src/lib/GradedAxes/src/gradedunitrange.jl
@@ -0,0 +1,23 @@
+using BlockArrays: BlockArrays, BlockedUnitRange, blockedrange
+
+struct GradedUnitRange{T,G,S} <: AbstractGradedUnitRange{T,G}
+  blockedrange::BlockedUnitRange{T}
+  sectors::Vector{G}
+  scale_factor::S
+end
+
+BlockArrays.blockedrange(s::GradedUnitRange) = s.blockedrange
+sectors(s::GradedUnitRange) = s.sectors
+scale_factor(s::GradedUnitRange) = s.scale_factor
+
+function gradedrange(sectors::Vector, blocklengths::Vector{Int}, scale_factor=1)
+  return GradedUnitRange(blockedrange(blocklengths), sectors, scale_factor)
+end
+
+function gradedrange(sectors_lengths::Vector{<:Pair{<:Any,Int}}, scale_factor=1)
+  return gradedrange(first.(sectors_lengths), last.(sectors_lengths), scale_factor)
+end
+
+function gradedrange(a::BlockedUnitRange, sectors::Vector, scale_factor=1)
+  return GradedUnitRange(a, sectors, scale_factor)
+end
diff --git a/NDTensors/src/lib/GradedAxes/src/groupsortperm.jl b/NDTensors/src/lib/GradedAxes/src/groupsortperm.jl
new file mode 100644
index 0000000000..4980fdb37b
--- /dev/null
+++ b/NDTensors/src/lib/GradedAxes/src/groupsortperm.jl
@@ -0,0 +1,15 @@
+using BlockArrays: BlockVector
+using SplitApplyCombine: groupcount
+
+function groupsorted(v)
+  return groupcount(identity, v)
+end
+
+# Get the permutation for sorting, then group by common elements.
+# groupsortperm([2, 1, 2, 3]) == [[2], [1, 3], [4]]
+function groupsortperm(v)
+  perm = sortperm(v)
+  v_sorted = @view v[perm]
+  group_lengths = groupsorted(v_sorted)
+  return blocks(BlockVector(perm, collect(group_lengths)))
+end
diff --git a/NDTensors/src/lib/GradedAxes/src/tensor_product.jl b/NDTensors/src/lib/GradedAxes/src/tensor_product.jl
new file mode 100644
index 0000000000..6d903cbeeb
--- /dev/null
+++ b/NDTensors/src/lib/GradedAxes/src/tensor_product.jl
@@ -0,0 +1,14 @@
+using BlockArrays: BlockedUnitRange, blocks
+
+# https://en.wikipedia.org/wiki/Tensor_product
+# https://github.com/KeitaNakamura/Tensorial.jl
+tensor_product(a1, a2, a3, as...) = foldl(tensor_product, (a1, a2, a3, as...))
+tensor_product(a1, a2) = error("Not implemented for $(typeof(a1)) and $(typeof(a2)).")
+
+function tensor_product(a1::Base.OneTo, a2::Base.OneTo)
+  return Base.OneTo(length(a1) * length(a2))
+end
+
+function tensor_product(a1::BlockedUnitRange, a2::BlockedUnitRange)
+  return blockedrange(prod.(length, vec(collect(Iterators.product(blocks.((a1, a2))...)))))
+end
diff --git a/NDTensors/src/lib/GradedAxes/test/Project.toml b/NDTensors/src/lib/GradedAxes/test/Project.toml
new file mode 100644
index 0000000000..9b1d5ccd25
--- /dev/null
+++ b/NDTensors/src/lib/GradedAxes/test/Project.toml
@@ -0,0 +1,2 @@
+[deps]
+NDTensors = "23ae76d9-e61a-49c4-8f12-3f1a16adf9cf"
diff --git a/NDTensors/src/lib/GradedAxes/test/runtests.jl b/NDTensors/src/lib/GradedAxes/test/runtests.jl
new file mode 100644
index 0000000000..354d54aa63
--- /dev/null
+++ b/NDTensors/src/lib/GradedAxes/test/runtests.jl
@@ -0,0 +1,77 @@
+@eval module $(gensym())
+using NDTensors.BlockArrays: Block, blocklength, blocklengths, findblock
+using NDTensors.GradedAxes:
+  GradedAxes,
+  blockmerge,
+  blockmergesortperm,
+  fuse,
+  gradedrange,
+  sector,
+  sectors,
+  tensor_product
+using Test: @test
+
+struct U1
+  dim::Int
+end
+Base.isless(l1::U1, l2::U1) = isless(l1.dim, l2.dim)
+Base.:*(c::Int, l::U1) = U1(c * l.dim)
+GradedAxes.fuse(l1::U1, l2::U1) = U1(l1.dim + l2.dim)
+
+a = gradedrange([U1(0), U1(1)], [2, 3])
+@test a isa GradedAxes.GradedUnitRange
+@test a == gradedrange([U1(0) => 2, U1(1) => 3])
+@test length(a) == 5
+@test a == 1:5
+@test a[Block(1)] == 1:2
+@test a[Block(2)] == 3:5
+@test blocklength(a) == 2 # Number of sectors
+@test blocklengths(a) == [2, 3]
+# TODO: Maybe rename to `labels`, `label`.
+@test sectors(a) == [U1(0), U1(1)]
+@test sector(a, Block(1)) == U1(0)
+@test sector(a, Block(2)) == U1(1)
+@test findblock(a, 1) == Block(1)
+@test findblock(a, 2) == Block(1)
+@test findblock(a, 3) == Block(2)
+@test findblock(a, 4) == Block(2)
+@test findblock(a, 5) == Block(2)
+@test sector(a, 1) == U1(0)
+@test sector(a, 2) == U1(0)
+@test sector(a, 3) == U1(1)
+@test sector(a, 4) == U1(1)
+@test sector(a, 5) == U1(1)
+
+# Naive tensor product, no sorting and merging
+a2 = tensor_product(a, a)
+@test a2 isa GradedAxes.GradedUnitRange
+@test a2 == gradedrange([U1(0) => 4, U1(1) => 6, U1(1) => 6, U1(2) => 9])
+@test length(a2) == 25
+@test a2 == 1:25
+@test blocklength(a2) == 4
+@test blocklengths(a2) == [4, 6, 6, 9]
+@test sectors(a2) == [U1(0), U1(1), U1(1), U1(2)]
+@test sector(a2, Block(1)) == U1(0)
+@test sector(a2, Block(2)) == U1(1)
+@test sector(a2, Block(3)) == U1(1)
+@test sector(a2, Block(4)) == U1(2)
+
+# Fusion tensor product, with sorting and merging
+a2 = fuse(a, a)
+@test a2 isa GradedAxes.GradedUnitRange
+@test a2 == gradedrange([U1(0) => 4, U1(1) => 12, U1(2) => 9])
+@test length(a2) == 25
+@test a2 == 1:25
+@test blocklength(a2) == 3
+@test blocklengths(a2) == [4, 12, 9]
+@test sectors(a2) == [U1(0), U1(1), U1(2)]
+@test sector(a2, Block(1)) == U1(0)
+@test sector(a2, Block(2)) == U1(1)
+@test sector(a2, Block(3)) == U1(2)
+
+# The partitioned permutation needed to sort
+# and merge an unsorted graded space
+perm_a = blockmergesortperm(tensor_product(a, a))
+@test perm_a == [[1], [2, 3], [4]]
+@test blockmerge(tensor_product(a, a), perm_a) == fuse(a, a)
+end
diff --git a/NDTensors/src/lib/SparseArrayInterface/src/SparseArrayInterface.jl b/NDTensors/src/lib/SparseArrayInterface/src/SparseArrayInterface.jl
index 5b76dbe4b4..9468eaf17c 100644
--- a/NDTensors/src/lib/SparseArrayInterface/src/SparseArrayInterface.jl
+++ b/NDTensors/src/lib/SparseArrayInterface/src/SparseArrayInterface.jl
@@ -1,4 +1,5 @@
 module SparseArrayInterface
+include("densearray.jl")
 include("interface.jl")
 include("interface_optional.jl")
 include("indexing.jl")
diff --git a/NDTensors/src/lib/SparseArrayInterface/src/densearray.jl b/NDTensors/src/lib/SparseArrayInterface/src/densearray.jl
new file mode 100644
index 0000000000..9b0f738432
--- /dev/null
+++ b/NDTensors/src/lib/SparseArrayInterface/src/densearray.jl
@@ -0,0 +1,12 @@
+# Generic functionality for converting to a
+# dense array, trying to preserve information
+# about the array (such as which device it is on).
+# TODO: Maybe call `densecopy`?
+# TODO: Make sure this actually preserves the device,
+# maybe use `NDTensors.Unwrap.unwrap_type`.
+function densearray(a::AbstractArray)
+  # TODO: `set_ndims(unwrap_type(a), ndims(a))(a)`
+  # Maybe define `densetype(a) = set_ndims(unwrap_type(a), ndims(a))`.
+  # Or could use `unspecify_parameters(unwrap_type(a))(a)`.
+  return Array(a)
+end
diff --git a/NDTensors/src/lib/SparseArrayInterface/test/runtests.jl b/NDTensors/src/lib/SparseArrayInterface/test/runtests.jl
index f00698db26..2ae116cf59 100644
--- a/NDTensors/src/lib/SparseArrayInterface/test/runtests.jl
+++ b/NDTensors/src/lib/SparseArrayInterface/test/runtests.jl
@@ -36,6 +36,10 @@ using Test: @test, @testset, @test_broken, @test_throws
     @test iszero(a)
     @test iszero(SparseArrayInterface.nstored(a))
 
+    a_dense = SparseArrayInterface.densearray(a)
+    @test a_dense == a
+    @test a_dense isa Array{elt,ndims(a)}
+
     a = SparseArray{elt}(2, 3)
     fill!(a, 2)
     @test size(a) == (2, 3)
@@ -288,6 +292,10 @@ using Test: @test, @testset, @test_broken, @test_throws
     @test a[1, 1] == 11
     @test a[2, 2] == 22
 
+    a_dense = SparseArrayInterface.densearray(a)
+    @test a_dense == a
+    @test a_dense isa Array{elt,ndims(a)}
+
     b = similar(a)
     @test b isa DiagonalArray
     @test size(b) == (2, 3)
diff --git a/NDTensors/test/arraytensor/Project.toml b/NDTensors/test/arraytensor/Project.toml
new file mode 100644
index 0000000000..2238172387
--- /dev/null
+++ b/NDTensors/test/arraytensor/Project.toml
@@ -0,0 +1,3 @@
+[deps]
+BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
+NDTensors = "23ae76d9-e61a-49c4-8f12-3f1a16adf9cf"
diff --git a/NDTensors/test/arraytensor/array.jl b/NDTensors/test/arraytensor/array.jl
index 98be7d2ed9..95c44925a3 100644
--- a/NDTensors/test/arraytensor/array.jl
+++ b/NDTensors/test/arraytensor/array.jl
@@ -1,9 +1,8 @@
-using NDTensors
-using LinearAlgebra
-using Test
-
-using NDTensors: storage, storagetype
-
+@eval module $(gensym())
+using LinearAlgebra: svd
+using NDTensors: array, contract, inds, storage, storagetype, tensor
+using Random: randn!
+using Test: @test, @testset, @test_broken
 @testset "Tensor wrapping Array" begin
   is1 = (2, 3)
   D1 = randn(is1)
@@ -49,3 +48,4 @@ using NDTensors: storage, storagetype
   D12 = contract(D1, (1, -1), D2, (-1, 2))
   @test D12 ≈ Array(T12)
 end
+end
diff --git a/NDTensors/test/arraytensor/arraytensor.jl b/NDTensors/test/arraytensor/arraytensor.jl
deleted file mode 100644
index a71944d01c..0000000000
--- a/NDTensors/test/arraytensor/arraytensor.jl
+++ /dev/null
@@ -1,8 +0,0 @@
-using NDTensors
-using LinearAlgebra
-using Test
-
-@testset "Tensor wrapping AbstractArrays" begin
-  include("array.jl")
-  include("blocksparsearray.jl")
-end
diff --git a/NDTensors/test/arraytensor/diagonalarray.jl b/NDTensors/test/arraytensor/diagonalarray.jl
index 70c0029dfc..d5d73b88db 100644
--- a/NDTensors/test/arraytensor/diagonalarray.jl
+++ b/NDTensors/test/arraytensor/diagonalarray.jl
@@ -1,10 +1,8 @@
-using NDTensors
-using NDTensors.DiagonalArrays
-using LinearAlgebra
-using Test
-
-using NDTensors: storage, storagetype
-
+@eval module $(gensym())
+using NDTensors: contract, tensor
+using NDTensors.SparseArrayInterface: densearray
+using NDTensors.DiagonalArrays: DiagonalArray
+using Test: @test, @testset
 @testset "Tensor wrapping DiagonalArray" begin
   D = DiagonalArray(randn(3), 3, 4, 5)
   Dᵈ = densearray(D)
@@ -24,3 +22,4 @@ using NDTensors: storage, storagetype
   @test contract(Dᵗ, (-1, -2, -3), Aᵗ, (-1, -2, -3)) ≈
     contract(Dᵈᵗ, (-1, -2, -3), Aᵗ, (-1, -2, -3))
 end
+end
diff --git a/NDTensors/test/arraytensor/runtests.jl b/NDTensors/test/arraytensor/runtests.jl
new file mode 100644
index 0000000000..a423524d05
--- /dev/null
+++ b/NDTensors/test/arraytensor/runtests.jl
@@ -0,0 +1,8 @@
+@eval module $(gensym())
+using Test: @testset
+@testset "Tensor wrapping AbstractArrays $(f)" for f in [
+  "array.jl", "blocksparsearray.jl", "diagonalarray.jl"
+]
+  include(f)
+end
+end
diff --git a/NDTensors/test/lib/Project.toml b/NDTensors/test/lib/Project.toml
index 6c95bd52e2..b86a01c764 100644
--- a/NDTensors/test/lib/Project.toml
+++ b/NDTensors/test/lib/Project.toml
@@ -6,3 +6,4 @@ Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
 NDTensors = "23ae76d9-e61a-49c4-8f12-3f1a16adf9cf"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
diff --git a/NDTensors/test/lib/runtests.jl b/NDTensors/test/lib/runtests.jl
index 8f0fbbdfc8..976f0e601d 100644
--- a/NDTensors/test/lib/runtests.jl
+++ b/NDTensors/test/lib/runtests.jl
@@ -6,6 +6,7 @@ using Test: @testset
   "BlockSparseArrays",
   "BroadcastMapConversion",
   "DiagonalArrays",
+  "GradedAxes",
   "NamedDimsArrays",
   "SetParameters",
   "SmallVectors",
diff --git a/NDTensors/test/runtests.jl b/NDTensors/test/runtests.jl
index fad15eb5f0..1e8beed07b 100644
--- a/NDTensors/test/runtests.jl
+++ b/NDTensors/test/runtests.jl
@@ -19,17 +19,17 @@ end
 
 @safetestset "NDTensors" begin
   @testset "$filename" for filename in [
+    "arraytensor/runtests.jl",
+    "ITensors/runtests.jl",
     "lib/runtests.jl",
-    "linearalgebra.jl",
-    "dense.jl",
     "blocksparse.jl",
-    "diagblocksparse.jl",
+    "combiner.jl",
+    "dense.jl",
     "diag.jl",
+    "diagblocksparse.jl",
     "emptynumber.jl",
     "emptystorage.jl",
-    "combiner.jl",
-    "arraytensor/arraytensor.jl",
-    "ITensors/runtests.jl",
+    "linearalgebra.jl",
   ]
     println("Running $filename")
     include(filename)