title: OpenMP parallelization of for loop : poor efficiency in my code
body: <p>I have a function which is apparently the bottleneck of my whole program. I thought that parallelization with OpenMP could be helpful.</p>
  
  <p>Here is an working example of my computation (sorry, the function is a little bit long). In my program, some of the work before the 5 nested loops is done somewhere else and is not a problem at all for the efficiency.</p>
  
  <pre><code>#include &lt;vector&gt;
  #include &lt;iostream&gt;
  #include &lt;cmath&gt;
  #include &lt;cstdio&gt;
  #include &lt;chrono&gt;
  #include "boost/dynamic_bitset.hpp"
  
  using namespace std::chrono;
  
  void compute_mddr(unsigned Ns, unsigned int block_size, unsigned int sector)
  {
    std::vector&lt;unsigned int&gt; basis;
    for (std::size_t s = 0; s != std::pow(2,Ns); s++) {
      boost::dynamic_bitset&lt;&gt; s_bin(Ns,s);
      if (s_bin.count() == Ns/2) {
        basis.push_back(s);
      }
    }
    std::vector&lt;double&gt; gs(basis.size());
    for (unsigned int i = 0; i != gs.size(); i++)
      gs[i] = double(std::rand())/RAND_MAX;
  
    unsigned int ns_A = block_size;
    unsigned int ns_B = Ns-ns_A;
    boost::dynamic_bitset&lt;&gt; mask_A(Ns,(1&lt;&lt;ns_A)-(1&lt;&lt;0));
    boost::dynamic_bitset&lt;&gt; mask_B(Ns,((1&lt;&lt;ns_B)-(1&lt;&lt;0))&lt;&lt;ns_A);
  
    // Find the basis of the A block
    unsigned int NAsec = sector;
    std::vector&lt;double&gt; basis_NAsec;
    for (unsigned int s = 0; s &lt; std::pow(2,ns_A); s++) {
      boost::dynamic_bitset&lt;&gt; s_bin(ns_A,s);
      if (s_bin.count() == NAsec)
        basis_NAsec.push_back(s);
    }
    unsigned int bs_A = basis_NAsec.size();
  
    // Find the basis of the B block
    unsigned int NBsec = (Ns/2)-sector;
    std::vector&lt;double&gt; basis_NBsec;
    for (unsigned int s = 0; s &lt; std::pow(2,ns_B); s++) {
      boost::dynamic_bitset&lt;&gt; s_bin(ns_B,s);
      if (s_bin.count() == NBsec)
        basis_NBsec.push_back(s);
    }
    unsigned int bs_B = basis_NBsec.size();
  
    std::vector&lt;std::vector&lt;double&gt; &gt; mddr(bs_A);
    for (unsigned int i = 0; i != mddr.size(); i++) {
      mddr[i].resize(bs_A);
      for (unsigned int j = 0; j != mddr[i].size(); j++) {
        mddr[i][j] = 0.0;
      }
    }
  
    // Main calculation part
    for (unsigned int mu_A = 0; mu_A != bs_A; mu_A++) { // loop 1
      boost::dynamic_bitset&lt;&gt; mu_A_bin(ns_A,basis_NAsec[mu_A]);
      for (unsigned int nu_A = mu_A; nu_A != bs_A; nu_A++) { // loop 2
        boost::dynamic_bitset&lt;&gt; nu_A_bin(ns_A,basis_NAsec[nu_A]);
  
        double sum = 0.0;
  #pragma omp parallel for reduction(+:sum)
        for (unsigned int mu_B = 0; mu_B &lt; bs_B; mu_B++) { // loop 3
          boost::dynamic_bitset&lt;&gt; mu_B_bin(ns_B,basis_NBsec[mu_B]);
  
          for (unsigned int si = 0; si != basis.size(); si++) { // loop 4
            boost::dynamic_bitset&lt;&gt; si_bin(Ns,basis[si]);
            boost::dynamic_bitset&lt;&gt; si_A_bin = si_bin &amp; mask_A;
            si_A_bin.resize(ns_A);
            if (si_A_bin != mu_A_bin)
              continue;
            boost::dynamic_bitset&lt;&gt; si_B_bin = (si_bin &amp; mask_B)&gt;&gt;ns_A;
            si_B_bin.resize(ns_B);
            if (si_B_bin != mu_B_bin)
              continue;
  
            for (unsigned int sj = 0; sj &lt; basis.size(); sj++) { // loop 5
              boost::dynamic_bitset&lt;&gt; sj_bin(Ns,basis[sj]);
              boost::dynamic_bitset&lt;&gt; sj_A_bin = sj_bin &amp; mask_A;
              sj_A_bin.resize(ns_A);
              if (sj_A_bin != nu_A_bin)
                continue;
              boost::dynamic_bitset&lt;&gt; sj_B_bin = (sj_bin &amp; mask_B)&gt;&gt;ns_A;
              sj_B_bin.resize(ns_B);
              if (sj_B_bin != mu_B_bin)
                continue;
              sum += gs[si]*gs[sj];
            }
          }
        }
        mddr[nu_A][mu_A] = mddr[mu_A][nu_A] = sum;
      }
    }
  }
  
  
  int main()
  {
    unsigned int l = 8;
    unsigned int Ns = 2*l;
    unsigned block_size = 6; // must be between 1 and l
    unsigned sector = (block_size%2 == 0) ? block_size/2 : (block_size+1)/2;
  
    high_resolution_clock::time_point t1 = high_resolution_clock::now();
    compute_mddr(Ns,block_size,sector);
    high_resolution_clock::time_point t2 = high_resolution_clock::now();
    duration&lt;double&gt; time_span = duration_cast&lt;duration&lt;double&gt;&gt;(t2 - t1);
    std::cout &lt;&lt; "Function took " &lt;&lt; time_span.count() &lt;&lt; " seconds.";
    std::cout &lt;&lt; std::endl;
  }
  </code></pre>
  
  <p>The <code>compute_mddr</code> function is basically entirely filling up the matrix <code>mddr</code>, and this corresponds to the outermost loops 1 and 2. 
  I decided to parallelize the loop 3, since it's essentially computing a sum. To give order of magnitudes, the loop 3 is over ~50-100 elements in the <code>basis_NBsec</code> vector, while the two innermost loops <code>si</code> and <code>sj</code> run over the ~10000 elements for the vector <code>basis</code>.</p>
  
  <p>However, when running the code (compiled with -O3 -fopenmp on gcc 5.4.0, ubuntu 16.0.4 and i5-4440 cpu) I see either no speed-up (2 threads) or a very limited gain (3 and 4 threads):</p>
  
  <pre><code>time OMP_NUM_THREADS=1 ./a.out
  Function took 230.435 seconds.
  real    3m50.439s
  user    3m50.428s
  sys 0m0.000s
  
  
  time OMP_NUM_THREADS=2 ./a.out 
  Function took 227.754 seconds.
  real    3m47.758s
  user    7m2.140s
  sys 0m0.048s
  
  
  time OMP_NUM_THREADS=3 ./a.out 
  Function took 181.492 seconds.
  real    3m1.495s
  user    7m36.056s
  sys 0m0.036s
  
  
  time OMP_NUM_THREADS=4 ./a.out 
  Function took 150.564 seconds.
  real    2m30.568s
  user    7m56.156s
  sys 0m0.096s
  </code></pre>
  
  <p>If I understand correctly the numbers from user, for 3 and 4 threads the cpu usage is not good (and indeed, when the code is running I get ~250% cpu usage for 3 threads and barely 300% for 4 threads).</p>
  
  <p>It's my first use of OpenMP, I just played with it very quickly on simple examples. Here, as far as I can see, I'm not modifying any of the shared vectors <code>basis_NAsec</code>, <code>basis_NBsec</code> and <code>basis</code> in the parallel part, only reading (this was an aspect pointed out in several related questions I read).</p>
  
  <p>So, what am I doing wrong ?</p>