From 723a0585d9c537f575cc1f35460569ec05f6058d Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Thu, 28 Sep 2023 10:38:29 -0400
Subject: [PATCH 1/5] [NDTensor] Allow general objects are storage backends for
 Tensor

---
 NDTensors/src/NDTensors.jl                |  7 ++
 NDTensors/src/arraytensor/array.jl        | 64 ++++++++++++++++
 NDTensors/src/arraytensor/arraytensor.jl  | 90 +++++++++++++++++++++++
 NDTensors/src/tensor/tensor.jl            | 13 ++--
 NDTensors/test/arraytensor/arraytensor.jl | 53 +++++++++++++
 NDTensors/test/runtests.jl                |  1 +
 src/itensor.jl                            |  4 +-
 7 files changed, 223 insertions(+), 9 deletions(-)
 create mode 100644 NDTensors/src/arraytensor/array.jl
 create mode 100644 NDTensors/src/arraytensor/arraytensor.jl
 create mode 100644 NDTensors/test/arraytensor/arraytensor.jl

diff --git a/NDTensors/src/NDTensors.jl b/NDTensors/src/NDTensors.jl
index ae138d11fc..d2ffcafb59 100644
--- a/NDTensors/src/NDTensors.jl
+++ b/NDTensors/src/NDTensors.jl
@@ -109,6 +109,13 @@ include("empty/EmptyTensor.jl")
 include("empty/tensoralgebra/contract.jl")
 include("empty/adapt.jl")
 
+#####################################
+# Array Tensor (experimental)
+# TODO: Move to `Experimental` module.
+#
+include("arraytensor/arraytensor.jl")
+include("arraytensor/array.jl")
+
 #####################################
 # Deprecations
 #
