title: OpenMP slower than serial codes
body: <p>I am trying to do the parallelization of a serial preconditioned conjugate gradient solver codes for 3D fire simulation using OpenMP (Intel compiler). But the performance seems not to be improved.</p>
  
  <p>The grid dimension is 79x81x79 and the solver can converge after 565 iterations. The serial codes cost 3.39 seconds and the OpenMP version needs 3.86 seconds on Intel i7 2600 <strong>(OS: openSUSE 13.1)</strong>.</p>
  
  <p>Please help me to check  the codes. Thanks a lot.</p>
  
  <pre><code>//   preconditioned conjugate gradient solver  ...
  void PCGSSolver::solveNew(const Array3D&lt;double&gt;&amp; sn, const Array3D&lt;double&gt;&amp; ae,  const Array3D&lt;double&gt;&amp;aw,
                            const Array3D&lt;double&gt;&amp; as,  const Array3D&lt;double&gt;&amp; an,  const Array3D&lt;double&gt;&amp;at,  const Array3D&lt;double&gt;&amp;ab,
                            const Array3D&lt;double&gt;&amp; ap, Array3D&lt;double&gt;&amp; ptmp){
      std::size_t dimX=sn.getDimI();
      std::size_t dimY=sn.getDimJ();
      std::size_t dimZ=sn.getDimK();
  
  
      Array3D&lt;double&gt; p1(dimX,dimY,dimZ,0.0);
      Array3D&lt;double&gt; res(dimX,dimY,dimZ,0.0);
      Array3D&lt;double&gt; d(dimX,dimY,dimZ,0.0);
      Array3D&lt;double&gt; ain(dimX,dimY,dimZ,0.0);
  
  
      double tiny=1.0e-30;
  #pragma omp parallel
  {
          //Jacobi preconditioner
  #pragma omp for nowait
      for(std::size_t k=1;k&lt;dimZ-1; k++){
          for(std::size_t j=1; j&lt;dimY-1; j++){
              for(std::size_t i=1; i&lt;dimX-1; i++){
                  d(i,j,k)=1./ap(i,j,k);
              }
          }
      }
  #pragma omp for nowait
      for(std::size_t k=1;k&lt;dimZ-1; k++){
          for(std::size_t j=1; j&lt;dimY-1; j++){
              for(std::size_t i=1; i&lt;dimX-1; i++){
                  res(i,j,k)=ae(i,j,k)*ptmp(i+1,j,k) + aw(i,j,k)*ptmp(i-1,j,k)+an(i,j,k)*ptmp(i,j+1,k)+as(i,j,k)*ptmp(i,j-1,k)+
                                  at(i,j,k)*ptmp(i,j,k+1)+ab(i,j,k)*ptmp(i,j,k-1)+sn(i,j,k)-ap(i,j,k)*ptmp(i,j,k);
              }
          }
      }
  
  }
  
  
      double big =1.0e+30;
      double s1old=big;
      //start iteration
      for(std::size_t intswp=0; intswp&lt;this-&gt;nswpvr; intswp++){
  
          double alpha=0.0;
          double bbeta=0.0;
          double s1=0.0;
          double s2=0.0;
          double testir=0.0;
  #pragma omp parallel
  {
  #pragma omp for reduction(+:s1)
          for(std::size_t k=1;k&lt;dimZ-1; k++){
              for(std::size_t j=1; j&lt;dimY-1; j++){
                  for(std::size_t i=1; i&lt;dimX-1; i++){
                      ain(i,j,k)=res(i,j,k)*d(i,j,k);
                      s1+=(res(i,j,k)*ain(i,j,k));
                  }
              }
          }
  
  #pragma omp single
  {
          bbeta=s1/(s1old+tiny);
  }
  #pragma omp for
          for(std::size_t k=1;k&lt;dimZ-1; k++){
              for(std::size_t j=1; j&lt;dimY-1; j++){
                  for(std::size_t i=1; i&lt;dimX-1; i++){
                      p1(i,j,k)=ain(i,j,k)+bbeta*p1(i,j,k);
                  }
              }
          }
  
  #pragma omp for reduction(+:s2)
          for(std::size_t k=1;k&lt;dimZ-1; k++){
              for(std::size_t j=1; j&lt;dimY-1; j++){
                  for(std::size_t i=1; i&lt;dimX-1; i++){
                      ain(i,j,k)=ap(i,j,k)*p1(i,j,k)-ae(i,j,k)*p1(i+1,j,k)-aw(i,j,k)*p1(i-1,j,k)-
                                      an(i,j,k)*p1(i,j+1,k)-as(i,j,k)*p1(i,j-1,k)-
                                      at(i,j,k)*p1(i,j,k+1)-ab(i,j,k)*p1(i,j,k-1);
                      s2+=(p1(i,j,k)*ain(i,j,k));
                  }
              }
          }
  
  #pragma omp single
  {
          alpha=s1/(s2+tiny);
  }
  #pragma omp for reduction(+:testir)
          for(std::size_t k=1;k&lt;dimZ-1; k++){
              for(std::size_t j=1; j&lt;dimY-1; j++){
                  for(std::size_t i=1; i&lt;dimX-1; i++){
                      ptmp(i,j,k)=ptmp(i,j,k)+alpha*p1(i,j,k);
                      res(i,j,k)=res(i,j,k)-alpha*ain(i,j,k);
                      testir+=fabs(res(i,j,k));
                  }
              }
          }
  
  }//==openmp region end
          s1old=s1;
          //test stop criteria
          if(testir &lt; ccvar){
              std::cout&lt;&lt;"PCGS solver coverage at "&lt;&lt;(intswp+1)&lt;&lt;" iterations!"&lt;&lt;std::scientific&lt;&lt;testir&lt;&lt;std::endl;
              return;
          }
  
      }
      std::cout&lt;&lt;"PCGS solver can not coverage "&lt;&lt;std::endl;
  }
  </code></pre>
  
  <p>The Array3D is a my 3 dimension array class.</p>
  
  <pre><code>#ifndef ARRAY3D_H
  #define ARRAY3D_H
  
  #include &lt;vector&gt;
  #include &lt;algorithm&gt;
  
  template&lt;typename T&gt; class Array3D
  {
  public:
      typedef T value_type;
      Array3D(){
          dim_i=dim_j=dim_k=0;
          dim_ij=0;
      }
  
      Array3D(std::size_t size_i, std::size_t size_j, std::size_t size_k){
          this-&gt;resize(size_i,size_j,size_k);
      }
  
      Array3D(std::size_t size_i, std::size_t size_j, std::size_t size_k,const value_type&amp; defaultValue){
          this-&gt;resize(size_i,size_j,size_k,defaultValue);
      }
  
      virtual ~Array3D(){}
  
      std::size_t getDimI()const{
          return this-&gt;dim_i;
      }
  
      std::size_t getDimJ()const{
          return this-&gt;dim_j;
      }
  
      std::size_t getDimK()const{
          return this-&gt;dim_k;
      }
  
      //check if valid indices
      bool checkIndices(std::size_t i, std::size_t j, std::size_t k){
          return (i&lt;this-&gt;dim_i ) &amp;&amp; (j&lt;this-&gt;dim_j) &amp;&amp; (k&lt;this-&gt;dim_k);
      }
      void resize(std::size_t size_i, std::size_t size_j, std::size_t size_k,const value_type&amp; defaultValue){
          this-&gt;resize(size_i,size_j,size_k);
          this-&gt;fillValue(defaultValue);
      }
      //resize the array. The data will be ereased.
      void resize(std::size_t size_i, std::size_t size_j, std::size_t size_k){
          this-&gt;dim_i=size_i;
          this-&gt;dim_j=size_j;
          this-&gt;dim_k=size_k;
          this-&gt;dim_ij=this-&gt;dim_i*this-&gt;dim_j;
  
          std::size_t totalSize=this-&gt;dim_i*this-&gt;dim_j*this-&gt;dim_k;
          this-&gt;data.resize(totalSize);
      }
      std::size_t size()const{
          return this-&gt;data.size();
      }
  
      void fillValue(const value_type&amp; defaultValue){
          std::fill(this-&gt;data.begin(),this-&gt;data.end(),defaultValue);
      }
      value_type minValue()const{
          return *(std::min_element(data.begin(),data.end()));
      }
      value_type maxValue()const{
          return *(std::max_element(data.begin(),data.end()));
      }
  
      //Fill the array value using the sum of two array
      void setValueSum(const Array3D&amp; array1, const Array3D&amp; array2){
          size_t minSize=std::min(std::min(array1.data.size(),array2.data.size()),this-&gt;data.size());
          for(size_t i=0; i&lt;minSize; i++)
              this-&gt;data[i]=array1.data[i]+array2.data[i];
      }
  
      void clear(){
          dim_i=dim_j=dim_k=0;
          dim_ij=0;
          this-&gt;data.clear();
      }
  
      //get value reference at (i,j,k) or (x,y,z) or (u,v,w)...
      const value_type&amp; operator () (std::size_t i, std::size_t j, std::size_t k )const{
          return this-&gt;data.at(this-&gt;calIndex(i,j,k));
      }
  
      value_type&amp; operator ()(std::size_t i, std::size_t j, std::size_t k ){
          return this-&gt;data.at(this-&gt;calIndex(i,j,k));
      }
      //access the raw data by 1D index
      const value_type&amp; operator [] (std::size_t i )const{
          return this-&gt;data.at(i);
      }
      value_type&amp; operator [](std::size_t i ){
          return this-&gt;data.at(i);
      }
      std::vector&lt;value_type&gt;* rawData(){
          return &amp;(data);
      }
  private:
      inline std::size_t calIndex(std::size_t i, std::size_t j, std::size_t k )const{
          return k*this-&gt;dim_ij+j*this-&gt;dim_i+i;
      }
  
  private:
      //dimension of array (i,j,k)(x,y,z)(u,v,w)...
      std::size_t dim_i, dim_j, dim_k;
      //raw data, order is I-J-K
      std::vector&lt;value_type&gt; data;
  
      //dim_i*dim_j
      std::size_t dim_ij;
  };
  
  #endif // ARRAY3D_H
  </code></pre>
  
  <p>I measure the time using a Timer class codes downloaded from internet.</p>
  
  <pre><code>    timer.start();
      PCGSSolver solver;
      solver.setTolerance(this-&gt;ccvar);
      solver.setIteNum(this-&gt;nswpp);
      solver.solveNew(sn,ae,aw,as,an,at,ab,ap,ptmp);
      timer.stop();
      std::cout&lt;&lt;"PCGS time:"&lt;&lt;timer.getElapsedTimeInSec()&lt;&lt;"sec"&lt;&lt;std::endl;
  </code></pre>
  
  <p>Timer.h</p>
  
  <pre><code>    //////////////////////////////////////////////////////////////////////////////
  // Timer.h
  // =======
  // High Resolution Timer.
  // This timer is able to measure the elapsed time with 1 micro-second accuracy
  // in both Windows, Linux and Unix system 
  //
  //  AUTHOR: Song Ho Ahn (song.ahn@gmail.com)
  // CREATED: 2003-01-13
  // UPDATED: 2006-01-13
  //
  // Copyright (c) 2003 Song Ho Ahn
  //////////////////////////////////////////////////////////////////////////////
  
  #ifndef TIMER_H_DEF
  #define TIMER_H_DEF
  
  #ifdef WIN32   // Windows system specific
  #include &lt;windows.h&gt;
  #else          // Unix based system specific
  #include &lt;sys/time.h&gt;
  #endif
  
  class Timer
  {
  public:
      Timer();                                    // default constructor
      ~Timer();                                   // default destructor
  
      void   start();                             // start timer
      void   stop();                              // stop the timer
      double getElapsedTime();                    // get elapsed time in second
      double getElapsedTimeInSec();               // get elapsed time in second (same as getElapsedTime)
      double getElapsedTimeInMilliSec();          // get elapsed time in milli-second
      double getElapsedTimeInMicroSec();          // get elapsed time in micro-second
  
  
  protected:
  
  
  private:
      double startTimeInMicroSec;                 // starting time in micro-second
      double endTimeInMicroSec;                   // ending time in micro-second
      int    stopped;                             // stop flag 
  #ifdef WIN32
      LARGE_INTEGER frequency;                    // ticks per second
      LARGE_INTEGER startCount;                   //
      LARGE_INTEGER endCount;                     //
  #else
      timeval startCount;                         //
      timeval endCount;                           //
  #endif
  };
  
  #endif // TIMER_H_DEF
  </code></pre>
  
  <p>Timer.cpp</p>
  
  <pre><code>    //////////////////////////////////////////////////////////////////////////////
  // Timer.cpp
  // =========
  // High Resolution Timer.
  // This timer is able to measure the elapsed time with 1 micro-second accuracy
  // in both Windows, Linux and Unix system 
  //
  //  AUTHOR: Song Ho Ahn (song.ahn@gmail.com)
  // CREATED: 2003-01-13
  // UPDATED: 2006-01-13
  //
  // Copyright (c) 2003 Song Ho Ahn
  //////////////////////////////////////////////////////////////////////////////
  
  #include  "Timer.h"
  #include &lt;stdlib.h&gt;
  
  ///////////////////////////////////////////////////////////////////////////////
  // constructor
  ///////////////////////////////////////////////////////////////////////////////
  Timer::Timer()
  {
  #ifdef WIN32
      QueryPerformanceFrequency(&amp;frequency);
      startCount.QuadPart = 0;
      endCount.QuadPart = 0;
  #else
      startCount.tv_sec = startCount.tv_usec = 0;
      endCount.tv_sec = endCount.tv_usec = 0;
  #endif
  
      stopped = 0;
      startTimeInMicroSec = 0;
      endTimeInMicroSec = 0;
  }
  
  
  
  ///////////////////////////////////////////////////////////////////////////////
  // distructor
  ///////////////////////////////////////////////////////////////////////////////
  Timer::~Timer()
  {
  }
  
  
  
  ///////////////////////////////////////////////////////////////////////////////
  // start timer.
  // startCount will be set at this point.
  ///////////////////////////////////////////////////////////////////////////////
  void Timer::start()
  {
      stopped = 0; // reset stop flag
  #ifdef WIN32
      QueryPerformanceCounter(&amp;startCount);
  #else
      gettimeofday(&amp;startCount, NULL);
  #endif
  }
  
  
  
  ///////////////////////////////////////////////////////////////////////////////
  // stop the timer.
  // endCount will be set at this point.
  ///////////////////////////////////////////////////////////////////////////////
  void Timer::stop()
  {
      stopped = 1; // set timer stopped flag
  
  #ifdef WIN32
      QueryPerformanceCounter(&amp;endCount);
  #else
      gettimeofday(&amp;endCount, NULL);
  #endif
  }
  
  
  
  ///////////////////////////////////////////////////////////////////////////////
  // compute elapsed time in micro-second resolution.
  // other getElapsedTime will call this first, then convert to correspond resolution.
  ///////////////////////////////////////////////////////////////////////////////
  double Timer::getElapsedTimeInMicroSec()
  {
  #ifdef WIN32
      if(!stopped)
          QueryPerformanceCounter(&amp;endCount);
  
      startTimeInMicroSec = startCount.QuadPart * (1000000.0 / frequency.QuadPart);
      endTimeInMicroSec = endCount.QuadPart * (1000000.0 / frequency.QuadPart);
  #else
      if(!stopped)
          gettimeofday(&amp;endCount, NULL);
  
      startTimeInMicroSec = (startCount.tv_sec * 1000000.0) + startCount.tv_usec;
      endTimeInMicroSec = (endCount.tv_sec * 1000000.0) + endCount.tv_usec;
  #endif
  
      return endTimeInMicroSec - startTimeInMicroSec;
  }
  
  
  
  ///////////////////////////////////////////////////////////////////////////////
  // divide elapsedTimeInMicroSec by 1000
  ///////////////////////////////////////////////////////////////////////////////
  double Timer::getElapsedTimeInMilliSec()
  {
      return this-&gt;getElapsedTimeInMicroSec() * 0.001;
  }
  
  
  
  ///////////////////////////////////////////////////////////////////////////////
  // divide elapsedTimeInMicroSec by 1000000
  ///////////////////////////////////////////////////////////////////////////////
  double Timer::getElapsedTimeInSec()
  {
      return this-&gt;getElapsedTimeInMicroSec() * 0.000001;
  }
  
  
  
  ///////////////////////////////////////////////////////////////////////////////
  // same as getElapsedTimeInSec()
  ///////////////////////////////////////////////////////////////////////////////
  double Timer::getElapsedTime()
  {
      return this-&gt;getElapsedTimeInSec();
  }
  </code></pre>