! File: order_stat_test.f90
! Public domain 2006 James Van Buskirk
module mykinds
   implicit none
   integer, parameter :: sp = selected_real_kind(6,30)
   integer, parameter :: dp = selected_real_kind(15,300)
end module mykinds

module type1_mod
   use mykinds, only: sp
   implicit none
   private
   public type1, assignment(=), operator(<=), get_opcounts_type1

   integer :: opcounts = 0

   type type1
      real(sp) x
   end type type1

   interface assignment(=)
      module procedure type1_assign
   end interface assignment(=)

   interface operator(<=)
      module procedure type1_LE
   end interface operator(<=)

   contains
      elemental subroutine type1_assign(lhs, rhs)
         type(type1), intent(out) :: lhs
         type(type1), intent(in) :: rhs

         lhs%x = rhs%x
      end subroutine type1_assign

      function type1_LE(x, y)
         logical type1_LE
         type(type1), intent(in) :: x, y

         type1_LE = x%x <= y%x
         opcounts = opcounts+1
      end function type1_LE

      function get_opcounts_type1()
         integer get_opcounts_type1

         get_opcounts_type1 = opcounts
      end function get_opcounts_type1
end module type1_mod

module type2_mod
   implicit none
   private
   public type2, assignment(=), operator(<=), get_opcounts_type2

   integer :: opcounts = 0

   type type2
      integer x
   end type type2

   interface assignment(=)
      module procedure type2_assign
   end interface assignment(=)

   interface operator(<=)
      module procedure type2_LE
   end interface operator(<=)

   contains
      elemental subroutine type2_assign(lhs, rhs)
         type(type2), intent(out) :: lhs
         type(type2), intent(in) :: rhs

         lhs%x = rhs%x
      end subroutine type2_assign

      function type2_LE(x, y)
         logical type2_LE
         type(type2), intent(in) :: x, y

         type2_LE = x%x <= y%x
         opcounts = opcounts+1
      end function type2_LE

      function get_opcounts_type2()
         integer get_opcounts_type2

         get_opcounts_type2 = opcounts
      end function get_opcounts_type2
end module type2_mod

module type1_order_stat
   use type1_mod, T1 =>type1
   implicit none
   private
   public select
   contains
include 'order_stat.i90'
end module type1_order_stat

module type2_order_stat
   use type2_mod, T1 =>type2
   implicit none
   private
   public select
   contains
include 'order_stat.i90'
end module type2_order_stat

module generic_recombination
   use type1_order_stat, select_type1 => select
   use type2_order_stat, select_type2 => select
   implicit none
   interface select
      module procedure select_type1, select_type2
   end interface select
end module generic_recombination

program test
   use type1_mod, only: type1, get_opcounts_type1
   use type2_mod, only: type2, get_opcounts_type2
   use mykinds, only: sp
   use generic_recombination
   implicit none
   integer N
   type(type1), allocatable :: array1(:)
   type(type2), allocatable :: array2(:)
   integer i
   integer opcountsi, opcountsf
   real(sp) harvest
   real(sp) median

   write(*,'(a)',advance='no') 'Enter the size of the arrays:> '
   read(*,*) N
   if(N < 1) then
      write(*,'(a,i0,a)') "Sorry, arrays of size = ", N, &
         " are not allowed"
      stop
   end if
   allocate(array1(N), array2(N))
   call random_seed()
   do i = 1, N
      call random_number(harvest)
      array1(i)%x = harvest
      array2(i)%x = harvest*(N/10)
   end do
   opcountsi = get_opcounts_type1()
   if(mod(N,2) == 0) then
      call select(array1, N, (/ N/2, N/2+1 /))
      median = 0.5*(array1(N/2)%x+array1(N/2+1)%x)
   else
      call select(array1, N, (/ (N+1)/2 /))
      median = array1((N+1)/2)%x
   end if
   opcountsf = get_opcounts_type1()
   write(*,'(a)') 'Results for array1:'
   write(*,'(a,i0)') 'Number of comparisons = ',opcountsf-opcountsi
   if(N < 60) then
      do i = 1, N
         write(*,*) i, array1(i)%x <= median, array1(i)
      end do
   else
      write(*,'(a,i0,a,i0)') 'From ',1,' to ',N/2
      write(*,'(i0,a,f0.6)') count(array1(1:N/2)%x <= median),' are <= ',median
      write(*,'(i0,a,f0.6)') count(array1(1:N/2)%x > median),' are > ',median
      write(*,'(a,i0,a,i0)') 'From ',(N+3)/2,' to ',N
      write(*,'(i0,a,f0.6)') count(array1((N+3)/2:N)%x < median),' are < ',median
      write(*,'(i0,a,f0.6)') count(array1((N+3)/2:N)%x >= median),' are >= ',median
   end if
   opcountsi = get_opcounts_type2()
   if(mod(N,2) == 0) then
      call select(array2, N, (/ N/2, N/2+1 /))
      median = 0.5*(array2(N/2)%x+array2(N/2+1)%x)
   else
      call select(array2, N, (/ (N+1)/2 /))
      median = array2((N+1)/2)%x
   end if
   opcountsf = get_opcounts_type2()
   write(*,'(a)') 'Results for array2:'
   write(*,'(a,i0)') 'Number of comparisons = ',opcountsf-opcountsi
   if(N < 60) then
      do i = 1, N
         write(*,*) i, array2(i)%x <= median, array2(i)
      end do
   else
      write(*,'(a,i0,a,i0)') 'From ',1,' to ',N/2
      write(*,'(i0,a,f0.6)') count(array2(1:N/2)%x <= median),' are <= ',median
      write(*,'(i0,a,f0.6)') count(array2(1:N/2)%x > median),' are > ',median
      write(*,'(a,i0,a,i0)') 'From ',(N+3)/2,' to ',N
      write(*,'(i0,a,f0.6)') count(array2((N+3)/2:N)%x < median),' are < ',median
      write(*,'(i0,a,f0.6)') count(array2((N+3)/2:N)%x >= median),' are >= ',median
   end if
end program test
! End of file: order_stat_test.f90