!  simple program to compare System_Clock to RDTSC for reliability
!
      integer*4, parameter :: million    = 1000000
      integer*4, parameter :: num_cycles = 10000
      real*8,    parameter :: one        = 1
!
      integer*8 calls(0:num_cycles)
      integer*8 system_tick, system_clock_ticks(0:num_cycles), system_rate
      integer*8 rdtsc_tick, rdtsc_clock_ticks(0:num_cycles)
      external  clock_tick, rdtsc_tick, system_rate
      integer*4 i, last_tick, next_tick
      real*8    seconds, dt, n_dt, s_dt, ss_dt, dx
      real*8    System_Clock_rate 
!
      write (*,*) 'Program to test variability of System_Clock'
      system_clock_ticks = 0
      rdtsc_clock_ticks  = 0
      calls              = 0
      last_tick          = 0
      next_tick          = last_tick+1
!
!  Timing loop
      do i = 1, 20*million
        system_clock_ticks(next_tick) = system_tick ()
        rdtsc_clock_ticks(next_tick)  = rdtsc_tick ()
        if ( system_clock_ticks(last_tick) == system_clock_ticks(next_tick) ) cycle
!
        calls(next_tick) = i
        last_tick        = next_tick
        next_tick        = last_tick+1
        if (last_tick == num_cycles) exit
      end do
!
!  Reporting
      n_dt  = 0
      s_dt  = 0
      ss_dt = 0
      do i = 3,last_tick
        dt = (rdtsc_clock_ticks(i)  - rdtsc_clock_ticks(i-1)) /   & 
             (system_clock_ticks(i) - system_clock_ticks(i-1))
!
        write (*,*) i,  calls(i)-calls(i-1),                          &
                    system_clock_ticks(i) - system_clock_ticks(i-1),  & 
                    rdtsc_clock_ticks(i)  - rdtsc_clock_ticks(i-1),   & 
                    nint (dt)
!
        n_dt  = n_dt  + one
        s_dt  = s_dt  + dt
        ss_dt = ss_dt + dt**2
      end do
!
      System_Clock_rate     = system_rate ()
      calls(0)              = calls(last_tick)              - calls(2)
      system_clock_ticks(0) = system_clock_ticks(last_tick) - system_clock_ticks(2)
      rdtsc_clock_ticks(0)  = rdtsc_clock_ticks(last_tick)  - rdtsc_clock_ticks(2)
      seconds               = system_clock_ticks(0) / System_Clock_rate
!
      write (*,*) 'Elapsed seconds          =', seconds
      write (*,*) 'Total calls              =', calls(0)
      write (*,*) 'Total ticks              =', last_tick-2
      write (*,*) 'Total System_Clock ticks =', system_clock_ticks(0)
      write (*,*) 'Total RDTSC ticks        =', rdtsc_clock_ticks(0)
      write (*,*) 'Calls per second         =', calls(0) / seconds
      write (*,*) 'Cycles per second        =', (last_tick-2) / seconds
      write (*,*) 'System_Clock accuracy    =', seconds / (last_tick-2),' seconds'
!
      s_dt  = s_dt/n_dt
      ss_dt = sqrt ( ( ss_dt - n_dt*s_dt**2 ) / max(n_dt-one,one) )
      write (*,*) ' '
      write (*,*) 'average RDTSC ticks per cycle', s_dt
      write (*,*) 'standard deviation           ', ss_dt
      write (*,*) 'variability                  ', ss_dt / s_dt
      write (*,*) ' '
      write (*,*) 'Extreme values'
      n_dt = 0
      do i = 3,last_tick
        dt = (rdtsc_clock_ticks(i)  - rdtsc_clock_ticks(i-1)) /   & 
             (system_clock_ticks(i) - system_clock_ticks(i-1))
        dx = (dt-s_dt)/ss_dt
        if (abs(dx) < 4.) cycle
        n_dt = n_dt+one
        write (*,2001) i,  calls(i)-calls(i-1),                       &
                    system_clock_ticks(i) - system_clock_ticks(i-1),  & 
                    rdtsc_clock_ticks(i)  - rdtsc_clock_ticks(i-1),   & 
                    nint (dt), dx
      end do
      write (*,*) nint(n_dt),' extreme values'
2001  format (i6,4i10,f8.2)      
      end
!
      integer*8 function system_tick ()
      integer*8 tick, clock_rate, clock_max
      logical :: first_call = .true.
!
      if (first_call) then
        call system_clock ( tick, clock_rate, clock_max )
        write (*,*) 'clock_tick =',tick
        write (*,*) 'clock_rate =',clock_rate
        write (*,*) 'clock_max  =',clock_max
        first_call = .false.
      end if
!
      call system_clock ( tick )
      system_tick = tick
      end

      integer*8 function system_rate ()
      integer*8 tick, clock_rate, clock_max
      logical :: first_call = .true.
!
        call system_clock ( tick, clock_rate, clock_max )
