Mersenne Twister not giving uniform random number distributi

Accelerated Computing HPC Compilers Legacy PGI Compilers
1

Hi.

I can’t work out what I’ve done wrong here. I have simply followed the recipe for generating random numbers from the pginsider Monte Carlo example. In my code below, all I am trying to do is generate 4096 ranom numbers and write them to a file to check they look OK. The results I get depend on the number I use for the seed. With the value of 777 from the example, I don’t get any numbers above 0.25. Can anyone spot a stupid mistake?

Also, the CUDA Fortran Programming Guide and Reference seems to suggest that the random_seed and random_number F90 intrinsics are available for device subprograms…

Rob.

pgfortran -c mersenne.cuf 
nvcc -O3 -c -Iinc MersenneTwister_kernel.cu
pgfortran -fast -Mcuda mersenne.o MersenneTwister_kernel.o




program random

use cudafor

interface

    ! C function to load the Mersenne Twister input data set
    subroutine loadMTGPU(dat_path) bind(c,name='loadMTGPU')
       use iso_c_binding
       character(len=80) :: dat_path
    end subroutine loadMTGPU

    ! C function to initialize the random seed on the GPU
    subroutine seedMTGPU(seed) bind(c,name='seedMTGPU')
       use iso_c_binding
        integer(c_int), value :: seed   ! Pass by value
    end subroutine seedMTGPU

    ! CUDA C Mersenne Twister algorithm to generate an array of random
    ! numbers in parallel.  Note that nvcc has added '__entry' to the
    ! end of the function name.
    subroutine RandomGPU(d_Rand, n_per_rng) bind(c,name='randomgpu__entry')
        use iso_c_binding
        real(c_float), dimension(:), device :: d_Rand
        integer, value :: n_per_rng  ! pass by value
    end subroutine RandomGPU

end interface

character(len=80)                     :: dpath
real, dimension(32*128*1),device      :: dXD
real, dimension(32*128*1)             :: dXH
integer                               :: iflag

iflag = cudaSetDevice(0)
dXH   = 0.0
dXD   = 0.0

write(dpath,'(a)') 'MersenneTwister.dat'//char(0)
call loadMTGPU(dpath)
call seedmtgpu(777)
call randomgpu<<<32,128>>>(dXD,1)

dXH = dXD

open(5,file='test.txt')
     write(5,*) dXH
close(5)

end
2

I have downloaded the Mersenne Twister example, and modified the code to write out the first 10e6 random numbers filled in dX. If I plot the values in the order in which they reside in dX (along the bottom), value on the y axis (as a scatter plot), I would hope to see a uniformly illuminated rectangle with a minimum of 0.0 and maximum of 1.0 on the y-axis. Unfortunately, I see a lot of structure. The structure gets less at the higher end of the array, but it’s still visible. The random number generator is “random enough” to give a sensible value of pi, but it’s not random enough for many apps.


I have submitted images via e-mail.

Rob.
[/img]

3

Hi Rob,

I took the Mersenne Twister code directly from the NVIDIA CUDA SDK examples (http://developer.download.nvidia.com/compute/cuda/sdk/website/samples.html#MersenneTwister) so unfortunately don’t have any insight into the poor distribution at smaller values of N. It could be a configuration error on my part or a problem with the algorithm. I can investigate this further, but it may be a while. Perhaps the author could help?

  • Mat
4

Thanks Mat. I have e-mailed Victor Podlozhnyuk to see if he can help. I’ll let you know what I find out…

Rob.