Keeping the same seeds for random numbers with OpenMP loop

Peter_F_5 · ‎02-18-2018

I have an OpenMP paralelized loop. I would like to be able to replicate experiments everytime I run the program. The program calls a lot of random numbers. The issue is given the nature of parallelization I cannot replicate the experiments from one run to the next run. Is there anyway of doing this? I can see two options:

1) Call all the random numbers before executing the paralelized loop (I really would like to avoid this option as potentially the number of random numbers to be generated could be unknown at the time one enters the loop);

2) Have an option to fix the seeds for each thread. This is something that is possible in matlab very easily by setting a seed that depends on the index of the parallelized loop. E.g:

%% Matlab example

parfor k = 1:nsim
   rng(k, 'twister');
end

Below a minimal working example of something similar to what I am running:

%% Fortran minimal working example

program mwe 
use tools
implicit none

integer :: ind1, ind2, ind3
integer :: n1 = 1000, n2 = 750, n3 = 300
real :: mu, mu1=0.01, mu2=0.02, sigma, sigma1=0.20, sigma2=0.40
real :: aux
real, allocatable :: z(:,:,:)
allocate(z(n1, n2, n3))


!$omp parallel do default(private) shared(n1, n2, n3, z, mu1, mu2, sigma1, sigma2)
do ind1 = 1,n1
    
    
    do ind2=1,n2
    do ind3 = 1,n3
    
         call RANDOM_NUMBER(aux)   
         if (aux .le. 0.5) then
                mu    = mu1
                sigma = sigma1
         else
                mu = mu2 
                sigma = sigma2
        end if
    
        z(ind1,ind2,ind3) = normal(mu,sigma)
        
        
    end do
    end do
    
end do
!$omp end parallel do
ind1 = 1

end program



MODULE tools

contains 
    
    
    function ran1()  !returns random number between 0 - 1 
        implicit none 
        real :: ran1,x 
        call random_number(x) ! built in fortran 90 random number function 
        ran1=x 
    end function ran1
    
    function normal(mean,sigma) !returns a normal distribution 

        implicit none 
        real :: normal,tmp 
        real :: mean,sigma 
        integer flag 
        real :: fac,gsave,rsq,r1,r2 
        save flag,gsave 
        data flag /0/ 
        if (flag.eq.0) then 
        rsq=2.0 
            do while(rsq.ge.1.0.or.rsq.eq.0.0) ! new from for do 
                r1=2.0*ran1()-1.0
                r2=2.0*ran1()-1.0 
                rsq=r1*r1+r2*r2 
            enddo 
            fac=sqrt(-2.0*log(rsq)/rsq) 
            gsave=r1*fac 
            tmp=r2*fac 
            flag=1 
        else 
            tmp=gsave 
            flag=0 
        endif 
        normal=tmp*sigma+mean 
        return 
    end function normal
    
    
end module

Steve_Lionel · ‎02-18-2018

There is no support in RANDOM_NUMBER for thread-specific sequences. Fortran 2018 sort-of has this but only at the coarray image level.

Other than your first option, all I can think of is to provide your own RNG that allows you to pass the seed in, and have that seed be a thread-private (or FIRSTPRIVATE) variable.

How much do you care about the "quality" of the RNG? You could use one of the portability library routines that have seed arguments.

Peter_F_5 · ‎02-18-2018

Steve Lionel (Ret.) wrote:

There is no support in RANDOM_NUMBER for thread-specific sequences. Fortran 2018 sort-of has this but only at the coarray image level.

Other than your first option, all I can think of is to provide your own RNG that allows you to pass the seed in, and have that seed be a thread-private (or FIRSTPRIVATE) variable.

How much do you care about the "quality" of the RNG? You could use one of the portability library routines that have seed arguments.

I do not really care about "quality" of the RNG (unless there is something that somehow bias the random numbers and consequently the overall results). Other than that I am happy to use any routine that generates random numbers. If they are thread safe and allow me to replicate results even with OpenMP perfect.

jimdempseyatthecove · ‎02-18-2018

Check the Intel Visual Fortran document via index for RANDOM_NUMBER. You will find an example using RANDOM_SEED and RANDOM_NUMBER. You can modify this such that you can:

