Copy module RANDOM from GODUNOV and make it working.

This commit is contained in:
Grzegorz Kowal 2012-07-22 19:37:58 -03:00
parent 77d9e511c1
commit 3d7a74f61d
4 changed files with 296 additions and 179 deletions

View File

@ -56,7 +56,7 @@ program amun
#ifdef MPI #ifdef MPI
use parameters , only : redistribute_parameters use parameters , only : redistribute_parameters
#endif /* MPI */ #endif /* MPI */
use random , only : init_generator use random , only : initialize_random, finalize_random
use timers , only : initialize_timers, start_timer, stop_timer & use timers , only : initialize_timers, start_timer, stop_timer &
, set_timer, get_timer, get_timer_total & , set_timer, get_timer, get_timer_total &
, timer_enabled, timer_description, ntimers , timer_enabled, timer_description, ntimers
@ -200,9 +200,9 @@ program amun
em = 59 em = 59
es = 59 es = 59
! initialize random number generator ! initialize the random number generator
! !
call init_generator() call initialize_random(nprocs, nproc)
! initialize block module ! initialize block module
! !
@ -431,6 +431,10 @@ program amun
! !
call clear_coords() call clear_coords()
! finalize the random number generator
!
call finalize_random()
! stop time accounting for the termination ! stop time accounting for the termination
! !
call stop_timer(itm) call stop_timer(itm)

View File

@ -1111,7 +1111,7 @@ module io
! set the seed values ! set the seed values
! !
call set_seeds(seeds(:)) call set_seeds(nseeds, seeds(:))
! deallocate seed array ! deallocate seed array
! !

View File

@ -227,7 +227,7 @@ parameters.o : parameters.F90 mpitools.o
problem.o : problem.F90 blocks.o constants.o coords.o mpitools.o \ problem.o : problem.F90 blocks.o constants.o coords.o mpitools.o \
random.o scheme.o variables.o random.o scheme.o variables.o
scheme.o : scheme.F90 blocks.o config.o interpolation.o variables.o scheme.o : scheme.F90 blocks.o config.o interpolation.o variables.o
random.o : random.F90 mpitools.o random.o : random.F90 mpitools.o parameters.o
timers.o : timers.F90 timers.o : timers.F90
variables.o : variables.F90 variables.o : variables.F90

View File

