module thrust

interface thrustmin

	function min_integer(input,N) bind(C,name="min_int_wrapper")
		use iso_c_binding
		integer(c_int),device:: input(*)
		integer(c_int),value:: N, min_integer
	end function

	function min_float(input,N) bind(C,name="min_float_wrapper")
		use iso_c_binding
		real(c_float),device:: input(*)
		integer(c_int),value:: N, min_float
	end function

	function min_double(input,N) bind(C,name="min_double_wrapper")
		use iso_c_binding
		real(c_double),device:: input(*)
		integer(c_int),value:: N, min_double
	end function

end interface

interface thrustmax

	function max_integer(input,N) bind(C,name="max_int_wrapper")
		use iso_c_binding
		integer(c_int),device:: input(*)
		integer(c_int),value:: N, max_integer
	end function

	function max_float(input,N) bind(C,name="max_float_wrapper")
		use iso_c_binding
		real(c_float),device:: input(*)
		integer(c_int),value:: N, max_float
	end function

	function max_double(input,N) bind(C,name="max_double_wrapper")
		use iso_c_binding
		real(c_double),device:: input(*)
		integer(c_int),value:: N, max_double
	end function

end interface

interface thrustsort

	subroutine sort_integer(input,N) bind(C,name="sort_int_wrapper")
		use iso_c_binding
		integer(c_int),device:: input(*)
		integer(c_int),value:: N
	end subroutine

	subroutine sort_float(input,N) bind(C,name="sort_float_wrapper")
		use iso_c_binding
		real(c_float),device:: input(*)
		integer(c_int),value:: N
	end subroutine

	subroutine sort_double(input,N) bind(C,name="sort_double_wrapper")
		use iso_c_binding
		real(c_double),device:: input(*)
		integer(c_int),value:: N
	end subroutine

end interface

interface thrustsum

	function reduce_integer(input,N) bind(C,name="reduce_int_wrapper")
		use iso_c_binding
		integer(c_int),device:: input(*)
		integer(c_int) :: reduce_integer
		integer(c_int),value:: N
	end function

	function reduce_float(input,N) bind(C,name="reduce_float_wrapper")
		use iso_c_binding
		real(c_float),device:: input(*)
		real(c_float) :: reduce_float
		integer(c_int),value:: N
	end function

	function reduce_double(input,N) bind(C,name="reduce_double_wrapper")
		use iso_c_binding
		real(c_double),device:: input(*)
		real(c_double) :: reduce_double
		integer(c_int),value:: N
	end function

end interface

end module thrust


program GPU
  use cudafor
  use thrust
  implicit none
  real, allocatable, device :: dev_array(:)
  real, allocatable :: host_array(:)
  real :: varsum
  integer N / 1024 /
  integer :: x
  
  allocate(dev_array(N))
  allocate(host_array(N))

  call random_seed
  call random_number(host_array)
  

  print *, "COPY!"
  host_array = host_array * 10.0
  dev_array  = host_array

  print *, "host_array:"
  print *, host_array
  
  print *, "elements:", N
  !call thrustsort(dev_array, size(dev_array))
  print *, "START CUDA! min:"
  x = thrustmin(dev_array, size(dev_array))
  print *, "Min:", x

  x = thrustmax(dev_array, size(dev_array))
  print *, "Max:", x

  varsum = thrustsum(dev_array, size(dev_array))
  print *, "Sum:", varsum

  print *, "DONE  CUDA!"
  
  
  deallocate(dev_array)
  deallocate(host_array)  

  end program GPU