a) create a once-only parallel region that determines the number of threads
b) allocate a seed array the size required for the number of threads (zero based)
c) populate the seed array with appropriate seed values, one per thread
d) in you ~~main parallel loop (do ind1=1,n1)~~ inner most loop, use a critical section that PUT's the seed value for the thread, obtains a harvest, then GET's the updated seed value before exiting the critical section.

Note, you may want to have each thread obtain a harvest of n3 number of random numbers. IOW place the d) acquisition between the do ind2 and do ind3 to collect n3 number of random numbers and thus reduce the number of critical sections from n1*n2*n3 to n1*n2.

Jim Dempsey

Peter_F_5 · ‎02-18-2018

Jim, thank you for the answer. I will give that a go and see if I can make it work.

Steve_Lionel · ‎02-18-2018

The problem with Jim's suggestion, though it will work, is that it introduces a critical section that blocks any other thread from getting a random number until the first one exits. This is why I suggested a simpler one that doesn't have global state the way RANDOM_NUMBER does.

gib · ‎02-18-2018

I modified Marsaglia & Tsang's ziggurat RNG to work with OpenMP. One initializes a collection of npar RNGs, and calls the various RV-returning functions with the argument kpar = 0,..,npar-1. I've been using this for several years now, and it seems to work well.

!Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating random variables', J. Statist. Software, v5(8).

I'd be happy to send you the code if you wish.

Gib

jimdempseyatthecove · ‎02-19-2018

>>The problem with Jim's suggestion, though it will work, is that it introduces a critical section that blocks any other thread from getting a random number until the first one exits.

And this is why I suggested you get batches of random numbers as opposed to one at a time.

Gib's suggestion of using a random number generator that is passed in and out the seed value is the better choice. Each thread though must be careful enough to obtain an initial (and different) seed value for use on each call.

Jim Dempsey

Peter_F_5 · ‎02-19-2018

gib wrote:

I modified Marsaglia & Tsang's ziggurat RNG to work with OpenMP. One initializes a collection of npar RNGs, and calls the various RV-returning functions with the argument kpar = 0,..,npar-1. I've been using this for several years now, and it seems to work well.

!Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating random variables', J. Statist. Software, v5(8).

I'd be happy to send you the code if you wish.

Gib

Hi Gib,
It would be great if you could share the code. Thank you very much.
Peter

gib · ‎02-19-2018

This is how the generator is initialised (in a module that USEs Par_Zig_mod):

subroutine RngInitialisation(npar, seed)
integer :: npar, seed
integer, allocatable :: zig_seed(:)
integer :: i
integer :: grainsize = 32

allocate(zig_seed(0:npar-1))
do i = 0,npar-1
    zig_seed(i) = seed*(i+1)
enddo
call par_zigset(npar,zig_seed,grainsize)
par_zig_init = .true.
end subroutine

This is the module:

! Marsaglia & Tsang generator for random normals & random exponentials.
! Translated from C by Alan Miller (amiller@bigpond.net.au)

! Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating
! random variables', J. Statist. Software, v5(8).

! This is an electronic journal which can be downloaded from:
! http://www.jstatsoft.org/v05/i08

! N.B. It is assumed that all integers are 32-bit.
! N.B. The value of M2 has been halved to compensate for the lack of
!      unsigned integers in Fortran.

! Latest version - 1 January 2001

! Parallel version - October 2006

! This version has been customised for parallel processing use,
! specifically with OpenMP.  Each thread uses its own pseudo-random
! sequence. (Gib Bogle)
! Note: thread index kpar is 0-based
!--------------------------------------------------------------------------

MODULE Par_Zig_mod

   IMPLICIT NONE

   PRIVATE

   INTEGER,  PARAMETER  ::  DP=SELECTED_REAL_KIND( 12, 60 )
   REAL(DP), PARAMETER  ::  m1=2147483648.0_DP,   m2=2147483648.0_DP,      &
                            half=0.5_DP
   REAL(DP)             ::  dn0=3.442619855899_DP, tn0=3.442619855899_DP,    &
                            vn=0.00991256303526217_DP,                     &
                            q,                    de0=7.697117470131487_DP, &
                            te0=7.697117470131487_DP,                       &
                            ve=0.003949659822581572_DP
