module square_kernel
 implicit none
contains
 ! 簡単な kernel プログラム
 ! Square the entries of 'a'
 attributes(global) subroutine s_kernel( a )
  real, dimension(*) :: a
  integer i
  i = (blockidx%x-1)*blockdim%x + threadidx%x
  a(i) = a(i) * a(i)
 end subroutine
end module

program main
 ! C言語関数との interface が必須となるため iso_c_binding 使用
 use, intrinsic :: iso_c_binding
 use cudafor
 use square_kernel
 implicit none

 ! CUDA Runtime API の中のC言語ベースのライブラリを使用するため、C ポインタ体系の世界
 ! そこで、Fortran上では、C ポインタ と Fortran ポインタの変換が必要となる
 ! 二つのポインタを用意しておく
 real, dimension(:), pointer :: ha      ! host pointer (Fortran pointer)
 type(c_ptr) :: local_ha                ! C pointer which we convert to ha
 real, dimension(:), device, allocatable :: da   ! device pointer (Fortran pointer)
 type(c_devptr) :: local_da             ! C pointer which we convert to da

 integer :: i, n, err, istat
 real :: e

 ! Properties check
 type(cudaDeviceProp) :: prop
 integer :: idevice, ilen

 ! Device properties
 idevice = 0
 istat = cudaSetDevice(idevice)
 istat = cudaGetDeviceProperties(prop,idevice)
 ilen = verify(prop%name, ' ', .true.)
 write (*,900) prop%name(1:ilen), &
               real(prop%clockRate)/1000.0, &
               real(prop%totalGlobalMem)/1024.0/1024.0, &
               prop%canMapHostMemory
900 format(/,' Device name: ',a,/,1x, f6.1,' MHz clock, ',f6.1,' MB memory.'/, &
           ' canMapHostMemory Flag : ', i2/)

     ! ホスト側で pinned メモリの使用が可能かどうか確認する
     ! UVA 下では、必ず mapped (Zero-copy memory) ホストメモリを割り付けることを指示する
     ! (cudaDeviceMapHost)
     ! This flag must be set in order to allocate pinned host memory
     ! that is accessible to the device.

 if (prop%canMapHostMemory == 1) then
  istat = cudasetdeviceflags(cudaDeviceMapHost)
 !istat = cudasetdeviceflags(cudaDeviceScheduleBlockingSync)
 end if

 n = 2048

     ! cudahostalloc を使用し、ホスト側に pinned and mapped memory 領域をアロケートする
     ! その際の CUDA API 管理下の（世界の）ホスト側領域のポインタを local_ha とする
     ! Allocate (mapped) pinned host memory.
     ! Maps the allocation into the CUDA address space. The device pointer to the memory
     ! may be obtained by calling cudaHostGetDevicePointer().
     ! local_ha : type(c_ptr) = C pointer

 istat = cudahostalloc( local_ha, sizeof(1.0)*n, cudaHostAllocMapped )

     ! CUDA API の世界のポインタ local_ha を Fortran の世界のポインタに変換する
     ! Convert C(local_ha) pointer to Fortran pointer

 call c_f_pointer( local_ha, ha, (/ n /) )

     ! ホスト側の配列 ha にデータを入れる＝ local_ha の場所にデータが入る
 do i = 1, n
  ha(i) = i+2
 enddo

     ! ホスト側の local_ha ポインタ領域を UVA 下の「一つのデバイス」領域として認識し、
     ! UVA上のポインタ(local_da) に付け替える
     ! convert the host pointer to a device address
     ! Passes back the device pointer corresponding to the mapped,
     ! pinned host buffer allocated by cudaHostAlloc() .
     ! local_da : Returned device pointer for mapped memory : type(C_DEVPTR)
     ! local_ha : Requested host pointer mapping : type(C_PTR)

 istat = cudahostgetdevicepointer( local_da, local_ha, 0 )

     ! Convert C pointer to Fortran pointer

 call c_f_pointer( local_da, da, (/ n /) )
 
      ! GPUデバイスには、UVA 管理下のポインタ(local_da)を渡すことが必要
      ! プログラムは Fortran の世界で coding しているため、s_kernel への引数には、
      ! Fortran の世界の da ポインタを渡す必要あり
      
 call s_kernel <<< 16,128 >>> ( da )
     
     ! Wait for the kernel to finish before testing the results
 istat = cudaThreadSynchronize();

     ! UVA 下では、デバイス上の演算結果は、 da ポインタ= ha ポインタが示すホスト側領域に
     ! ストアされているため、デバイス側からホストへのデータ転送は必要ない
     ! At thist point, da(i) ( = ha(i) ) should be (i+2)**2
 err = 0
 do i = 1, n
  e = float(i+2)
  e = e * e
  if( ha(i) - e .ne. 0 )then
   if( err <= 10 )then
    write(*,10) i, ha(i), e
10  format( 'a(',i4,') = ', e14.8, ', expecting ', e14.8 )
   endif
  endif
 enddo
 print *, err, ' errors found'

 ! shutdown ... host pinned memory
 istat = cudaFreehost( local_ha )
!print *, 'error code = ', cudaGetErrorString( cudaGetLastError() )

end program