Introduction

Monte Carlo algorithms have been heavily used in transport applications since the earliest days of electronic computing. In spite of this very long history, little has changed in how random number generators are used in these applications. Random numbers have typically been generated with a linear congruential recurrence formula of the form

x[n] = a x[n - 1] + b   (modm)   , (1)

where x is the hidden random number state stored in statically allocated memory. The constants a and b are chosen to produce a maximal period and good randomness properties. The modulus of the arithmetic, m, is typically chosen to be a power-of-two in order to minimize the cost of performing modular arithmetic using binary operations on a computer. The code developer is usually provided with routines that allow one to save and restore the internal state x of the generator. These routines are used for job check pointing purposes and can also be used to set the starting state for a statistically independent instance of the program run.

The increasing popularity of shared memory parallel computing has exposed a hazard in the traditional approach when the internal state, x, is shared among the threads or processors executing the application. The internal state can be protected with a lock, removing the read/modify/write hazard, but such an arrangement still results in non-deterministic (but correct) behavior because the processors can access the random number generator in any order influenced by machine time sharing effects. One must localize the internal random number state, on a per-thread or per-processor basis, in order to remove this hazard. The coding required to localize memory in this manner is specific to a given parallel programming environment, leading to a multiplicity of random number library implementations (possibly on a given shared memory multiprocessor) and the attendant risk of loading against the wrong one.

Further complications arise in a distributed memory programming environment where the number of processors employed can be large. The straight forward extension of the per-thread or per-processor approach has led to the development of random number generators with parameterized constants that provide independent full period sequences [1], assuring the independence of random numbers generated in each processor. Although this solution provides independent generators for each processor, it can lead to problems with determinism when a particle must be moved between processors for load balancing purposes. When the particle is moved, it encounters a different random number stream and its history is changed as a result. More significantly, its presence affects the history of all of the particles being processed by the processor that it is moved to, leading to undesired avalanche effects in application behavior.

In the RNG library described in this paper, we take an alternative approach of efficiently supporting a per-particle random number state. The random number generator library is stateless, leading to a single implementation for serial, shared memory parallel, and distributed memory applications. The user may use the library to implement a single random number state, a per-processor or per-thread random number state for parallel applications, or a per-particle random number state in the limit of fine granularity. A pointer to the desired random number state is supplied as an argument to the random number routines provided in the library.