parallel C code for distance computation

Question

I've a C code which compute the distance between two sets of nodes (three coordinate each), even though my code has been fast enough yet, I want to boost it up a bit more using parallel computing. I've already found some info about openMP and I'm trying to use it right now, but there's something a bit weird. Without omp the code cpu time is 20s, adding the two pragma lines it takes 160s! How could it happen?

I append my code down here

float computedist(float **vG1, float **vG2, int ncft, int ntri2, int jump, float *dist){
    int k = 0, i, j;
    float min = 0;
    float max = 0;
    float avg = 0;
    float *d = malloc(3*sizeof(float));
    float diff;

    #pragma omp parallel
    for(i=0;i<ncft;i+=jump){
        #pragma omp parallel
        for(j=0;j<ntri2;j++){
            d[0] = vG1[i][0] - vG2[j][0];
            d[1] = vG1[i][1] - vG2[j][1];
            d[2] = vG1[i][2] - vG2[j][2];
            diff = sqrt(pow(d[0],2) + pow(d[1],2) + pow(d[2],2));
            if(j==0)
                dist[k] = diff;
            else
                if(diff<dist[k])
                    dist[k] = diff;

        }
        avg += dist[k];
        if(dist[k]>max)
            max = dist[k];
        k++;
    }

    printf("max distance: %f\n",max);
    printf("average distance: %f\n",avg/(int)(ncft/jump));

    free(d);

    return max;
}

Thank you so much for any help

"How could it happen?" - usual cause is inappropriate parallelisation scheme, either through locality of reference or too much synchronisation (or both). — Mitch Wheat
If you set an environment variable OMP_NUM_THREADS to 1, and run your OpenMP program with a single thread, how much time will it take? — Alexey Kukanov
@AlexeyKukanov is it ok putting void omp_set_num_threads(int num_threads) before the parallel loop? — Nicholas
@AlexeyKukanov I put omp_set_num_threads(1) before the parallel loop and it takes 23s using omp_set_num_threads(2) it gives 50s! — Nicholas
The new code looks fine to me. Do you call computedist only once, or in a loop (and how many times)? How big are ntri1 and ntri2? — Alexey Kukanov

Alexey Kukanov Alexey Kukanov · Accepted Answer · 2012-01-30T09:58:43

(The answer below refers to the initial code in the question, which since then was improved with applying these suggestions)

You need to read more about how to use OpenMP. The specification is available at http://www.openmp.org; and there are links to tutorials and other resources.

I will point out some issues in your code and give recommendations how to fix those.

    float *d = malloc(3*sizeof(float));
    float diff;

d is used as a temporary variable, so should be marked as private in #pragma omp parallel for (see below) to avoid data races. Meanwhile, instead of dynamic allocation, I would just use 3 separate floats. diff also holds a temporary value, so should also be private.

    #pragma omp parallel
    for(i=0;i<ncft;i+=jump){
        #pragma omp parallel
        for(j=0;j<ntri2;j++){

You created a parallel region where each thread executes the whole loop (because the region does not contain any work sharing constructs), and inside it you created a nested region with a new(!) set of threads, also each executing the whole inner loop. It adds lots of overhead and unnecessary computation to your program. What you need is #pragma omp parallel for, and only applied to the outer loop.

            d[0] = vG1[i][0] - vG2[j][0];
            d[1] = vG1[i][1] - vG2[j][1];
            d[2] = vG1[i][2] - vG2[j][2];
            diff = sqrt(pow(d[0],2) + pow(d[1],2) + pow(d[2],2));

Not related to parallelism, but why calling pow just to compute squares? A good old multiplication would likely be both simpler to read and faster.

            if(j==0)
                dist[k] = diff;
            else
                if(diff<dist[k])
                    dist[k] = diff;

Since the action is the same (dist[k]=diff;), the code can be simplified by combining two conditions with || (logical or).

        }
        avg += dist[k];
        if(dist[k]>max)
            max = dist[k];

Here you compute aggregate values across the outer loop. In OpenMP, this is done with reduction clause of #pragma omp for.

        k++;
    }

Currently, you increment k at each iteration, thus creating an unnecessary dependency between iterations that leads to a data race in parallel code. According to your code, k is just a convenient "alias" for i/jump - so just assign it as such at the beginning of the iteration, and make private.

parallel C code for distance computation

3 Answers