diff --git a/NDTensors/src/arraytensor/array.jl b/NDTensors/src/arraytensor/array.jl
new file mode 100644
index 0000000000..ba0dbab601
--- /dev/null
+++ b/NDTensors/src/arraytensor/array.jl
@@ -0,0 +1,64 @@
+# Combiner
+promote_rule(::Type{<:Combiner}, arraytype::Type{<:MatrixOrArrayLike}) = arraytype
+
+# Generic AbstractArray code
+function contract(
+  array1::MatrixOrArrayLike,
+  labels1,
+  array2::MatrixOrArrayLike,
+  labels2,
+  labelsR=contract_labels(labels1, labels2),
+)
+  output_array = contraction_output(array1, labels1, array2, labels2, labelsR)
+  contract!(output_array, labelsR, array1, labels1, array2, labels2)
+  return output_array
+end
+
+function contraction_output(array1::MatrixOrArrayLike, array2::MatrixOrArrayLike, indsR)
+  arraytypeR = contraction_output_type(typeof(array1), typeof(array2), indsR)
+  return NDTensors.similar(arraytypeR, indsR)
+end
+
+function contraction_output_type(
+  arraytype1::Type{<:MatrixOrArrayLike}, arraytype2::Type{<:MatrixOrArrayLike}, inds
+)
+  return similartype(promote_type(arraytype1, arraytype2), inds)
+end
+
+function contraction_output(
+  array1::MatrixOrArrayLike,
+  labelsarray1,
+  array2::MatrixOrArrayLike,
+  labelsarray2,
+  labelsoutput_array,
+)
+  # TODO: Maybe use `axes` here to be more generic, for example for BlockArrays?
+  indsoutput_array = contract_inds(
+    size(array1), labelsarray1, size(array2), labelsarray2, labelsoutput_array
+  )
+  output_array = contraction_output(array1, array2, indsoutput_array)
+  return output_array
+end
+
+# Required interface for specific AbstractArray types
+function contract!(
+  arrayR::MatrixOrArrayLike,
+  labelsR,
+  array1::MatrixOrArrayLike,
+  labels1,
+  array2::MatrixOrArrayLike,
+  labels2,
+)
+  props = ContractionProperties(labels1, labels2, labelsR)
+  compute_contraction_properties!(props, array1, array2, arrayR)
+  # TODO: Change this to just `contract!`, or maybe `contract_ttgt!`?
+  _contract!(arrayR, array1, array2, props)
+  return arrayR
+end
+
+function permutedims!(
+  output_array::MatrixOrArrayLike, array::MatrixOrArrayLike, perm, f::Function
+)
+  @strided output_array .= f.(output_array, permutedims(array, perm))
+  return output_array
+end
diff --git a/NDTensors/src/arraytensor/arraytensor.jl b/NDTensors/src/arraytensor/arraytensor.jl
new file mode 100644
index 0000000000..76d7aa7dae
--- /dev/null
+++ b/NDTensors/src/arraytensor/arraytensor.jl
@@ -0,0 +1,90 @@
+# Used for dispatch to distinguish from Tensors wrapping TensorStorage.
+# Remove once TensorStorage is removed.
+const ArrayLike{T,N} = Union{Array{T,N},ReshapedArray{T,N},SubArray{T,N}}
+const MatrixLike{T} = Union{
+  Matrix{T},
+  Transpose{T},
+  Adjoint{T},
+  Symmetric{T},
+  Hermitian{T},
+  UpperTriangular{T},
+  LowerTriangular{T},
+  UnitUpperTriangular{T},
+  UnitLowerTriangular{T},
+  Diagonal{T},
+}
+const MatrixOrArrayLike{T} = Union{MatrixLike{T},ArrayLike{T}}
+
+const ArrayLikeTensor{T,N,S,I} = Tensor{T,N,S,I} where {S<:ArrayLike{T,N}}
+const MatrixLikeTensor{T,S,I} = Tensor{T,2,S,I} where {S<:MatrixLike{T}}
+const MatrixOrArrayLikeTensor{T,S,I} = Tensor{T,N,S,I} where {N,S<:MatrixOrArrayLike{T}}
+
+function getindex(tensor::MatrixOrArrayLikeTensor, I::Integer...)
+  return storage(tensor)[I...]
+end
+
+function setindex!(tensor::MatrixOrArrayLikeTensor, v, I::Integer...)
+  storage(tensor)[I...] = v
+  return tensor
+end
+
+function contraction_output(
+  tensor1::MatrixOrArrayLikeTensor, tensor2::MatrixOrArrayLikeTensor, indsR
+)
+  tensortypeR = contraction_output_type(typeof(tensor1), typeof(tensor2), indsR)
+  return NDTensors.similar(tensortypeR, indsR)
+end
+
+function contract!(
+  tensorR::MatrixOrArrayLikeTensor,
+  labelsR,
+  tensor1::MatrixOrArrayLikeTensor,
+  labels1,
+  tensor2::MatrixOrArrayLikeTensor,
+  labels2,
+)
+  contract!(storage(tensorR), labelsR, storage(tensor1), labels1, storage(tensor2), labels2)
+  return tensorR
+end
+
+function permutedims!(
+  output_tensor::MatrixOrArrayLikeTensor, tensor::MatrixOrArrayLikeTensor, perm, f::Function
+)
+  permutedims!(storage(output_tensor), storage(tensor), perm, f)
+  return output_tensor
+end
+
+# Linear algebra (matrix algebra)
+Base.adjoint(tens::MatrixLikeTensor) = tensor(adjoint(storage(tens)), reverse(inds(tens)))
+
+function LinearAlgebra.mul!(C::MatrixLikeTensor, A::MatrixLikeTensor, B::MatrixLikeTensor)
+  mul!(storage(C), storage(A), storage(B))
+  return C
+end
+
+function LinearAlgebra.svd(tens::MatrixLikeTensor)
+  F = svd(storage(tens))
+  U, S, V = F.U, F.S, F.Vt
+  i, j = inds(tens)
+  # TODO: Make this more general with a `similar_ind` function,
+  # so the dimension can be determined from the length of `S`.
+  min_ij = dim(i) ≤ dim(j) ? i : j
+  α = sim(min_ij) # similar_ind(i, space(S))
+  β = sim(min_ij) # similar_ind(i, space(S))
+  Utensor = tensor(U, (i, α))
+  # TODO: Remove conversion to `Matrix` to make more general.
+  # Used for now to avoid introducing wrapper types.
+  Stensor = tensor(Diagonal(S), (α, β))
+  Vtensor = tensor(V, (β, j))
+  return Utensor, Stensor, Vtensor, Spectrum(nothing, 0.0)
+end
+
+array(tensor::MatrixOrArrayLikeTensor) = storage(tensor)
+
+# Combiner
+function contraction_output(
+  tensor1::MatrixOrArrayLikeTensor, tensor2::CombinerTensor, indsR
+)
+  tensortypeR = contraction_output_type(typeof(tensor1), typeof(tensor2), indsR)
+  return NDTensors.similar(tensortypeR, indsR)
+end
diff --git a/NDTensors/src/tensor/tensor.jl b/NDTensors/src/tensor/tensor.jl
index 84699222f1..4db6714cd4 100644
--- a/NDTensors/src/tensor/tensor.jl
+++ b/NDTensors/src/tensor/tensor.jl
@@ -1,11 +1,10 @@
-
 """
 Tensor{StoreT,IndsT}
 
 A plain old tensor (with order independent
 interface and no assumption of labels)
 """