@ -21,314 +21,427 @@
!! !!
!!****************************************************************************** !!******************************************************************************
!! !!
!! module: RANDOM - handles random number generators by Marsaglia & Tsang !! module: RANDOM
!! !!
!! references: !! This module provides subroutines to random number generators.
!! Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating
!! random variables', J. Statist. Software, v5(8).
!! !!
!! Copyright (C) 2000 George Marsaglia, Wai Wan Tsang !! References:
!! Copyright (C) 2008 John Burkardt !!
!! [1] Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating
!! random variables', J. Statist. Software, v5(8)
!! !!
!!****************************************************************************** !!******************************************************************************
! !
module random module random
#ifdef PROFILE
! include external subroutines
!
use timers, only : set_timer, start_timer, stop_timer
#endif /* PROFILE */
! declare all module variables as implicit
!
implicit none implicit none
#ifdef PROFILE
! timer indices
!
integer , save :: iri, irc
#endif /* PROFILE */
! random generator type
!
character(len = 32), save :: gentype = "same"
! variables to store seeds for random generator
!
integer , save :: kp = 0
integer , save :: nseeds, lseed
integer(kind=4), dimension(:), allocatable, save :: seeds
! by default everything is private
!
private private
integer , save :: nseeds ! declare public subroutines
integer(kind=4), dimension(:), allocatable, save :: seeds !
integer(kind=4), dimension(128) , save :: kn public :: initialize_random, finalize_random
real (kind=4), dimension(128) , save :: fn, wn public :: nseeds, set_seeds, get_seeds, randomu, randomz, randomn
public :: init_generator, nseeds, get_seeds, set_seeds &
, randomu, randomz, randomn
contains contains
! !
!=============================================================================== !===============================================================================
! !
! init_generator: subroutine initializes the random number generator ! subroutine INITIALIZE_RANDOM:
! ----------------------------
!
! subroutine initializes random number generator;
! !
!=============================================================================== !===============================================================================
! !
subroutine init_generator() subroutine initialize_random(nprocs, nproc)
! obtain required variables from other modules
!
use parameters, only : get_parameter_string
! declare all variables as implicit
!
implicit none implicit none
! subroutine arguments
!
integer, intent(in) :: nprocs, nproc
! local variables ! local variables
! !
integer(kind=4) :: i integer :: i
real(kind=8) :: dn, tn, q real :: r
! parameters
!
real(kind=8), parameter :: m1 = 2.14748364800000d+09
real(kind=8), parameter :: vn = 9.91256303526217d-03
! OpenMP functions
!
!$ integer :: omp_get_num_threads, omp_get_thread_num
! !
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
! initialize the seeds required by all subroutines #ifdef PROFILE
! set timer descriptions
! !
nseeds = 1 call set_timer('random generator initialization', iri)
!$omp parallel call set_timer('random number generation' , irc)
!$ nseeds = omp_get_num_threads()
!$omp end parallel
! allocate array for seeds ! start accounting time for the random number generator
! !
allocate(seeds(0:nseeds-1)) call start_timer(iri)
#endif /* PROFILE */
! fill seeds with random numbers ! set the processor number
! !
do i = 0, nseeds - 1 kp = nproc
call random_number(q)
seeds(i) = 123456789 * q ! obtain the generator type
!
call get_parameter_string("gentype", gentype)
! calculate the number of seeds
!
nseeds = 2 * nprocs
lseed = nseeds - 1
! allocate seeds for random number generator
!
allocate(seeds(0:lseed))
! prepare the seeds depending on the type of generator
!
select case(gentype)
case('random')
do i = 0, lseed
call random_number(r)
seeds(i) = 123456789 * r
end do end do
case default
! prepare the arrays used by nonunifor distribution generators call random_number(r)
! do i = 0, nprocs - 1
dn = 3.442619855899d+00 seeds(i) = 123456789 * r
tn = 3.442619855899d+00
q = vn / exp(-0.5d+00 * dn * dn)
kn(1) = int((dn / q) * m1)
kn(2) = 0
wn(1) = real(q / m1, kind=4)
wn(128) = real(dn / m1, kind=4)
fn(1) = 1.0e+00
fn(128) = real(exp(-0.5d+00 * dn * dn), kind=4)
do i = 127, 2, -1
dn = sqrt(-2.0d+00 * log(vn / dn + exp(-0.5d+00 * dn * dn)))
kn(i+1) = int((dn / tn) * m1)
tn = dn
fn(i) = real(exp(-0.5d+00*dn*dn), kind=4)
wn(i) = real(dn / m1, kind=4)
end do end do
call random_number(r)
do i = nprocs, lseed
seeds(i) = 123456789 * r
end do
end select
#ifdef PROFILE
! stop accounting time for the random number generator
!
call stop_timer(iri)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine initialize_random
!
!===============================================================================
!
! subroutine FINALIZE_RANDOM:
! --------------------------
!
! subroutine releases memory allocated by random number generator variables;
!
!===============================================================================
!
subroutine finalize_random()
! declare all variables as implicit
!
implicit none
! !
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
end subroutine init_generator #ifdef PROFILE
! start accounting time for the random number generator
! !
!=============================================================================== call start_timer(iri)
! #endif /* PROFILE */
! get_seeds: subroutine returns the generator seeds
!
!===============================================================================
!
subroutine get_seeds(seed)
! deallocate seeds if they are allocated
!
if (allocated(seeds)) deallocate(seeds)
#ifdef PROFILE
! stop accounting time for the random number generator
!
call stop_timer(iri)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine finalize_random
!
!===============================================================================
!
! subroutine SET_SEEDS:
! --------------------
!
! subroutine sets the seeds from the given array;
!
!===============================================================================
!
subroutine set_seeds(np, seed)
! declare all variables as implicit
!
implicit none implicit none
! subroutine arguments ! input arguments
! !
integer(kind=4), dimension(0:nseeds-1), intent(out) :: seed integer , intent(in) :: np
integer(kind=4), dimension(0:np-1), intent(in) :: seed
! !
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
seed(:) = seeds(:) #ifdef PROFILE
! start accounting time for the random number generator
! !
!------------------------------------------------------------------------------- call start_timer(iri)
! #endif /* PROFILE */
end subroutine get_seeds
!
!===============================================================================
!
! set_seeds: subroutine sets the generator seeds
!
!===============================================================================
!
subroutine set_seeds(seed)
implicit none ! set the seeds only if the input array and seeds have the same sizes
!
if (np .eq. nseeds) then
! subroutine arguments seeds(0:lseed) = seed(0:lseed)
!
integer(kind=4), dimension(0:nseeds-1), intent(in) :: seed end if
!
!------------------------------------------------------------------------------- #ifdef PROFILE
! ! stop accounting time for the random number generator
seeds(:) = seed(:)
! !
call stop_timer(iri)
#endif /* PROFILE */
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
end subroutine set_seeds end subroutine set_seeds
! !
!=============================================================================== !===============================================================================
! !
! shr3: subroutine evaluates the SHR3 generator for integers ! subroutine GET_SEEDS:
! --------------------
! !
! input arguments: ! subroutine returns the seeds through an array;
!
! np - the process number (for parallel operations)
! !
!=============================================================================== !===============================================================================
! !
function shr3(np) result(jr) subroutine get_seeds(seed)
! declare all variables as implicit
!
implicit none implicit none
integer(kind=4) :: np, jp, jt, jr ! output arguments
!
integer(kind=4), dimension(0:lseed), intent(out) :: seed
! !
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
jp = seeds(np) #ifdef PROFILE
jt = jp ! start accounting time for the random number generator
jt = ieor(jt, ishft(jt, 13))
jt = ieor(jt, ishft(jt, -17))
jt = ieor(jt, ishft(jt, 5))
seeds(np) = jt
jr = jp + jt
return
! !
call start_timer(iri)
#endif /* PROFILE */
seed(0:lseed) = seeds(0:lseed)
#ifdef PROFILE
! stop accounting time for the random number generator
!
call stop_timer(iri)
#endif /* PROFILE */
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
end function shr3 end subroutine get_seeds
! !
!=============================================================================== !===============================================================================
! !
! randomu: function generates uniformly distributed random numbers from the ! function RANDOMU:
! range 0.0...1.0 ! ----------------
! !
! input arguments: ! function generates uniformly distributed random numbers in range 0.0..1.0;
!
! np - the process number (for parallel operations)
! !
!=============================================================================== !===============================================================================
! !
function randomu(np) result(val) function randomu() result(val)
integer(kind=4) :: np ! declare all variables as implicit
!
implicit none
! output variables
!
real :: val real :: val
! local variables
!
integer(kind=4) :: jz, jsr
! !
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
val = 0.5 + 0.23283064365e-9 * shr3(np) #ifdef PROFILE
! start accounting time for the random number generation
!
call start_timer(irc)
#endif /* PROFILE */
jsr = seeds(kp)
jz = jsr
jsr = ieor( jsr, ishft( jsr, 13 ) )
jsr = ieor( jsr, ishft( jsr, -17 ) )
jsr = ieor( jsr, ishft( jsr, 5 ) )
seeds(kp) = jsr
val = 0.5 + 0.23283064365e-9 * (jz + jsr)
#ifdef PROFILE
! stop accounting time for the random number generation
!
call stop_timer(irc)
#endif /* PROFILE */
return return
!
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
end function randomu end function randomu
! !
!=============================================================================== !===============================================================================
! !
! randomz: function generates uniformly distributed random numbers from the ! function RANDOMZ:
! range -1.0...1.0 ! ----------------
! !
! input arguments: ! function generates uniformly distributed random numbers in range -0.5..0.5;
!
! np - the process number (for parallel operations)
! !
!=============================================================================== !===============================================================================
! !
function randomz(np) result(val) function randomz() result(val)
integer(kind=4) :: np ! declare all variables as implicit
!
implicit none
! output variables
!
real :: val real :: val
! local variables
!
integer(kind=4) :: jz, jsr
! !
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
val = 0.46566128730e-9 * shr3(np) #ifdef PROFILE
! start accounting time for the random number generation
!
call start_timer(irc)
#endif /* PROFILE */
jsr = seeds(kp)
jz = jsr
jsr = ieor( jsr, ishft( jsr, 13 ) )
jsr = ieor( jsr, ishft( jsr, -17 ) )
jsr = ieor( jsr, ishft( jsr, 5 ) )
seeds(kp) = jsr
val = 0.23283064365e-9 * (jz + jsr)
#ifdef PROFILE
! stop accounting time for the random number generation
!
call stop_timer(irc)
#endif /* PROFILE */
return return
!
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
end function randomz end function randomz
! !
!=============================================================================== !===============================================================================
! !
! randomn: function generates normally distributed random numbers ! function RANDOMN:
! ----------------
! !
! input arguments: ! function generates uniformly distributed random numbers in range -1.0..1.0;
!
! np - the process number (for parallel operations)
! !
!=============================================================================== !===============================================================================
! !
function randomn(np) result(val) function randomn() result(val)
! declare all variables as implicit
!
implicit none implicit none
! input/output arguments ! output variables
! !
integer(kind=4) :: np
real :: val real :: val
! local variables ! local variables
! !
integer(kind=4) :: hz, iz integer(kind=4) :: jz, jsr
real(kind=4) :: x, y, z, t
! parameters
!
real(kind=4), parameter :: r = 3.442620e+00
! !
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
hz = shr3(np) #ifdef PROFILE
iz = iand(hz, 127) ! start accounting time for the random number generation
!
call start_timer(irc)
#endif /* PROFILE */
if (abs(hz) .lt. kn(iz+1)) then jsr = seeds(kp)
val = real(hz, kind=4) * wn(iz+1) jz = jsr
else
do
if (iz .eq. 0) then
do
x = - 0.2904764e+00 * log(randomu(np))
y = - log(randomu(np))
if ((x * x) .le. (y + y)) then
exit
end if
end do
if (hz .le. 0) then jsr = ieor( jsr, ishft( jsr, 13 ) )
val = - r - x jsr = ieor( jsr, ishft( jsr, -17 ) )
else jsr = ieor( jsr, ishft( jsr, 5 ) )
val = + r + x
end if
exit seeds(kp) = jsr
end if
x = real(hz, kind=4) * wn(iz+1) val = 0.46566128730e-9 * (jz + jsr)
z = fn(iz+1) + randomu(np) * (fn(iz) - fn(iz+1))
t = exp(-0.5e+00 * x * x)
if (z .lt. t) then #ifdef PROFILE
val = x ! stop accounting time for the random number generation
exit !
end if call stop_timer(irc)
#endif /* PROFILE */
hz = shr3(np)
iz = iand(hz, 127)
if (abs(hz) .lt. kn(iz+1)) then
val = real(hz, kind=4) * wn(iz+1)
exit
end if
end do
end if
return return
!
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
end function randomn end function randomn
!===============================================================================
! !
end module random end module random