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

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, Ziggurat/ziggurat.f90 at main · Beliavsky/Ziggurat · GitHub
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

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

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.

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