! File: order_stat.i90
! Public domain 2006 James Van Buskirk
! Include file for generic order statistics
      recursive subroutine select(array, N, k)
         integer, intent(in) :: N
         type(T1) array(N)
         integer, intent(in) :: k(:)
         type(T1) med_loc((N+4)/5)
         integer med_last
         integer i
         integer N_last
         type(T1) mark
         integer mark_count
         integer lo_insert, hi_insert
         type(T1) lo_hold(10), hi_hold(10)
         integer lo_temp, hi_temp
         integer i_hi, i_lo
         integer j
         type(T1), allocatable :: merge_temp(:)
         character(80) fmt

! Handle wierd cases
         if(any(k <= 0 .OR. k > N)) then
            j = count(k <= 0 .OR. k > N)
            if(j == 1) then
               write(fmt,'(a)') '(a,i0,a)'
            else
               write(fmt,'(a,i0,a)') '(a,',j,'(i0,","),i0,a)'
            end if
            write(*,fmt) "Invalid k = (/", pack(k,k <= 0 .OR. k > N), &
               "/) detected in subroutine select"
            stop
         else if(all(k == 1)) then
            i_lo = sum(minloc(array))
            lo_hold(1) = array(i_lo)
            array(i_lo) = array(1)
            array(1) = lo_hold(1)
            return
         else if(any(k < N)) then
            ! Normal case follows IF block
         else ! all(k == N)
            i_hi = sum(maxloc(array))
            hi_hold(1) = array(i_hi)
            array(i_hi) = array(N)
            array(N) = hi_hold(1)
            return
         end if

! Normal case
! Handle small arrays via sorting
         if(N <= 20) then
            allocate(merge_temp((N+1)/2))
            call merge_sort(array, N, merge_temp)
            return
         end if
! Find medians of subsets of 5
         do i = 1, N/5
            call median_5(array(5*i-4:5*i))
            med_loc(i) = array(5*i-2)
         end do
         N_last = mod(N,5)
         med_last = N+1-(N_last+1)/2
         if(N_last > 0) then
            call sort_short(array(N-N_last+1:N),N_last)
            med_loc(i) = array(med_last)
         end if
! Find median of medians
         call select(med_loc,size(med_loc),(/(size(med_loc)+1)/2/))
! Partition array around mark = median of medians
         mark = med_loc((size(med_loc)+1)/2)
! Initialize partition code variables
         mark_count = 0
         lo_insert = 1
         hi_insert = N
         lo_temp = 1
         hi_temp = 1
! Handle possible deficient subset
         if(N_last > 0) then
            if(mark <= array(med_last)) then
               if(array(med_last) <= mark) then
                  mark_count = mark_count+1
                  med_loc(mark_count) = array(med_last)
                  hi_insert = med_last
               else
                  hi_insert = med_last-1
               end if
               do j = med_last-1, N-N_last+1, -1
                  if(array(j) <= mark) then
                     lo_hold(lo_temp) = array(j)
                     lo_temp = lo_temp+1
                  else
                     array(hi_insert) = array(j)
                     hi_insert = hi_insert-1
                  end if
               end do
            else
               do j = N, med_last+1, -1
                  if(array(j) <= mark) then
                     lo_hold(lo_temp) = array(j)
                     lo_temp = lo_temp+1
                  else
                     array(hi_insert) = array(j)
                     hi_insert = hi_insert-1
                  end if
               end do
               do j = med_last, N-N_last+1, -1
                  lo_hold(lo_temp) = array(j)
                  lo_temp = lo_temp+1
               end do
            end if
         end if
         i_hi = i-1
         i_lo = 1
