# Is there a gaussian random number generator which can generate n gaussian variables very fast? **URL:** \ **Category:** Uncategorized\ **Created:** [October 29, 2021, 7:46am UTC](https://fortran-lang.discourse.group/t/is-there-a-gaussian-random-number-generator-which-can-generate-n-gaussian-variables-very-fast/2190 "2021-10-29T07:46:05Z")\ **Posts on this page:** 4\ **Page:** 1
**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\ **Post date:** [October 29, 2021, 7:46am UTC](https://fortran-lang.discourse.group/t/is-there-a-gaussian-random-number-generator-which-can-generate-n-gaussian-variables-very-fast/2190/1 "2021-10-29T07:46:05Z")
Deal all, I am shocked by the length of my thread title 🙂 Inspired from @Beliavsky thread > [@Pure random number generators](https://fortran-lang.discourse.group/t/pure-random-number-generators/2170): > > Are there pure random number generators in Fortran? They could be used within pure procedures. You would pass the function a seed and it would return a vector of deviates that depends on the seed. The main program program xzig ! 10/26/2021 01:38 PM demonstrate random number generator that is effectively pure use ziggurat\_mod, only: uni implicit none integer, parameter :: n = 5, nseeds = 2, seeds(nseeds) = [654321,123456] integer :: iter, iseed do iter=1,3 do iseed=1,nseeds write (\*,"(i… 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](https://github.com/Beliavsky/Ziggurat/blob/main/ziggurat.f90) I saw there is a function ```auto 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! ---
**Author:** ![Arjen](https://avatars.discourse-cdn.com/v4/letter/a/b9bd4f/32.png) [@Arjen](https://fortran-lang.discourse.group/u/Arjen)\ **Post date:** [October 29, 2021, 8:29am UTC](https://fortran-lang.discourse.group/t/is-there-a-gaussian-random-number-generator-which-can-generate-n-gaussian-variables-very-fast/2190/2 "2021-10-29T08:29:10Z")
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? ---
**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\ **Post date:** [October 29, 2021, 8:42am UTC](https://fortran-lang.discourse.group/t/is-there-a-gaussian-random-number-generator-which-can-generate-n-gaussian-variables-very-fast/2190/3 "2021-10-29T08:42:41Z")
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 ```auto 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 ```auto 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^6_1000_4\*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. ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [October 29, 2021, 9:34am UTC](https://fortran-lang.discourse.group/t/is-there-a-gaussian-random-number-generator-which-can-generate-n-gaussian-variables-very-fast/2190/4 "2021-10-29T09:34:34Z")
Intel oneAPI MKL has routines for random number generation including discrete and continuous distributions. To explore what it offers, start from here: > **[Random Number Generators](https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-fortran/top/statistical-functions/random-number-generators.html)** > > Intel® oneAPI Math Kernel Library Developer Reference for Fortran If you already downloaded the MKL libraries, the Intel Fortran compiler has a flag [`/Qmkl:[lib]`](https://www.intel.com/content/www/us/en/develop/documentation/oneapi-dpcpp-cpp-compiler-dev-guide-and-reference/top/compiler-reference/compiler-options/compiler-option-details/advanced-optimization-options/qmkl-qmkl.html) (Windows) that will add the necessary linkage flags.