-struct Tensor{ElT,N,StoreT<:TensorStorage,IndsT} <: AbstractArray{ElT,N}
+struct Tensor{ElT,N,StoreT,IndsT} <: AbstractArray{ElT,N}
   storage::StoreT
   inds::IndsT
 
@@ -21,8 +20,8 @@ struct Tensor{ElT,N,StoreT<:TensorStorage,IndsT} <: AbstractArray{ElT,N}
   and tensor(store::TensorStorage, inds) constructors.
   """
   function Tensor{ElT,N,StoreT,IndsT}(
-    ::AllowAlias, storage::TensorStorage, inds::Tuple
-  ) where {ElT,N,StoreT<:TensorStorage,IndsT}
+    ::AllowAlias, storage, inds::Tuple
+  ) where {ElT,N,StoreT,IndsT}
     @assert ElT == eltype(StoreT)
     @assert length(inds) == N
     return new{ElT,N,StoreT,IndsT}(storage, inds)
@@ -64,7 +63,7 @@ The Tensor holds a copy of the storage data.
 
 The indices `inds` will be converted to a `Tuple`.
 """
-function Tensor(as::AliasStyle, storage::TensorStorage, inds::Tuple)
+function Tensor(as::AliasStyle, storage, inds::Tuple)
   return Tensor{eltype(storage),length(inds),typeof(storage),typeof(inds)}(
     as, storage, inds
   )
@@ -74,12 +73,12 @@ end
 # already (like a Vector). In the future this may be lifted
 # to allow for very large tensor orders in which case Tuple
 # operations may become too slow.
-function Tensor(as::AliasStyle, storage::TensorStorage, inds)
+function Tensor(as::AliasStyle, storage, inds)
   return Tensor(as, storage, Tuple(inds))
 end
 
 tensor(args...; kwargs...) = Tensor(AllowAlias(), args...; kwargs...)
-Tensor(storage::TensorStorage, inds::Tuple) = Tensor(NeverAlias(), storage, inds)
+Tensor(storage, inds::Tuple) = Tensor(NeverAlias(), storage, inds)
 
 function Tensor(eltype::Type, inds::Tuple)
   return Tensor(AllowAlias(), default_storagetype(eltype, inds)(dim(inds)), inds)