!
      system_rate = clock_rate
      end function system_rate

      integer*8 function clock_tick ()
      integer*8 tick
      real*8    time
      logical :: first_call = .true.
!
      if (first_call) then
        call CPU_Time ( time )
        tick = (time*1.d6)
        write (*,*) 'clock_tick =',tick
        first_call = .false.
      end if
!
      call CPU_Time ( time )
      tick = (time*1.d6)
      clock_tick = tick
      end


module setup_rdtsc

   use ifwin

   implicit none
   private
   public initialize_rdtsc, rdtsc, cp_rdtsc, rdtsc_rate

   integer, parameter :: code32(1)  = [-1866256113]
   integer, parameter :: code64(3)  = [-1052233457,155721954,-1869560880]
!
   integer*8          :: rdtsc_rate = -1    ! ticks per second ~ processor clock rate
!
   interface

      function rdtsc()
         integer(8) rdtsc
      end function rdtsc

   end interface

   pointer (cp_rdtsc,rdtsc)

   contains

      subroutine initialize_rdtsc

         integer code(*)
         pointer (ap,code)

         integer, parameter :: nbits = bit_size(ap)
         ap = VirtualAlloc(NULL,12_HANDLE,MEM_COMMIT,PAGE_EXECUTE_READWRITE)

         if (nbits == 32) then
            code(1:size(code32)) = code32
         else
            code(1:size(code64)) = code64
         end if
         cp_rdtsc = ap

      end subroutine initialize_rdtsc

end module setup_rdtsc

      integer*8 function rdtsc_tick ()
use setup_rdtsc
!
      integer*8 :: tick_start  = -1
      integer*8 :: tick
!
      if (tick_start < 0) then
         call initialize_rdtsc
         tick_start = rdtsc()     ! API
      end if
!
      tick = rdtsc() - tick_start
!
      rdtsc_tick = tick
      end


      subroutine use_RDTSC (time, num)   ! API
use setup_rdtsc
!
      real*8,    intent (out) :: time
      integer*8, intent (out) :: num
!
      integer*8 :: num_start  = -1
!
      if (num_start < 0) then
         call get_rdtsc_rate
         num_start = rdtsc()
      end if
!
      num  = rdtsc() - num_start
      time = dble(num) / dble(rdtsc_rate)
!
      end subroutine use_RDTSC      

   subroutine get_rdtsc_rate
use setup_rdtsc
!
! initialises rdtsc pointer and estimates rdtsc tick rate
!
!    integer*8 tick_rate   ! ticks per second ~ processor clock rate
!
    integer*8  rd, rd_0, num
    real*8     elapse, elapse_0
    integer*4  i
    integer*4, parameter :: million   = 1000000
    integer*4, parameter :: num_cycle = 10 * million
!
    integer*8 :: known_rate = -1   !   2666701126    !  -1 or 2666700003
!
    write (*,11) 'Rdtsc Initialise       ', known_rate
    call initialize_rdtsc
    call use_QueryPerform (elapse, num)
!
    if ( known_rate <= 0) then
      do i = 0, num_cycle
         rd = rdtsc ()
         call use_QueryPerform (elapse, num)
!
         if (i==0) then
            rd_0     = rd
            elapse_0 = elapse
         end if
      end do
!
      known_rate = dble (rd-rd_0) / (elapse-elapse_0)
      i          = dble (i)       / (elapse-elapse_0)
!
      write (*,12) ' initialise duration  =', elapse-elapse_0, ' seconds'
      write (*,11) ' rdtsc duration       =', rd - rd_0,       ' ticks'
      write (*,11) ' rdtsc cycles         =', num_cycle,       ' calls'
      write (*,11) ' rdtsc rate           =', i,               ' calls per second'
    end if
!
    rdtsc_rate = known_rate
      write (*,11) ' rdtsc rate           =', rdtsc_rate, ' ticks per second'
      write (*,10) ' rdtsc initialised'
      write (*,10) ' '
10 format (a)
11 format (a,i12,a)
12 format (a,f12.6,a)
!
end subroutine get_rdtsc_rate

      SUBROUTINE use_QueryPerform (time, num)     !     Win API Routine
!
use ifwin, only: T_LARGE_INTEGER,QueryPerformanceCounter, QueryPerformanceFrequency
!
      real*8,    intent (out) :: time
      integer*8, intent (out) :: num
!
      type(T_LARGE_INTEGER) :: arg
      real*8    :: freq = -1.
      logical*4 :: ll
!
      if (freq < 0) then
         ll   = QueryPerformanceFrequency (arg)     ! cycle rate
         num  = transfer (arg,num)
         freq = 1.0d0 / dble (num)
         write (*,*) 'QueryPerformance frequency =',num,' cycles per second'
      end if
!
      ll   = QueryPerformanceCounter (arg)          ! cycle count; assume ll is ok
      num  = transfer (arg,num)
!
      time = dble (num) * freq
!
      END SUBROUTINE use_QueryPerform
              