module device_code

contains

! CUDA Kernel !

attributes(global) subroutine VectorSum(a,b,c,size)
  implicit none
  real, dimension(0:size-1) :: a,b,c
  integer, value :: size
  integer :: i
  do i=(threadIdx%x-1)+(blockIdx%x-1)*blockDim%x,size-1,blockDim%x*gridDim%x
    c(i)=a(i)+b(i)
  end do
end subroutine VectorSum

end module device_code

program exercise03

use cudafor
use device_code

implicit none

! define vector size !
integer, parameter :: vector_size = 100000000
real, dimension(:),pinned,allocatable :: vector_a_cpu
real, dimension(:),pinned,allocatable :: vector_b_cpu
real, dimension(:),pinned,allocatable :: vector_c_cpu
real, dimension(0:vector_size-1), device :: vector_a_gpu
real, dimension(0:vector_size-1), device :: vector_b_gpu
real, dimension(0:vector_size-1), device :: vector_c_gpu
real(kind=8) :: timer1,timer2,durationCopyIn,durationKernel,durationCopyOut
integer :: i
logical :: success
integer :: error

! CPU code !

call random_seed()

allocate(vector_a_cpu(0:vector_size-1))
allocate(vector_b_cpu(0:vector_size-1))
allocate(vector_c_cpu(0:vector_size-1))

! initialize cpu input vector to 0 !
do i=0,vector_size-1
  call random_number(vector_a_cpu(i))
  call random_number(vector_b_cpu(i))
end do
! initialize cpu output vector to 0 !
vector_c_cpu=0

! initialize gpu output vector to 0 !
vector_c_gpu=vector_c_cpu

call cpu_time(timer1)

! copy input vectors a,b to gpu !
vector_a_gpu=vector_a_cpu
vector_b_gpu=vector_b_cpu

call cpu_time(timer2)
durationCopyIn=timer2-timer1
timer1=timer2

! cuda kernel call !
! Note: good values highly depend on the concrete GPU.
call VectorSum<<<2048,512>>>(vector_a_gpu,vector_b_gpu,vector_c_gpu,vector_size)
error = cudaDeviceSynchronize()

call cpu_time(timer2)
durationKernel=timer2-timer1
timer1=timer2

! copy output vector from gpu to cpu !
vector_c_cpu=vector_c_gpu

call cpu_time(timer2)
durationCopyOut=timer2-timer1

! verify results !
success=.true.
do i=0,vector_size-1
  if (abs(vector_a_cpu(i)+vector_b_cpu(i)-vector_c_cpu(i))>1e-5) then
    print *,'computation result is wrong: index=',i,         &
        ' expected result=',vector_a_cpu(i)*vector_b_cpu(i), &
        'cuda result=', vector_c_cpu(i)
    success=.false.
    exit
  end if
end do
if (success) then
  print *,'computation result is correct.'
end if

print *,'duration (copy in)          =',durationCopyIn,'seconds'
print *,'memory bandwidth (copy in)  =',vector_size*8./durationCopyIn*1e-6,'MB/s'
print *,'duration (kernel)           =',durationKernel,'seconds' 
print *,'performance                 =',vector_size/durationKernel*1e-6,'MFLOPS'
print *,'memory bandwidth (kernel)   =',vector_size*12./durationKernel*1e-6,'MB/s'
print *,'duration (copy out)         =',durationCopyOut,'seconds'
print *,'memory bandwidth (copy out) =',vector_size*4./durationCopyOut*1e-6,'MB/s'

end program exercise03