program perf_alloc

use OMP_LIB
use ISO_C_BINDING
implicit none

include 'mkl_service.fi'

real(8) :: tstart
real(8), dimension(1:2) :: telaps

integer, parameter :: cnt = 100000000
integer :: i

integer, dimension(:), pointer :: a

print *,'num threads = ', omp_get_max_threads()

! non-parallel
tstart = getTime()
do i = 1,cnt
   allocate(a(1:40))
   deallocate(a)
end do
telaps(1) = elapsedTime(tstart)

write(*,*) 'Starting parallel'; flush(6)
! parallel
tstart = getTime()
!$omp parallel do schedule(dynamic, 100) default(shared) private(i, a)
!!$omp parallel do schedule(static) default(shared) private(i) firstprivate(a)
do i = 1,cnt
   call allocate_memory(40,a)
   call deallocate_memory(a)
   !allocate(a(1:40))
   !deallocate(a)
end do
!$omp end parallel do
telaps(2) = elapsedTime(tstart)

print '(a,2f15.5,"s")', 'non/parallel: ', telaps(1:2)
print '(a,f15.5)', 'ratio = ', telaps(1)/telaps(2)


contains

function getTime() result(tstamp)
   real(8) :: tstamp

   integer(8) :: cnt, cntRate
   real(8) :: tdouble

   call system_clock(cnt, cntRate)
   tdouble = real(cnt,8) / real(cntRate,8)
   tstamp = real(tdouble,8)
end function getTime


function elapsedTime(tstart) result(telapsed)
   real(8) :: telapsed
   real(8), intent(in) :: tstart

   integer(8) :: cnt, cntRate
   real(8) :: tdouble

   call system_clock(cnt, cntRate)
   tdouble = real(cnt,8) / real(cntRate,8)
   telapsed = real(tdouble,8) - tstart
end function elapsedTime

!*************************************************************************************************
subroutine allocate_memory(n,ptr)
   implicit none
   integer, intent(in) :: n
   integer, pointer :: ptr(:)


   type (C_PTR) :: cptr

   integer (kind=C_INTPTR_T) :: iptr
   integer(kind=C_SIZE_T) :: size_in_bytes

   integer, parameter :: I8B = selected_int_kind(12)  ! 8-byte integer
   integer, parameter :: I4B = selected_int_kind(9)   ! 4-byte integer

   integer, parameter :: bytes = sizeof(1)
   integer, parameter :: bit_align = storage_size(1,I4B)

   integer :: bsize

   bsize = n*bytes
   size_in_bytes = bsize

   iptr = mkl_malloc(size_in_bytes,bit_align)

   ! \todo - error checking

   cptr = transfer(iptr,cptr)
   call c_f_pointer(cptr,ptr,shape=(/n/))

end subroutine allocate_memory
!*************************************************************************************************
subroutine deallocate_memory(ptr)
   implicit none

   integer, pointer :: ptr(:)

   type (C_PTR) :: cptr

   integer (kind=C_INTPTR_T) :: iptr

   cptr = c_loc(ptr)
   iptr = transfer(cptr,iptr)
   call mkl_free(iptr)

end subroutine deallocate_memory

end program perf_alloc