diff --git a/NDTensors/test/arraytensor/arraytensor.jl b/NDTensors/test/arraytensor/arraytensor.jl
new file mode 100644
index 0000000000..50f5e316ea
--- /dev/null
+++ b/NDTensors/test/arraytensor/arraytensor.jl
@@ -0,0 +1,53 @@
+using NDTensors
+using LinearAlgebra
+using Test
+
+using NDTensors: storage, storagetype
+
+@testset "Tensor wrapping Array" begin
+  is1 = (2, 3)
+  D1 = randn(is1)
+
+  is2 = (3, 4)
+  D2 = randn(is2)
+
+  T1 = tensor(D1, is1)
+  T2 = tensor(D2, is2)
+
+  @test T1[1, 1] == D1[1, 1]
+
+  x = rand()
+  T1[1, 1] = x
+
+  @test T1[1, 1] == x
+  @test array(T1) == D1
+  @test storagetype(T1) <: Matrix{Float64}
+  @test storage(T1) == D1
+  @test eltype(T1) == eltype(D1)
+  @test inds(T1) == is1
+
+  R = T1 * T2
+  @test storagetype(R) <: Matrix{Float64}
+  @test Array(R) ≈ Array(T1) * Array(T2)
+
+  T1r = randn!(similar(T1))
+  @test Array(T1r + T1) ≈ Array(T1r) + Array(T1)
+  @test Array(permutedims(T1, (2, 1))) ≈ permutedims(Array(T1), (2, 1))
+
+  U, S, V = svd(T1)
+  @test U * S * V ≈ T1
+
+  T12 = contract(T1, (1, -1), T2, (-1, 2))
+  @test T12 ≈ T1 * T2
+
+  D12 = contract(D1, (1, -1), D2, (-1, 2))
+  @test D12 ≈ Array(T12)
+
+  ## IT1 = itensor(T1)
+  ## IT2 = itensor(T2)
+  ## IR = IT1 * IT2
+  ## IX = IT1 + IT1
+
+  ## IU, IS, IV = svd(IT1, i)
+  ## @test IU * IS * IV ≈ IT1
+end
diff --git a/NDTensors/test/runtests.jl b/NDTensors/test/runtests.jl
index 90aeeea118..cdd0d6cf35 100644
--- a/NDTensors/test/runtests.jl
+++ b/NDTensors/test/runtests.jl
@@ -28,6 +28,7 @@ end
     "emptynumber.jl",
     "emptystorage.jl",
     "combiner.jl",
+    "arraytensor/arraytensor.jl",
   ]
     println("Running $filename")
     include(filename)
diff --git a/src/itensor.jl b/src/itensor.jl
index 19286f065c..c26786e5e2 100644
--- a/src/itensor.jl
+++ b/src/itensor.jl
@@ -82,7 +82,7 @@ NDTensors.Dense{Float64,Array{Float64,1}}
 ## Accessor Functions, Index Functions and Operations
 mutable struct ITensor
   tensor::Tensor
-  function ITensor(::AllowAlias, T::Tensor{<:Any,<:Any,<:TensorStorage,<:Tuple})
+  function ITensor(::AllowAlias, T::Tensor{<:Any,<:Any,<:Any,<:Tuple})
     @debug_check begin
       is = inds(T)
       if !allunique(is)
