gslib_acorni (DEPRECATED) - opengeostat/gslib_new_code GitHub Wiki
Acorni is part of the core GSLIB library. It contains the function acorni that generates and generates random numbers.
This is considered not parallel safe and it is considered deprecated. Use instead RANDOM_NUMBER, which is a built-in Fortran 90 subroutine and is documented at https://gcc.gnu.org/onlinedocs/gcc-10.3.0/gfortran/RANDOM_005fNUMBER.html#RANDOM_005fNUMBER. To use with seed number you need to call RANDOM_SEED, documented at https://gcc.gnu.org/onlinedocs/gcc-10.3.0/gfortran/RANDOM_005fSEED.html#RANDOM_005fSEED.
To use in OpenMP you need to keep in mind that each tread has its own random number state.
see also RANDOM_INIT at https://gcc.gnu.org/onlinedocs/gcc-10.3.0/gfortran/RANDOM_005fINIT.html#RANDOM_005fINIT. Only works for gfortran 9 or above (2018 std) and it is intended for coarrays.
This example always produces the same sequence of random numbers in each run, for each tread. The trick is using the same seed.
program test_random_number
use OMP_LIB
real :: r(5)
integer, allocatable :: seed(:)
integer :: n
call random_seed(size = n)
allocate(seed(n))
seed = 878768768
!$OMP PARALLEL
call random_seed(put=seed)
call random_number(r)
print *, 'Hello from process: ', OMP_GET_THREAD_NUM(), r
!$OMP END PARALLEL
end program
The output is
Hello from process: 5 0.334234238 0.642744005 0.283221900 0.998922229 0.924292505
Hello from process: 1 0.334234238 0.642744005 0.283221900 0.998922229 0.924292505
Hello from process: 6 0.334234238 0.642744005 0.283221900 0.998922229 0.924292505
Hello from process: 7 0.334234238 0.642744005 0.283221900 0.998922229 0.924292505
Hello from process: 4 0.334234238 0.642744005 0.283221900 0.998922229 0.924292505
Hello from process: 0 0.334234238 0.642744005 0.283221900 0.998922229 0.924292505
Hello from process: 2 0.334234238 0.642744005 0.283221900 0.998922229 0.924292505
Hello from process: 3 0.334234238 0.642744005 0.283221900 0.998922229 0.924292505
This hack seems to work, but I am not sure if this is the correct approach:
program test_random_number
use OMP_LIB
real :: r(10000), rr
integer, allocatable :: seed(:)
integer :: n, i
call random_seed(size = n)
allocate(seed(n*3))
! call random_seed(get = seed)
seed = 878768768
call random_seed(put=seed)
!$OMP PARALLEL private(seed, r, i, rr)
do i = 1,10000
call random_number(rr)
r(i) = rr
end do
print *, 'Hello from process: ', OMP_GET_THREAD_NUM(), sum(r)/10000
!$OMP END PARALLEL
end program
The output:
Hello from process: 0 0.334234238 0.642744005 0.283221900
Hello from process: 4 0.162345409 0.840058029 9.94555354E-02
Hello from process: 3 7.85911679E-02 0.241505861 0.400358200
Hello from process: 6 1.73936486E-02 0.396930397 0.280517757
Hello from process: 1 0.186948657 0.335481524 0.356720865
Hello from process: 5 0.844142258 0.983546615 0.513327420
Hello from process: 7 0.766133487 0.591100812 0.382381201
Hello from process: 2 0.252353966 0.576505065 2.22165585E-02
To test that the mean is correct (for uniform in ]0,1[ interval)
program test_random_number
use OMP_LIB
integer, parameter :: l = 10000000
real :: r, rr
integer, allocatable :: seed(:)
integer :: n, i
call random_seed(size = n)
allocate(seed(n*3))
! call random_seed(get = seed)
seed = 878768768
call random_seed(put=seed)
!$OMP PARALLEL private(seed, r, i, rr)
r = 0
do i = 1,l
call random_number(rr)
r = r+ rr
end do
print *, 'Hello from process: ', OMP_GET_THREAD_NUM(), r/real(l)
!$OMP END PARALLEL
end program
Hello from process: 0 0.499775350
Hello from process: 6 0.499944240
Hello from process: 3 0.500050247
Hello from process: 2 0.499839336
Hello from process: 4 0.500093520
Hello from process: 7 0.500090837
Hello from process: 5 0.499973357
Hello from process: 1 0.499919236
with l = 100000000 the program fails (the mean is 0.167772159), due to precision in summation. Using real*8 fix the problem.