Is there a straightforward way to generate random numbers inside an OpenACC parallel loop? I want the seed within each loop/thread to be different. The reference code that demonstrates the required task to be done is as follows:
PROGRAM genRandNums
implicit none
integer :: i
double precision :: a(100)
!$acc parallel loop
DO i = 1,100
call random_seed()
call random_number(a(i))
ENDDO
END PROGRAM
This code works on the gnu’s 'gfortran ’ but throws an error with the 'nvfortran ’ compiler as the intersic functions random_seed() and random_number() are not supported. How can I modify this code to work with the nvfortran compiler?
In general, running a RNG in parallel is problematic. RNGs contain state which is often shared, making the RNG unsafe to parallelize. Instead, each iteration of the parallel loop would need to maintain it’s own state in order to avoid race conditions.
You can do this with the device side cuRAND (examples include with the compilers under the “examples/CUDA-Libraries/cuRAND” directory), the cost of maintaining state for each iteration is high. Plus, you need to pass in a set of randomly generated seeds so each instance of the RNG is unique. Since your only using one random number per iteration, you’re better off calling cuRAND from the host to generate an array of random values and then using this array in the parallel loop.
I had the opportunity to work with Johan Carlsson on a pure OpenACC device side RNG implementation. Like cuRAND, you’d want to use it if your generating many random number per loop iteration. For cases where you’re using one random number per iteration, its still better to precompute the random values. Though unlike cuRAND Johan’s DES PRNG implementation is much lighter weight so has less overhead. For full details on DES PRNG see: Pseudo Random Number Generation by Lightweight Threads | OpenACC
I’ve only used DES PRNG with C/C++ programs, but presume if you add a Fortran interface module, should be able to be used within Fortran as well.
Hi Mat. I took note of your reply regarding using DESPRNG. I am using nvfortran, so would have to use a Fortran interface module. I have been looking into this, and found the following statement in the book by Metcalf, Reid and Cohen on “Modern Fortran”, p 321, - “No Fortran type is interoperable with…a C structure type that contains a bit field…” Unfortunately, it looks like DESPRNG does contain such a field, “bytebit”. Am I correct? If so, is there any way around this? Thanks. Malcolm.
I believe you’re correct, but I’m not sure it’s relevant. The interfaces themselves are passing C pointers to the array of structures, not the structures themselves so you shouldn’t need to create an equivalent Fortran type.
You’d likely need to add a routine to the DESPRNG library to perform the allocation, but again only the C pointer needs to be part of the interface.
I wish I had time to write the interface myself (it’s something I’ve wanted to do), but please let me know if you encounter any issue and it would be great if you can share your work.
Mat, thanks for the fast response. I was intending to use ISO_C_BINDING to call DESPRNG. Is this what you are thinking about when you say I shouldn’t need to create an equivalent Fortran type? So how would I call the routine - where would the pointer point to? Also, when you suggest that I’d likely need to add a routine to the DESPRNG library to perform the allocation, what allocation are you talking about? You obviously have a certain approach in mind, do you care to elaborate? If I am successful, I will certainly share the results. Malcolm.
You’d use the ISO_C_BINDING’s “C_PTR” type in the interface. A C pointer is just an unsigned long variable whose value is an address so it can be passed from Fortran to C and back.
/* Make some workspace on the stack for the DES PRNGs */
nident = alloca(8 * Npart);
thread_data = alloca(sizeof(desprng_individual_t) * Npart);
process_data = alloca(sizeof(desprng_common_t));
xi = alloca(8 * Npart);
Since the Fortran won’t know the size of the structs, you’ll want to move this to a routine in DESPNG. Though change alloca to malloc so the allocation is on the heap and not the stack.
Of course, I haven’t actually done the port so there could be issues, but it’s the approach I’d take.
Effectively you want to move all the allocation and data management to the DESPNG code so it’s self contained and all the Fortran side needs to do is pass the pointer handles to the structs between the calls to DESPNG.
Hi Mat, I am making good(!) progress with a Fortran interface to DESPRNG.c. I need some clarification on the intent of some variables. In particular, in the call to initialize_desprng what are “seed, key, length” looking for? Similarly, what are " output, length" in generate_random_numbers? Are the two “lengths” the same? I have my own ideas but don’t want to waste time guessing! Malcolm.
Sorry but I don’t see an “initialize_desprng” or “generate_random_numbers” routine in Johan’s code. Maybe I’m missing them, but wanted to double check we’re looking at the same code:
The routines I see are “initialize_common”, “initialize_individual”, and “get_uniform_prn” but none have the args you mention.