@@ -761,7 +761,7 @@ end
 
 Return a view of the TensorStorage of the ITensor.
 """
-storage(T::ITensor)::TensorStorage = storage(tensor(T))
+storage(T::ITensor) = storage(tensor(T))
 
 storagetype(x::ITensor) = storagetype(tensor(x))
 

From b6efb313d7e9392b4f0b46e92f87c87567c54752 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Thu, 28 Sep 2023 10:54:17 -0400
Subject: [PATCH 2/5] Rename type union defining non-TensorStorage storage

---
 NDTensors/src/arraytensor/array.jl       | 24 ++++++------
 NDTensors/src/arraytensor/arraytensor.jl | 48 ++++++++++++++----------
 2 files changed, 42 insertions(+), 30 deletions(-)

diff --git a/NDTensors/src/arraytensor/array.jl b/NDTensors/src/arraytensor/array.jl
index ba0dbab601..40e89c5181 100644
--- a/NDTensors/src/arraytensor/array.jl
+++ b/NDTensors/src/arraytensor/array.jl
@@ -1,11 +1,11 @@
 # Combiner
-promote_rule(::Type{<:Combiner}, arraytype::Type{<:MatrixOrArrayLike}) = arraytype
+promote_rule(::Type{<:Combiner}, arraytype::Type{<:MatrixOrArrayStorage}) = arraytype
 
 # Generic AbstractArray code
 function contract(
-  array1::MatrixOrArrayLike,
+  array1::MatrixOrArrayStorage,
   labels1,
-  array2::MatrixOrArrayLike,
+  array2::MatrixOrArrayStorage,
   labels2,
   labelsR=contract_labels(labels1, labels2),
 )
@@ -14,21 +14,23 @@ function contract(
   return output_array
 end
 
-function contraction_output(array1::MatrixOrArrayLike, array2::MatrixOrArrayLike, indsR)
+function contraction_output(
+  array1::MatrixOrArrayStorage, array2::MatrixOrArrayStorage, indsR
+)
   arraytypeR = contraction_output_type(typeof(array1), typeof(array2), indsR)
   return NDTensors.similar(arraytypeR, indsR)
 end
 
 function contraction_output_type(
-  arraytype1::Type{<:MatrixOrArrayLike}, arraytype2::Type{<:MatrixOrArrayLike}, inds
+  arraytype1::Type{<:MatrixOrArrayStorage}, arraytype2::Type{<:MatrixOrArrayStorage}, inds
 )
   return similartype(promote_type(arraytype1, arraytype2), inds)
 end
 
 function contraction_output(
-  array1::MatrixOrArrayLike,
+  array1::MatrixOrArrayStorage,
   labelsarray1,
-  array2::MatrixOrArrayLike,
+  array2::MatrixOrArrayStorage,
   labelsarray2,
   labelsoutput_array,
 )
@@ -42,11 +44,11 @@ end
 
 # Required interface for specific AbstractArray types
 function contract!(
-  arrayR::MatrixOrArrayLike,
+  arrayR::MatrixOrArrayStorage,
   labelsR,
-  array1::MatrixOrArrayLike,
+  array1::MatrixOrArrayStorage,
   labels1,
-  array2::MatrixOrArrayLike,
+  array2::MatrixOrArrayStorage,
   labels2,
 )
   props = ContractionProperties(labels1, labels2, labelsR)
@@ -57,7 +59,7 @@ function contract!(
 end
 
 function permutedims!(
-  output_array::MatrixOrArrayLike, array::MatrixOrArrayLike, perm, f::Function
+  output_array::MatrixOrArrayStorage, array::MatrixOrArrayStorage, perm, f::Function
 )
   @strided output_array .= f.(output_array, permutedims(array, perm))
   return output_array
diff --git a/NDTensors/src/arraytensor/arraytensor.jl b/NDTensors/src/arraytensor/arraytensor.jl
index 76d7aa7dae..c88a5b8c98 100644
--- a/NDTensors/src/arraytensor/arraytensor.jl
+++ b/NDTensors/src/arraytensor/arraytensor.jl
@@ -1,8 +1,10 @@
 # Used for dispatch to distinguish from Tensors wrapping TensorStorage.
 # Remove once TensorStorage is removed.
-const ArrayLike{T,N} = Union{Array{T,N},ReshapedArray{T,N},SubArray{T,N}}
-const MatrixLike{T} = Union{
-  Matrix{T},
+const ArrayStorage{T,N} = Union{
+  Array{T,N},ReshapedArray{T,N},SubArray{T,N},PermutedDimsArray{T,N},StridedView{T,N}
+}
+const MatrixStorage{T} = Union{
+  ArrayStorage{T,2},
   Transpose{T},
   Adjoint{T},
   Symmetric{T},
@@ -13,34 +15,35 @@ const MatrixLike{T} = Union{
   UnitLowerTriangular{T},
   Diagonal{T},
 }
-const MatrixOrArrayLike{T} = Union{MatrixLike{T},ArrayLike{T}}
+const MatrixOrArrayStorage{T} = Union{MatrixStorage{T},ArrayStorage{T}}
 
-const ArrayLikeTensor{T,N,S,I} = Tensor{T,N,S,I} where {S<:ArrayLike{T,N}}
-const MatrixLikeTensor{T,S,I} = Tensor{T,2,S,I} where {S<:MatrixLike{T}}
-const MatrixOrArrayLikeTensor{T,S,I} = Tensor{T,N,S,I} where {N,S<:MatrixOrArrayLike{T}}
+const ArrayStorageTensor{T,N,S,I} = Tensor{T,N,S,I} where {S<:ArrayStorage{T,N}}
+const MatrixStorageTensor{T,S,I} = Tensor{T,2,S,I} where {S<:MatrixStorage{T}}
+const MatrixOrArrayStorageTensor{T,S,I} =
+  Tensor{T,N,S,I} where {N,S<:MatrixOrArrayStorage{T}}
 
-function getindex(tensor::MatrixOrArrayLikeTensor, I::Integer...)
+function getindex(tensor::MatrixOrArrayStorageTensor, I::Integer...)
   return storage(tensor)[I...]
 end
 
-function setindex!(tensor::MatrixOrArrayLikeTensor, v, I::Integer...)
+function setindex!(tensor::MatrixOrArrayStorageTensor, v, I::Integer...)
   storage(tensor)[I...] = v
   return tensor
 end
 
 function contraction_output(
-  tensor1::MatrixOrArrayLikeTensor, tensor2::MatrixOrArrayLikeTensor, indsR
+  tensor1::MatrixOrArrayStorageTensor, tensor2::MatrixOrArrayStorageTensor, indsR
 )
   tensortypeR = contraction_output_type(typeof(tensor1), typeof(tensor2), indsR)
   return NDTensors.similar(tensortypeR, indsR)
 end
 
 function contract!(
-  tensorR::MatrixOrArrayLikeTensor,
+  tensorR::MatrixOrArrayStorageTensor,
   labelsR,
-  tensor1::MatrixOrArrayLikeTensor,
+  tensor1::MatrixOrArrayStorageTensor,
   labels1,
-  tensor2::MatrixOrArrayLikeTensor,
+  tensor2::MatrixOrArrayStorageTensor,
   labels2,
 )
   contract!(storage(tensorR), labelsR, storage(tensor1), labels1, storage(tensor2), labels2)
@@ -48,21 +51,28 @@ function contract!(
 end
 
 function permutedims!(
-  output_tensor::MatrixOrArrayLikeTensor, tensor::MatrixOrArrayLikeTensor, perm, f::Function
+  output_tensor::MatrixOrArrayStorageTensor,
+  tensor::MatrixOrArrayStorageTensor,
+  perm,
+  f::Function,
 )
   permutedims!(storage(output_tensor), storage(tensor), perm, f)
   return output_tensor
 end
 
 # Linear algebra (matrix algebra)
-Base.adjoint(tens::MatrixLikeTensor) = tensor(adjoint(storage(tens)), reverse(inds(tens)))
+function Base.adjoint(tens::MatrixStorageTensor)
+  return tensor(adjoint(storage(tens)), reverse(inds(tens)))
+end
 
-function LinearAlgebra.mul!(C::MatrixLikeTensor, A::MatrixLikeTensor, B::MatrixLikeTensor)
+function LinearAlgebra.mul!(
+  C::MatrixStorageTensor, A::MatrixStorageTensor, B::MatrixStorageTensor
+)
   mul!(storage(C), storage(A), storage(B))
   return C
 end
 
-function LinearAlgebra.svd(tens::MatrixLikeTensor)
+function LinearAlgebra.svd(tens::MatrixStorageTensor)
   F = svd(storage(tens))
   U, S, V = F.U, F.S, F.Vt
   i, j = inds(tens)
@@ -79,11 +89,11 @@ function LinearAlgebra.svd(tens::MatrixLikeTensor)
   return Utensor, Stensor, Vtensor, Spectrum(nothing, 0.0)
 end
 
-array(tensor::MatrixOrArrayLikeTensor) = storage(tensor)
+array(tensor::MatrixOrArrayStorageTensor) = storage(tensor)
 
 # Combiner
 function contraction_output(
-  tensor1::MatrixOrArrayLikeTensor, tensor2::CombinerTensor, indsR
+  tensor1::MatrixOrArrayStorageTensor, tensor2::CombinerTensor, indsR
 )
   tensortypeR = contraction_output_type(typeof(tensor1), typeof(tensor2), indsR)
   return NDTensors.similar(tensortypeR, indsR)

From 7d256ace92976250257a41a254cfb6a3b643e557 Mon Sep 17 00:00:00 2001
From: Matt Fishman <mtfishman@users.noreply.github.com>
Date: Thu, 28 Sep 2023 11:15:12 -0400
Subject: [PATCH 3/5] Update comment

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

diff --git a/NDTensors/src/arraytensor/arraytensor.jl b/NDTensors/src/arraytensor/arraytensor.jl
index c88a5b8c98..67b51cef93 100644
--- a/NDTensors/src/arraytensor/arraytensor.jl
+++ b/NDTensors/src/arraytensor/arraytensor.jl
@@ -82,7 +82,7 @@ function LinearAlgebra.svd(tens::MatrixStorageTensor)
   α = sim(min_ij) # similar_ind(i, space(S))
   β = sim(min_ij) # similar_ind(i, space(S))
   Utensor = tensor(U, (i, α))
-  # TODO: Remove conversion to `Matrix` to make more general.
+  # TODO: Remove conversion to `Diagonal` to make more general, or make a generic `Diagonal` concept that works for `BlockSparseArray`.
   # Used for now to avoid introducing wrapper types.
   Stensor = tensor(Diagonal(S), (α, β))
   Vtensor = tensor(V, (β, j))

From d5a84095d287d6b4f9122eaa5a190ef04b8856bd Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Thu, 28 Sep 2023 11:43:57 -0400
Subject: [PATCH 4/5] Add back some type constraints in Tensor constructors

---
 NDTensors/src/arraytensor/arraytensor.jl                 | 8 ++++++++
 NDTensors/src/tensor/tensor.jl                           | 4 ++--
 NDTensors/src/tensoralgebra/generic_tensor_operations.jl | 2 +-
 3 files changed, 11 insertions(+), 3 deletions(-)

diff --git a/NDTensors/src/arraytensor/arraytensor.jl b/NDTensors/src/arraytensor/arraytensor.jl
index c88a5b8c98..26dc85f437 100644
--- a/NDTensors/src/arraytensor/arraytensor.jl
+++ b/NDTensors/src/arraytensor/arraytensor.jl
@@ -22,6 +22,14 @@ const MatrixStorageTensor{T,S,I} = Tensor{T,2,S,I} where {S<:MatrixStorage{T}}
 const MatrixOrArrayStorageTensor{T,S,I} =
   Tensor{T,N,S,I} where {N,S<:MatrixOrArrayStorage{T}}
 
+Tensor(storage::MatrixOrArrayStorageTensor, inds::Tuple) = Tensor(NeverAlias(), storage, inds)
+
+function Tensor(as::AliasStyle, storage::MatrixOrArrayStorage, inds::Tuple)
+  return Tensor{eltype(storage),length(inds),typeof(storage),typeof(inds)}(
+    as, storage, inds
+  )
+end
+
 function getindex(tensor::MatrixOrArrayStorageTensor, I::Integer...)
   return storage(tensor)[I...]
 end
diff --git a/NDTensors/src/tensor/tensor.jl b/NDTensors/src/tensor/tensor.jl
index 4db6714cd4..7e4af3f403 100644
--- a/NDTensors/src/tensor/tensor.jl
+++ b/NDTensors/src/tensor/tensor.jl
@@ -63,7 +63,7 @@ The Tensor holds a copy of the storage data.
 
 The indices `inds` will be converted to a `Tuple`.
 """
