Strange behavior in a parallelized loop

[yurytsoy]

Hi, everybody!

There is a code, which works very strange to me (this is to model some situation elsewhere, so it doesn't have much sense, but it works in an unexpected manner):

Code: Select all

#include <iostream>
#include <vector>

using namespace std;
const int SIZE = 5000;

class Fitness {
public:
   float value;

   Fitness (float val) {
      value = val;
   }
};

Fitness* func2(float* mas) {
   float res = 0;
   for (int i=0; i<SIZE; ++i) {
      res += mas[i];
   }
   return new Fitness (res);
}

void func(float**);

int main () {
   float** mas = new float*[SIZE];
   for (int i=0; i<SIZE; ++i) {
      mas[i] = new float[SIZE];
      for (int j=0; j<SIZE; ++j) {
         mas[i][j] = i;   // each row simply contains all elements equal to the row's index.
      }
   }

   func (mas);

   for (int i=0; i<SIZE; ++i)
      delete[] mas[i];
   delete[] mas;

   return 0;
}

void func (float** mas) {
   int size = 10;
   Fitness** fits = new Fitness*[SIZE];

#pragma omp parallel for shared(fits, SIZE, mas)
   for (int i = 0; i < SIZE; i++) {
      fits[i] = func2(mas[i]);
#pragma omp critical
      {
         if (i != 0) {
            cout<<i<<"\t"<<fits[i]->value<<"\t"<<fits[i]->value / i<<"\n";    // should always print SIZE in the last column
         } else {
            cout<<i<<"\t"<<fits[i]->value<<"\n";
         }
      }
   }

   for (int i=0; i<SIZE; ++i) delete fits[i];
   delete[] fits;
}


The problem is that the line for checking correctness of the func function doesn't always prints SIZE, but sometimes (~15-20%) something very close to it, like 4999.9 or 5000.1 (if SIZE is 5000). The problem doesn't emerge when SIZE is rather small (10, 100, 1000), but it's there for SIZE > 5000. I can't understand it at all neither can find any potential explanation.

Could anyone please help me out what I'm missing here? Thanks a lot!

Yury

[MarkB]

This is just a loss of precision problem in floating point arithmetic: by the time the row sums are large enough (around 10^7), they can no longer be represented at full accuracy in a float type. If you use doubles instead of floats you should see the problem go away (at least for SIZE = 5000). This is nothing to do with OpenMP or parallelism, by the way: the problem is just the same in the sequential code.

[yurytsoy]

Oops, that was really a lame mistake of mine. Thanks!