uniformize into similarstoragetype

fix ambiguity

more careful with storagetypes

more careful with tensoroperations

even more careful

the carefulest!

Author: lkdvos

src/tensors/abstracttensor.jl

@@ -45,11 +45,57 @@ end
 Return the type of vector that stores the data of a tensor.
 """ storagetype
 
-similarstoragetype(TT::Type{<:AbstractTensorMap}) = similarstoragetype(TT, scalartype(TT))
+# storage type determination and promotion - hooks for specializing
+# the default implementation tries to leverarge inference and `similar`
+@doc """
+    similarstoragetype(t, [T = scalartype(t)]) -> Type{<:DenseVector{T}}
+    similarstoragetype(TT, [T = scalartype(t)]) -> Type{<:DenseVector{T}}
+    similarstoragetype(A, [T = scalartype(t)]) -> Type{<:DenseVector{T}}
+    similarstoragetype(D, [T = scalartype(t)]) -> Type{<:DenseVector{T}}
+
+    similarstoragetype(T::Type{<:Number}) -> Vector{T}
+
+For a given tensor `t`, tensor type `TT <: AbstractTensorMap`, array type `A <: AbstractArray`,
+or sector dictionary type `D <: AbstractDict{<:Sector, <:AbstractMatrix}`, compute an appropriate
+storage type for tensors. Optionally, a different scalar type `T` can be supplied as well.
+
+This function determines the type of newly allocated `TensorMap`s throughout TensorKit.jl.
+It does so by leveraging type inference and calls to `Base.similar` for automatically determining
+appropriate storage types. Additionally this registers the default storage type when only a type
+`T <: Number` is provided, which is `Vector{T}`.
+""" similarstoragetype
+
+# implement in type domain
+similarstoragetype(t) = similarstoragetype(typeof(t))
+similarstoragetype(t, ::Type{T}) where {T <: Number} = similarstoragetype(typeof(t), T)
+
+# avoid infinite recursion
+similarstoragetype(X::Type) =
+    throw(ArgumentError("Cannot determine a storagetype for tensor / array type `$X`"))
+similarstoragetype(X::Type, ::Type{T}) where {T <: Number} =
+    throw(ArgumentError("Cannot determine a storagetype for tensor / array type `$X` and/or scalar type `$T`"))
+
+# implement on tensors
+similarstoragetype(::Type{TT}) where {TT <: AbstractTensorMap} = similarstoragetype(storagetype(TT))
+similarstoragetype(::Type{TT}, ::Type{T}) where {TT <: AbstractTensorMap, T <: Number} =
+    similarstoragetype(storagetype(TT), T)
+
+# implement on arrays
+similarstoragetype(::Type{A}) where {A <: DenseVector{<:Number}} = A
+Base.@assume_effects :foldable similarstoragetype(::Type{A}) where {A <: AbstractArray{<:Number}} =
+    Core.Compiler.return_type(similar, Tuple{A, Int})
+Base.@assume_effects :foldable similarstoragetype(::Type{A}, ::Type{T}) where {A <: AbstractArray, T <: Number} =
+    Core.Compiler.return_type(similar, Tuple{A, Type{T}, Int})
+
+# implement on sectordicts
+similarstoragetype(::Type{D}) where {D <: AbstractDict{<:Sector, <:AbstractMatrix}} =
+    similarstoragetype(valtype(D))
+similarstoragetype(::Type{D}, ::Type{T}) where {D <: AbstractDict{<:Sector, <:AbstractMatrix}, T <: Number} =
+    similarstoragetype(valtype(D), T)
+
+# default storage type for numbers
+similarstoragetype(::Type{T}) where {T <: Number} = Vector{T}
 
-function similarstoragetype(TT::Type{<:AbstractTensorMap}, ::Type{T}) where {T}
-    return Core.Compiler.return_type(similar, Tuple{storagetype(TT), Type{T}})
-end
 
 # tensor characteristics: space and index information
 #-----------------------------------------------------
@@ -175,7 +221,6 @@ end
 InnerProductStyle(t::AbstractTensorMap) = InnerProductStyle(typeof(t))
 storagetype(t::AbstractTensorMap) = storagetype(typeof(t))
 blocktype(t::AbstractTensorMap) = blocktype(typeof(t))
-similarstoragetype(t::AbstractTensorMap, T = scalartype(t)) = similarstoragetype(typeof(t), T)
 
 numout(t::AbstractTensorMap) = numout(typeof(t))
 numin(t::AbstractTensorMap) = numin(typeof(t))
@@ -503,17 +548,17 @@ end
 
 # 3 arguments
 Base.similar(t::AbstractTensorMap, codomain::TensorSpace, domain::TensorSpace) =
-    similar(t, similarstoragetype(t), codomain ← domain)
+    similar(t, similarstoragetype(t, scalartype(t)), codomain ← domain)
 Base.similar(t::AbstractTensorMap, ::Type{T}, codomain::TensorSpace) where {T} =
     similar(t, T, codomain ← one(codomain))
 
 # 2 arguments
 Base.similar(t::AbstractTensorMap, codomain::TensorSpace) =
-    similar(t, similarstoragetype(t), codomain ← one(codomain))
-Base.similar(t::AbstractTensorMap, V::TensorMapSpace) = similar(t, similarstoragetype(t), V)
+    similar(t, codomain ← one(codomain))
+Base.similar(t::AbstractTensorMap, V::TensorMapSpace) = similar(t, scalartype(t), V)
 Base.similar(t::AbstractTensorMap, ::Type{T}) where {T} = similar(t, T, space(t))
 # 1 argument
-Base.similar(t::AbstractTensorMap) = similar(t, similarstoragetype(t), space(t))
+Base.similar(t::AbstractTensorMap) = similar(t, scalartype(t), space(t))
 
 # generic implementation for AbstractTensorMap -> returns `TensorMap`
 function Base.similar(t::AbstractTensorMap, ::Type{TorA}, V::TensorMapSpace) where {TorA}
@@ -524,7 +569,7 @@ end
 
 # implementation in type-domain
 function Base.similar(::Type{TT}, V::TensorMapSpace) where {TT <: AbstractTensorMap}
-    TT′ = tensormaptype(spacetype(V), numout(V), numin(V), similarstoragetype(TT))
+    TT′ = tensormaptype(spacetype(V), numout(V), numin(V), similarstoragetype(TT, scalartype(TT)))
     return TT′(undef, V)
 end
 Base.similar(::Type{TT}, cod::TensorSpace, dom::TensorSpace) where {TT <: AbstractTensorMap} =

src/tensors/tensor.jl

@@ -48,22 +48,11 @@ i.e. a tensor map with only a non-trivial output space.
 """
 const Tensor{T, S, N, A} = TensorMap{T, S, N, 0, A}
 
