3

I am trying to use fortran intrinsic PRNG in a MPI code.

I understand from this link that GFortran implement the PRNG using xorshift1024* which has a period of 2^1024 - 1. It also says:

Note that in a multi-threaded program (e.g. using OpenMP directives), each thread will have its own random number state.

Then reading this I found:

When a new thread uses RANDOM_NUMBER for the first time, the seed is copied from the master seed, and forwarded N * 2^512 steps to guarantee that the random stream does not alias any other stream in the system, where N is the number of threads that have used RANDOM_NUMBER so far during the program execution

If this is an automatic feature of GFortran, it works only in OpenMP? What if I want to have parallel PRNG using MPI? How can I ensure the portability of the code to other compilers?

In other words: Is there any way to do what GFortran says it does (i.e. guarantee real parallel PRNG) in a portable way using the fortran intrinsic instructions?

NOTE: I was using the PRNG of Numerical Recipes in MPI. That worked well for some years, but now I am getting some errors in some assumptions on the integer model that Numerical Recipes says goes beyond fortran... so I don't see how to solve that and that is way I want to use the intrinsic PRNG if is possible.

alexis
  • 410
  • 4
  • 18
  • If you care about portability/consistency across compilers (even compiler versions) then `random_number` (and `random_seed`) will present many problems. – francescalus Sep 08 '16 at 15:12
  • *Is there any way to do what GFortran says it does ... in a portable way using the fortran intrinsic instructions?* No. This -- http://stackoverflow.com/questions/8920411/possible-sources-for-random-number-seeds -- may be useful. – High Performance Mark Sep 08 '16 at 16:34
  • Thank you for your comments. If the intrinsic are not portable and the Numerical Recipes seems to be not portable too (in my experience). What other libraries, subroutines, or methods do you suggest? – alexis Sep 08 '16 at 18:04

2 Answers2

3

If you want to have a portable version of a multi-stream random number generator for a Fortran program, there is a multi-stream Fortran version of the Mersenne Twister. See http://theo.phys.sci.hiroshima-u.ac.jp/~ishikawa/PRNG/mt_stream_en.html . It uses the concept of advancing the PRNG by a very large number of steps for the different threads. It's setup and configured by subroutine calls so you should be able to use it from various multi-threading environments.

M. S. B.
  • 28,968
  • 2
  • 46
  • 73
3

Note that the use of xorshoft1024* is a very new feature in GFortran, it's only available in the development trunk version, no released version has it yet at the time of writing this. It will be released as part of GCC 7, probably in spring 2017.

So when you're using MPI, each MPI rank is a separate process and the random number generators in each process is completely separate with no communication between the PRNG's in different processes (unless you handle it yourself with MPI, of course). The thing with forwarding the PRNG stream 2^512 steps happens only when using the PRNG from multiple threads within the same process.

That being said, xorshift1024* has a fairly long period (2^1024-1), and the first time the PRNG is used in a process (again, think MPI rank) it is initialized with random data from the OS (/dev/urandom on POSIX systems), unless it has already been explicitly initialized with RANDOM_SEED. So in practice I think you'll be fine, it's exceedingly unlikely that the PRNG streams for different MPI ranks will alias.

And no, the above describes the PRNG in GFortran version 7. If you want something portable you cannot rely on anything beyond what the standard guarantees. Beyond the parallel aspects, for portable high quality random numbers you're probably better of with using a known good PRNG rather than relying on the one provided by the compiler (I have personal experience of at least one compiler producing poor quality random numbers with the RANDOM_NUMBER intrinsic, but I'll refrain from naming the vendor since it was many years ago and they might have fixed it since, if they are even in business anymore, I don't know).

(If you find the semantics of the new xorshift1024* implementation difficult, blame a) me, since I devised it and implemented it b) the Fortran standard which makes it impossible to have a parallel PRNG with simple semantics)

janneb
  • 36,249
  • 2
  • 81
  • 97
  • Thank you a lot @janneb. Could you suggest a known good and portable PRNG in fortran? – alexis Sep 09 '16 at 14:10
  • Well, I think xorshift1024* ia a fairly good choice. The reference implementation is in C though, and uses unsigned 64 bit arithmetic which is hard to emulate in Fortran. I suggest using the C implementation via ISO_C_BINDING. – janneb Sep 10 '16 at 18:20
  • 1
    This is an old thread, but I wanted to update it to note that the xorshift1024* was the default RNG in Gfortran version 8. I'm not sure when it changed, but in version 10, the RNG has been changed to xoshiro256** – NuclearFission Mar 04 '23 at 01:51
  • @NuclearFission: I did that change in August 2019: https://gcc.gnu.org/git/?p=gcc.git&a=commit;h=0e99e0933984e0c30fda1d089bfbd6857fc9273f – janneb Mar 04 '23 at 12:00