! Loop for general subsets
         outer_partition: do
            inner_partition_lo: do
            if(i_lo > i_hi) exit outer_partition
               if(array(5*i_lo-2) <= mark) then
                  if(mark <= array(5*i_lo-2)) then
                     mark_count = mark_count+1
                     med_loc(mark_count) = array(5*i_lo-2)
                     do j = 5*i_lo-4, 5*i_lo-3
                        array(lo_insert) = array(j)
                        lo_insert = lo_insert+1
                     end do
                  else
                     do j = 5*i_lo-4, 5*i_lo-2
                        array(lo_insert) = array(j)
                        lo_insert = lo_insert+1
                     end do
                  end if
                  do j = 5*i_lo-1, 5*i_lo
                     if(array(j) <= mark) then
                        array(lo_insert) = array(j)
                        lo_insert = lo_insert+1
                     else
                        hi_hold(hi_temp) = array(j)
                        hi_temp = hi_temp+1
                     end if
                  end do
               else
                  do j = 5*i_lo-4, 5*i_lo-3
                     if(array(j) <= mark) then
                        array(lo_insert) = array(j)
                        lo_insert = lo_insert+1
                     else
                        hi_hold(hi_temp) = array(j)
                        hi_temp = hi_temp+1
                     end if
                  end do
                  do j = 5*i_lo-2, 5*i_lo
                     hi_hold(hi_temp) = array(j)
                     hi_temp = hi_temp+1
                  end do
               end if
               i_lo = i_lo+1
            if(hi_temp > 5) exit inner_partition_lo
            end do inner_partition_lo
            inner_partition_hi: do
            if(i_lo > i_hi) exit outer_partition
               if(mark <= array(5*i_hi-2)) then
                  if(array(5*i_hi-2) <= mark) then
                     mark_count = mark_count+1
                     med_loc(mark_count) = array(5*i_hi-2)
                     do j = 5*i_hi, 5*i_hi-1, -1
                        array(hi_insert) = array(j)
                        hi_insert = hi_insert-1
                     end do
                  else
                     do j = 5*i_hi, 5*i_hi-2, -1
                        array(hi_insert) = array(j)
                        hi_insert = hi_insert-1
                     end do
                  end if
                  do j = 5*i_hi-3, 5*i_hi-4, -1
                     if(mark <= array(j)) then
                        array(hi_insert) = array(j)
                        hi_insert = hi_insert-1
                     else
                        lo_hold(lo_temp) = array(j)
                        lo_temp = lo_temp+1
                     end if
                  end do
               else
                  do j = 5*i_hi, 5*i_hi-1, -1
                     if(mark <= array(j)) then
                        array(hi_insert) = array(j)
                        hi_insert = hi_insert-1
                     else
                        lo_hold(lo_temp) = array(j)
                        lo_temp = lo_temp+1
                     end if
                  end do
                  do j = 5*i_hi-2, 5*i_hi-4, -1
                     lo_hold(lo_temp) = array(j)
                     lo_temp = lo_temp+1
                  end do
               end if
               i_hi = i_hi-1
            if(lo_temp > 5) exit inner_partition_hi
            end do inner_partition_hi
            do j = 1, 5
               lo_temp = lo_temp-1
               array(lo_insert) = lo_hold(lo_temp)
               lo_insert = lo_insert+1
               hi_temp = hi_temp-1
               array(hi_insert) = hi_hold(hi_temp)
               hi_insert = hi_insert-1
            end do
         end do outer_partition
         do lo_temp = lo_temp-1, 1, -1
            array(lo_insert) = lo_hold(lo_temp)
            lo_insert = lo_insert+1
         end do
         do hi_temp = hi_temp-1, 1, -1
            array(hi_insert) = hi_hold(hi_temp)
            hi_insert = hi_insert-1
         end do
         array(lo_insert:hi_insert) = med_loc(1:mark_count)
         if(any(k < lo_insert)) then
            call select(array(1:lo_insert-1), lo_insert-1, &
               pack(k, k < lo_insert))
         end if
         if(any(k > hi_insert)) then
            call select(array(hi_insert+1:N), N-hi_insert, &
               pack(k, k > hi_insert)-hi_insert)
         end if

         contains
            subroutine median_5(A5)
               type(T1) A5(5)

               if(.NOT. A5(1) <= A5(2)) then
                  A5(1:2) = A5((/2,1/))
               end if
               if(.NOT. A5(2) <= A5(3)) then
                  if(.NOT. A5(1) <= A5(3)) then
                     A5(1:3) = A5((/3,1,2/))
                  else
                     A5(2:3) = A5((/3,2/))
                  end if
               end if
               if(.NOT. A5(4) <= A5(5)) then
                  A5(4:5) = A5((/5,4/))
               end if
               if(A5(4) <= A5(2)) then
                  if(A5(2) <= A5(5)) then
                     A5(2:4) = A5((/4,2,3/))
                  else if(A5(5) <= A5(1)) then
                     A5 = A5((/4,5,1,2,3/))
                  else
                     A5(2:5) = A5((/4,5,2,3/))
                  end if
               else
                  if(.NOT. A5(3) <= A5(4)) then
                     A5(3:4) = A5((/4,3/))
                  end if
               end if
            end subroutine median_5

            subroutine sort_short(A4,N)
               integer, intent(in) :: N
               type(T1) :: A4(N)
               integer i, j

               do j = N, 2, -1
                  do i = 1, j-1
                     if(.NOT. A4(i) <= A4(i+1)) then
                        A4(i:i+1) = A4((/i+1,i/))
                     end if
                  end do
               end do
            end subroutine sort_short
      end subroutine select

      recursive subroutine merge_sort(array, N, temp)
         integer, intent(in) :: N
         type(T1) array(N)
         type(T1) temp(*)
         integer lo_remove, hi_remove
         integer j

