! Benchmark generating random numbers
! Janne Blomqvist 2015
program randbench
#ifdef _OPENMP
  use omp_lib
#endif
  implicit none
  integer, parameter :: dp=kind(0.d0)          ! double precision
  integer, parameter :: i64 = selected_int_kind(18) ! At least 64-bit integer

#ifdef _OPENMP
  print *, "Using up to ", omp_get_max_threads(), " threads."
#endif
  call genr4
  call genr8
contains
  subroutine genr4
    integer, parameter :: n = int(1e7)
    real, save :: r(n)
    integer :: i
    integer(i64) :: t1, t2, td
#ifdef _OPENMP
    integer :: blocks, blocksize, l, h
#endif
    Print *, "Generate default real random variables"
    call system_clock (t1)
    !$omp parallel do private(i)
    do i = 1, n
       call random_number(r(i))
    end do
    !$omp end parallel do
    call system_clock (t2)
    td = t2 - t1
    print *, "Generating ", n, " default reals individually took ", td, " ticks."
    call system_clock (t1)
#ifdef _OPENMP
    blocks = omp_get_max_threads()
    blocksize = n / blocks
    !$omp parallel do private(l,h,i)
    do i = 0, blocks - 1
       l = i * blocksize + 1
       h = l + blocksize - 1
       !print *, "Low: ", l, " High: ", h
       call random_number(r(l:h))
    end do
#else
    call random_number(r)
#endif
    Call system_clock (t2)
    print *, "Generating ", n, " default reals as an array took  ", t2-t1, &
         " ticks. => ind/arr = ", real(td, dp) / (t2-t1)
  end subroutine genr4

  subroutine genr8
    integer, parameter :: n = int(1e7)
    real(dp), save :: r(n)
    integer :: i
    integer(i64) :: t1, t2, td
#ifdef _OPENMP
    integer :: blocks, blocksize, l, h
#endif
    print *, "Generate double real random variables"
    call system_clock (t1)
    !$omp parallel do
    do i = 1, n
       call random_number(r(i))
    end do
    call system_clock (t2)
    td = t2 - t1
    print *, "Generating ", n, " double reals individually took  ", td, " ticks."
    call system_clock (t1)
#ifdef _OPENMP
    blocks = omp_get_max_threads()
    blocksize = n / blocks
    !$omp parallel do private(l,h,i)
    do i = 0, blocks - 1
       l = i * blocksize + 1
       h = l + blocksize - 1
       !print *, "Low: ", l, " High: ", h
       call random_number(r(l:h))
    end do
#else
    call random_number(r)
#endif
    call system_clock (t2)
    print *, "Generating ", n, " double reals as an array took   ", t2-t1, &
         " ticks. => ind/arr = ", real(td, dp) / (t2 -t1)
  end subroutine genr8

end program