!   INTEGER,  SAVE       ::  iz, jz, jsr=123456789, kn(0:127),              &
!                            ke(0:255), hz
!   REAL(DP), SAVE       ::  wn(0:127), fn(0:127), we(0:255), fe(0:255)
!   LOGICAL,  SAVE       ::  initialized=.FALSE.

	integer, save :: par_n = 0, par_step
	integer, allocatable, save ::  par_jsr(:), par_kn(:,:), par_ke(:,:)
	real(DP), allocatable, save :: par_wn(:,:), par_fn(:,:), par_we(:,:), par_fe(:,:)

   PUBLIC  :: par_zigset, par_zigfree, par_shr3, par_uni, par_rnor, par_rexp, par_test

CONTAINS

!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
SUBROUTINE par_zigset( npar, par_jsrseed, grainsize)

   INTEGER, INTENT(IN)  :: npar, grainsize, par_jsrseed(0:*)

   INTEGER  :: i, kpar
   REAL(DP) :: dn, tn, de, te

   par_n = npar
   par_step = grainsize
! First we need to allocate all the non-volatile arrays with the size npar
   allocate(par_jsr(0:npar*par_step))
   allocate(par_kn(0:127,0:npar-1))
   allocate(par_ke(0:255,0:npar-1))
   allocate(par_wn(0:127,0:npar-1))
   allocate(par_fn(0:127,0:npar-1))
   allocate(par_we(0:255,0:npar-1))
   allocate(par_fe(0:255,0:npar-1))

! Now treat each instance separately
do kpar = 0,npar-1
   !  Set the seed
   par_jsr(kpar*par_step) = par_jsrseed(kpar)

   !  Tables for RNOR
   dn = dn0
   tn = tn0
   q = vn*EXP(half*dn*dn)
   par_kn(0,kpar) = (dn/q)*m1
   par_kn(1,kpar) = 0
   par_wn(0,kpar) = q/m1
   par_wn(127,kpar) = dn/m1
   par_fn(0,kpar) = 1.0_DP
   par_fn(127,kpar) = EXP( -half*dn*dn )
   DO  i = 126, 1, -1
      dn = SQRT( -2.0_DP * LOG( vn/dn + EXP( -half*dn*dn ) ) )		! dn
      par_kn(i+1,kpar) = (dn/tn)*m1
      tn = dn														! tn
      par_fn(i,kpar) = EXP(-half*dn*dn)
      par_wn(i,kpar) = dn/m1
   END DO

   !  Tables for REXP
   de = de0
   te = te0
   q = ve*EXP( de )
   par_ke(0,kpar) = (de/q)*m2
   par_ke(1,kpar) = 0
   par_we(0,kpar) = q/m2
   par_we(255,kpar) = de/m2
   par_fe(0,kpar) = 1.0_DP
   par_fe(255,kpar) = EXP( -de )
   DO  i = 254, 1, -1
      de = -LOG( ve/de + EXP( -de ) )								! de
      par_ke(i+1,kpar) = m2 * (de/te)
      te = de														! te
      par_fe(i,kpar) = EXP( -de )
      par_we(i,kpar) = de/m2
   END DO
enddo
RETURN
END SUBROUTINE par_zigset

!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
subroutine par_zigfree
deallocate(par_jsr)
deallocate(par_kn)
deallocate(par_ke)
deallocate(par_wn)
deallocate(par_fn)
deallocate(par_we)
deallocate(par_fe)
end subroutine

!--------------------------------------------------------------------------
!  Generate random 32-bit integers
!--------------------------------------------------------------------------
FUNCTION par_shr3(kpar) RESULT( ival )
   INTEGER  ::  ival, kpar
   integer :: jz, jsr

   jsr = par_jsr(kpar*par_step)
   jz = jsr
   jsr = IEOR( jsr, ISHFT( jsr,  13 ) )
   jsr = IEOR( jsr, ISHFT( jsr, -17 ) )
   jsr = IEOR( jsr, ISHFT( jsr,   5 ) )
   par_jsr(kpar*par_step) = jsr
   ival = jz + jsr
   RETURN
END FUNCTION par_shr3

