fix fermionic tensor trace

Author: Jutho

src/tensors/tensoroperations.jl

@@ -163,7 +163,10 @@ const _add_kernels = (_add_trivial_kernel!, _add_abelian_kernel!, _add_general_k
 function trace!(α, tsrc::AbstractTensorMap{S}, β, tdst::AbstractTensorMap{S, N₁, N₂},
                 p1::IndexTuple{N₁}, p2::IndexTuple{N₂},
                 q1::IndexTuple{N₃}, q2::IndexTuple{N₃}) where {S, N₁, N₂, N₃}
-    # TODO: check Frobenius-Schur indicators!, and  add fermions!
+
+    if !(BraidingStyle(sectortype(S)) isa SymmetricBraiding)
+        throw(SectorMismatch("only tensors with symmetric braiding rules can be contracted; try `@planar` instead"))
+    end
     @boundscheck begin
         all(i->space(tsrc, p1[i]) == space(tdst, i), 1:N₁) ||
             throw(SpaceMismatch("trace: tsrc = $(codomain(tsrc))←$(domain(tsrc)),
@@ -203,6 +206,11 @@ function trace!(α, tsrc::AbstractTensorMap{S}, β, tdst::AbstractTensorMap{S, N
                 f2′′, g2 = split(f2′, N₂)
                 if g1 == g2
                     coeff *= dim(g1.coupled)/dim(g1.uncoupled[1])
+                    for i = 2:length(g1.uncoupled)
+                        if !(g1.isdual[i])
+                            coeff *= twist(g1.uncoupled[i])
+                        end
+                    end
                     TO._trace!(α*coeff, tsrc[f1, f2], true, tdst[f1′′, f2′′], pdata, q1, q2)
                 end
             end
@@ -250,21 +258,21 @@ function contract!(α, A::AbstractTensorMap{S}, B::AbstractTensorMap{S},
     memcost4 += dB*(!hsp(B, oindB, cindB′′)) +
                 dA*(!hsp(A, cindA′′, oindA))
 
-    # if min(memcost1, memcost2) <= min(memcost3, memcost4)
+    if min(memcost1, memcost2) <= min(memcost3, memcost4)
         if memcost1 <= memcost2
             return _contract!(α, A, B, β, C, oindA, cindA′, oindB, cindB′, p1, p2, syms)
         else
             return _contract!(α, A, B, β, C, oindA, cindA′′, oindB, cindB′′, p1, p2, syms)
         end
-    # else
-    #     p1′ = map(n->ifelse(n>N₁, n-N₁, n+N₂), p1)
-    #     p2′ = map(n->ifelse(n>N₁, n-N₁, n+N₂), p2)
-    #     if memcost3 <= memcost4
-    #         return _contract!(α, B, A, β, C, oindB, cindB′, oindA, cindA′, p1′, p2′, syms)
-    #     else
-    #         return _contract!(α, B, A, β, C, oindB, cindB′′, oindA, cindA′′, p1′, p2′, syms)
-    #     end
-    # end
+    else
+        p1′ = map(n->ifelse(n>N₁, n-N₁, n+N₂), p1)
+        p2′ = map(n->ifelse(n>N₁, n-N₁, n+N₂), p2)
+        if memcost3 <= memcost4
+            return _contract!(α, B, A, β, C, oindB, cindB′, oindA, cindA′, p1′, p2′, syms)
+        else
+            return _contract!(α, B, A, β, C, oindB, cindB′′, oindA, cindA′′, p1′, p2′, syms)
+        end
+    end
 end
 
 function _contract!(α, A::AbstractTensorMap{S}, B::AbstractTensorMap{S},

test/tensors.jl

@@ -49,13 +49,13 @@ VSU₃ = (ℂ[SU3Irrep]((0,0,0)=>3, (1,0,0)=>1),
                ℂ[SU3Irrep]((1,0,0)=>1, (2,0,0)=>1),
                ℂ[SU3Irrep]((0,0,0)=>1, (1,0,0)=>1, (1,1,0)=>1)')
 
-for V in (Vtr, Vℤ₂, Vfℤ₂, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂, VfSU₂)#, VSU₃)
+for V in (Vtr, Vℤ₂, Vfℤ₂, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂, VfSU₂, VSU₃)
     V1, V2, V3, V4, V5 = V
     @assert V3 * V4 * V2 ≿ V1' * V5' # necessary for leftorth tests
     @assert V3 * V4 ≾ V1' * V2' * V5' # necessary for rightorth tests
 end
 
-for V in (Vtr, Vℤ₂, Vℤ₃, VU₁, VCU₁, VSU₂, VSU₃)
+for V in (Vtr, Vℤ₂, Vfℤ₂, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂, VfSU₂, VSU₃)
     I = sectortype(first(V))
     Istr = TensorKit.type_repr(I)
     println("---------------------------------------")
@@ -205,11 +205,18 @@ for V in (Vtr, Vℤ₂, Vℤ₃, VU₁, VCU₁, VSU₂, VSU₃)
         t2 = permute(t, (1,2), (4,3))
         s = @constinferred tr(t2)
         @test conj(s) ≈ tr(t2')
+        if !isdual(V1)
+            t2 = twist!(t2, 1)
+        end
+        if isdual(V2)
+            t2 = twist!(t2, 2)
+        end
+        ss = tr(t2)
         @tensor s2 = t[a,b,b,a]
         @tensor t3[a,b] := t[a,c,c,b]
         @tensor s3 = t3[a,a]
-        @test s ≈ s2
-        @test s ≈ s3
+        @test ss ≈ s2
+        @test ss ≈ s3
     end
     @timedtestset "Partial trace: test self-consistency" begin
         t = Tensor(rand, ComplexF64, V1 ⊗ V2' ⊗ V3 ⊗ V2 ⊗ V1' ⊗ V3')