[test] add block-sparse Green's functions tests using all symmetries

Author: paco-ri

c++/triqs_xca/block_sparse_backbone.cpp

@@ -501,21 +501,20 @@ void DiagramEvaluator::find_path_correlator(CorrelatorBackbone &backbone, std::v
           }
           ++w;
         }
-        if (f_ix == 3 && kap == 2 && b_ix == 1) {}
         int ip_vct0 = ip;
 
         // now consider mu_ops[mu]
         if (path_all_nonzero) {
           for (int mu = 0; mu < mu_ops.size(); ++mu) {
+            path_all_nonzero  = true; // reset for each mu
             const auto &mu_op = mu_ops[mu];
             ip_old            = ip_vct0;
             ip                = mu_op.get_block_index(ip_vct0);
             if (ip == -1 || Gt.get_zero_block_index(ip) == -1) { // only continue if mu_ops[mu] and Gt have nonzero blocks in right place
               path_all_nonzero = false;
             } else {
               block_dims(vct0 + 1) = mu_op.get_block_size(ip_old, 0);
-              // replace m!!!!
-              ind_path(vct0) = ip;
+              ind_path(vct0)       = ip;
 
               // traverse factors on left-hand side of diagram
               w = vct0 + 1; // reset w to the vertex connected to vertex 0
@@ -532,10 +531,9 @@ void DiagramEvaluator::find_path_correlator(CorrelatorBackbone &backbone, std::v
                 }
                 ++w;
               }
-              if (path_all_nonzero) {
+              if (path_all_nonzero && b_ix == ind_path(2 * m - 1)) {
                 // evaluate the diagram with these directions, poles, and orbital indices
                 // b_ix is the block index for the first edge
-                // DOUBLE-CHECK ALL INDEXING
                 eval_correlator_fixed_indices(backbone, ind_path, block_dims);
                 for (int t = 0; t < r; ++t) {
                   Tmuop(t, range(0, block_dims(2 * m)), range(0, block_dims(1))) =

examples/block_sparse_backbone_third_order.cpp

@@ -5,8 +5,6 @@
 #include <triqs_xca/block_sparse.hpp>
 #include <triqs_xca/block_sparse_manual.hpp>
 #include <triqs_xca/block_sparse_backbone.hpp>
-// temporary include
-#include <triqs_xca/dense_backbone.hpp>
 
 using namespace nda;
 using namespace cppdlr;
@@ -267,10 +265,6 @@ int main() {
 
   auto D = DiagramEvaluator(beta, Lambda, eps, Deltat, nda::make_regular(dlr_rf / beta), Gt, Fq); // create DiagramEvaluator object
 
-  // temporary dense calculation
-  auto Fset = DenseFSet(Fs_dense, F_dags_dense, hyb_coeffs);
-  auto DD = DenseDiagramEvaluator(beta, itops, Deltat, Deltat, dlr_rf, Gt_dense, Fset);
-
   auto Deltadlr                            = itops.vals2coefs(Deltat); //obtain dlr coefficient of Delta(t)
   nda::vector<double> dlr_rf_reflect       = -dlr_rf;
   nda::array<dcomplex, 3> Deltadlr_reflect = Deltadlr * 1.0;
@@ -293,12 +287,6 @@ int main() {
     third_order_result  = BlockDiagOpFun(Sigma_third_gf);
     D.reset(); // reset the DiagramEvaluator for the next topology
 
-    // temp for verification
-    auto BB = Backbone(topologies(i, _, _), n);
-    DD.eval_self_energy(BB);
-    auto dense_result = DD.Sigma;
-    DD.reset();
-
     // Compute third-order contribution using old code
     nda::array<dcomplex, 3> TCA_old(r, N, N);
     TCA_old = 0;
@@ -318,9 +306,6 @@ int main() {
 
     std::cout << "max error in block 0: " << nda::max_element(nda::abs(TCA_old(_, range(0, 4), range(0, 4)) - third_order_result.get_block(0)))
               << std::endl;
-    std::cout << "TCA_old block 0:\n" << TCA_old(5, range(0, 4), range(0, 4)) << std::endl;
-    std::cout << "third_order_result block 0:\n" << third_order_result.get_block(0)(5, range(0, 4), range(0, 4)) << std::endl;
-    std::cout << "dense result block 0:\n" << dense_result(5, range(0, 4), range(0, 4)) << std::endl;
     std::cout << "max error in block 1: " << nda::max_element(nda::abs(TCA_old(_, range(4, 10), range(4, 10)) - third_order_result.get_block(1)))
               << std::endl;
     std::cout << "max error in block 2: " << nda::max_element(nda::abs(TCA_old(_, range(10, 11), range(10, 11)) - third_order_result.get_block(2)))

test/c++/test_block_sparse_OCA_gf_manual.cpp

@@ -213,7 +213,6 @@ TEST(DenseOCAGF, two_band_discrete_bath_dense) {
   hyb_F Delta_F(16, r, dim), Delta_F_reflect(16, r, dim);
   Delta_F.update_inplace(Delta_decomp, dlr_it, Fs_dense, F_dags_dense); // Compression of Delta(t) and F, F_dag matrices
   Delta_F_reflect.update_inplace(Delta_decomp_reflect, dlr_it, F_dags_dense, Fs_dense);
-  nda::vector<int> fb      = {1, 0};
   nda::array<int, 2> D_OCA = {{0, 2}, {1, 3}};
   auto OCA_gf_old          = G_Diagram_calc_sum_all(Delta_F, Delta_F_reflect, D_OCA, Gt_dense, itops, beta, Fs_dense, F_dags_dense);
   auto OCA_gf_old_eq       = eval_eq(itops, OCA_gf_old, n_quad);

test/c++/test_block_sparse_backbone_eval.cpp

@@ -380,7 +380,6 @@ TEST(Backbone, OCA_semicircle_bath_aaa) {
   auto OCA_result_gf            = D.compute_self_energy(topology);
   auto OCA_result               = BlockDiagOpFun(OCA_result_gf);
 
-  // Use slightly relaxed tolerance to account for accumulated truncation errors from SVD in hyb_decomp
   ASSERT_LE(nda::max_element(nda::abs(OCA_result.get_block(0) - OCA_old(_, range(0, 4), range(0, 4)))), eps);
   ASSERT_LE(nda::max_element(nda::abs(OCA_result.get_block(1) - OCA_old(_, range(4, 10), range(4, 10)))), eps);
   ASSERT_LE(nda::max_element(nda::abs(OCA_result.get_block(2) - OCA_old(_, range(10, 11), range(10, 11)))), eps);