-function tensormaptype(S::Type{<:IndexSpace}, N₁, N₂, TorA::Type)
-    A = _tensormap_storagetype(TorA)
-    A <: DenseVector || throw(ArgumentError("Cannot determine a valid storage type from argument $TorA"))
+function tensormaptype(::Type{S}, N₁, N₂, ::Type{TorA}) where {S <: IndexSpace, TorA}
+    A = similarstoragetype(TorA)
     return TensorMap{scalartype(A), S, N₁, N₂, A}
 end
 
-# hook for mapping input types to storage types -- to be implemented in extensions
-_tensormap_storagetype(::Type{A}) where {A <: DenseVector{<:Number}} = A
-_tensormap_storagetype(::Type{A}) where {A <: Array} = _tensormap_storagetype(scalartype(A))
-_tensormap_storagetype(::Type{T}) where {T <: Number} = Vector{T}
-function _tensormap_storagetype(::Type{A}) where {A <: AbstractArray}
-    PA = parenttype(A)
-    PA === A && throw(MethodError(_tensormap_storagetype, A)) # avoid infinite recursion
-    return _tensormap_storagetype(PA)
-end
-
 # Basic methods for characterising a tensor:
 #--------------------------------------------
 space(t::TensorMap) = t.space
@@ -201,9 +190,8 @@ cases.
     the specified symmetry structure, up to a tolerance `tol`.
 """
 function TensorMap(data::AbstractArray, V::TensorMapSpace; tol = sqrt(eps(real(float(eltype(data))))))
-    T = eltype(data)
-    A = _tensormap_storagetype(typeof(data))
-    return TensorMapWithStorage{T, A}(data, V; tol)
+    A = similarstoragetype(data)
+    return TensorMapWithStorage{scalartype(A), A}(data, V; tol)
 end
 TensorMap(data::AbstractArray, codom::TensorSpace, dom::TensorSpace; kwargs...) =
     TensorMap(data, codom ← dom; kwargs...)
@@ -259,7 +247,7 @@ Construct a `TensorMap` by explicitly specifying its block data.
 - `domain::ProductSpace{S, N₂}`: the domain as a `ProductSpace` of `N₂` spaces of type `S <: ElementarySpace`.
 """
 function TensorMap(data::_BlockData, V::TensorMapSpace)
-    A = _tensormap_storagetype(valtype(data))
+    A = similarstoragetype(data)
     return TensorMapWithStorage{scalartype(A), A}(data, V)
 end
 TensorMap(data::_BlockData, codom::TensorSpace, dom::TensorSpace) =

src/tensors/tensoroperations.jl

@@ -155,12 +155,16 @@ function TO.tensorcontract_type(
     ) where {N₁, N₂}
     spacetype(A) == spacetype(B) || throw(SpaceMismatch("incompatible space types"))
     I = sectortype(A)
-    M = similarstoragetype(A, sectorscalartype(I) <: Real ? TC : complex(TC))
-    MB = similarstoragetype(B, sectorscalartype(I) <: Real ? TC : complex(TC))
-    M == MB || throw(ArgumentError("incompatible storage types:\n$(M) ≠ $(MB)"))
+    TC′ = isreal(I) ? TC : complex(TC)
+    M = promote_storagetype(similarstoragetype(A, TC′), similarstoragetype(B, TC′))
     return tensormaptype(spacetype(A), N₁, N₂, M)
 end
 
+# TODO: handle actual promotion rule system
+function promote_storagetype(::Type{M₁}, ::Type{M₂}) where {M₁, M₂}
+    return M₁ === M₂ ? M₁ : throw(ArgumentError("Cannot determine storage type for combining `$M₁` and `$M₂`"))
+end
+
 function TO.tensorcontract_structure(
         A::AbstractTensorMap, pA::Index2Tuple, conjA::Bool,
         B::AbstractTensorMap, pB::Index2Tuple, conjB::Bool,