openmp bug with radial computation

Question

#pragma omp parallel for schedule(static) default(none)
for(int row = 0; row < m_height; row++)
{
    for(int col = 0; col < m_width; col++)
    {
       int RySqr, RxSqr;

       SettingSigmaN(eta, m_RxInitial + col, m_RyInitial + row , RxSqr, RySqr);

       FunctionUsing(RySqr,RxSqr); 
    }
} 

void CImagePro::SettingSigmaN(int Eta, int x, int y, int &RxSqr, int &RySqr, int &returnValue)
{
    int rSqr = GetRadius(x,y,RxSqr,RySqr);
    returnValue = GetNumberFromTable(rsqr);
}


int CImagePro::GetRadius(int x, int y, int &RxSqr, int &RySqr) 
{   
    if (x == m_RxInitial)
    {       
        RxSqr = m_RxSqrInitial;
        if (y == m_RyInitial) 
        {
            RySqr = m_RySqrInitial;
        }
        else if ( abs(y) % 2 == abs(m_RyInitial) % 2)
        {
            RySqr = RySqr + (y<<2) + 4; //(y+2)^2
        }
    } 
    else 
    {
        RxSqr = RxSqr + ( x << 1) + 1; //(x+1)^2
    }   
    return clamp(((RxSqr+RySqr)>>RAD_RES_REDUCTION),0,(1<<(RAD_RES-RAD_RES_REDUCTION))-1);
}

ok so here is my code and my problem is inside GetRadius Function. since i have many threads each threads starts at a different place of x,y. however i don't really understand where is the bug inside GetRadius().

i thought maybe is it the RySqr computation. can you suggest a way to debug? or can you see my problem?

UPDATE:

this has fixed most of my code: i still don't really understand, why the there are jumps between the different threads.

  int CImagePro::GetRadius(int x, int y, int &RxSqr, int &RySqr) 
  { 
     if (x == m_RxInitial)
     {      
        RxSqr = m_RxSqrInitial;
     }
     else
     {
         RxSqr = x * x;
     }

    if (y == m_RyInitial)
    {
        RySqr = m_RySqrInitial;
    }
    else if (abs(y) % 2 == abs(m_RyInitial) % 2)
    {
        RySqr = y * y;
    }
    return clamp(( (RxSqr + RySqr) >> RAD_RES_REDUCTION), 0, ( 1 << (RAD_RES - RAD_RES_REDUCTION) ) - 1);   
}

Answer 1

I really wonder if this thing compiles? You specify default(none), but consistently use data members of your class. Are they all static?

What you could do is either i) leave default(none) away, which means default(shared), ii) have a shared access to the values by explicitly sharing them, or iii) initialise the variables you use inside the parallel region so that each thread has it's own private copy of, say, m_RxInitial called p_RxInitial etc. The first option is almost guaranteed to get you into trouble.

Following illustrates option ii):

1) Make a helper class containing everything you need to pass, for you this could be

struct ShareData{
    int s_RxInitial
    /* ... */
}

2) In the member function containing parallel section, before parallel loop define

ShareData SD;
SD.s_RxInitial = m_RxInitial;
/* ... */

3) Give it to the parallel section

#pragma omp parallel for schedule(static), default(none), shared(SD)

4) Use the SD datamembers in function calls.

I hope this was clear enough. I would appreciate it if someone had a more elegant solution to offer.

If you wanted private variables of option iii), you could say firstprivate(SD) instead of shared(SD). This would give each thread initialized (to the original values) private copy of SD. It may or may not give some performance advantage by avoiding serial access. I had a similar problem few days ago and there was no difference.

Answer 2

You cannot guarantee the order in which the threads execute, if you need to guarantee that either make if statements as you did or simply don't parallelize it because it is a critical section.

http://bisqwit.iki.fi/story/howto/openmp/