Usage of vsrnggaussian() - very easy question

[johansoderstrom]

Dear FORTRAN experts

I am stuck at a very beginner problem with FORTRAN and I hope that you can help me. I want to get a zero-mean, unit-varience (sigma=1) Gaussian-distributed number. I found the function vsrnggaussian() which should do the trick, but I have no idea of how to implement it in my FORTAN code. I found an example (http://software.intel.com/en-us/forums/showthread.php?t=48992) but that seems to be not for FORTRAN.

So, in essance, how do I go about to get such a random (non-integer) number?

As you probably see I am very new to FORTRAN so I have no idea of how to initialize streams, and how to de-initialize them... If someone whould be kind enough to provide me with a short example I would be very happy and grateful

Thank you so much,

Johan S

[TimP]

Quoting - johansoderstrom

Dear FORTRAN experts

I am stuck at a very beginner problem with FORTRAN and I hope that you can help me. I want to get a zero-mean, unit-varience (sigma=1) Gaussian-distributed number. I found the function vsrnggaussian() which should do the trick, but I have no idea of how to implement it in my FORTAN code. I found an example (http://software.intel.com/en-us/forums/showthread.php?t=48992) but that seems to be not for FORTRAN.

The Fortran interface to VSL is described in the mklman.pdf in the MKL documentation. As this is part of the MKL library, you would be more likely to get expert answers on the MKL forum.

[Steven_L_Intel1]

I'll move this to the MKL forum (now that I can!)

[anthonyrichards]

Quoting - Steve Lionel (Intel)

I'll move this to the MKL forum (now that I can!)

Here is a simple program to generate real numbers having a Guassian distribution with standard deviation SIGMA. It uses standard routines to generate a pseudorandom sequence of numbers uniform in the range 0 to 1.

!****************************************************************************
!
! PROGRAM: gaussianrand
!
! PURPOSE: generate gaussian random numbers with given sigma
! Uses FUNCTION RAND(0) to generate uniform random numbers
! between 0 and 1, and ROUTINE SEED to set the seed
! Uses the Box-Muller algorithm to generate numbers with
! a Gaussian distribution of probability.
!
!****************************************************************************
program gaussianrandimplicit none

integer,parameter:: nrand=1000000
REAL*8 UNIF1(nrand), UNIF2(nrand)
!
REAL*8 Y1,Y2,SIGMA,binsize, binhalf
INTEGER I,nbins, bins(201),ix

! supply seed to get the same pseudorandom sequence each time
! program is run. To get a different sequence, change the seed
! by calling SEED with a different integer value.
call seed(1)
! set the sigma value you want
SIGMA=2
!
! get the random numbers. RGAUSS returns two gaussian random numbers
! for each call.
!
DO I=1,NRAND
CALL RGAUSS(SIGMA,Y1,Y2)
UNIF1(I)=Y1
UNIF2(I)=Y2
END DO
!
! there are NRAND values in each array UNIF1 and UNIF2, each with
! a Gaussian distribution with standard deviation SIGMA
!
! bin the results, bin size=0.05*sigma
nbins=201
binsize=0.05*sigma
binhalf=0.5*binsize

do i=1,nrand
! compute which bin the number goes into
ix=101+floor((unif1(i)+binhalf)/binsize)
if(ix.gt.0.and.ix.lt.201) bins(ix)=bins(ix)+1
end do

! save the numbers in each bin to display later
open(unit=10,file='rgauss.dat',form='formatted',status='unknown')
write(10,1000) (bins(i),i=1,201)
1000 format(10i8)
close(10)

end program gaussianrand

subroutine rgauss(sigma, y1,y2)
real*8 x1, x2, w, y1, y2, sigma

do while ( (w .ge. 1.0).or.(w.eq.0) )
x1 = 2.0 * rand(0) - 1.0
x2 = 2.0 * rand(0) - 1.0
w = x1 * x1 + x2 * x2
end do

w = sigma*sqrt( (-2.0 * log( w ) ) / w )
y1 = x1 * w
y2 = x2 * w

end subroutine rgauss

[anthonyrichards]

I wanted to edit my post to correct pasting error - but I cannot find an edit button

Anyway, please note the error in

gaussianrandimplicit none

which should of course read

Program gaussianrand

implicit none

[davidspurr]

Johan - following is a snipped of code from a test program I used to check different random methods. I have cut down to code to make it as simple as possible. Hope I have not introduced any errors. Yes - the descriptions in the documentation are confusing, but actual application is very simple. (Steve - not certain if this is the right forum, but have posted here in to increase chance of Johan finding it)

USE mkl_vsl
IMPLICIT NONE

integer status !Error / status
integer brng !Index of generator to initialize stream
integer method !Method = 0 in case of uniform distribution
integer dimen !Initial condition of the stream
TYPE (VSL_STREAM_STATE) stream !Random stream
integer, PARAMETER :: nr = 1000 !Number of values required (whatever you want)
real gauss(nr)
real mean, sigma

mean = 0.0
sigma = 1.0

!SOBEL
method = VSL_METHOD_SGAUSSIAN_ICDF !NB - are different alternatives, but I found this worked well
brng = VSL_BRNG_SOBOL
dimen = 1 ! <== different values may be needed depending on application
gauss = 0.0
status = vslnewstream( stream, brng, dimen ) !Create new stream (SOBOL in this case
status = vsrnggaussian( method, stream, nr, gauss, mean, sigma )
status = vsldeletestream( stream ) !Delete the stream


!For alternate random types, substitute
! brng = VSL_BRNG_MT19937
! brng = VSL_BRNG_MRG32K3A
!NB: I found VSL_BRNG_MT19937 overall best for my application

OT: "Insert Code using SyntaxHighlighter" facility has no option for Fortran !!! (in a Fortran forum)

[johansoderstrom]

I found all the answers (I think) and will try to implement them in my code as soon as I can. Thanks a lot everyone! Once I get it working (or not) I will let you know :)

~Johan

[Andrey_N_Intel]

Quoting - johansoderstrom

I found all the answers (I think) and will try to implement them in my code as soon as I can. Thanks a lot everyone! Once I get it working (or not) I will let you know :)

~Johan

Please, have a look at Fortran VSLexamples available in examplesvslfsource folder of MKL installation directory. In particular, vdrnggaussian.f and vsrnggaussian.f files show how to use Gaussian generator.