Monte-Carlo simulations

Random numbers generators

The Linear Congruential random number generator (LC RND) kernel

__kernel void kRand(global uint* rand ){
  int x = get_global_id(0);
  rand[x] = 1664525u*rand[x] + 1013904223u;
}
stores its state in a global buffer array. You see below random number generator based on the Tiny Encryption Algorithm (TEA) made by Tomi Aarnio.
uint2 rand(uint2 seed, int iterations) {
  uint sum = 0;
  uint delta = 0x9E3779B9;
  uint k[4] = { 0xA341316C, 0xC8013EA4, 0xAD90777D, 0x7E95761E };
  for (int j=0; j < iterations; j++) {
    sum += delta;
    seed.x += ((seed.y << 4) + k[0]) & (seed.y + sum) & ((seed.y >> 5) + k[1]);
    seed.y += ((seed.x << 4) + k[2]) & (seed.x + sum) & ((seed.x >> 5) + k[3]);
  }
  return seed;
}
It generates a pair of 32-bit random numbers from the given pair of 32-bit seed values seed and t (see below)
  int ix = get_global_id(0),  iy = get_global_id(1),  t = ix + n*iy;
  uint2 rnd = (seed, seed << 3);
  rnd.x += t + (t << 11) + (t << 19);
  rnd.y += t + (t << 9) + (t << 21);
  rnd = rand(rnd, iter);
The number of iterations iter can be adjusted depending on the required degree of randomness: 16 iterations is typically enough for any use.

For iter = 8 XY model TEA RNG is ~3 times slower than XY model LC RNG on AMD E-350.

Look at [2] for comparison of LC and TEA pseudo-random number generators.

Modified chess-board algorithm

Direct Chess-Board algorithm (CBA) is used in Ising model CBA. At first all white independant spins are calculated in parallel then all grey ones in the second pass (see Fig.1a and [3]). In CBA only 1/2 of shaders ("grey" or "white") are used in every pass.
constant int n = 512;
constant float h = 1.0f/512.0f;

__kernel void kSpin(read_only image2d_t s, global uint* rand,
  write_only image2d_t s1, int bw, constant float* w ) {
  const sampler_t samp =
    CLK_NORMALIZED_COORDS_TRUE|CLK_ADDRESS_REPEAT|CLK_FILTER_NEAREST;
  int ix = get_global_id(0),  iy = get_global_id(1);
  if ( ((ix + iy) & 1) == bw) return;
  int t = ix + n*iy;
  uint r = 1664525u * rand[t] + 1013904223u;
  rand[t] = r;
  float x = h*(ix + .5f),  y = h*(iy + .5f);
  int sum = read_imagei(s, samp, (float2)(x + h, y)).x;
  sum += read_imagei(s, samp, (float2)(x - h, y)).x;
  sum += read_imagei(s, samp, (float2)(x, y + h)).x;
  sum += read_imagei(s, samp, (float2)(x, y - h)).x;
  int c = (r < w[sum]) ? 0 : 1;
  write_imagei(s1, (int2)(ix, iy), c );
}
To avoid spins processing in the chess-board order in the Modified Chess-Board algorithm (MCBA) the odd columns are displaced as shown in Fig.1 (so you have to calculate neighbours with accurancy). Furthermore "grey" and "white" spins (texels) are stored in two different textures. Therefore MCBA is 2 timed faster than CBA. See Ising model MCBA.

[1] GPU Gems 3. Efficient Random Number Generation and Application Using CUDA.
[2] Pseudo-random number generation for Brownian Dynamics and Dissipative Particle Dynamics simulations on GPU devices
Carolyn L. Phillips, Joshua A. Anderson, Sharon C. Glotzer
Journal of Computational Physics 230 (2011) 7191-7201
[3] CPS 713 Monte Carlo Simulation for Statistical Physics
Paul Coddington, Northeast Parallel Architectures Center at Syracuse University,   Jan 1996


WebCL Demos     updated 18 Sep 2011