title: OpenMP parallelize but not accelerate code
body: <p>I've been trying to accelerate the following code</p>
  <pre><code>
  void calcstiffTestOMP(TPZCompMesh *cmesh,int nthread){
      int64_t nelem = cmesh-&gt;NElements();
      omp_set_num_threads(nthread);
  
      auto beginCalcStiff = std::chrono::high_resolution_clock::now();
      int s=0;
  #pragma omp parallel for reduction(+:s)
      for (int64_t iel = 0; iel &lt; nelem; iel++)
      {
          TPZCompEl *el = cmesh-&gt;Element(iel);
          if (!el) continue;
          TPZElementMatrix ek(cmesh, TPZElementMatrix::EK), ef(cmesh, TPZElementMatrix::EF);
          el-&gt;CalcStiff(ek, ef);
      }
      auto endCalcStiff = std::chrono::high_resolution_clock::now();
      auto elapsed = std::chrono::duration_cast&lt;std::chrono::nanoseconds&gt;(endCalcStiff - beginCalcStiff);
      unsigned long int duration = static_cast&lt;unsigned long int&gt;(elapsed.count());
      std::cout &lt;&lt; &quot; CalcStiff parallel duration omp = &quot; &lt;&lt; duration*1E-9 &lt;&lt; &quot; seconds, with nthreads= &quot; &lt;&lt; nthread &lt;&lt; std::endl;
  }
  </code></pre>
  <p>but the simulation time hasn't changed, independently of the number of threads.
  The same code has been accelerated when using TBB directives, as follows</p>
  <pre><code>void calcstiffTestTBB(TPZCompMesh *cmesh,int nthread){
      auto beginCalcStiff = std::chrono::high_resolution_clock::now();
  
      int64_t nelem = cmesh-&gt;NElements();
      tbb::task_scheduler_init init(nthread);
      tbb::parallel_for( tbb::blocked_range&lt;int64_t&gt;(0,nelem),
                        [&amp;](tbb::blocked_range&lt;int64_t&gt; r){
      for (int64_t iel = r.begin(); iel &lt; r.end(); iel++)
      {
          TPZCompEl *el = cmesh-&gt;Element(iel);
          if (!el) continue;
          TPZElementMatrix ek(cmesh, TPZElementMatrix::EK), ef(cmesh, TPZElementMatrix::EF);
  
          el-&gt;CalcStiff(ek, ef);
      }
      });
      auto endCalcStiff = std::chrono::high_resolution_clock::now();
      auto elapsed = std::chrono::duration_cast&lt;std::chrono::nanoseconds&gt;(endCalcStiff - beginCalcStiff);
      unsigned long int duration = static_cast&lt;unsigned long int&gt;(elapsed.count());
      std::cout &lt;&lt; &quot; CalcStiff duration tbb= &quot; &lt;&lt; duration*1E-9 &lt;&lt; &quot; seconds, with nthreads= &quot;&lt;&lt;nthread &lt;&lt; std::endl;
  }
  
  </code></pre>
  <p>The calcstiff function is given by:</p>
  <pre><code>void TPZMultiphysicsCompEl&lt;TGeometry&gt;::CalcStiff(TPZElementMatrix &amp;ek, TPZElementMatrix &amp;ef)
  {
      TPZMaterial * material = Material();
      if(!material){
          PZError &lt;&lt; &quot;Error at &quot; &lt;&lt; __PRETTY_FUNCTION__ &lt;&lt; &quot; this-&gt;Material() == NULL\n&quot;;
          ek.Reset();
          ef.Reset();
          return;
      }
      
      TPZNullMaterial *nullmat = dynamic_cast&lt;TPZNullMaterial *&gt;(material);
      if(nullmat)
      {
          ek.Reset();
          ef.Reset();
          ek.fType = TPZElementMatrix::EK;
          ef.fType = TPZElementMatrix::EF;
          return;
      }
      InitializeElementMatrix(ek,ef);
      
      if (this-&gt;NConnects() == 0) return;//boundary discontinuous elements have this characteristic
      
      TPZManVector&lt;TPZMaterialData,6&gt; datavec;
      const int64_t nref = fElementVec.size();
      datavec.resize(nref);
      InitMaterialData(datavec);
      
      TPZManVector&lt;TPZTransform&lt;&gt; &gt; trvec;
      AffineTransform(trvec);
      
      int dim = Dimension();
      TPZAutoPointer&lt;TPZIntPoints&gt; intrule;
      
      TPZManVector&lt;REAL,4&gt; intpointtemp(TGeometry::Dimension,0.);
      REAL weight = 0.;
      
      TPZManVector&lt;int,4&gt; ordervec;
      //ordervec.resize(nref);
      for (int64_t iref=0;  iref&lt;nref; iref++)
      {
          TPZInterpolationSpace *msp  = dynamic_cast &lt;TPZInterpolationSpace *&gt;(fElementVec[iref].Element());
          int svec;
          if(msp)
          {
              ordervec.Resize(ordervec.size()+1);
              svec = ordervec.size();
          }
          else
          {
              continue;
          }
          datavec[iref].p =1; //msp-&gt;MaxOrder();
          ordervec[svec-1] = datavec[iref].p;
      }
      int order = material-&gt;IntegrationRuleOrder(ordervec);
      
      TPZGeoEl *ref = this-&gt;Reference();
      intrule = ref-&gt;CreateSideIntegrationRule(ref-&gt;NSides()-1, order);
      
      TPZManVector&lt;int,4&gt; intorder(dim,order);
      intrule-&gt;SetOrder(intorder);
      int intrulepoints = intrule-&gt;NPoints();
      if(intrulepoints &gt; 1000) {
          DebugStop();
      }
      
      TPZFMatrix&lt;REAL&gt; jac, axe, jacInv;
      REAL detJac;
      for(int int_ind = 0; int_ind &lt; intrulepoints; ++int_ind)
      {
          intrule-&gt;Point(int_ind,intpointtemp,weight);
          ref-&gt;Jacobian(intpointtemp, jac, axe, detJac , jacInv);
          weight *= fabs(detJac);
          for (int i = 0; i &lt; fElementVec.size(); i++) {
              TPZInterpolationSpace *msp  = dynamic_cast &lt;TPZInterpolationSpace *&gt;(fElementVec[i].Element());
              if (!msp) {
                  continue;
              }
              datavec[i].intLocPtIndex = int_ind;
          }
          
          this-&gt;ComputeRequiredData(intpointtemp,trvec,datavec);
          
  
      
         material-&gt;Contribute(datavec,weight,ek.fMat,ef.fMat);
          
  
          }//loop over integration points
      
      CleanupMaterialData(datavec);
  }
  
  </code></pre>
  <p>The code has been parallelized but not accelerated. Are we using omp directives correctly?
  Does OMP create mutex sections on its own? Is there a better way to parallelize this code using OMP?</p>
  <p>PZ is a library which serve as a template for finite element code, and can be found in <a href="https://github.com/labmec/neopz" rel="nofollow noreferrer">https://github.com/labmec/neopz</a></p>
  <p>SOLUTION: Adding the 'schedule(dynamic,1)'  flag to '#pragma' directive solved this problem. The code had different thread workloads
  and OMP was using a single thread to handle the slow section of the loop. The 'schedule(dynamic,1)' directive solved this.</p>