openMP for loop increment statment handling

Question

for (uint i = 0; i < x; i++) {
   for (uint j = 0; j < z; j++) {
           if (inFunc(p, index)) {
                XY[2*nind] = i;
                XY[2*nind + 1] = j;
                nind++;
           }
   }
}

here x = 512 and z = 512 and nind = 0 initially and XY[2*x*y].

I want to optimize this for loops with openMP but 'nind' variable is closely binded serially to for loop. I have no clue because I am also checking a condition and so some of the time it will not enter in if and will skip increment or it will enter increment nind. openMP threads will increment nind variable as first come will increment nind firstly. Is there any way to unbind it. ('binding' I mean only can be implemented serially).

Answer 1

A typical cache-friendly solution in that case is to collect the (i,j) pairs in private arrays, then concatenate those private arrays at the end, and finally sort the result if needed:

#pragma omp parallel
{
  uint myXY[2*z*x];
  uint mynind = 0;

  #pragma omp for collapse(2) schedule(dynamic,N)
  for (uint i = 0; i < x; i++) {
    for (uint j = 0; j < z; j++) {
      if (inFunc(p, index)) {
        myXY[2*mynind] = i;
        myXY[2*mynind + 1] = j;
        mynind++;
      }
    }
  }

  #pragma omp critical(concat_arrays)
  {
    memcpy(&XY[2*nind], myXY, 2*mynind*sizeof(uint));
    nind += mynind;
  }
}

// Sort the pairs if needed
qsort(XY, nind, 2*sizeof(uint), compar);

int compar(const uint *p1, const uint *p2)
{
   if (p1[0] < p2[0])
     return -1;
   else if (p1[0] > p2[0])
     return 1;
   else
   {
     if (p1[1] < p2[1])
       return -1;
     else if (p1[1] > p2[1])
       return 1;
   }
   return 0;
}

You should experiment with different values of N in the schedule(dynamic,N) clause in order to achieve the best trade-off between overhead (for small values of N) and load imbalance (for large values of N). The comparison function compar could probably be written in a more optimal way.

The assumption here is that the overhead from merging and sorting the array is small. Whether that will be the case depends on many factors.

Answer 2

Here is a variation on Hristo Iliev's good answer.

The important parameter to act on here is the index of the pairs rather than the pairs themselves.

We can fill private arrays of the pair indices in parallel for each thread. The arrays for each thread will be sorted (irrespective of the scheduling).

The following function merges two sorted arrays

void merge(int *a, int *b, int*c, int na, int nb) {
    int i=0, j=0, k=0;
    while(i<na && j<nb) c[k++] = a[i] < b[j] ? a[i++] : b[j++];
    while(i<na) c[k++] = a[i++];
    while(j<nb) c[k++] = b[j++];
}

Here is the remaining code

uint nind = 0;
uint *P;
#pragma omp parallel
{
    uint myP[x*z];
    uint mynind = 0;
    #pragma omp for schedule(dynamic) nowait
    for(uint k = 0 ; k < x*z; k++) {
        if (inFunc(p, index)) myP[mynind++] = k;
    }
    #pragma omp critical
    {
        uint *t = (uint*)malloc(sizeof *P * (nind+mynind));
        merge(P, myP, t, nind, mynind);
        free(P);
        P = t;
        nind += mynind;
    }
}

Then given an index k in P the pair is (k/z, k%z).

The merging can be improved. Right now it goes at O(omp_get_num_threads()) but it could be done in O(log2(omp_get_num_threads())). I did not bother with this.

---

Hristo Iliev's pointed out that dynamic scheduling does not guarantee that the iterations per thread increase monotonically. I think in practice they are but it's not guaranteed in principle.

If you want to be 100% sure that the iterations increase monotonically you can implement dynamic scheduling by hand.

Answer 3

The code you provide looks like you are trying to fill the XY data in sequential order. In this case OMP multithreading is probably not the tool for the job as threads (in a best case) should avoid communication as much as possible. You could introduce an atomic counter, but then again, it is probably going to be faster just doing it sequentially.

Also what do you want to achieve by optimizing it? The x and z are not too big, so I doubt that you will get a substantial speed increase even if you reformulate your problem in a parallel fashion.

If you do want parallel execution - map your indexes to the array, e.g. (not tested, but should do)

#pragma omp parallel for shared(XY)
for (uint i = 0; i < x; i++) {
   for (uint j = 0; j < z; j++) {
           if (inFunc(p, index)) {
                uint idx = (2 * i) * x + 2 * j; 
                XY[idx] = i;
                XY[idx + 1] = j;
           }
   }
}

However, you will have gaps in your array XY then. Which may or may not be a problem for you.