-function Tensor(as::AliasStyle, storage, inds::Tuple)
+function Tensor(as::AliasStyle, storage::TensorStorage, inds::Tuple)
   return Tensor{eltype(storage),length(inds),typeof(storage),typeof(inds)}(
     as, storage, inds
   )
@@ -78,7 +78,7 @@ function Tensor(as::AliasStyle, storage, inds)
 end
 
 tensor(args...; kwargs...) = Tensor(AllowAlias(), args...; kwargs...)
-Tensor(storage, inds::Tuple) = Tensor(NeverAlias(), storage, inds)
+Tensor(storage::TensorStorage, inds::Tuple) = Tensor(NeverAlias(), storage, inds)
 
 function Tensor(eltype::Type, inds::Tuple)
   return Tensor(AllowAlias(), default_storagetype(eltype, inds)(dim(inds)), inds)
diff --git a/NDTensors/src/tensoralgebra/generic_tensor_operations.jl b/NDTensors/src/tensoralgebra/generic_tensor_operations.jl
index 170b472bd8..db0a740e22 100644
--- a/NDTensors/src/tensoralgebra/generic_tensor_operations.jl
+++ b/NDTensors/src/tensoralgebra/generic_tensor_operations.jl
@@ -18,7 +18,7 @@ end
 function permutedims!(output_tensor::Tensor, tensor::Tensor, perm, f::Function=(r, t) -> t)
   Base.checkdims_perm(output_tensor, tensor, perm)
   error(
-    "`perutedims!(output_tensor::Tensor, tensor::Tensor, perm, f::Function=(r, t) -> t)` not implemented for `typeof(output_tensor) = $(typeof(output_tensor))`, `typeof(tensor) = $(typeof(tensor))`, `perm = $perm`, and `f = $f`.",
+    "`permutedims!(output_tensor::Tensor, tensor::Tensor, perm, f::Function=(r, t) -> t)` not implemented for `typeof(output_tensor) = $(typeof(output_tensor))`, `typeof(tensor) = $(typeof(tensor))`, `perm = $perm`, and `f = $f`.",
   )
   return output_tensor
 end

