PGI User Forum
 SearchSearch   MemberlistMemberlist     RegisterRegister   ProfileProfile    Log inLog in 

Free OpenACC Webinar

Mersenne Twister not giving uniform random number distributi

 
Post new topic   Reply to topic    PGI User Forum Forum Index -> Accelerator Programming
View previous topic :: View next topic  
Author Message
alfvenwave



Joined: 08 Apr 2010
Posts: 79

PostPosted: Tue May 04, 2010 8:27 am    Post subject: Mersenne Twister not giving uniform random number distributi Reply with quote

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.

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


Code:
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


Back to top
View user's profile
alfvenwave



Joined: 08 Apr 2010
Posts: 79

PostPosted: Wed May 05, 2010 1:25 am    Post subject: Reply with quote

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]
Back to top
View user's profile
mkcolg



Joined: 30 Jun 2004
Posts: 6211
Location: The Portland Group Inc.

PostPosted: Wed May 05, 2010 10:44 am    Post subject: Reply with quote

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
Back to top
View user's profile
alfvenwave



Joined: 08 Apr 2010
Posts: 79

PostPosted: Wed May 05, 2010 12:42 pm    Post subject: Reply with quote

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

Rob.
Back to top
View user's profile
Display posts from previous:   
Post new topic   Reply to topic    PGI User Forum Forum Index -> Accelerator Programming All times are GMT - 7 Hours
Page 1 of 1

 
Jump to:  
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum


Powered by phpBB © phpBB Group