title: Simulation result differing between parallel and serial implementations (And Compiler settings)
body: <p>I'm conducting some Monte Carlo simulations for uni and I've noticed that I get slightly different results depending on whether I run my code serially or in parallel. I generate data points (x[i]) from the Weibull distribution, compute the (theoretical) CDF at x[i] and then compare that to the empirical CDF (1/n) through the KS test. I do this 10,000 times and then print the 95% value to a csv file. This whole procedure is repeated 100 times, so the resulting csv is one column with 100 rows. I take the mean of these (manually in Libre Math) and get the following results for the different scenarios</p>
  
  <p>Serially:
  Mean = 1.340 +/- 0.001 which is the correct value (I have a benchmark to compare against)</p>
  
  <p>Parallel:
  Mean = 1.345 +/- 0.002 </p>
  
  <p>I'm Compiling with g++ 5.4.0 and the following flags </p>
  
  <pre><code>    -fopenmp -fno-inline -std=c++17 -g -Wfatal-errors -Wextra -Wall
  </code></pre>
  
  <p>The error margins are just what I've observed and aren't statistically rigorous. The difference is enough to effect my conclusion. This is a drastically simplified version of the real code and all of my initial thoughts (e.g. objects being shared) didn't fix the problem. </p>
  
  <p>Additionally I've noticed that if I include the following optimisation flag I get the same ( and correct ) for both serial and parallel cases.</p>
  
  <pre><code> -Ofast
  </code></pre>
  
  <p>My question is </p>
  
  <p>1) Why do I get different results?</p>
  
  <p>2) What am I doing to upset the compiler in one case but not the other?</p>
  
  <p>For interests sake, I'm running this on a dual core machine, that's been hyper-threaded to have 4 logical cores.
  Thank you very much.</p>
  
  <pre><code>#include &lt;iostream&gt;
  #include &lt;cmath&gt;
  #include &lt;thread&gt;
  #include &lt;stdio.h&gt;
  #include &lt;stdlib.h&gt;
  #include &lt;omp.h&gt; 
  #include &lt;random&gt;
  #include &lt;algorithm&gt;
  #include &lt;fstream&gt;
  
  
  // INITIAL
  const float scale = 400.0; const float shape = 0.577; 
  const int n=100;    
  const int numOfCV=100;const int numOfSims=10000;
  std::random_device rd;
  std::mt19937_64 e2(rd());
  std::uniform_real_distribution&lt;double&gt; uniform_dist(0, 1);
  
  int main(){
  
     // DETERMINE NUMBER OF CORES
      const unsigned Cores = std::thread::hardware_concurrency();
      std::cout&lt;&lt;"Cores="&lt;&lt;Cores&lt;&lt;std::endl;
      omp_set_num_threads(Cores);
      //OPEN OUTPUT FILE
      std::ofstream csvfile; csvfile.open("KS_D95.csv");
  
      #pragma omp parallel for // I comment this line out to run serially
      for(int CV_inc = 0; CV_inc &lt; numOfCV; ++CV_inc){
  
        double x[n];
        double CDF_emp[n];
        double CDF_theory[n];
        for(int i=0; i&lt;n ;++i){
            CDF_emp[i]=(1.0+i)/n;
        }
          double DnP[numOfSims];
          double D95;
          double dt;
          double dt0;
  
          int simNum = -1;
          while(simNum &lt; numOfSims-1 ){simNum++;
              // GENERATE DATA POINTS
              for(int i=0; i&lt;n ;++i){
                  x[i] = scale*pow(-log(uniform_dist(e2)),1/shape);
              }
              // CALCULATE CDF        
              for(int i=0; i&lt;n ;++i){
                  CDF_theory[i]=1-exp(-pow(x[i]/scale,shape));
                  if ( std::isinf(CDF_theory[i])) {
                      printf("CDF_theory[i]=%f \n",CDF_theory[i]);
                  }
              }
              std::sort(CDF_theory,CDF_theory+n);
  
              // KS TEST 
              dt=0.0;dt0=0.0;
              for(int i=1; i&lt;n ;++i){
              if (dt&lt;std::abs(CDF_theory[i]-CDF_emp[i])){
              dt=std::abs(CDF_theory[i]-CDF_emp[i]);
                  }
              if (dt0&lt;std::abs(CDF_theory[i]-CDF_emp[i-1])){
                  dt0=std::abs(CDF_theory[i]-CDF_emp[i-1]);
              }
          }
              DnP[simNum]= sqrt(n)*std::max(dt,dt0);
          }
          // PRINT RESULT
          std::sort(DnP,DnP+numOfSims);
          D95=DnP[(int)(numOfSims*0.95)-1];
          std::cout&lt;&lt;"line complete, KS_D95="&lt;&lt; D95 &lt;&lt;std::endl;
          csvfile&lt;&lt;D95&lt;&lt;std::endl;
  
      }
      return 0;
  }
  </code></pre>