Fix handling of real tensors with complex scalartype

Project.toml

@@ -52,7 +52,7 @@ Random = "1"
 SafeTestsets = "0.1"
 ScopedValues = "1.3.0"
 Strided = "2"
-TensorKitSectors = "0.3.3"
+TensorKitSectors = "0.3.5"
 TensorOperations = "5.1"
 Test = "1"
 TestExtras = "0.2,0.3"

docs/src/lib/tensors.md

@@ -184,7 +184,7 @@ type `add_transform!`, for additional expert-mode options that allows for additi
 scaling, as well as the selection of a custom backend.
 
 ```@docs
-permute(::AbstractTensorMap, ::Index2Tuple{N₁,N₂}) where {N₁,N₂}
+permute(::AbstractTensorMap, ::Index2Tuple)
 braid(::AbstractTensorMap, ::Index2Tuple, ::IndexTuple)
 transpose(::AbstractTensorMap, ::Index2Tuple)
 repartition(::AbstractTensorMap, ::Int, ::Int)

src/planar/planaroperations.jl

@@ -69,8 +69,7 @@ function planartrace!(
         α::Number, β::Number,
         backend, allocator
     )
-    (S = spacetype(C)) == spacetype(A) ||
-        throw(SpaceMismatch("incompatible spacetypes"))
+    S = check_spacetype(C, A)
     if BraidingStyle(sectortype(S)) == Bosonic()
         return trace_permute!(C, A, (p₁, p₂), (q₁, q₂), α, β, backend)
     end

src/spaces/vectorspaces.jl

@@ -376,19 +376,31 @@ abstract type CompositeSpace{S <: ElementarySpace} <: VectorSpace end
 InnerProductStyle(::Type{<:CompositeSpace{S}}) where {S} = InnerProductStyle(S)
 
 """
-    spacetype(a) -> Type{S<:IndexSpace}
-    spacetype(::Type) -> Type{S<:IndexSpace}
+    spacetype(a) -> Type{S <: IndexSpace}
+    spacetype(::Type) -> Type{S <: IndexSpace}
 
-Return the type of the elementary space `S` of object `a` (e.g. a tensor). Also works in
-type domain.
+Return the type of the elementary space `S` of object `a` (e.g. a tensor).
+Also works in type domain.
 """
 spacetype(x) = spacetype(typeof(x))
-function spacetype(::Type{T}) where {T}
-    throw(MethodError(spacetype, (T,)))
-end
+spacetype(::Type{T}) where {T} = throw(MethodError(spacetype, (T,)))
 spacetype(::Type{E}) where {E <: ElementarySpace} = E
 spacetype(::Type{S}) where {E, S <: CompositeSpace{E}} = E
 
+"""
+    check_spacetype(Bool, x, y, z...) -> Bool
+    check_spacetype(x, y, z...) -> Type{<:IndexSpace}
+
+Check whether the given inputs have matching spacetypes.
+The first signature returns a `Bool` to indicate whether all spacetypes are equal,
+while the second will return the spacetype if all types are equal, and throw a [`SpaceMismatch`](@ref) if not.
+"""
+check_spacetype(::Type{Bool}, x, y, z...) = _allequal(spacetype, (x, y, z...))
+@noinline function check_spacetype(x, y, z...)
+    check_spacetype(Bool, x, y, z...) || throw(SpaceMismatch("incompatible space types"))
+    return spacetype(x)
+end
+
 # make ElementarySpace instances behave similar to ProductSpace instances
 blocksectors(V::ElementarySpace) = collect(sectors(V))
 blockdim(V::ElementarySpace, c::Sector) = dim(V, c)

src/tensors/abstracttensor.jl

@@ -102,7 +102,6 @@ similarstoragetype(::Type{D}, ::Type{T}) where {D <: AbstractDict{<:Sector, <:Ab
 # default storage type for numbers
 similarstoragetype(::Type{T}) where {T <: Number} = Vector{T}
 
-
 # tensor characteristics: space and index information
 #-----------------------------------------------------
 """

src/tensors/diagonal.jl

@@ -260,11 +260,10 @@ function TO.tensorcontract_type(
         B::DiagonalTensorMap, ::Index2Tuple{1, 1}, ::Bool,
         ::Index2Tuple{1, 1}
     )
-    M = similarstoragetype(A, TC)
-    M == similarstoragetype(B, TC) ||
-        throw(ArgumentError("incompatible storage types:\n$(M) ≠ $(similarstoragetype(B, TC))"))
-    spacetype(A) == spacetype(B) || throw(SpaceMismatch("incompatible space types"))
-    return DiagonalTensorMap{TC, spacetype(A), M}
+    S = check_spacetype(A, B)
+    TC′ = promote_permute(TC, sectortype(S))
+    M = promote_storagetype(similarstoragetype(A, TC′), similarstoragetype(B, TC′))
+    return DiagonalTensorMap{TC, S, M}
 end
 
 function TO.tensoralloc(
@@ -291,6 +290,15 @@ function Base.zero(d::DiagonalTensorMap)
     return DiagonalTensorMap(zero.(d.data), d.domain)
 end
 
+function compose_dest(A::DiagonalTensorMap, B::DiagonalTensorMap)
+    S = check_spacetype(A, B)
+    TC = TO.promote_contract(scalartype(A), scalartype(B), One)
+    M = promote_storagetype(similarstoragetype(A, TC), similarstoragetype(B, TC))
+    TTC = DiagonalTensorMap{TC, S, M}
+    structure = codomain(A) ← domain(B)
+    return TO.tensoralloc(TTC, structure, Val(false))
+end
+
 function LinearAlgebra.mul!(
         dC::DiagonalTensorMap, dA::DiagonalTensorMap, dB::DiagonalTensorMap, α::Number, β::Number
     )

src/tensors/indexmanipulations.jl

@@ -1,5 +1,28 @@
 # Index manipulations
 #---------------------
+
+# find the scalartype after applying operations: take into account fusion and/or braiding
+# might need to become Float or Complex to capture complex recoupling coefficients but don't alter precision
+for (operation, manipulation) in (
+        :flip => :sector, :twist => :braiding,
+        :transpose => :fusion, :permute => :sector, :braid => :sector,
+    )
+    promote_op = Symbol(:promote_, operation)
+    manipulation_scalartype = Symbol(manipulation, :scalartype)
+
+    @eval begin
+        $promote_op(t::AbstractTensorMap) = $promote_op(typeof(t))
+        $promote_op(::Type{T}) where {T <: AbstractTensorMap} =
+            $promote_op(scalartype(T), sectortype(T))
+        $promote_op(::Type{T}, ::Type{I}) where {T <: Number, I <: Sector} =
+            sectorscalartype(I) <: Integer ? T :
+            sectorscalartype(I) <: Real ? float(T) : complex(T)
+        # TODO: currently the manipulations all use sectorscalartype, change to:
+        # $manipulation_scalartype(I) <: Integer ? T :
+        # $manipulation_scalartype(I) <: Real ? float(T) : complex(T)
+    end
+end
+
 """
     flip(t::AbstractTensorMap, I) -> t′::AbstractTensorMap
 
@@ -15,7 +38,7 @@ Return a new tensor that is isomorphic to `t` but where the arrows on the indice
 """
 function flip(t::AbstractTensorMap, I; inv::Bool = false)
     P = flip(space(t), I)
-    t′ = similar(t, P)
+    t′ = similar(t, promote_flip(t), P)
     for (f₁, f₂) in fusiontrees(t)
         (f₁′, f₂′), factor = only(flip(f₁, f₂, I; inv))
         scale!(t′[f₁′, f₂′], t[f₁, f₂], factor)
@@ -39,53 +62,36 @@ See [`permute`](@ref) for creating a new tensor and [`add_permute!`](@ref) for a
 end
 
 """
-    permute(tsrc::AbstractTensorMap, (p₁, p₂)::Index2Tuple;
-            copy::Bool=false)
-        -> tdst::TensorMap
+    permute(tsrc::AbstractTensorMap, (p₁, p₂)::Index2Tuple; copy::Bool = false) -> tdst::TensorMap
 
 Return tensor `tdst` obtained by permuting the indices of `tsrc`.
 The codomain and domain of `tdst` correspond to the indices in `p₁` and `p₂` of `tsrc` respectively.
 
-If `copy=false`, `tdst` might share data with `tsrc` whenever possible. Otherwise, a copy is always made.
+If `copy = false`, `tdst` might share data with `tsrc` whenever possible.
+Otherwise, a copy is always made.
 
 To permute into an existing destination, see [permute!](@ref) and [`add_permute!`](@ref)
 """
-function permute(
-        t::AbstractTensorMap, (p₁, p₂)::Index2Tuple{N₁, N₂}; copy::Bool = false
-    ) where {N₁, N₂}
-    space′ = permute(space(t), (p₁, p₂))
-    # share data if possible
-    if !copy && p₁ === codomainind(t) && p₂ === domainind(t)
-        return t
-    end
-    # general case
-    @inbounds begin
-        return permute!(similar(t, space′), t, (p₁, p₂))
-    end
-end
-function permute(t::TensorMap, (p₁, p₂)::Index2Tuple{N₁, N₂}; copy::Bool = false) where {N₁, N₂}
-    space′ = permute(space(t), (p₁, p₂))
+function permute(t::AbstractTensorMap, p::Index2Tuple; copy::Bool = false)
     # share data if possible
     if !copy
-        if p₁ === codomainind(t) && p₂ === domainind(t)
+        if p == (codomainind(t), domainind(t))
             return t
-        elseif has_shared_permute(t, (p₁, p₂))
-            return TensorMap(t.data, space′)
+        elseif t isa TensorMap && has_shared_permute(t, p)
+            return TensorMap(t.data, permute(space(t), p))
         end
     end
+
     # general case
-    @inbounds begin
-        return permute!(similar(t, space′), t, (p₁, p₂))
-    end
+    tdst = similar(t, promote_permute(t), permute(space(t), p))
+    return @inbounds permute!(tdst, t, p)
 end
 function permute(t::AdjointTensorMap, (p₁, p₂)::Index2Tuple; copy::Bool = false)
     p₁′ = adjointtensorindices(t, p₂)
     p₂′ = adjointtensorindices(t, p₁)
     return adjoint(permute(adjoint(t), (p₁′, p₂′); copy))
 end
-function permute(t::AbstractTensorMap, p::IndexTuple; copy::Bool = false)
-    return permute(t, (p, ()); copy)
-end
+permute(t::AbstractTensorMap, p::IndexTuple; copy::Bool = false) = permute(t, (p, ()); copy)
 
 function has_shared_permute(t::AbstractTensorMap, (p₁, p₂)::Index2Tuple)
     return (p₁ === codomainind(t) && p₂ === domainind(t))
@@ -143,20 +149,16 @@ If `copy=false`, `tdst` might share data with `tsrc` whenever possible. Otherwis
 To braid into an existing destination, see [braid!](@ref) and [`add_braid!`](@ref)
 """
 function braid(
-        t::AbstractTensorMap, (p₁, p₂)::Index2Tuple, levels::IndexTuple; copy::Bool = false
+        t::AbstractTensorMap, p::Index2Tuple, levels::IndexTuple; copy::Bool = false
     )
-    @assert length(levels) == numind(t)
-    if BraidingStyle(sectortype(t)) isa SymmetricBraiding
-        return permute(t, (p₁, p₂); copy = copy)
-    end
-    if !copy && p₁ == codomainind(t) && p₂ == domainind(t)
-        return t
-    end
+    length(levels) == numind(t) || throw(ArgumentError("invalid levels"))
+
+    BraidingStyle(sectortype(t)) isa SymmetricBraiding && return permute(t, p; copy)
+    (!copy && p == (codomainind(t), domainind(t))) && return t
+
     # general case
-    space′ = permute(space(t), (p₁, p₂))
-    @inbounds begin
-        return braid!(similar(t, space′), t, (p₁, p₂), levels)
-    end
+    tdst = similar(t, promote_braid(t), permute(space(t), p))
+    return @inbounds braid!(tdst, t, p, levels)
 end
 # TODO: braid for `AdjointTensorMap`; think about how to map the `levels` argument.
 
@@ -196,20 +198,15 @@ If `copy=false`, `tdst` might share data with `tsrc` whenever possible. Otherwis
 To permute into an existing destination, see [permute!](@ref) and [`add_permute!`](@ref)
 """
 function LinearAlgebra.transpose(
-        t::AbstractTensorMap, (p₁, p₂)::Index2Tuple = _transpose_indices(t);
+        t::AbstractTensorMap, p::Index2Tuple = _transpose_indices(t);
         copy::Bool = false
     )
-    if sectortype(t) === Trivial
-        return permute(t, (p₁, p₂); copy = copy)
-    end
-    if !copy && p₁ == codomainind(t) && p₂ == domainind(t)
-        return t
-    end
+    sectortype(t) === Trivial && return permute(t, p; copy)
+    (!copy && p == (codomainind(t), domainind(t))) && return t
+
     # general case
-    space′ = permute(space(t), (p₁, p₂))
-    @inbounds begin
-        return transpose!(similar(t, space′), t, (p₁, p₂))
-    end
+    tdst = similar(t, promote_transpose(t), permute(space(t), p))
+    return @inbounds transpose!(tdst, t, p)
 end
 
 function LinearAlgebra.transpose(
@@ -296,6 +293,10 @@ function twist!(t::AbstractTensorMap, inds; inv::Bool = false)
         throw(ArgumentError(msg))
     end
     has_shared_twist(t, inds) && return t
+
+    (scalartype(t) <: Real && !(sectorscalartype(sectortype(t)) <: Real)) &&
+        throw(ArgumentError("No in-place `twist!` for a real tensor with complex sector type"))
+
     N₁ = numout(t)
     for (f₁, f₂) in fusiontrees(t)
         θ = prod(i -> i <= N₁ ? twist(f₁.uncoupled[i]) : twist(f₂.uncoupled[i - N₁]), inds)
@@ -317,7 +318,9 @@ See [`twist!`](@ref) for storing the result in place.
 """
 function twist(t::AbstractTensorMap, inds; inv::Bool = false, copy::Bool = false)
     !copy && has_shared_twist(t, inds) && return t
-    return twist!(Base.copy(t), inds; inv)
+    tdst = similar(t, promote_twist(t))
+    copy!(tdst, t)
+    return twist!(tdst, inds; inv)
 end
 
 # Methods which change the number of indices, implement using `Val(i)` for type inference
@@ -413,7 +416,7 @@ end
 spacecheck_transform(f, tdst::AbstractTensorMap, tsrc::AbstractTensorMap, args...) =
     spacecheck_transform(f, space(tdst), space(tsrc), args...)
 @noinline function spacecheck_transform(f, Vdst::TensorMapSpace, Vsrc::TensorMapSpace, p::Index2Tuple)
-    spacetype(Vdst) == spacetype(Vsrc) || throw(SectorMismatch("incompatible sector types"))
+    check_spacetype(Vdst, Vsrc)
     f(Vsrc, p) == Vdst ||
         throw(
         SpaceMismatch(
@@ -427,7 +430,7 @@ spacecheck_transform(f, tdst::AbstractTensorMap, tsrc::AbstractTensorMap, args..
     return nothing
 end
 @noinline function spacecheck_transform(::typeof(braid), Vdst::TensorMapSpace, Vsrc::TensorMapSpace, p::Index2Tuple, levels::IndexTuple)
-    spacetype(Vdst) == spacetype(Vsrc) || throw(SectorMismatch("incompatible sector types"))
+    check_spacetype(Vdst, Vsrc)
     braid(Vsrc, p, levels) == Vdst ||
         throw(
         SpaceMismatch(

src/tensors/linalg.jl

@@ -19,17 +19,15 @@ LinearAlgebra.normalize!(t::AbstractTensorMap, p::Real = 2) = scale!(t, inv(norm
 LinearAlgebra.normalize(t::AbstractTensorMap, p::Real = 2) = scale(t, inv(norm(t, p)))
 
 # destination allocation for matrix multiplication
+# note that we don't fall back to `tensoralloc_contract` since that needs to account for
+# permutations, which might require complex scalartypes even if the inputs are real.
 function compose_dest(A::AbstractTensorMap, B::AbstractTensorMap)
+    S = check_spacetype(A, B)
     TC = TO.promote_contract(scalartype(A), scalartype(B), One)
-    pA = (codomainind(A), domainind(A))
-    pB = (codomainind(B), domainind(B))
-    pAB = (codomainind(A), ntuple(i -> i + numout(A), numin(B)))
-    return TO.tensoralloc_contract(
-        TC,
-        A, pA, false,
-        B, pB, false,
-        pAB, Val(false)
-    )
+    M = promote_storagetype(similarstoragetype(A, TC), similarstoragetype(B, TC))
+    TTC = tensormaptype(S, numout(A), numin(B), M)
+    structure = codomain(A) ← domain(B)
+    return TO.tensoralloc(TTC, structure, Val(false))
 end
 
 """
@@ -292,8 +290,9 @@ function LinearAlgebra.rank(
         t::AbstractTensorMap;
         atol::Real = 0, rtol::Real = atol > 0 ? 0 : _default_rtol(t)
     )
-    r = 0 * dim(first(allunits(sectortype(t))))
-    dim(t) == 0 && return r
+    r = dim(t)
+    iszero(r) && return r
+    r = zero(r)
     S = MatrixAlgebraKit.svd_vals(t)
     tol = max(atol, rtol * maximum(parent(S)))
     for (c, b) in pairs(S)
@@ -325,7 +324,7 @@ end
 function LinearAlgebra.tr(t::AbstractTensorMap)
     domain(t) == codomain(t) ||
         throw(SpaceMismatch("Trace of a tensor only exist when domain == codomain"))
-    s = zero(scalartype(t))
+    s = zero(scalartype(t)) * zero(dimscalartype(sectortype(t)))
     for (c, b) in blocks(t)
         s += dim(c) * tr(b)
     end
@@ -538,8 +537,7 @@ absorb(tdst::AbstractTensorMap, tsrc::AbstractTensorMap) = absorb!(copy(tdst), t
 function absorb!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap)
     numin(tdst) == numin(tsrc) && numout(tdst) == numout(tsrc) ||
         throw(DimensionError("Incompatible number of indices for source and destination"))
-    S = spacetype(tdst)
-    S == spacetype(tsrc) || throw(SpaceMismatch("incompatible spacetypes"))
+    S = check_spacetype(tdst, tsrc)
     dom = mapreduce(infimum, ⊗, domain(tdst), domain(tsrc); init = one(S))
     cod = mapreduce(infimum, ⊗, codomain(tdst), codomain(tsrc); init = one(S))
     for (f1, f2) in fusiontrees(cod ← dom)
@@ -561,7 +559,7 @@ new `TensorMap` instance whose codomain is `codomain(t1) ⊗ codomain(t2)` and w
 is `domain(t1) ⊗ domain(t2)`.
 """
 function ⊗(A::AbstractTensorMap, B::AbstractTensorMap)
-    (S = spacetype(A)) === spacetype(B) || throw(SpaceMismatch("incompatible space types"))
+    check_spacetype(A, B)
 
     # allocate destination with correct scalartype
     pA = ((codomainind(A)..., domainind(A)...), ())

src/tensors/tensor.jl

@@ -105,6 +105,11 @@ TensorMapWithStorage{T, A}(::UndefInitializer, codomain::TensorSpace, domain::Te
     TensorMapWithStorage{T, A}(undef, codomain ← domain)
 TensorWithStorage{T, A}(::UndefInitializer, V::TensorSpace) where {T, A} = TensorMapWithStorage{T, A}(undef, V ← one(V))
 
+# Utility constructors
+# --------------------
+TensorMap(t::TensorMap) = copy(t)
+
+
 # raw data constructors
 # ---------------------
 # - dispatch starts with TensorMap{T}(::DenseVector{T}, ...)