Large overheads for product range style loops on Green's functions

[HugoStrand]

Dear all,

I am doing timings for the two-particle accumulation in cthyb where the computational bottleneck is the product of two M-matrices in frequency space into a contribution to the two-particle Green's function.

In clef syntax I want to do

G2(w,n1,n2)(i,j,k,l) << G2(w,n1,n2)(i j,k,l) - s * M_il(n1,n2)(i,l) * M_kj(n2+w,n1+w)(k,j);

However this is 10x slower than the direct naive 7-level for loop

    for (const auto &w : b_mesh)
      for (const auto &n1 : f_mesh)
        for (const auto &n2 : f_mesh)
	  for (const auto i : range(G2.target_shape()[0]))
	  for (const auto j : range(G2.target_shape()[1]))
	  for (const auto k : range(G2.target_shape()[2]))
	  for (const auto l : range(G2.target_shape()[3]))
	  G2[w, n1, n2](i, j, k, l) -= s * M_il[n1, n2](i, l) * M_kj[n2 + w, n1 + w](k, j);

Using the c++17 style product range for loops looks nicer

for (const auto &[w, n1, n2] : G2.mesh())
      for (const auto [i, j, k, l] : G2.target_indices())
	G2[w, n1, n2](i, j, k, l) -= s * M_il[n1, n2](i, l) * M_kj[n2 + w, n1 + w](k, j);

however, it comes with a factor of 2 runtime overhead, where 75% of the time is spent in the first loop over the mesh. This is rather far from zero cost abstraction...

Can something be done to get us closer to "zero cost" ?

Best, Hugo

[Wentzell]

Hey @HugoStrand
Thanks for the information!
This is indeed too much overhead for use in performance critical parts.

Could you possibly wrap that up into a small google benchmark?

[Wentzell]

To be rechecked with TRIQS 3.1 (using nda).