I am working with Monte Carlo simulations to find the decimal places of PI. So far so good but OpenMP came in and I realize that ran2, arguably the best RGN, is not threadsafe! The implementation is here. Since I have not worked with OpenMP and neither a lot on multi-threading I am stuck at making this thread safe using OpenMP.
So far what I know is that a function is already thread-safe if it doesn't modify non-local memory and it doesn't call any function that does. In this case, there are 3 variables which are static and thus will be modified if gets used by different threads.
One possible solution is to call this in a thread safe way by enclosing the calling of ran2 in a critical section but that makes no sense as I get no speedup.
Can somebody give me pointers on how to proceed with this or if somebody has any reference that will be great too!
What is generally done to render a procedure thread safe is to associate previously static data to a thread local data. Look for instance at the man page of rand_r() that is a thread safe version of rand().
So in your version of L'Ecuyer:
define a struct (say state) that holds the static data
redefine procedure ran2() to have an additional parameter that is a pointer to this struct state and modify the code accordingly. Let ran2_r() be the new name.
define in every thread a local struct state to hold the state
probably state needs to be seeded. You can use get_thread_num() to provide a per thread seed to initialize properly the state when entering the thread.
Now you just need to call your new ran2_r() with a pointer to this state. It will be modified by the procedure, but modifications will be stored in the thread local state var.