Openmp results are not robust?

Question

I am new to openmp, when I add openmp into my code, I found the results are not the same in different run. is this the inherent problem of openmp or my code problem? Thank you!

#include "stdafx.h"
#include <iostream>
#include <vector>
#include <fstream>
#include <math.h>
#include<sstream>
#include <omp.h>
using namespace std;

int main()
{
    double a[1000];
    for (int i = 0; i < 1000; i++)
    {
        a[i] = 0;
    }
    for (int m = 0; m < 1000; m++)
    {
    #pragma omp parallel for shared(a)
        for (int i = 0; i < 1000000; i++)
        {
            int b = pow(i, 0.5);
            if (b < 1000)
            {
                //cout << i <<" "<<sin(i)<< endl;

                a[b] += sin(i);
            }
        }
    }

    fstream all_temp;
    all_temp.open("temperatureabcd.dat", fstream::out);
    for (int aaa = 0; aaa < 1000; aaa++)
    {
        all_temp << a[aaa] << endl;
    }
    all_temp.close();
    return 0;
}

Answer 1

Your code is doing an array reduction. The simple solution would be to do

#pragma omp parallel for reduction(+:a[:1000])

However MSVC which you are using (which I infer from the precompiled header stdafx.h) does not support OpenMP array reduction. You can do the array reduction by hand by changing the code a bit like this

double a[1000] = {0};
for (int m = 0; m < 1000; m++) {
  #pragma omp parallel
  {
    double a2[1000] = {0};
    #pragma omp for nowait
    for (int i = 0; i < 1000000; i++) {
      int b = pow(i, 0.5);
      if (b < 1000) a2[b] += sin(i);
    }
    #pragma omp critical
    for(int i = 0; i < 1000; i++) a[i] += a2[i];
  }
}

The other problem is that floating point addition is not associative so the order that the reduction is done matters. This can be fixed with a bit more work but it's probably not your main problem.