[NDTensors] In the effort to replace EmptyStorage with an empty DataType

NDTensors/Project.toml

@@ -8,6 +8,7 @@ Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
 FLoops = "cc61a311-1640-44b5-9fba-1b764f453329"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 Folds = "41a02a25-b8f0-4f67-bc48-60067656b558"
 Functors = "d9f16b24-f501-4c13-a1f2-28368ffc5196"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
@@ -41,11 +42,11 @@ julia = "1.6"
 [weakdeps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 Metal = "dde4c033-4e86-420c-a63e-0dd931031962"
-TBLIS = "48530278-0828-4a49-9772-0f3830dfa1e9"
 Octavian = "6fd5a793-0b7e-452c-907f-f8bfe9c57db4"
+TBLIS = "48530278-0828-4a49-9772-0f3830dfa1e9"
 
 [extensions]
 NDTensorCUDA = "CUDA"
 NDTensorMetal = "Metal"
-NDTensorTBLIS = "TBLIS"
 NDTensorOctavian = "Octavian"
+NDTensorTBLIS = "TBLIS"

NDTensors/src/NDTensors.jl

@@ -4,6 +4,7 @@ using Adapt
 using Base.Threads
 using Compat
 using Dictionaries
+using FillArrays
 using FLoops
 using Folds
 using Random
@@ -99,12 +100,10 @@ include("blocksparse/combiner.jl")
 include("blocksparse/linearalgebra.jl")
 
 #####################################
-# Empty
+# Zeros
 #
-include("empty/empty.jl")
-include("empty/EmptyTensor.jl")
-include("empty/tensoralgebra/contract.jl")
-include("empty/adapt.jl")
+include("zeros/zeros.jl")
+include("zeros/set_types.jl")
 
 #####################################
 # Deprecations

NDTensors/src/exports.jl

@@ -59,11 +59,6 @@ export
   Diag,
   DiagTensor,
 
-  # empty.jl
-  EmptyStorage,
-  EmptyTensor,
-  EmptyBlockSparseTensor,
-
   # tensorstorage.jl
   data,
   TensorStorage,

NDTensors/src/tensor/tensor.jl

@@ -107,6 +107,10 @@ function Tensor(datatype::Type{<:AbstractArray}, inds::Tuple)
   return Tensor(generic_zeros(datatype, dim(inds)), inds)
 end
 
+function Tensor()
+  return Tensor(Zeros{default_eltype(),1,default_datatype(default_eltype())}(()), ())
+end
+
 ## End Tensor constructors
 
 ## Random Tensor
@@ -184,6 +188,8 @@ setstorage(T, nstore) = tensor(nstore, inds(T))
 
 setinds(T, ninds) = tensor(storage(T), ninds)
 
+isemptystorage(t::Tensor) = datatype(t) isa NDTensors.Zeros
+
 #
 # Generic Tensor functions
 #

NDTensors/src/zeros/set_types.jl

@@ -0,0 +1,44 @@
+import .SetParameters: set_parameter, nparameters, default_parameter
+
+# `SetParameters.jl` overloads.
+get_parameter(::Type{<:Zeros{P1}}, ::Position{1}) where {P1} = P1
+get_parameter(::Type{<:Zeros{<:Any,P2}}, ::Position{2}) where {P2} = P2
+get_parameter(::Type{<:Zeros{<:Any,<:Any,P3}}, ::Position{3}) where {P3} = P3
+
+# Set parameter 1
+set_parameter(::Type{<:Zeros}, ::Position{1}, P1) = Zeros{P1}
+set_parameter(::Type{<:Zeros{<:Any,P2}}, ::Position{1}, P1) where {P2} = Zeros{P1,P2}
+function set_parameter(::Type{<:Zeros{<:Any,<:Any,P3}}, ::Position{1}, P1) where {P3}
+  @assert eltype(P3) == P1
+  return Zeros{P1,<:Any,P3}
+end
+function set_parameter(::Type{<:Zeros{<:Any,P2,P3}}, ::Position{1}, P1) where {P2,P3}
+  @assert eltype(P3) == P1
+  return Zeros{P1,P2,P3}
+end
+
+# Set parameter 2
+set_parameter(::Type{<:Zeros}, ::Position{2}, P2) = Zeros{<:Any,P2}
+set_parameter(::Type{<:Zeros{P1}}, ::Position{2}, P2) where {P1} = Zeros{P1,P2}
+function set_parameter(::Type{<:Zeros{<:Any,<:Any,P3}}, ::Position{2}, P2) where {P3}
+  P1 = eltype(P3)
+  return Zeros{P1,P2,P3}
+end
+function set_parameter(::Type{<:Zeros{P1,<:Any,P3}}, ::Position{2}, P2) where {P1,P3}
+  return Zeros{P1,P2,P3}
+end
+
+# Set parameter 3
+set_parameter(::Type{<:Zeros}, ::Position{3}, P3) = Zeros{<:Any,<:Any,P3}
+set_parameter(::Type{<:Zeros{P1}}, ::Position{3}, P3) where {P1} = Zeros{P1,<:Any,P3}
+function set_parameter(::Type{<:Zeros{<:Any,P2}}, ::Position{3}, P3) where {P2}
+  P1 = eltype(P3)
+  return Zeros{P1,P2,P3}
+end
+set_parameter(::Type{<:Zeros{P1,P2}}, ::Position{3}, P3) where {P1,P2} = Zeros{P1,P2,P3}
+
+default_parameter(::Type{<:Zeros}, ::Position{1}) = Float64
+default_parameter(::Type{<:Zeros}, ::Position{2}) = 1
+default_parameter(::Type{<:Zeros}, ::Position{3}) = Vector{Float64}
+
+nparameters(::Type{<:Zeros}) = Val(3)

NDTensors/src/zeros/zeros.jl

@@ -0,0 +1,35 @@
+struct Zeros{ElT,N,DataT} <: AbstractArray{ElT,N}
+  z::FillArrays.Zeros
+  is::Tuple
+  function NDTensors.Zeros{ElT,N,DataT}(inds::Tuple) where {ElT,N,DataT}
+    @assert eltype(DataT) == ElT
+    @assert ndims(DataT) == N
+    z = FillArrays.Zeros(ElT, dim(inds))
+    return new{ElT,N,DataT}(z, inds)
+  end
+end
+
+Base.ndims(::NDTensors.Zeros{ElT,N}) where {ElT,N} = N
+ndims(::NDTensors.Zeros{ElT,N}) where {ElT,N} = N
+Base.eltype(::Zeros{ElT}) where {ElT} = ElT
+datatype(::NDTensors.Zeros{ElT,N,DataT}) where {ElT,N,DataT} = DataT
+datatype(::Type{<:NDTensors.Zeros{ElT,N,DataT}}) where {ElT,N,DataT} = DataT
+
+Base.size(zero::Zeros) = Base.size(zero.z)
+
+Base.print_array(io::IO, X::Zeros) = Base.print_array(io, X.z)
+
+data(zero::Zeros) = zero.z
+getindex(zero::Zeros) = getindex(zero.z)
+
+array(zero::Zeros) = datatype(zero)(zero.z)
+Array(zero::Zeros) = array(zero)
+dims(z::Zeros) = z.is
+copy(z::Zeros) = Zeros{eltype(z),1,datatype(z)}(dims(z))
+
+Base.convert(x::Type{T}, z::NDTensors.Zeros) where {T<:Array} = Base.convert(x, z.z)
+
+Base.getindex(a::Zeros, i) = Base.getindex(a.z, i)
+Base.sum(z::Zeros) = sum(z.z)
+LinearAlgebra.norm(z::Zeros) = norm(z.z)
+setindex!(A::NDTensors.Zeros, v, I) = setindex!(A.z, v, I)

NDTensors/test/runtests.jl

@@ -25,8 +25,6 @@ end
     "blocksparse.jl",
     "diagblocksparse.jl",
     "diag.jl",
-    "emptynumber.jl",
-    "emptystorage.jl",
     "combiner.jl",
   ]
     println("Running $filename")