!--------------------------------------------------------------------------
!  Generate uniformly distributed random numbers, sequence kpar
!--------------------------------------------------------------------------
FUNCTION par_uni(kpar) RESULT( fn_val )
   integer :: kpar
   REAL(DP)  ::  fn_val

   if (kpar >= par_n) then
       write(*,*) 'thread number: ',kpar,' exceeds initialized max: ',par_n-1
       write(21,*) 'thread number: ',kpar,' exceeds initialized max: ',par_n-1
       stop
   endif
   fn_val = half + 0.2328306e-9_DP * par_shr3(kpar)
   RETURN
END FUNCTION par_uni

FUNCTION par_test(kpar) RESULT( fn_val )
    integer :: kpar
    integer :: fn_val
    if (kpar >= par_n) then
        write(*,*) 'thread number: ',kpar,' exceeds initialized max: ',par_n-1
        stop
    endif
    fn_val = kpar
    return
end function

!--------------------------------------------------------------------------
!  Generate random normals, sequence kpar
!--------------------------------------------------------------------------
FUNCTION par_rnor(kpar) RESULT( fn_val )
   REAL(DP) ::  fn_val
   integer :: kpar

   REAL(DP), PARAMETER  ::  r = 3.442620_DP
   REAL(DP)             ::  x, y
   integer :: iz, hz

!   IF( .NOT. initialized ) CALL zigset( jsr )
   if (kpar >= par_n) then
       write(*,*) 'thread number exceeds initialized max: ',kpar,par_n-1
       stop
   endif
   hz = par_shr3(kpar)
   iz = IAND( hz, 127 )
   IF( ABS( hz ) < par_kn(iz,kpar) ) THEN
      fn_val = hz * par_wn(iz,kpar)
   ELSE
      DO
         IF( iz == 0 ) THEN
            DO
               x = -0.2904764_DP* LOG( par_uni(kpar) )
               y = -LOG( par_uni(kpar) )
               IF( y+y >= x*x ) EXIT
            END DO
            fn_val = r+x
            IF( hz <= 0 ) fn_val = -fn_val
            RETURN
         END IF
         x = hz * par_wn(iz,kpar)
         IF( par_fn(iz,kpar) + par_uni(kpar)*(par_fn(iz-1,kpar)-par_fn(iz,kpar)) < EXP(-half*x*x) ) THEN
            fn_val = x
            RETURN
         END IF
         hz = par_shr3(kpar)
         iz = IAND( hz, 127 )
         IF( ABS( hz ) < par_kn(iz,kpar) ) THEN
            fn_val = hz * par_wn(iz,kpar)
            RETURN
         END IF
      END DO
   END IF
   RETURN
END FUNCTION par_rnor

!--------------------------------------------------------------------------
!  Generate random exponentials, sequence kpar
!--------------------------------------------------------------------------
FUNCTION par_rexp(kpar) RESULT( fn_val )
   REAL(DP)  ::  fn_val
   integer :: kpar

   REAL(DP)  ::  x
   integer :: iz, jz

!   IF( .NOT. initialized ) CALL Zigset( jsr )
   if (kpar >= par_n) then
       write(*,*) 'thread number exceeds initialized max: ',kpar,par_n-1
       stop
   endif
   jz = par_shr3(kpar)
   iz = IAND( jz, 255 )
   IF( ABS( jz ) < par_ke(iz,kpar) ) THEN
      fn_val = ABS(jz) * par_we(iz,kpar)
      RETURN
   END IF
   DO
      IF( iz == 0 ) THEN
         fn_val = 7.69711 - LOG(par_uni(kpar) )
         RETURN
      END IF
      x = ABS( jz ) * par_we(iz,kpar)
      IF( par_fe(iz,kpar) + par_uni(kpar)*(par_fe(iz-1,kpar) - par_fe(iz,kpar)) < EXP( -x ) ) THEN
         fn_val = x
         RETURN
      END IF
      jz = par_shr3(kpar)
      iz = IAND( jz, 255 )
      IF( ABS( jz ) < par_ke(iz,kpar) ) THEN
         fn_val = ABS( jz ) * par_we(iz,kpar)
         RETURN
      END IF
   END DO
   RETURN
END FUNCTION par_rexp

END MODULE par_zig_mod