randn() in Fortran

[yingwu]

Hi, I am back again. Is there any function to implement the randn() in Fortran? Or I need to use MKL for this. Thanks very much.

Ying

[Jugoslav_Dujic · Accepted solution]

Quoting - yingwu

Hi, I am back again. Is there any function to implement the randn() in Fortran? Or I need to use MKL for this. Thanks very much.

Ying

There are several routines at:

http://users.bigpond.net.au/amiller/random.html

I wasn't aware that MKL contains an appropriate routine -- does it? If it does, what's wrong with using it?

[yingwu]

vsrnggaussian

Quoting - Jugoslav Dujic

There are several routines at:

http://users.bigpond.net.au/amiller/random.html

I wasn't aware that MKL contains an appropriate routine -- does it? If it does, what's wrong with using it?

Hi, yes, The MKL supports this. In matlab, randn() generates Normally distributed pseudorandom numbers, which is standard normal distribution. In fact, it is a special Gaussian distribution. So using vsrnggaussian() from MKL is OK, but we need to specify the mean is 0 and deviation is 1.

I know why there is no randn() in IVF. It is reasonable.....

[Steven_L_Intel1]

IVF itself provides uniform random number generators only. The one in the standard intrinsic RANDOM_NUMBER is the one I would recommend. If you have other distribution requirements, you'll need to find another library, such as MKL, that has what you want.

[yingwu]

Yes, Steve, I agree with this. The MKL is good enough to do this.

I just switch from matlab to fortran, so I need some time and help to be used to Fortran. But anyway, it is not difficult. The matlab provides a rich set of functions to use, which is the same as IVF with MKL or IMSL.

[anthonyrichards]

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 gaussianrand
implicit 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
! from the same distribution for each call.
!
DO I=1,NRAND
CALL RGAUSS(SIGMA,Y1,Y2)
UNIF1(I)=Y1
UNIF2(I)=Y2
END DO
!
! After this loop, there are NRAND values in each array UNIF1 and UNIF2, each with
! a Gaussian distribution with standard deviation SIGMA
!
! To display the resultsfor UNIF1, bin the results usingbin size=0.05*sigma
nbins=201
binsize=0.05*sigma
binhalf=0.5*binsize

do i=1,nrand
! compute which bineach 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