Is there a gaussian random number generator which can generate n gaussian variables very fast?

1

Deal all,

I am shocked by the length of my thread title :slight_smile:
Inspired from @Beliavsky thread

I just curious, is there some RNG that is particularly good/fast at generating n gaussian random number instead of just generating one by one?

I mean I looked at @Beliavsky his ziggurat rng, https://github.com/Beliavsky/Ziggurat/blob/main/ziggurat.f90
I saw there is a function

function rnor_vec(n) result(ran)
    ! generate n random normal variates
    integer, intent(in) :: n
    real(kind=dp)       :: ran(n)
    integer             :: i
    do i=1,n
       ran(i) = rnor()
    end do
    end function rnor_vec

which can generate a n-element vector gaussian random number, which is exactly what I want. The function rnor_vec(n) is great, but it still generator random number one by one.
I wonder is there is a more `vectorized’ way or very fast way to generate n such random number?

Thank you very much in advance!

1 Like
2

Well, you could merge the rnor_vec and rnor functions, so that you get a do-loop over n and a do-loop inside that to determine the random number A smart compiler might be able to vectorise such a construction, but my guess is that the code involves too many variables for that to be feasible.

However, can you indicate why the current implementation is too slow?

1 Like
3

Thank you very much! @Arjen
It is loop over rnor() is not bad. Actually the rnor() seems faster than many other RNGs. I am just curious if there are faster ways to directly generator n gaussian random numbers instead of looping over rnor() for n times.

I have a subroutine

	subroutine rk4_ti_fullvec_test ( x0, np, nstep, q, h, fi_gi_in, nd, x )
    ! https://stackoverflow.com/questions/69147944/is-there-room-to-further-optimize-the-stochastic-rk-fortran-90-code 
	! https://stackoverflow.com/questions/32809769/how-to-pass-subroutine-names-as-arguments-in-fortran
	use random
	use fg
	implicit none 
    integer(kind = i8), intent(in) :: np
	integer(kind = i8), intent(in) :: nstep
	integer(kind = i8), intent(in) :: nd
	procedure(fi_gi_fullvec_03) :: fi_gi_in
	real(kind = r8), intent(in) :: q,h
	real(kind = r8), intent(in) :: x0(np,nd)
	real ( kind = r8 ), parameter :: alphas(4) = [ 0.47012396888046_r8 &
    ,0.36597075368373_r8 &
	,0.08906615686702_r8 &
	,0.07483912056879_r8 ]
	real(kind = r8), parameter :: as(4,4) = reshape([ &
		&  0.0_r8,              0.0_r8,              0.0_r8,              0.0_r8, &
		&  2.71644396264860_r8, 0.0_r8,              0.0_r8,              0.0_r8, &
		& -6.95653259006152_r8, 0.78313689457981_r8, 0.0_r8,              0.0_r8, &
		&  0.0_r8,              0.48257353309214_r8, 0.26171080165848_r8, 0.0_r8 ], [4,4])
	real(kind = r8), parameter :: qs(4) = [  2.12709852335625_r8, &
		&  2.73245878238737_r8, &
		& 11.22760917474960_r8, &
		& 13.36199560336697_r8 ]
	real(kind = r8) :: ks(np,nd,4)
	real(kind = r8) :: xs(np,nd,4)
	integer(kind = i8) :: i,j,k,l,m,n   
	real(kind = r8), intent(out) :: x(np,nd,0:nstep)
    real(kind = r8) :: xstar(np,nd)
    real ( kind = r8 ) :: f(np,nd), g(np,nd)
    real(kind = r8) :: warray(np,nd,4)
    real(kind = r8) :: sigma(4)  
    !real (kind = r8), allocatable :: normal(:,:,:,:)  
    !allocate(normal(np,nd,4,nstep))
    !normal = reshape(rnor_vec(int(np*nd*4*nstep,kind=i4)),shape(normal))    
    sigma=sqrt(qs*q/h)  
	x(:,:,0) = x0	
	do k = 1, nstep
        xstar = x(:,:,k-1)
		do j = 1,4
            do l = 1,nd 
                 xs(:,l,j) = x(:,l,k-1) + matmul(ks(:,l,:j-1),as(:j-1,j))
                 do n = 1,np 
                      warray(n,l,j) = rnor()*sigma(j)   !gaussian(2)*sqrt(qs(j)*qoh)                      
                 enddo
            enddo  
            call fi_gi_in( np, nd, xs(:,:,j), x0(1,:), f, g ) 
			ks(:,:,j) = h * ( f + g*warray(:,:,j) )
            !ks(:,:,j) = h * ( f + g*normal(:,:,j,k)*sigma(j) )
			xstar = xstar + alphas(j)*ks(:,:,j)     
        enddo    
        x(:,:,k) = xstar   
	enddo
	return
    end subroutine rk4_ti_fullvec_test

	pure subroutine fi_gi_fullvec_03 (np,nd, x, x0, f, g ) 
	use constants
	implicit none
    integer (kind = i8), intent(in) :: np,nd
    real ( kind = r8 ), intent(in) :: x(np,nd)
    real ( kind = r8 ), intent(in) :: x0(nd)
	real ( kind = r8 ), intent(out) :: f(np,nd)
    real ( kind = r8 ), intent(out) :: g(np,nd) 
    integer (kind = i8) :: i
    real ( kind = r8 ), parameter :: gp(2) = [one,zero] 
	f(:,1) = -x(:,2)*x(:,1) 
	f(:,2) = -x(:,2) + x0(2)  
    g(:,1) = gp(1)
    g(:,2) = gp(2)
	return
	end subroutine fi_gi_fullvec_03

Inside the loop of subroutine rk4_ti_fullvec_test, I need to generator many gaussian random number such as

                 do n = 1,np 
                      warray(n,l,j) = rnor()*sigma(j)   !gaussian(2)*sqrt(qs(j)*qoh)                      
                 enddo

where np can be 10^6, nstep can be 1000, nd=2, j=1:4. So the speed of rnor() becomes a little bit crucial since need to generate like 10^610004*2 random gaussian numbers.
If there are ways to directly generator many gaussian random number and store them in an array outside the loops, this may be efficient.

4

Intel oneAPI MKL has routines for random number generation including discrete and continuous distributions.

To explore what it offers, start from here:

If you already downloaded the MKL libraries, the Intel Fortran compiler has a flag /Qmkl:[lib] (Windows) that will add the necessary linkage flags.

2 Likes