! timer.f90
!
! (c)2007 Jan H. Meinke
!
#ifdef BGL
module rts
! interface definition of rts_get_timebase, Fortran2003, available in XLF Version 9
  interface
    function rts_get_timebase() bind(c)
      use, intrinsic       :: iso_c_binding
      integer(c_long_long) :: rts_get_timebase
    end function rts_get_timebase
  end interface
end module rts
#endif

!! This module contains a minimal set of routines to time sections of code.
!  It can keep track of several timers at once. Each timer is associated with
!  an integer ID. Individual timing events are not stored.
!  
!  Currently only timerIDs between 1 and 10 are valid.
module timer
#ifdef BGL
    use rts
    use, intrinsic :: iso_c_binding
#endif
    implicit none
    
    !! Contains the minimum time ever measured, the maximum time ever measured,
    !  the average time, the standard deviation, and the number of measurements.
    type TimingStructure
        integer :: counter
#ifdef BGL
        integer(c_long_long) :: t0, t1
#else
        real(kind=8) :: t0, t1
#endif
        real(kind=8) :: minTime, maxTime
        real(kind=8) :: sumOfTimes, sumOfSquares
        real(kind=8) :: average, sigma
    end type TimingStructure
            
    !! Contains the timing data for all timers.
    type(TimingStructure), allocatable :: timingData(:)
            
    integer :: initialTimerCount = 10
#ifdef BGL
    real(kind=8) :: clockspeed = 700.0d6
#endif
            
    contains    
        !! Initializes timer
        !  This routine must be called before starting any timer.
        !  
        !  @author Jan H. Meinke
        subroutine init_timer()
            integer :: i
            
            ! Initialize timers
            allocate (timingData(initialTimerCount))
!            forall(i = 1:initialTimerCount)
            do i = 1, initialTimerCount
                timingData(i) = TimingStructure(0, -1.0, -1.0, -1.0, -1.0, 0.0, &
                        0.0, -1.0, -1.0)
!            end forall
            end do
        end subroutine init_timer
    
        !! Start a new measurement. 
        !
        !  @param id ID of the timer
        !
        !  @author Jan H. Meinke
        subroutine start_timer(id)
            integer, intent(in) :: id
            timingData(id)%t0 = getCurrentTime()
        end subroutine start_timer
                
        !! Take a measurement and stop the timer with ID id
        !
        !  @param id ID of the timer
        !
        !  @author Jan H. Meinke
        subroutine stop_timer(id)
            integer, intent(in) :: id
            real(kind=8) :: dt
                    
            timingData(id)%t1 = getCurrentTime()
#ifdef BGL
            dt = (timingData(id)%t1 - timingData(id)%t0) / clockspeed
#else 
            dt = (timingData(id)%t1 - timingData(id)%t0)
#endif
            timingData(id)%counter = timingData(id)%counter + 1
            timingData(id)%maxTime = max(dt, timingData(id)%maxTime)
            if (timingData(id)%minTime > 0) then
                timingData(id)%minTime = min(dt, timingData(id)%minTime)
            else 
                timingData(id)%minTime = dt
            end if
            timingData(id)%sumOfTimes = timingData(id)%sumOfTimes + dt
            timingData(id)%sumOfSquares = timingData(id)%sumOfSquares + dt ** 2
                    
        end subroutine stop_timer
              
        subroutine evaluate(id)
            integer, intent(in) :: id
            if (timingData(id)%counter > 0) then
                timingData(id)%average = timingData(id)%sumOfTimes / timingData(id)%counter
            end if
            if (timingData(id)%counter > 1) then
!                write (*,*) timingData(id)%average, timingData(id)%sumOfSquares,&
!                        timingData(id)%counter
                timingData(id)%sigma = sqrt((abs(timingData(id)%average ** 2 &
                        - timingData(id)%sumOfSquares / timingData(id)%counter)) &
                        / (timingData(id)%counter - 1))
            end if
        end subroutine evaluate
        
        !! Wrapper around the timing function used.
        !  If we are using MPI, we can use MPI_WTime, otherwise we use cpu_time(time)
#ifdef BGL
        integer(c_long_long) function getCurrentTime()
            integer(c_long_long) :: time
            time = rts_get_timebase()
            getCurrentTime = time
        end function getCurrentTime
#elif MPI
        real(kind=8) function getCurrentTime()
            real(kind=8) :: time
            include 'mpif.h'
            time = MPI_Wtime()
            getCurrentTime = time
        end function getCurrentTime
#else
        real(kind=8) function getCurrentTime()
            real(kind=8) :: time
            call cpu_time(time)
            getCurrentTime = time
        end function getCurrentTime
#endif 
        
end module timer