diff --git a/include/Symmetry.hpp b/include/Symmetry.hpp index 702954b73..d31fb9ccf 100644 --- a/include/Symmetry.hpp +++ b/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 }; diff --git a/src/BlockFermionicUniTensor.cpp b/src/BlockFermionicUniTensor.cpp index 11a354441..4b24c22a1 100644 --- a/src/BlockFermionicUniTensor.cpp +++ b/src/BlockFermionicUniTensor.cpp @@ -1356,6 +1356,8 @@ namespace cytnx { std::vector lidx(lhs_rank); std::vector 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 perm(this->rank()); diff --git a/src/Bond.cpp b/src/Bond.cpp index aac1769cd..e241f39fc 100644 --- a/src/Bond.cpp +++ b/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( - out != this->_syms[i].get_fermion_parity( - qnum[i])); // false (ODD) if the symmetries are not equal + out = static_cast(out != this->_syms[i].get_fermion_parity(qnum[i])); } return out; diff --git a/tests/BlockFermionicUniTensor_test.cpp b/tests/BlockFermionicUniTensor_test.cpp index 9c206535d..87d22c477 100644 --- a/tests/BlockFermionicUniTensor_test.cpp +++ b/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& loc) { + auto p = ut.at(loc); + return cytnx_complex128(double(p.real()), is_complex ? double(p.imag()) : 0.0); + }; + const std::vector ashape = A.shape(); + const std::vector bshape = B.shape(); + std::vector 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 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());