Advanced Topic: Seeds For Different Threads

Adding OpenMP pragmas on the ‘workhorse’ for loops where most of the computation is being done is often a helpful way to make your code run faster. In the case of Monte Carlo simulations, there is one issue that should be addressed to ensure the best random distribution of numbers from the random number generator functions. We must start each thread with a different seed.

Recall that random number generators start from a ‘seed’ large integer and create a sequence of integers by permuting the seed and each successive integer in a manner that ensures they are distributed across the range of all integers. The key point is this: the sequence of numbers from a random generator is always the same when it starts with the same seed. In code where we fork threads to do the work of generating random numbers, we lose the desired random distribution if each thread begins generating random numbers from the same seed.

The solution to this issue for threaded code, which you can download as coinFlip_omp_seeds.cpp, is to ensure that each thread has its own seed from which it begins generating its sequence of integers. Let’s revisit the coin flipping example. Instead of generating one seed in main using time(), we can save a seed for each thread in an array and devise a function to create all of the seeds, based on the number of threads to run. We can add this code at the beginning of our original file:

/***  OMP ***/
const int nThreads = 4;  // number of threads to use
unsigned int seeds[nThreads];

void seedThreads() {
    int my_thread_id;
    unsigned int seed;
    #pragma omp parallel private (seed, my_thread_id)
    {
        my_thread_id = omp_get_thread_num();
        
        //create seed on thread using current time
        unsigned int seed = (unsigned) time(NULL);
        
        //munge the seed using our thread number so that each thread has its
        //own unique seed, therefore ensuring it will generate a different set of numbers
        seeds[my_thread_id] = (seed & 0xFFFFFFF0) | (my_thread_id + 1);
        
        printf("Thread %d has seed %u\n", my_thread_id, seeds[my_thread_id]);
    }
    
}
/***  OMP ***/

Not how we change the seed value for each thread by using the thread’s id to manipulate the original integer obtained from time().

Then later in the main function, we add a call to this function:

    /***  OMP ***/ 
    omp_set_num_threads(nThreads);  
    seedThreads();
    /***  OMP ***/ 

For each trial, we still parallelize the workhorse for loop, while also ensuring that each thread running concurrently has its own seed as the starting point for later numbers.

#pragma omp parallel num_threads(nThreads) default(none) \
        private(numFlips, tid, seed) \
        shared(trialFlips, seeds) \
        reduction(+:numHeads, numTails)
    {
        tid = omp_get_thread_num();   // my thread id
        seed = seeds[tid];            // it is much faster to keep a private copy of our seed
		srand(seed);	              //seed rand_r or rand
        
        #pragma omp for
        for (numFlips=0; numFlips<trialFlips; numFlips++) {
//          in Windows, can use rand()
//            if (rand()%2 == 0) // if random number is even, call it heads
            // linux: rand_r() is thread safe, to be run on separate threads concurrently
            if (rand_r(&seed)%2 == 0) // if random number is even, call it heads
                numHeads++;       
            else
                numTails++;
        }
        
    }

Study the above code carefully and compare it to our first version below. The pragma omp directive above is forking the new set of threads, which do a bit of work to set up their own seeds. Then the pragma omp for directive is indicating that those same threads should now split up the work of the for loop, just as in our previous example using the OpenMP pragma. The first OpenMP version we showed you looked like this:

#pragma omp parallel for num_threads(nThreads) default(none) \
        private(numFlips, seed) \
        shared(trialFlips) \
        reduction(+:numHeads, numTails)
        for (numFlips=0; numFlips<trialFlips; numFlips++) {
            // rand() is not thread safe in linux
            // rand_r() is available in linux and thread safe,
            // to be run on separate threads concurrently.
            // On windows in visual studio, use rand(), which is thread safe.
            if (rand_r(&seed)%2 == 0) // if random number is even, call it heads
                numHeads++;       
            else
                numTails++;
        }

Note

A common ‘gotcha’ that can cause trouble is if you accidentally use the original pragma omp parallel for directive near the for loop in the new version. This causes incorrect unintended behavior. Remember to remove the parallel keyword in the inner block when nesting bloaks as shown in the new version where we set up seeds first before splitting the loop work.

Note that as before, in linux we need to use the rand_r() function for thread-safe generation of the numbers. However, in Windows, the rand() function is thread-safe.

Try this yourself

Try creating versions of the Roulette wheel simulation or drawing four suits that ensure that each thread is generating numbers from its own seed.