! 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