! Handle small arrays
         if(N < 0) then
            write(*,'(a,i0,a)') "Invalid array size = ", N, &
               "detected in subroutine merge_sort"
            stop
         else if(N < 4) then
            select case(N)
            case(0:1)
            case(2)
               if(.NOT. array(1) <= array(2)) then
                  array(1:2) = array((/2,1/))
               end if
            case(3)
               if(.NOT. array(1) <= array(2)) then
                  array(1:2) = array((/2,1/))
               end if
               if(.NOT. array(2) <= array(3)) then
                  if(.NOT. array(1) <= array(3)) then
                     array(1:3) = array((/3,1,2/))
                  else
                     array(2:3) = array((/3,2/))
                  end if
               end if
            end select
            return
         end if
! General case
! Merge sort the two halves of the array
         temp(1:N/2) = array(1:N/2)
         call merge_sort(temp(1:N/2), N/2, array(1:N/2))
         call merge_sort(array(N/2+1:N), (N+1)/2, array(1:N/2))
! Merge the two sorted halves
         lo_remove = 1
         hi_remove = N/2+1
         do j = 1, N-1
            if(temp(lo_remove) <= array(hi_remove)) then
               array(j) = temp(lo_remove)
               lo_remove = lo_remove+1
               if(lo_remove > N/2) then
                  return
               end if
            else
               array(j) = array(hi_remove)
               hi_remove = hi_remove+1
               if(hi_remove > N) then
                  array(j+1:N) = temp(lo_remove:N/2)
                  return
               end if
            end if
         end do
      end subroutine merge_sort

      function minloc(array)
         integer minloc(1)
         type(T1), intent(in) :: array(:)
         integer i
         integer curr_min

         if(size(array) == 0) then
            minloc = 0
         else
            curr_min = 1
            do i = 2, size(array)
               if(.NOT.array(curr_min) <= array(i)) then
                  curr_min = i
               end if
            end do
            minloc = curr_min
         end if
      end function minloc

      function maxloc(array)
         integer maxloc(1)
         type(T1), intent(in) :: array(:)
         integer i
         integer curr_max

         if(size(array) == 0) then
            maxloc = 0
         else
            curr_max = 1
            do i = 2, size(array)
               if(.NOT.array(i) <= array(curr_max)) then
                  curr_max = i
               end if
            end do
            maxloc = curr_max
         end if
      end function maxloc
! End of file: order_stat.i90