From ecf10a27cfd0351b9667b4d1c75252de3744151e Mon Sep 17 00:00:00 2001
From: Junya Otsuki <otsuki@cmpt.phys.tohoku.ac.jp>
Date: Tue, 27 Mar 2018 21:30:30 +0900
Subject: [PATCH] Use MPI for frequency loop in G2 calc

---
 c++/pomerol_ed.cpp           | 33 +++++++++++++++++++++++++++++++--
 c++/pomerol_ed.hpp           |  2 +-
 python/pomerol2triqs_desc.py |  2 +-
 3 files changed, 33 insertions(+), 4 deletions(-)

diff --git a/c++/pomerol_ed.cpp b/c++/pomerol_ed.cpp
index b89bccc..fe8414c 100644
--- a/c++/pomerol_ed.cpp
+++ b/c++/pomerol_ed.cpp
@@ -428,12 +428,41 @@ namespace pomerol2triqs {
     std::iota(index_wb.begin(), index_wb.end(), 0);
     // for( auto i : iw_b )  std::cout << i << std::endl;
 
-    g2_t g2(p.n_b, 2*p.n_f, 2*p.n_f);
+    // std::cout << "Start freq loop: rank" << comm.rank() << std::endl;
+    // g2_t g2(p.n_b, 2*p.n_f, 2*p.n_f);
+    // // std::cout << typeid(g2).name() << std::endl;
+    // for(int ib=0; ib<index_wb.size(); ib++)
+    //   for(int if1=0; if1<index_wf.size(); if1++)
+    //     for(int if2=0; if2<index_wf.size(); if2++)
+    //       g2(ib, if1, if2) = -pom_g2(index_wb[ib]+index_wf[if1], index_wf[if2], index_wf[if1]);
+
+    // create a list of three frequency-indices
+    std::vector< std::tuple<int, int, int> > three_freqs;
     for(int ib=0; ib<index_wb.size(); ib++)
       for(int if1=0; if1<index_wf.size(); if1++)
         for(int if2=0; if2<index_wf.size(); if2++)
-          g2(ib, if1, if2) = -pom_g2(index_wb[ib]+index_wf[if1], index_wf[if2], index_wf[if1]);
+          three_freqs.push_back( std::make_tuple(ib, if1, if2) );
+    // std::cout << three_freqs.size() << std::endl;
+
+    // compute g2 value using MPI
+    g2_t g2(p.n_b, 2*p.n_f, 2*p.n_f);
+    for(int i=comm.rank(); i<three_freqs.size(); i+=comm.size()){
+      int ib = std::get<0>(three_freqs[i]);
+      int if1 = std::get<1>(three_freqs[i]);
+      int if2 = std::get<2>(three_freqs[i]);
+      g2(ib, if1, if2) = -pom_g2(index_wb[ib]+index_wf[if1], index_wf[if2], index_wf[if1]);
+    }
+
+    // broadcast results
+    for(int i=0; i<three_freqs.size(); i++){
+      int ib = std::get<0>(three_freqs[i]);
+      int if1 = std::get<1>(three_freqs[i]);
+      int if2 = std::get<2>(three_freqs[i]);
+      int sender = i % comm.size();
+      boost::mpi::broadcast(comm, g2(ib, if1, if2), sender);
+    }
 
+    // std::cout << "End freq loop: rank" << comm.rank() << std::endl;
     return g2;
   }
 
diff --git a/c++/pomerol_ed.hpp b/c++/pomerol_ed.hpp
index fa33505..6dd546c 100644
--- a/c++/pomerol_ed.hpp
+++ b/c++/pomerol_ed.hpp
@@ -82,7 +82,7 @@ namespace pomerol2triqs {
 
     // using w_nu_nup_t = cartesian_product<imfreq, imfreq, imfreq>;
     // using w_l_lp_t   = cartesian_product<imfreq, legendre, legendre>;
-    using g2_t = array<std::complex<double>, 3>;
+    using g2_t = triqs::arrays::array<std::complex<double>, 3>;
     // template <typename Mesh, typename Filler>
     // block2_gf<Mesh, tensor_valued<4>> compute_g2(gf_struct_t const &gf_struct, gf_mesh<Mesh> const &mesh, block_order_t block_order,
     //                                              g2_blocks_t const &g2_blocks, Filler filler) const;
diff --git a/python/pomerol2triqs_desc.py b/python/pomerol2triqs_desc.py
index 0e41211..511deff 100644
--- a/python/pomerol2triqs_desc.py
+++ b/python/pomerol2triqs_desc.py
@@ -67,7 +67,7 @@
 c.add_method("""block_gf<refreq> G_w (gf_struct_t gf_struct, double beta, std::pair<double,double> energy_window, int n_w, double im_shift = 0)""",
              doc = r"""Retarded Green's function on real energy axis""")
 
-c.add_method("""array<std::complex<double>, 3> G2_iw (**pomerol2triqs::g2_iw_inu_inup_params_t)""",
+c.add_method("""triqs::arrays::array<std::complex<double>, 3> G2_iw (**pomerol2triqs::g2_iw_inu_inup_params_t)""",
              doc = r"""Two-particle Green's function, Matsubara frequencies""")
 
 module.add_class(c)