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