CODEWITHSUNDEEP
c++gccopenmp
sarema
sarema

Reputation: 725

C++ openmp parallel computation calculates wrong results

I have an algorithm that I want to execute in parallel using openmp. I could validate the results for single-threaded execution (OMP_NUM_THREADS=1), but the results are slightly different as soon as I set the number of threads to 2 or higher. I also tried parallelizing the inner for-loops, but that doesn't yield correct results, either.

I'm fairly new to openmp and multi-threading in general. I suspect that my implementation shares variables between threads improperly somehow, but I can't figure it out.

extern "C" double *lomb_scargle(double *x, double *y, double *f, int NX, int NF) {
    double *result = (double*) malloc(2 * NF * sizeof(double));
    double w, tau, SS, SC, SST1, SST2, SCT1, SCT2, Ai, Bi;
    int i, j;

    #pragma omp parallel for
    for (i=0; i<NF; i++) {
        w = 2 * M_PI * f[i];
        SS = 0.;
        SC = 0.;
        for (j=0; j<NX; j++) {
            SS += sin(2*w*x[j]);
            SC += cos(2*w*x[j]);
        }
        tau = atan2(SS, SC) / (2 * w);
        SCT1 = 0.;
        SCT2 = 0.;
        SST1 = 0.;
        SST2 = 0.;
        for (j=0; j<NX; j++) {
           SCT1 += y[j] * cos(w * (x[j] - tau)); 
           SCT2 +=    pow(cos(w * (x[j] - tau)),2); 
           SST1 += y[j] * sin(w * (x[j] - tau)); 
           SST2 +=    pow(sin(w * (x[j] - tau)),2); 
        }
        Ai = SCT1 / SCT2;
        Bi = SST1 / SST2;
        // result contains the amplitude first, and then the phase 
        result[i] = sqrt(Ai*Ai + Bi*Bi);
        result[i+NF] = - (atan2(Bi, Ai) + w * tau);
    }
    return result;
}

EDIT: typo

Upvotes: 1

Views: 418

Answers (1)

Florian Weimer
Florian Weimer

Reputation: 33757

By default, OpenMP shares all variables declared in an outer scope between all workers. You need to move the temporary variables into the inner block (or declare them private).

Upvotes: 3

Related Questions