From fefae9d4aa604c5ee1ec3339047250b1906e0df1 Mon Sep 17 00:00:00 2001
From: Matt Fishman <mtfishman@users.noreply.github.com>
Date: Thu, 28 Sep 2023 12:01:48 -0400
Subject: [PATCH 5/5] Format

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
---
 NDTensors/src/arraytensor/arraytensor.jl | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/NDTensors/src/arraytensor/arraytensor.jl b/NDTensors/src/arraytensor/arraytensor.jl
index abcc2a7bbf..d91b4a6d01 100644
--- a/NDTensors/src/arraytensor/arraytensor.jl
+++ b/NDTensors/src/arraytensor/arraytensor.jl
@@ -22,7 +22,9 @@ const MatrixStorageTensor{T,S,I} = Tensor{T,2,S,I} where {S<:MatrixStorage{T}}
 const MatrixOrArrayStorageTensor{T,S,I} =
   Tensor{T,N,S,I} where {N,S<:MatrixOrArrayStorage{T}}
 
-Tensor(storage::MatrixOrArrayStorageTensor, inds::Tuple) = Tensor(NeverAlias(), storage, inds)
+function Tensor(storage::MatrixOrArrayStorageTensor, inds::Tuple)
+  return Tensor(NeverAlias(), storage, inds)
+end
 
 function Tensor(as::AliasStyle, storage::MatrixOrArrayStorage, inds::Tuple)
   return Tensor{eltype(storage),length(inds),typeof(storage),typeof(inds)}(