[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,8 +42,8 @@ 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]
 NDTensorsCUDAExt = "CUDA"

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
@@ -101,12 +102,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)
 
+iszerodata(t::Tensor) = data(t) isa NDTensors.Zeros
+
 #
 # Generic Tensor functions
 #

NDTensors/src/zeros/set_types.jl

@@ -0,0 +1,42 @@
+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}
+  return Zeros{P1,<:Any,P3}
+end
+function set_parameter(::Type{<:Zeros{<:Any,P2,P3}}, ::Position{1}, P1) where {P2,P3}
+  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,46 @@
+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)
+
+Base.iszero(t::Tensor) = iszero(storage(t))
+Base.iszero(st::TensorStorage) = data(st) isa Zeros
+
+function (arraytype::Type{<:Zeros})(::AllowAlias, A::Zeros)
+  return A
+end
+
+function (arraytype::Type{<:Zeros})(::NeverAlias, A::Zeros)
+  return copy(A)
+end

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/ITensors.jl

@@ -129,8 +129,13 @@ include("indexset.jl")
 #####################################
 # ITensor
 #
-include("itensor.jl")
-include("oneitensor.jl")
+include("itensor/itensor.jl")
+include("itensor/indexops.jl")
+include("itensor/oneitensor.jl")
+include("itensor/emptyitensor.jl")
+include("itensor/diagitensor.jl")
+include("itensor/randomitensor.jl")
+include("itensor/specialitensors.jl")
 include("tensor_operations/tensor_algebra.jl")
 include("tensor_operations/matrix_algebra.jl")
 include("tensor_operations/permutations.jl")
@@ -152,7 +157,7 @@ include("qn/flux.jl")
 include("qn/qn.jl")
 include("qn/qnindex.jl")
 include("qn/qnindexset.jl")
-include("qn/qnitensor.jl")
+include("itensor/qnitensor.jl")
 include("nullspace.jl")
 
 #####################################

src/imports.jl

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

src/itensor/diagitensor.jl

@@ -0,0 +1,107 @@
+#
+# Diag ITensor constructors
+#
+
+"""
+    diagITensor([ElT::Type, ]v::Vector, inds...)
+    diagitensor([ElT::Type, ]v::Vector, inds...)
+
+Make a sparse ITensor with non-zero elements only along the diagonal.
+In general, the diagonal elements will be those stored in `v` and
+the ITensor will have element type `eltype(v)`, unless specified explicitly
+by `ElT`. The storage will have `NDTensors.Diag` type.
+
+In the case when `eltype(v) isa Union{Int, Complex{Int}}`, by default it will
+be converted to `float(v)`. Note that this behavior is subject to change
+in the future.
+
+The version `diagITensor` will never output an ITensor whose storage data
+is an alias of the input vector data.
+
+The version `diagitensor` might output an ITensor whose storage data
+is an alias of the input vector data in order to minimize operations.
+"""
+function diagITensor(
+  as::AliasStyle, ElT::Type{<:Number}, v::AbstractVector{<:Number}, is::Indices
+)
+  length(v) ≠ mindim(is) && error(
+    "Length of vector for diagonal must equal minimum of the dimension of the input indices",
+  )
+  data = set_eltype(typeof(v), ElT)(as, v)
+  return itensor(Diag(data), is)
+end
+
+function diagITensor(as::AliasStyle, ElT::Type{<:Number}, v::AbstractVector{<:Number}, is...)
+  return diagITensor(as, ElT, v, indices(is...))
+end
+
+function diagITensor(as::AliasStyle, v::AbstractVector, is...)
+  return diagITensor(as, eltype(v), v, is...)
+end
+
+function diagITensor(as::AliasStyle, v::AbstractVector{<:RealOrComplex{Int}}, is...)
+  return diagITensor(as, float(eltype(v)), v, is...)
+end
+
+diagITensor(v::AbstractVector{<:Number}, is...) = diagITensor(NeverAlias(), v, is...)
+function diagITensor(ElT::Type{<:Number}, v::AbstractVector{<:Number}, is...)
+  return diagITensor(NeverAlias(), ElT, v, is...)
+end
+
+diagitensor(args...; kwargs...) = diagITensor(AllowAlias(), args...; kwargs...)
+
+# XXX TODO: explain conversion from Int
+# XXX TODO: proper conversion
+"""
+    diagITensor([ElT::Type, ]x::Number, inds...)
+    diagitensor([ElT::Type, ]x::Number, inds...)
+
+Make a sparse ITensor with non-zero elements only along the diagonal.
+In general, the diagonal elements will be set to the value `x` and
+the ITensor will have element type `eltype(x)`, unless specified explicitly
+by `ElT`. The storage will have `NDTensors.Diag` type.
+
+In the case when `x isa Union{Int, Complex{Int}}`, by default it will
+be converted to `float(x)`. Note that this behavior is subject to change
+in the future.
+"""
+function diagITensor(::AliasStyle, ElT::Type{<:Number}, x::Number, is::Indices)
+  return diagITensor(AllowAlias(), ElT, fill(eltype(x), mindim(is)), is...)
+end
+
+function diagITensor(as::AliasStyle, ElT::Type{<:Number}, x::Number, is...)
+  return diagITensor(as, ElT, x, indices(is...))
+end
+
+function diagITensor(as::AliasStyle, x::Number, is...)
+  return diagITensor(as, typeof(x), x, is...)
+end
+
+function diagITensor(as::AliasStyle, x::RealOrComplex{Int}, is...)
+  return diagITensor(as, float(typeof(x)), x, is...)
+end
+
+function diagITensor(ElT::Type{<:Number}, x::Number, is...)
+  return diagITensor(NeverAlias(), ElT, x, is...)
+end
+
+"""
+    diagITensor([::Type{ElT} = Float64, ]inds)
+    diagITensor([::Type{ElT} = Float64, ]inds::Index...)
+
+Make a sparse ITensor of element type `ElT` with only elements
+along the diagonal stored. Defaults to having `zero(T)` along
+the diagonal.
+
+The storage will have `NDTensors.Diag` type.
+"""
+function diagITensor(::Type{ElT}, is::Indices) where {ElT<:Number}
+  return diagITensor(NeverAlias(), ElT, 0, is)
+end
+
+diagITensor(::Type{ElT}, is...) where {ElT<:Number} = diagITensor(ElT, indices(is...))
+
+diagITensor(is::Indices) = diagITensor(Float64, is)
+diagITensor(is...) = diagITensor(indices(is...))
+
+diagITensor(x::Number, is...) = diagITensor(NeverAlias(), x, is...)

src/itensor/emptyitensor.jl

@@ -0,0 +1,26 @@
+#
+# EmptyStorage ITensor constructors
+#
+
+# TODO: Deprecated!
+"""
+    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.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(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(NDTensors.default_eltype(), is)
+
+emptyITensor(is...) = emptyITensor(NDTensors.default_eltype(), indices(is...))
+
+function emptyITensor(::Type{ElT}=NDTensors.default_eltype()) where {ElT<:Number}
+  return itensor(NDTensors.Zeros{ElT,1,NDTensors.default_datatype(ElT)}(()), ())
+end