diff --git a/Project.toml b/Project.toml index ddcbe4141..22a4915d7 100644 --- a/Project.toml +++ b/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" diff --git a/docs/src/lib/tensors.md b/docs/src/lib/tensors.md index 3e9520cbc..6bb19dee7 100644 --- a/docs/src/lib/tensors.md +++ b/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) diff --git a/src/planar/planaroperations.jl b/src/planar/planaroperations.jl index a06f4431d..982a848f8 100644 --- a/src/planar/planaroperations.jl +++ b/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 diff --git a/src/spaces/vectorspaces.jl b/src/spaces/vectorspaces.jl index 5ef171c31..867d18959 100644 --- a/src/spaces/vectorspaces.jl +++ b/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) diff --git a/src/tensors/abstracttensor.jl b/src/tensors/abstracttensor.jl index 4b3149a4c..823280e07 100644 --- a/src/tensors/abstracttensor.jl +++ b/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 #----------------------------------------------------- """ diff --git a/src/tensors/diagonal.jl b/src/tensors/diagonal.jl index c191bb6b5..447a4b8eb 100644 --- a/src/tensors/diagonal.jl +++ b/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 ) diff --git a/src/tensors/indexmanipulations.jl b/src/tensors/indexmanipulations.jl index 0c918e85e..906ad6379 100644 --- a/src/tensors/indexmanipulations.jl +++ b/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( diff --git a/src/tensors/linalg.jl b/src/tensors/linalg.jl index d174bef71..b1fb26a09 100644 --- a/src/tensors/linalg.jl +++ b/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)...), ()) diff --git a/src/tensors/tensor.jl b/src/tensors/tensor.jl index 526f4e489..6e6cce626 100644 --- a/src/tensors/tensor.jl +++ b/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}, ...) diff --git a/src/tensors/tensoroperations.jl b/src/tensors/tensoroperations.jl index 6b01d1b2c..81894d140 100644 --- a/src/tensors/tensoroperations.jl +++ b/src/tensors/tensoroperations.jl @@ -54,8 +54,7 @@ end function TO.tensoradd_type( TC, A::AbstractTensorMap, ::Index2Tuple{N₁, N₂}, ::Bool ) where {N₁, N₂} - I = sectortype(A) - M = similarstoragetype(A, sectorscalartype(I) <: Real ? TC : complex(TC)) + M = similarstoragetype(A, promote_permute(TC, sectortype(A))) return tensormaptype(spacetype(A), N₁, N₂, M) end @@ -103,7 +102,7 @@ end VB::TensorMapSpace, pB::Index2Tuple, conjB::Bool, pAB::Index2Tuple ) - spacetype(VC) == spacetype(VA) == spacetype(VB) || throw(SectorMismatch("incompatible sector types")) + check_spacetype(VC, VA, VB) TO.tensorcontract(VA, pA, conjA, VB, pB, conjB, pAB) == VC || throw( SpaceMismatch( @@ -153,11 +152,10 @@ function TO.tensorcontract_type( B::AbstractTensorMap, ::Index2Tuple, ::Bool, ::Index2Tuple{N₁, N₂} ) where {N₁, N₂} - spacetype(A) == spacetype(B) || throw(SpaceMismatch("incompatible space types")) - I = sectortype(A) - TC′ = isreal(I) ? TC : complex(TC) + S = check_spacetype(A, B) + TC′ = promote_permute(TC, sectortype(S)) M = promote_storagetype(similarstoragetype(A, TC′), similarstoragetype(B, TC′)) - return tensormaptype(spacetype(A), N₁, N₂, M) + return tensormaptype(S, N₁, N₂, M) end # TODO: handle actual promotion rule system @@ -213,8 +211,7 @@ function trace_permute!( backend = TO.DefaultBackend() ) # some input checks - (S = spacetype(tdst)) == spacetype(tsrc) || - throw(SpaceMismatch("incompatible spacetypes")) + S = check_spacetype(tdst, tsrc) if !(BraidingStyle(sectortype(S)) isa SymmetricBraiding) throw(SectorMismatch("only tensors with symmetric braiding rules can be contracted; try `@planar` instead")) end