scale_random.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
#include "scalelib.h"
module scale_random
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
#ifdef _OPENACC
  use openacc
#endif
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: random_setup
  public :: random_uniform
  public :: random_normal
  public :: random_finalize
  public :: random_knuth_shuffle
 
  interface random_uniform
     module procedure random_uniform_1d
     module procedure random_uniform_2d
     module procedure random_uniform_3d
  end interface
 
  interface random_normal
     module procedure random_normal_1d
     module procedure random_normal_2d
     module procedure random_normal_3d
  end interface random_normal
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  private :: random_reset
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  logical, private :: RANDOM_FIX = .false.
 
  integer, private, allocatable :: RANDOM_seedvar(:)
 
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine random_setup
    use scale_prc, only: &
       prc_abort
    implicit none
 
    namelist / param_random / &
       random_fix
 
    integer :: nseeds, ierr
    !---------------------------------------------------------------------------
 
    log_newline
    log_info("RANDOM_setup",*) 'Setup'
 
    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_random,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       log_info("RANDOM_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       log_error("RANDOM_setup",*) 'Not appropriate names in namelist PARAM_RANDOM. Check!'
       call prc_abort
    endif
    log_nml(param_random)
 
    call random_seed
    call random_seed(size=nseeds)
 
    allocate( random_seedvar(nseeds))
 
    log_newline
    log_info("RANDOM_setup",*) 'Array size for random seed:', nseeds
    if ( random_fix ) then
       log_info("RANDOM_setup",*) 'random seed is fixed.'
    endif
 
    call random_reset
 
    return
  end subroutine random_setup
 
  !-----------------------------------------------------------------------------
  subroutine random_reset
    use scale_prc, only: &
       prc_universal_myrank
    implicit none
 
    integer  :: time1(8) ! date and time information
    real(rp) :: time2    ! CPU time
    !---------------------------------------------------------------------------
 
    if ( random_fix ) then
       ! birthday of SCALE
       time1(1) = 2011 ! The year
       time1(2) = 12   ! The month
       time1(3) = 5    ! The day of the month
       time1(4) = 9    ! Time difference with UTC in minutes
       time1(5) = 10   ! The hour of the day
       time1(6) = 20   ! The minutes of the hour
       time1(7) = 41   ! The seconds of the minute
       time1(8) = 0    ! The milliseconds of the second
       time2    = 0.0_rp
    else
       call date_and_time(values=time1)
       call cpu_time(time2)
    endif
 
    random_seedvar(:) = prc_universal_myrank           &
                      + ( time1(1) - 1970 ) * 32140800 &
                      + time1(2) * 2678400             &
                      + time1(3) * 86400               &
                      + time1(5) * 60                  &
                      + time1(6) * 3600                &
                      + time1(7) * 60                  &
                      + time1(8)                       &
                      + int(time2*1.e6_rp)
 
    call random_seed(put=random_seedvar)
 
    return
  end subroutine random_reset
 
  !-----------------------------------------------------------------------------
  subroutine random_uniform_1d( var )
    implicit none
    real(RP), intent(out) :: var(:)
    !---------------------------------------------------------------------------
 
    call random_number(var)
    !$acc update device(var) if(acc_is_present(var))
 
    return
  end subroutine random_uniform_1d
 
  !-----------------------------------------------------------------------------
  subroutine random_uniform_2d( var )
    implicit none
    real(RP), intent(out) :: var(:,:)
    !---------------------------------------------------------------------------
 
    call random_number(var)
    !$acc update device(var) if(acc_is_present(var))
 
    return
  end subroutine random_uniform_2d
 
  !-----------------------------------------------------------------------------
  subroutine random_uniform_3d( var )
    implicit none
    real(RP), intent(out) :: var(:,:,:)
    !---------------------------------------------------------------------------
 
    call random_number(var)
    !$acc update device(var) if(acc_is_present(var))
 
    return
  end subroutine random_uniform_3d
 
  !-----------------------------------------------------------------------------
  subroutine random_normal_1d( var )
    implicit none
    real(RP), intent(out) :: var(:)
    integer :: n
 
    n = size(var)
    call get_normal( n, var(:) )
 
    return
  end subroutine random_normal_1d
 
  subroutine random_normal_2d( var )
    implicit none
    real(RP), intent(out) :: var(:,:)
    integer :: n
 
    n = size(var)
    call get_normal( n, var(:,:) )
 
    return
  end subroutine random_normal_2d
 
  subroutine random_normal_3d( var )
    implicit none
    real(RP), intent(out) :: var(:,:,:)
    integer :: n
 
    n = size(var)
    call get_normal( n, var(:,:,:) )
 
    return
  end subroutine random_normal_3d
 
  subroutine random_finalize
    implicit none
    !---------------------------------------------------------------------------
 
    random_fix = .false.
 
    deallocate( random_seedvar )
 
    return
  end subroutine random_finalize
 
  subroutine random_knuth_shuffle( num, a )
    implicit none
    integer, intent(in) :: num
    integer, intent(out) :: a(num)
    integer :: i, randpos, tmp
    real(dp) :: r
 
    do i = 1, num
      a(i) = i
    end do
 
    do i = num, 2, -1
      call random_number(r)
      randpos = int(r * i) + 1
      tmp = a(randpos)
      a(randpos) = a(i)
      a(i) = tmp
    end do
 
    return
  end subroutine random_knuth_shuffle
 
  ! private
 
  subroutine get_normal( n, var )
    use scale_const, only: &
       pi => const_pi
    implicit none
    integer,  intent(in)  :: n
    real(rp), intent(out) :: var(n)
 
    real(rp) :: rnd(n+1)
    real(rp) :: fact
    real(rp) :: theta
    integer :: i
    !---------------------------------------------------------------------------
 
    call random_number(rnd)
 
    !$omp parallel do &
    !$omp private(fact,theta)
    do i = 1, n/2
       fact = sqrt(-2.0_rp * log( 1.0_rp - rnd(i*2-1) ) ) ! 0 <= rnd < 1
       theta = 2.0_rp * pi * rnd(i*2)
       var(i*2-1) = fact * cos(theta)
       var(i*2  ) = fact * sin(theta)
    end do
    if ( mod(n,2) == 1 ) then
       fact = sqrt(-2.0_rp * log( 1.0_rp - rnd(n) ) )
       theta = 2.0_rp * pi * rnd(n+1)
       var(n) = fact * cos(theta)
    end if
 
    !$acc update device(var) if(acc_is_present(var))
 
    return
  end subroutine get_normal
 
end module scale_random