!-----------------------------------------------------------------------
!
!  diehard_data.f90
!
!  Utility to generate data file from random_pl package to be
!  tested by Marsaglia's 'diehard' package; output is to 
!  'diehard.dat'. Note that diehard expects 32-bit numbers, but
!  random_pl only produces 30- or 31-bit numbers; therefore we pick
!  twice as many numbers, and use half the numbers to replace the
!  two ls bits.
!
!-----------------------------------------------------------------------

program diehard_data

  use globals
  use random_pl

  implicit none

  integer(rpk), parameter :: m = 2147483647
  integer(rpk) :: ai, aa
  integer :: iol
  integer(rpk) :: j
  integer(rpk), parameter :: num=100000000
  real(wp) :: a
  integer :: timer_init, ntime, time_unit

  ai = 0
  inquire(iolength=iol) ai
  open(unit=30,file='diehard.dat',action='write',status='replace',&
               form='unformatted',access='direct',recl=iol)

  do j = 1, 3*1000**2
    ai = rand_pl() * m
    ai = ishft(ai, -1)
    ai = ishft(ai, 2)
    aa = rand_pl() * m
    ai = ai + iand(aa, 3_rpk)
    write(30, rec=j) ai
  end do

  close(30)

  call system_clock(count=timer_init, count_rate=time_unit)
  do j = 1, num
    a = rand_pl()
  end do
  call system_clock(count=ntime, count_rate=time_unit)
  write(0,*) 'time to fetch ', num, ' random numbers:'
  write(0,*) real(ntime - timer_init) / real(time_unit), ' seconds'
  write(0,*) 'final deviate: ', a


end program diehard_data