title: About the combination of OpenMP and -Ofast
body: <p>I implemented <strong>OpenMP parallelization in a for loop</strong> where I have a <strong>sum</strong> that is the principal cause of slowing down my code. When I did so, I found out that the final results were not the same that I obtained for the non-parallelize code (which is written in C). So first, one might think "well, I just didn't implemented well the parallelization" but the curious thing is that when I run the parallelized code with the <strong>-Ofast</strong> optimization suddenly the results are correct.</p>
  
  <p>That would be:</p>
  
  <ul>
  <li>-O0 correct</li>
  <li>-Ofast correct</li>
  <li>OMP -O0 wrong </li>
  <li>OMP -O1 wrong</li>
  <li>OMP -O2 wrong</li>
  <li>OMP -O3 wrong</li>
  <li>OMP -Ofast correct!</li>
  </ul>
  
  <p>What could -Ofast be doing that solves an error that only appears when I implement openmp? 
  Any recommendation of what could I check or test?
  Thanks!</p>
  
  <p><strong>EDIT</strong>
  Here I include the smallest version of my code that still reproduces the problem.</p>
  
  <pre><code>#include &lt;stdlib.h&gt;
  #include &lt;stdio.h&gt;
  #include &lt;math.h&gt;
  #include &lt;gsl/gsl_rng.h&gt;
  #include &lt;gsl/gsl_randist.h&gt;
  
  #define LENGTH 100
  #define R 50.0
  #define URD 1.0/sqrt(2.0)
  #define PI (4.0*atan(1.0)) //pi
  
  const gsl_rng_type * Type;
  gsl_rng * item;
  
  double CalcDeltaEnergy(double **M,int sx,int sy){
      double DEnergy,r,zz;
      int k,j;
      double rrx,rry;
      int rx,ry;
      double Energy, Cpm, Cmm, Cmp, Cpp;
      DEnergy = 0;
  
      //OpenMP parallelization:
      #pragma omp parallel for reduction (+:DEnergy)
      for (int index = 0; index &lt; LENGTH*LENGTH; index++){
          k = index % LENGTH;
          j = index / LENGTH;
  
      zz = 0.5*(1.0 - pow(-1.0, k + j + sx + sy));
      for (rx = -1; rx &lt;= 1; rx++){
          for (ry = -1; ry &lt;= 1; ry++){
              rrx = (sx - k - rx*LENGTH)*URD;
              rry = (sy - j - ry*LENGTH)*URD;
  
              r = sqrt(rrx*rrx + rry*rry + zz);
              if(r != 0 &amp;&amp; r &lt;= R){
                  Cpm = sqrt((rrx+0.5*(0.702*cos(M[k][j])-0.702*cos(M[sx][sy])))*(rrx+0.5*(0.702*cos(M[k][j])-0.702*cos(M[sx][sy]))) + (rry+0.5*(0.702*sin(M[k][j])-0.702*sin(M[sx][sy])))*(rry+0.5*(0.702*sin(M[k][j])-0.702*sin(M[sx][sy]))) + zz);
                  Cmm = sqrt((rrx-0.5*(0.702*cos(M[k][j])-0.702*cos(M[sx][sy])))*(rrx-0.5*(0.702*cos(M[k][j])-0.702*cos(M[sx][sy]))) + (rry-0.5*(0.702*sin(M[k][j])-0.702*sin(M[sx][sy])))*(rry-0.5*(0.702*sin(M[k][j])-0.702*sin(M[sx][sy]))) + zz);
                  Cpp = sqrt((rrx+0.5*(0.702*cos(M[k][j])+0.702*cos(M[sx][sy])))*(rrx+0.5*(0.702*cos(M[k][j])+0.702*cos(M[sx][sy]))) + (rry+0.5*(0.702*sin(M[k][j])+0.702*sin(M[sx][sy])))*(rry+0.5*(0.702*sin(M[k][j])+0.702*sin(M[sx][sy]))) + zz);
                  Cmp = sqrt((rrx-0.5*(0.702*cos(M[k][j])+0.702*cos(M[sx][sy])))*(rrx-0.5*(0.702*cos(M[k][j])+0.702*cos(M[sx][sy]))) + (rry-0.5*(0.702*sin(M[k][j])+0.702*sin(M[sx][sy])))*(rry-0.5*(0.702*sin(M[k][j])+0.702*sin(M[sx][sy]))) + zz);
                  Cpm = 1.0/Cpm;
                  Cmm = 1.0/Cmm;
                  Cpp = 1.0/Cpp;
                  Cmp = 1.0/Cmp;
                  Energy = (Cpm + Cmm - Cpp - Cmp)/(0.702*0.702); // S=cte=1
  
                  DEnergy -= 2.0*Energy;
              }
          }
      }
      }
  return DEnergy;
  }
  
  void Initialize(double **M){
  double random;
      for(int i=0;i&lt;(LENGTH-1);i=i+2){
            for(int j=0;j&lt;(LENGTH-1);j=j+2) {
                random=gsl_rng_uniform(item);
                if (random&lt;0.5) M[i][j]=PI/4.0;
                else M[i][j]=5.0*PI/4.0;
  
                random=gsl_rng_uniform(item);
                if (random&lt;0.5) M[i][j+1]=3.0*PI/4.0;
                else M[i][j+1]=7.0*PI/4.0;
  
                random=gsl_rng_uniform(item);
                if (random&lt;0.5) M[i+1][j]=3.0*PI/4.0;
                else M[i+1][j]=7.0*PI/4.0;
  
                random=gsl_rng_uniform(item);
                if (random&lt;0.5) M[i+1][j+1]=PI/4.0;
                else M[i+1][j+1]=5.0*PI/4.0;
            }
      }
  }
  
  int main(){
      //Choose and initiaze the random number generator
      gsl_rng_env_setup();
      Type = gsl_rng_default; //default=mt19937, ran2, lxs0
      item = gsl_rng_alloc (Type);
  
      double **S; //site matrix
      S = (double **) malloc(LENGTH*sizeof(double *));
      for (int i = 0; i &lt; LENGTH; i++)
          S[i] = (double *) malloc(LENGTH*sizeof(double ));
  
      //Initialization
      Initialize(S);
  
      int l,m;
      for (int cl = 0; cl &lt; LENGTH*LENGTH; cl++) {
          l = gsl_rng_uniform_int(item, LENGTH); // RNG[0, LENGTH-1]
          m = gsl_rng_uniform_int(item, LENGTH); // RNG[0, LENGTH-1]
          printf("%lf\n", CalcDeltaEnergy(S, l, m));
      }
  
  
      //Free memory
      for (int i = 0; i &lt; LENGTH; i++)
          free(S[i]);
      free(S);
      return 0;
  } 
  </code></pre>
  
  <p>I compile with: </p>
  
  <pre><code>g++ [optimization] -lm test.c -o test.x -lgsl -lgslcblas -fopenmp
  </code></pre>
  
  <p>and run with:</p>
  
  <pre><code>GSL_RNG_SEED=123; ./test.x &gt; test.dat
  </code></pre>
  
  <p>Comparing the outputs for different optimizations one can see what I stated before.</p>