Account for possible rounding error in beta-scaling with complex arit…

…hmetic. Fixes #30.
devinamatthews · Jun 8, 2020 · 3e4c4b8 · 3e4c4b8
1 parent 5ccacba
commit 3e4c4b8
Show file tree

Hide file tree

Showing 2 changed files with 24 additions and 6 deletions.
diff --git a/test/3t/contract.cxx b/test/3t/contract.cxx
@@ -50,7 +50,7 @@ REPLICATED_TEMPLATED_TEST_CASE(contract, R, T, all_types)
     TENSOR_INFO(C);
 
     auto idx_AB = intersection(idx_A, idx_B);
-    auto neps = prod(select_from(A.lengths(), idx_A, idx_AB))*prod(C.lengths());
+    auto neps = (prod(select_from(A.lengths(), idx_A, idx_AB))+1)*prod(C.lengths());
 
     impl = BLAS_BASED;
     D.reset(C);
@@ -103,7 +103,7 @@ REPLICATED_TEMPLATED_TEST_CASE(dpd_contract, R, T, all_types)
         unsigned irrep_AC = A.irrep()^irrep_AB;
         unsigned irrep_BC = B.irrep()^irrep_AB;
 
-        neps += size_AB[irrep_AB]*
+        neps += (size_AB[irrep_AB]+1)*
                 size_AC[irrep_AC]*
                 size_BC[irrep_BC];
     }
@@ -145,7 +145,7 @@ REPLICATED_TEMPLATED_TEST_CASE(indexed_contract, R, T, all_types)
     INDEXED_TENSOR_INFO(C);
 
     auto idx_AB = intersection(idx_A, idx_B);
-    auto neps = prod(select_from(A.lengths(), idx_A, idx_AB))*prod(C.lengths());
+    auto neps = (prod(select_from(A.lengths(), idx_A, idx_AB))+1)*prod(C.lengths());
 
     dpd_impl = dpd_impl_t::BLOCKED;
     D.reset(C);
@@ -189,7 +189,7 @@ REPLICATED_TEMPLATED_TEST_CASE(indexed_dpd_contract, R, T, all_types)
         unsigned irrep_AC = A.irrep()^irrep_AB;
         unsigned irrep_BC = B.irrep()^irrep_AB;
 
-        neps += size_AB[irrep_AB]*
+        neps += (size_AB[irrep_AB]+1)*
                 size_AC[irrep_AC]*
                 size_BC[irrep_BC];
     }

diff --git a/test/3t/mult.cxx b/test/3t/mult.cxx
@@ -55,7 +55,7 @@ REPLICATED_TEMPLATED_TEST_CASE(mult, R, T, all_types)
     TENSOR_INFO(C);
 
     auto idx_AB = exclusion(intersection(idx_A, idx_B), idx_C);
-    auto neps = prod(select_from(A.lengths(), idx_A, idx_AB))*prod(C.lengths());
+    auto neps = (prod(select_from(A.lengths(), idx_A, idx_AB))+1)*prod(C.lengths());
 
     impl = REFERENCE;
     D.reset(C);
@@ -116,6 +116,15 @@ REPLICATED_TEMPLATED_TEST_CASE(dpd_mult, R, T, all_types)
                 size_AC[irrep_AC]*
                 size_BC[irrep_BC];
     }
+    for (unsigned irrep_AC = 0;irrep_AC < nirrep;irrep_AC++)
+    for (unsigned irrep_BC = 0;irrep_BC < nirrep;irrep_BC++)
+    {
+        unsigned irrep_ABC = irrep_AC^irrep_BC^C.irrep();
+
+        neps += size_ABC[irrep_ABC]*
+                size_AC[irrep_AC]*
+                size_BC[irrep_BC];
+    }
 
     dpd_impl = dpd_impl_t::BLOCKED;
     D.reset(C);
@@ -148,7 +157,7 @@ REPLICATED_TEMPLATED_TEST_CASE(indexed_mult, R, T, all_types)
     INDEXED_TENSOR_INFO(C);
 
     auto idx_AB = exclusion(intersection(idx_A, idx_B), idx_C);
-    auto neps = prod(select_from(A.lengths(), idx_A, idx_AB))*prod(C.lengths());
+    auto neps = (prod(select_from(A.lengths(), idx_A, idx_AB))+1)*prod(C.lengths());
 
     dpd_impl = dpd_impl_t::BLOCKED;
     D.reset(C);
@@ -202,6 +211,15 @@ REPLICATED_TEMPLATED_TEST_CASE(indexed_dpd_mult, R, T, all_types)
                 size_AC[irrep_AC]*
                 size_BC[irrep_BC];
     }
+    for (unsigned irrep_AC = 0;irrep_AC < nirrep;irrep_AC++)
+    for (unsigned irrep_BC = 0;irrep_BC < nirrep;irrep_BC++)
+    {
+        unsigned irrep_ABC = irrep_AC^irrep_BC^C.irrep();
+
+        neps += size_ABC[irrep_ABC]*
+                size_AC[irrep_AC]*
+                size_BC[irrep_BC];
+    }
 
     dpd_impl = dpd_impl_t::BLOCKED;
     D.reset(C);