src/imports.jl

@@ -114,7 +114,6 @@ import LinearAlgebra:
 using ITensors.NDTensors:
   Algorithm,
   @Algorithm_str,
-  EmptyNumber,
   _Tuple,
   _NTuple,
   blas_get_num_threads,

src/itensor.jl

@@ -173,21 +173,24 @@ A = ITensor(i,j)
 B = ITensor(ComplexF64,k,j)
 ```
 """
-function ITensor(eltype::Type{<:Number}, is::Indices)
-  return itensor(EmptyStorage(eltype), is)
+function ITensor(ElT::Type{<:Number}, is::Indices)
+  z = NDTensors.Zeros{ElT,1,NDTensors.default_datatype(ElT)}(is)
+  return ITensor(AllowAlias(), NDTensors.default_storagetype(typeof(z), is)(z), is)
 end
 
-ITensor(eltype::Type{<:Number}, is...) = ITensor(eltype, indices(is...))
+ITensor(ElT::Type{<:Number}, is...) = ITensor(ElT, indices(is...))
 
-ITensor(is...) = ITensor(EmptyNumber, is...)
+ITensor(is...) = ITensor(NDTensors.default_eltype(), is...)
 
 # To fix ambiguity with QN Index version
 # TODO: define as `emptyITensor(ElT)`
-ITensor(eltype::Type{<:Number}=EmptyNumber) = ITensor(eltype, ())
+ITensor(ElT::Type{<:Number}=NDTensors.default_eltype()) = ITensor(ElT, ())
 
 # TODO: define as `emptyITensor(ElT)`
 function ITensor(::Type{ElT}, inds::Tuple{}) where {ElT<:Number}
-  return ITensor(EmptyStorage(ElT), inds)
+  z = NDTensors.Zeros{ElT,1,NDTensors.default_datatype(ElT)}(inds)
+  store = NDTensors.default_storagetype(typeof(z), inds)(z)
+  return ITensor(store, inds)
 end
 
 """
@@ -271,25 +274,25 @@ ITensor(x::RealOrComplex{Int}, is...) = ITensor(float(x), is...)
 
 # TODO: Deprecated!
 """
-    emptyITensor([::Type{ElT} = NDTensors.EmptyNumber, ]inds)
-    emptyITensor([::Type{ElT} = NDTensors.EmptyNumber, ]inds::Index...)
+    emptyITensor([::Type{ElT} = NDTensors.default_eltype(), ]inds)
+    emptyITensor([::Type{ElT} = NDTensors.default_eltype(), ]inds::Index...)
 
-Construct an ITensor with storage type `NDTensors.EmptyStorage`, indices `inds`, and element type `ElT`. If the element type is not specified, it defaults to `NDTensors.EmptyNumber`, which represents a number type that can take on any value (for example, the type of the first value it is set to).
+Construct an ITensor with storage type `NDTensors.EmptyStorage`, indices `inds`, and element type `ElT`. If the element type is not specified, it defaults to `NDTensors.default_eltype()`, which represents a number type that can take on any value (for example, the type of the first value it is set to).
 """
 function emptyITensor(::Type{ElT}, is::Indices) where {ElT<:Number}
-  return itensor(EmptyTensor(ElT, is))
+  return itensor(NDTensors.Zeros{ElT,1,NDTensors.default_datatype(ElT)}(is), is)
 end
 
 function emptyITensor(::Type{ElT}, is...) where {ElT<:Number}
   return emptyITensor(ElT, indices(is...))
 end
 
-emptyITensor(is::Indices) = emptyITensor(EmptyNumber, is)
+emptyITensor(is::Indices) = emptyITensor(NDTensors.default_eltype(), is)
 
-emptyITensor(is...) = emptyITensor(EmptyNumber, indices(is...))
+emptyITensor(is...) = emptyITensor(NDTensors.default_eltype(), indices(is...))
 
-function emptyITensor(::Type{ElT}=EmptyNumber) where {ElT<:Number}
-  return itensor(EmptyTensor(ElT, ()))
+function emptyITensor(::Type{ElT}=NDTensors.default_eltype()) where {ElT<:Number}
+  return itensor(NDTensors.Zeros{ElT,1,NDTensors.default_datatype(ElT)}(()), ())
 end
 
 """
@@ -392,7 +395,7 @@ end
 function ITensor(eltype::Type{<:Number}, A::AbstractArray{<:Number}; kwargs...)
   return ITensor(NeverAlias(), eltype, A; kwargs...)
 end
-function ITensor(A::AbstractArray{<:Number}; kwargs...)
+function ITensor(A::AbstractArray{ElT}; kwargs...) where {ElT}
   return ITensor(NeverAlias(), eltype(A), A; kwargs...)
 end
 
@@ -842,7 +845,7 @@ size(T::ITensor) = dims(T)
 size(A::ITensor, d::Int) = size(tensor(A), d)
 
 _isemptyscalar(A::ITensor) = _isemptyscalar(tensor(A))
-_isemptyscalar(A::Tensor) = ndims(A) == 0 && isemptystorage(A) && eltype(A) === EmptyNumber
+_isemptyscalar(A::Tensor) = ndims(A) == 0 && isemptystorage(A)
 
 """
     dir(A::ITensor, i::Index)