fixes a bug in contract when the two fermionic tensors have no common labels

include/Symmetry.hpp

@@ -31,9 +31,9 @@ namespace cytnx {
 
   /**
    * @brief fermionParity
-   * @details the parity of fermionis
+   * @details the parity of fermions
    * EVEN For bosons or an even number of fermions
-   * ODD  For an even number of fermions
+   * ODD  For an odd number of fermions
    */
   enum fermionParity : bool { EVEN = false, ODD = true };
 

src/BlockFermionicUniTensor.cpp

@@ -1356,6 +1356,8 @@ namespace cytnx {
       std::vector<cytnx_uint64> lidx(lhs_rank);
       std::vector<cytnx_uint64> ridx(rhs_rank);
       for (cytnx_int32 b = 0; b < tmp->_blocks.size(); b++) {
+        // b enumerates the ouptut blocks;
+        // idl[0] and idr[0] are the block indices in the left and right tensors
         const auto &outer_idx = tmp->_inner_to_outer_idx[b];
         std::copy_n(outer_idx.begin(), lhs_rank, lidx.begin());
         std::copy_n(outer_idx.begin() + lhs_rank, rhs_rank, ridx.begin());
@@ -1391,9 +1393,8 @@ namespace cytnx {
           cytnx_error_msg(out_block.shape() != tmp->_blocks[b].shape(),
                           "[ERROR][BlockFermionicUniTensors][contract] Mismatching shape!%s", "\n");
           tmp->_blocks[b] = out_block;
-          const bool lhs_signflip = this->_signflip[b];
-          const bool rhs_signflip = signflip_rhs[b];
-          tmp->_signflip[b] = (lhs_signflip == rhs_signflip) ? EVEN : ODD;
+          // the output sign is EVEN if the two source-block signs are EVEN-EVEN or ODD-ODD
+          tmp->_signflip[b] = (this->_signflip[idl[0]] == signflip_rhs[idr[0]]) ? EVEN : ODD;
         }
       }
 
@@ -1896,7 +1897,7 @@ namespace cytnx {
     if (ida < tmpRk) this->_rowrank--;
     if (idb < tmpRk) this->_rowrank--;
 
-    // make sure BRA_KET comes before BD_BRA to avoid supertrace
+    // make sure the BD_BRA bond comes before the BD_KET bond to avoid a supertrace
     if (this->_bonds[ida].type() == BD_KET) std::swap(ida, idb);
     // permute such that idb comes right after ida
     std::vector<cytnx_int64> perm(this->rank());

src/Bond.cpp

@@ -433,12 +433,10 @@ namespace cytnx {
         "[ERROR][get_fermion_parity] the qnum specify does not match the number of symmetries.%s",
         "\n");
 
+      // the total parity is the XOR over the parities of all symmetries of the bond
       fermionParity out = EVEN;
-      fermionParity curr = EVEN;
       for (cytnx_uint64 i = 0; i < qnum.size(); i++) {
-        out = static_cast<fermionParity>(
-          out != this->_syms[i].get_fermion_parity(
-                   qnum[i]));  // false (ODD) if the symmetries are not equal
+        out = static_cast<fermionParity>(out != this->_syms[i].get_fermion_parity(qnum[i]));
       }
 
       return out;

tests/BlockFermionicUniTensor_test.cpp

@@ -81,6 +81,85 @@ namespace cytnx {
         EXPECT_TRUE(std::abs(double(out.item().real()) - 32.0) < 1e-5);
       }
 
+      /*=====test info=====
+      describe:contraction without any common label (outer product): the output block
+               signs are the products of the two source-block signs, paired by the
+               source blocks that supply the data; contract-then-permute agrees with
+               permute-then-contract, and symmetry-forbidden blocks stay exactly zero
+      ====================*/
+      TEST_F(BlockFermionicUniTensorTest, NoCommonLabelContractCombinesSignflips) {
+        for (auto dtype : {Type.ComplexDouble, Type.Double}) {
+          const bool is_complex = (dtype == Type.ComplexDouble);
+          // permutations put pending signflips on some (but not all) blocks
+          UniTensor A0 =
+            UniTensor({B5Li, B5Lo, B5Ri, B5Ro}, {"v", "w", "x", "y"}, 2, dtype, Device.cpu);
+          random::uniform_(A0, -10.0, 10.0, 0);
+          UniTensor A = A0.permute({3, 1, 0, 2});
+          UniTensor B0 = BFUT1.astype(dtype);
+          random::uniform_(B0, -10.0, 10.0, 1);
+          UniTensor B = B0.permute({2, 0, 1});
+
+          // make sure A and B have some signflips but not on all blocks
+          bool a_flip = false, a_noflip = false, b_flip = false, b_noflip = false;
+          for (bool s : A.signflip()) (s ? a_flip : a_noflip) = true;
+          for (bool s : B.signflip()) (s ? b_flip : b_noflip) = true;
+          ASSERT_TRUE(a_flip);
+          ASSERT_TRUE(a_noflip);
+          ASSERT_TRUE(b_flip);
+          ASSERT_TRUE(b_noflip);
+
+          UniTensor out = A.contract(B);
+          EXPECT_EQ(out.uten_type(), UTenType.BlockFermionic);
+          ASSERT_EQ(out.rank(), A.rank() + B.rank());
+
+          // independent oracle: each physical (sign-applied) element of the outer
+          // product is the product of the corresponding physical elements
+          UniTensor outp = out.apply();
+          UniTensor Ap = A.apply();
+          UniTensor Bp = B.apply();
+          // the same physical tensor through the other operation order: contract the
+          // unpermuted operands first, then permute the labels of the result
+          UniTensor out2p = A0.contract(B0).permute({"y", "w", "v", "x", "c", "a", "b"}).apply();
+          ASSERT_EQ(out2p.labels(), outp.labels());
+          auto val = [is_complex](const UniTensor& ut, const std::vector<cytnx_uint64>& loc) {
+            auto p = ut.at(loc);
+            return cytnx_complex128(double(p.real()), is_complex ? double(p.imag()) : 0.0);
+          };
+          const std::vector<cytnx_uint64> ashape = A.shape();
+          const std::vector<cytnx_uint64> bshape = B.shape();
+          std::vector<cytnx_uint64> la(4), loc(7);
+          for (la[0] = 0; la[0] < ashape[0]; la[0]++)
+            for (la[1] = 0; la[1] < ashape[1]; la[1]++)
+              for (la[2] = 0; la[2] < ashape[2]; la[2]++)
+                for (la[3] = 0; la[3] < ashape[3]; la[3]++) {
+                  std::vector<cytnx_uint64> lb(3);
+                  for (lb[0] = 0; lb[0] < bshape[0]; lb[0]++)
+                    for (lb[1] = 0; lb[1] < bshape[1]; lb[1]++)
+                      for (lb[2] = 0; lb[2] < bshape[2]; lb[2]++) {
+                        std::copy(la.begin(), la.end(), loc.begin());
+                        std::copy(lb.begin(), lb.end(), loc.begin() + 4);
+                        EXPECT_EQ(out2p.at(loc).exists(), outp.at(loc).exists());
+                        if (outp.at(loc).exists()) {
+                          EXPECT_EQ(val(out2p, loc), val(outp, loc))
+                            << "dtype " << dtype << " contract/permute order dependence";
+                        }
+                        if (Ap.at(la).exists() && Bp.at(lb).exists()) {
+                          ASSERT_TRUE(outp.at(loc).exists());
+                          const cytnx_complex128 expected = val(Ap, la) * val(Bp, lb);
+                          EXPECT_LT(std::abs(val(outp, loc) - expected), 1e-14)
+                            << "dtype " << dtype << " element mismatch";
+                        } else if (outp.at(loc).exists()) {
+                          // the output tensor has all allowed blocks initialized, including those
+                          // where the flux of A equals the negative flux of B; therefore, these
+                          // blocks exist and are initialized to zero
+                          EXPECT_EQ(val(outp, loc), cytnx_complex128(0.0, 0.0))
+                            << "dtype " << dtype << " expected structural zero";
+                        }
+                      }
+                }
+        }
+      }
+
       TEST_F(BlockFermionicUniTensorTest, NormReturnsScalarTensor) {
         Tensor norm = BFUT1.Norm();
         EXPECT_TRUE(norm.is_scalar());