You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
855 lines
20 KiB
855 lines
20 KiB
! These various sort routines are provided as
|
|
! examples, rather than being fully optimised
|
|
! for the best possible performance. MJR 9/2019
|
|
!
|
|
! So far the collection is:
|
|
!
|
|
! bubble_sort (order N**2)
|
|
! selection_sort (order N**2)
|
|
! inline_selection_sort (order N**2)
|
|
! insertion_sort (order N**2)
|
|
! binary_insertion_sort (order N**2)
|
|
! pair_insertion_sort (order N**2)
|
|
! double_insertion_sort (order N**2)
|
|
! shell_sort (order N**(4/3)) (possibly)
|
|
! merge_sort (order N log N)
|
|
! alt_merge_sort (order N log N)
|
|
! heap_sort (order N log N)
|
|
! heap_sort_nr (order N log N)
|
|
! smoothsort (order N log N)
|
|
! quicksort (order N log N) (usually)
|
|
! quicksort_nr (order N log N) (usually)
|
|
! dpquicksort (order N log N) (usually)
|
|
|
|
! 14/10/19 and 18/10/19 indexing errors in heap sorts corrected
|
|
|
|
! 4/4/20 smooth sort and non-recursive quicksort added, and some
|
|
! minor code tidying
|
|
|
|
|
|
module sorts
|
|
use parameters
|
|
implicit none
|
|
integer, parameter :: prec=dp, i64=selected_int_kind(15)
|
|
|
|
! array L used by smoothsort
|
|
|
|
integer, parameter :: L(0:43) = (/ 1, 1, 3, 5, 9, 15, 25, 41, 67, 109, &
|
|
177, 287, 465, 753, 1219, 1973, 3193, 5167, 8361, &
|
|
13529, 21891, 35421, 57313, 92735, 150049, &
|
|
242785, 392835, 635621, 1028457, 1664079, 2692537, &
|
|
4356617, 7049155, 11405773, 18454929, 29860703, &
|
|
48315633, 78176337, 126491971, 204668309, 331160281, &
|
|
535828591, 866988873, 1402817465 /)
|
|
|
|
|
|
private :: binary_search, merge_arrays, heapify, heapify_nr, &
|
|
smooth_heapify, smooth_extract, sift_in, interheap_sift, L
|
|
|
|
contains
|
|
|
|
recursive subroutine quicksort(array)
|
|
real(prec), intent(inout)::array(:)
|
|
real(prec) :: temp,pivot
|
|
integer :: i,j,last,left,right
|
|
|
|
last=size(array)
|
|
|
|
if (last.lt.50) then ! use insertion sort on small arrays
|
|
do i=2,last
|
|
temp=array(i)
|
|
do j=i-1,1,-1
|
|
if (array(j).le.temp) exit
|
|
array(j+1)=array(j)
|
|
enddo
|
|
array(j+1)=temp
|
|
enddo
|
|
return
|
|
endif
|
|
! find median of three pivot
|
|
! and place sentinels at first and last elements
|
|
temp=array(last/2)
|
|
array(last/2)=array(2)
|
|
if (temp.gt.array(last)) then
|
|
array(2)=array(last)
|
|
array(last)=temp
|
|
else
|
|
array(2)=temp
|
|
endif
|
|
if (array(1).gt.array(last)) then
|
|
temp=array(1)
|
|
array(1)=array(last)
|
|
array(last)=temp
|
|
endif
|
|
if (array(1).gt.array(2)) then
|
|
temp=array(1)
|
|
array(1)=array(2)
|
|
array(2)=temp
|
|
endif
|
|
pivot=array(2)
|
|
|
|
left=3
|
|
right=last-1
|
|
do
|
|
do while(array(left).lt.pivot)
|
|
left=left+1
|
|
enddo
|
|
do while(array(right).gt.pivot)
|
|
right=right-1
|
|
enddo
|
|
if (left.ge.right) exit
|
|
temp=array(left)
|
|
array(left)=array(right)
|
|
array(right)=temp
|
|
left=left+1
|
|
right=right-1
|
|
enddo
|
|
if (left.eq.right) left=left+1
|
|
call quicksort(array(1:left-1))
|
|
call quicksort(array(left:))
|
|
|
|
end subroutine quicksort
|
|
|
|
! This version maintains its own stack, to avoid needing to call
|
|
! itself recursively. By always pushing the larger "half" to the
|
|
! stack, and moving directly to calculate the smaller "half",
|
|
! it can guarantee that the stack needs no more than log_2(N)
|
|
! entries
|
|
subroutine quicksort_nr(array)
|
|
real(prec), intent(inout)::array(:)
|
|
real(prec) :: temp,pivot
|
|
integer :: i,j,left,right,low,high
|
|
! If your compiler lacks storage_size(), replace
|
|
! storage_size(i) by 64
|
|
integer :: stack(2,storage_size(i)),stack_ptr
|
|
|
|
low=1
|
|
high=size(array)
|
|
stack_ptr=1
|
|
|
|
do
|
|
|
|
if (high-low.lt.50) then ! use insertion sort on small arrays
|
|
do i=low+1,high
|
|
temp=array(i)
|
|
do j=i-1,low,-1
|
|
if (array(j).le.temp) exit
|
|
array(j+1)=array(j)
|
|
enddo
|
|
array(j+1)=temp
|
|
enddo
|
|
! now pop from stack
|
|
if (stack_ptr.eq.1) return
|
|
stack_ptr=stack_ptr-1
|
|
low=stack(1,stack_ptr)
|
|
high=stack(2,stack_ptr)
|
|
cycle
|
|
endif
|
|
|
|
! find median of three pivot
|
|
! and place sentinels at first and last elements
|
|
temp=array((low+high)/2)
|
|
array((low+high)/2)=array(low+1)
|
|
if (temp.gt.array(high)) then
|
|
array(low+1)=array(high)
|
|
array(high)=temp
|
|
else
|
|
array(low+1)=temp
|
|
endif
|
|
if (array(low).gt.array(high)) then
|
|
temp=array(low)
|
|
array(low)=array(high)
|
|
array(high)=temp
|
|
endif
|
|
if (array(low).gt.array(low+1)) then
|
|
temp=array(low)
|
|
array(low)=array(low+1)
|
|
array(low+1)=temp
|
|
endif
|
|
pivot=array(low+1)
|
|
|
|
left=low+2
|
|
right=high-1
|
|
do
|
|
do while(array(left).lt.pivot)
|
|
left=left+1
|
|
enddo
|
|
do while(array(right).gt.pivot)
|
|
right=right-1
|
|
enddo
|
|
if (left.ge.right) exit
|
|
temp=array(left)
|
|
array(left)=array(right)
|
|
array(right)=temp
|
|
left=left+1
|
|
right=right-1
|
|
enddo
|
|
if (left.eq.right) left=left+1
|
|
! call quicksort(array(1:left-1))
|
|
! call quicksort(array(left:))
|
|
if (left.lt.(low+high)/2) then
|
|
stack(1,stack_ptr)=left
|
|
stack(2,stack_ptr)=high
|
|
stack_ptr=stack_ptr+1
|
|
high=left-1
|
|
else
|
|
stack(1,stack_ptr)=low
|
|
stack(2,stack_ptr)=left-1
|
|
stack_ptr=stack_ptr+1
|
|
low=left
|
|
endif
|
|
|
|
enddo
|
|
end subroutine quicksort_nr
|
|
|
|
|
|
! dual pivot quicksort
|
|
recursive subroutine dpquicksort(array)
|
|
real(prec), intent(inout)::array(:)
|
|
real(prec) :: temp,p1,p2
|
|
integer :: i,j,last,l,k,g
|
|
|
|
last=size(array)
|
|
|
|
if (last.lt.40) then ! use insertion sort on small arrays
|
|
do i=2,last
|
|
temp=array(i)
|
|
do j=i-1,1,-1
|
|
if (array(j).le.temp) exit
|
|
array(j+1)=array(j)
|
|
enddo
|
|
array(j+1)=temp
|
|
enddo
|
|
return
|
|
endif
|
|
p1=array(last/3)
|
|
p2=array(2*last/3)
|
|
if (p2.lt.p1) then
|
|
temp=p1
|
|
p1=p2
|
|
p2=temp
|
|
endif
|
|
array(last/3)=array(1)
|
|
array(1)=p1
|
|
array(2*last/3)=array(last)
|
|
array(last)=p2
|
|
|
|
g=last
|
|
l=2
|
|
do while (array(l).lt.p1)
|
|
l=l+1
|
|
enddo
|
|
k=l
|
|
|
|
do while(k.lt.g)
|
|
temp=array(k)
|
|
if (temp.lt.p1) then
|
|
array(k)=array(l)
|
|
array(l)=temp
|
|
l=l+1
|
|
else if (temp.gt.p2) then
|
|
do while(array(g-1).gt.p2)
|
|
g=g-1
|
|
enddo
|
|
if (k.ge.g) exit
|
|
g=g-1
|
|
if (array(g).lt.p1) then
|
|
array(k)=array(l)
|
|
array(l)=array(g)
|
|
array(g)=temp
|
|
l=l+1
|
|
else
|
|
array(k)=array(g)
|
|
array(g)=temp
|
|
endif
|
|
endif
|
|
k=k+1
|
|
enddo
|
|
if (l.gt.2) then
|
|
array(1)=array(l-1)
|
|
array(l-1)=p1
|
|
call dpquicksort(array(1:l-2))
|
|
endif
|
|
call dpquicksort(array(l:g-1))
|
|
if (g.lt.last) then
|
|
array(last)=array(g)
|
|
array(g)=p2
|
|
call dpquicksort(array(g+1:last))
|
|
endif
|
|
|
|
end subroutine dpquicksort
|
|
|
|
function binary_search(array, x, start, end)
|
|
integer binary_search
|
|
integer, intent(in) :: start,end
|
|
real(prec), intent(in) :: x,array(:)
|
|
integer :: a,b,mid
|
|
|
|
a=start
|
|
b=end
|
|
|
|
do
|
|
if (a.eq.b) then
|
|
if (array(a).gt.x) then
|
|
binary_search=a
|
|
return
|
|
else
|
|
binary_search=a+1
|
|
return
|
|
endif
|
|
endif
|
|
if (a.gt.b) then
|
|
binary_search=a
|
|
return
|
|
end if
|
|
|
|
|
|
mid=(a+b)/2
|
|
if (array(mid).lt.x) then
|
|
a=mid+1
|
|
else if (array(mid).gt.x) then
|
|
b=mid-1
|
|
else
|
|
binary_search=mid
|
|
return
|
|
endif
|
|
|
|
end do
|
|
end function binary_search
|
|
|
|
subroutine binary_insertion_sort(array)
|
|
real(prec), intent(inout) :: array(:)
|
|
integer :: i,j,pos
|
|
real(prec) :: x
|
|
|
|
do i=2,size(array)
|
|
x=array(i)
|
|
pos=binary_search(array,x,1,i-1)
|
|
do j=i,pos+1,-1
|
|
array(j)=array(j-1)
|
|
enddo
|
|
array(pos)=x
|
|
enddo
|
|
end subroutine binary_insertion_sort
|
|
|
|
subroutine insertion_sort(array)
|
|
real(prec), intent(inout) :: array(:)
|
|
integer :: i,j
|
|
real(prec) :: temp
|
|
|
|
do i=2,size(array)
|
|
temp=array(i)
|
|
do j=i-1,1,-1
|
|
if (array(j).le.temp) exit
|
|
array(j+1)=array(j)
|
|
enddo
|
|
array(j+1)=temp
|
|
enddo
|
|
end subroutine insertion_sort
|
|
|
|
subroutine double_insertion_sort(array)
|
|
real(prec), intent(inout) :: array(:)
|
|
integer :: i,j,left,right,last
|
|
real(prec) :: temp,p,next
|
|
|
|
last=size(array)
|
|
|
|
p=0.5*(array(1)+array(last))
|
|
if (array(1).gt.array(last)) then
|
|
temp=array(last)
|
|
array(last)=array(1)
|
|
array(1)=temp
|
|
endif
|
|
|
|
left=1
|
|
right=last
|
|
temp=array(2)
|
|
|
|
do i=2,last-1
|
|
if (temp.lt.p) then
|
|
do j=left,1,-1
|
|
if (array(j).le.temp) exit
|
|
array(j+1)=array(j)
|
|
enddo
|
|
array(j+1)=temp
|
|
temp=array(left+2)
|
|
left=left+1
|
|
else
|
|
next=array(right-1)
|
|
do j=right,last
|
|
if (array(j).ge.temp) exit
|
|
array(j-1)=array(j)
|
|
enddo
|
|
array(j-1)=temp
|
|
temp=next
|
|
right=right-1
|
|
endif
|
|
enddo
|
|
|
|
end subroutine double_insertion_sort
|
|
|
|
|
|
subroutine pair_insertion_sort(array)
|
|
real(prec), intent(inout) :: array(:)
|
|
integer :: i,j,last
|
|
real(prec) :: t1,t2
|
|
|
|
last=size(array)
|
|
do i=2,last-1,2
|
|
t1=min(array(i),array(i+1))
|
|
t2=max(array(i),array(i+1))
|
|
j=i-1
|
|
do while((j.ge.1).and.(array(j).gt.t2))
|
|
array(j+2)=array(j)
|
|
j=j-1
|
|
enddo
|
|
array(j+2)=t2
|
|
do while((j.ge.1).and.(array(j).gt.t1))
|
|
array(j+1)=array(j)
|
|
j=j-1
|
|
enddo
|
|
array(j+1)=t1
|
|
end do
|
|
|
|
if(mod(last,2).eq.0)then
|
|
t1=array(last)
|
|
do j=last-1,1,-1
|
|
if (array(j).le.t1) exit
|
|
array(j+1)=array(j)
|
|
end do
|
|
array(j+1)=t1
|
|
endif
|
|
|
|
end subroutine pair_insertion_sort
|
|
|
|
|
|
subroutine merge_arrays(a,b,out)
|
|
! routine used by merge_sort()
|
|
real(prec), intent(in) :: a(:),b(:)
|
|
real(prec), intent(out) :: out(:)
|
|
integer :: ai,bi,oi,i
|
|
|
|
ai=1
|
|
bi=1
|
|
oi=1
|
|
|
|
do
|
|
if (ai.gt.size(a)) then
|
|
do i=bi,size(b)
|
|
out(oi)=b(i)
|
|
oi=oi+1
|
|
enddo
|
|
return
|
|
endif
|
|
if (bi.gt.size(b)) then
|
|
do i=ai,size(a)
|
|
out(oi)=a(i)
|
|
oi=oi+1
|
|
enddo
|
|
return
|
|
endif
|
|
|
|
if (a(ai).lt.b(bi)) then
|
|
out(oi)=a(ai)
|
|
ai=ai+1
|
|
else
|
|
out(oi)=b(bi)
|
|
bi=bi+1
|
|
endif
|
|
oi=oi+1
|
|
enddo
|
|
end subroutine merge_arrays
|
|
|
|
recursive subroutine merge_sort(array)
|
|
! With a little skill, the extra copy of "array=temp"
|
|
! can be avoided
|
|
real(prec), intent(inout) :: array(:)
|
|
real(prec), allocatable :: temp(:)
|
|
integer :: last
|
|
|
|
last=size(array)
|
|
|
|
if (last.lt.40) then
|
|
call insertion_sort(array)
|
|
return
|
|
endif
|
|
|
|
call merge_sort(array(1:last/2))
|
|
call merge_sort(array(last/2+1:last))
|
|
allocate(temp(last))
|
|
call merge_arrays(array(1:last/2),array(last/2+1:last),temp)
|
|
array=temp
|
|
deallocate(temp)
|
|
end subroutine merge_sort
|
|
|
|
recursive subroutine alt_merge_sort(array)
|
|
! Alternate merges between array and temp
|
|
! to avoid copy
|
|
real(prec), intent(inout) :: array(:)
|
|
real(prec), allocatable :: temp(:)
|
|
integer :: last
|
|
|
|
last=size(array)
|
|
|
|
if (last.lt.40) then
|
|
call insertion_sort(array)
|
|
return
|
|
endif
|
|
|
|
call alt_merge_sort(array(1:last/4))
|
|
call alt_merge_sort(array(last/4+1:last/2))
|
|
call alt_merge_sort(array(last/2+1:last/2+last/4))
|
|
call alt_merge_sort(array(last/2+last/4+1:last))
|
|
|
|
allocate(temp(last))
|
|
call merge_arrays(array(1:last/4),array(last/4+1:last/2),temp(1:last/2))
|
|
call merge_arrays(array(last/2+1:last/2+last/4), &
|
|
array(last/2+last/4+1:last),temp(last/2+1:last))
|
|
call merge_arrays(temp(1:last/2),temp(last/2+1:last),array)
|
|
deallocate(temp)
|
|
end subroutine alt_merge_sort
|
|
|
|
subroutine selection_sort(array)
|
|
! This version uses Fortran's minloc intrinsic
|
|
real(prec), intent(inout) :: array(:)
|
|
real(prec) :: temp
|
|
integer :: i,j
|
|
|
|
do i=1,size(array)-1
|
|
j=minloc(array(i:),1)+i-1
|
|
temp=array(i)
|
|
array(i)=array(j)
|
|
array(j)=temp
|
|
enddo
|
|
|
|
end subroutine selection_sort
|
|
|
|
subroutine inline_selection_sort(array)
|
|
! This version does not use Fortran's minloc intrinsic
|
|
real(prec), intent(inout) :: array(:)
|
|
real(prec) :: temp
|
|
integer :: i,j,k
|
|
|
|
do i=1,size(array)-1
|
|
temp=array(i)
|
|
j=i
|
|
do k=i+1,size(array)
|
|
if (array(k).lt.temp) then
|
|
temp=array(k)
|
|
j=k
|
|
endif
|
|
enddo
|
|
array(j)=array(i)
|
|
array(i)=temp
|
|
enddo
|
|
|
|
end subroutine inline_selection_sort
|
|
|
|
subroutine shell_sort(array)
|
|
! There are many choices of intervals
|
|
! The sequence (3**k-1)/2 used below is given by Knuth
|
|
! The sort relies on the insertion sort of nearly sorted
|
|
! data being fast (i.e. not order N**2)
|
|
real(prec), intent(inout) :: array(:)
|
|
integer :: i,interval,last
|
|
|
|
interval=1
|
|
last=size(array)
|
|
do while(interval<last/3)
|
|
interval=3*interval+1
|
|
enddo
|
|
|
|
do while(interval>0)
|
|
do i=1,interval
|
|
! The following syntax passes a subarray with elements
|
|
! array(i),array(i+interval),array(i+2*interval),...
|
|
call insertion_sort(array(i::interval))
|
|
enddo
|
|
interval=(interval-1)/3
|
|
enddo
|
|
|
|
end subroutine shell_sort
|
|
|
|
subroutine bubble_sort(array)
|
|
|
|
real(prec), intent(inout) :: array(:)
|
|
real(prec) :: temp
|
|
integer :: i,j,last
|
|
|
|
last=size(array)
|
|
do i=last-1,1,-1
|
|
do j=1,i
|
|
if (array(j+1).lt.array(j)) then
|
|
temp=array(j+1)
|
|
array(j+1)=array(j)
|
|
array(j)=temp
|
|
endif
|
|
enddo
|
|
enddo
|
|
|
|
end subroutine bubble_sort
|
|
|
|
subroutine heap_sort(array)
|
|
real(prec), intent(inout) :: array(:)
|
|
real(prec) :: temp
|
|
integer :: i,last
|
|
|
|
last=size(array)
|
|
|
|
! Build heap
|
|
do i=last/2,1,-1
|
|
call heapify(array,i)
|
|
enddo
|
|
|
|
! Unpick heap
|
|
do i=last,2,-1
|
|
temp=array(1)
|
|
array(1)=array(i)
|
|
array(i)=temp
|
|
call heapify(array(1:i-1),1)
|
|
end do
|
|
|
|
end subroutine heap_sort
|
|
|
|
recursive subroutine heapify(array,root)
|
|
real(prec), intent(inout) :: array(:)
|
|
real(prec) :: temp
|
|
integer :: left,right,root,last,largest
|
|
|
|
last=size(array)
|
|
left=2*root
|
|
right=left+1
|
|
largest=root
|
|
|
|
if (left.le.last) then
|
|
if(array(left).gt.array(largest)) largest=left
|
|
endif
|
|
|
|
if (right.le.last) then
|
|
if(array(right).gt.array(largest)) largest=right
|
|
endif
|
|
|
|
if (largest.ne.root) then
|
|
temp=array(root)
|
|
array(root)=array(largest)
|
|
array(largest)=temp
|
|
call heapify(array,largest)
|
|
end if
|
|
|
|
end subroutine heapify
|
|
|
|
|
|
subroutine heap_sort_nr(array)
|
|
! The same but with a non-recursive heapify
|
|
real(prec), intent(inout) :: array(:)
|
|
real(prec) :: temp
|
|
integer :: i,last
|
|
|
|
last=size(array)
|
|
|
|
! Build heap
|
|
do i=last/2,1,-1
|
|
call heapify_nr(array,i)
|
|
enddo
|
|
|
|
! Unpick heap
|
|
do i=last,2,-1
|
|
temp=array(1)
|
|
array(1)=array(i)
|
|
array(i)=temp
|
|
call heapify_nr(array(1:i-1),1)
|
|
end do
|
|
|
|
end subroutine heap_sort_nr
|
|
|
|
subroutine heapify_nr(array,i)
|
|
! A non-recursive heapify
|
|
real(prec), intent(inout) :: array(:)
|
|
real(prec) :: root_val
|
|
integer :: i,left,right,root,last,largest
|
|
|
|
last=size(array)
|
|
root=i
|
|
left=2*root
|
|
root_val=array(root)
|
|
largest=root
|
|
|
|
do while(left.le.last)
|
|
right=left+1
|
|
|
|
if (left.le.last) then
|
|
if(array(left).gt.array(largest)) largest=left
|
|
endif
|
|
|
|
if (right.le.last) then
|
|
if(array(right).gt.array(largest)) largest=right
|
|
endif
|
|
|
|
if (largest.eq.root) exit
|
|
|
|
array(root)=array(largest)
|
|
array(largest)=root_val
|
|
root=largest
|
|
left=2*root
|
|
enddo
|
|
|
|
! array(largest)=root_val
|
|
|
|
end subroutine heapify_nr
|
|
|
|
! Smoothsort is introduced by Edsger Dijkstra
|
|
! https://doi.org/10.1016/0167-6423(82)90016-8
|
|
!
|
|
! This Fortran version owes much to https://www.keithschwarz.com/smoothsort/
|
|
! and the C version given by Martin Knoblauch Revuelta at
|
|
! http://code.google.com/archive/p/combsortcs2p-and-other-sorting-algorithms
|
|
|
|
subroutine smoothsort(array)
|
|
real(prec), intent(inout) :: array(:)
|
|
integer(i64) :: mask
|
|
integer :: offset
|
|
|
|
if (size(array).le.1) return
|
|
|
|
call smooth_heapify(array,mask,offset)
|
|
|
|
call smooth_extract(array,mask,offset)
|
|
end subroutine smoothsort
|
|
|
|
subroutine smooth_heapify(array,mask,offset)
|
|
real(prec), intent(inout) :: array(0:)
|
|
integer(i64), intent(out) :: mask
|
|
integer, intent(out) :: offset
|
|
integer :: i,last
|
|
logical :: wbf
|
|
|
|
mask=1
|
|
offset=1
|
|
last=size(array)
|
|
|
|
do i=1,last-1
|
|
if (iand(mask,2_i64).eq.2) then
|
|
mask=ior(ishft(mask,-2),1_i64)
|
|
offset=offset+2
|
|
else if (offset.eq.1) then
|
|
mask=ior(ishft(mask,1),1_i64)
|
|
offset=0
|
|
else
|
|
mask=ior(ishft(mask,offset-1),1_i64)
|
|
offset=1
|
|
endif
|
|
wbf=(((iand(mask,2_i64).eq.2).and.(i+1.lt.last)).or. &
|
|
((offset.gt.0).and.(1+i+L(offset-1).lt.last)))
|
|
if (wbf) then
|
|
call sift_in(array,i,offset)
|
|
else
|
|
call interheap_sift(array,i,mask,offset)
|
|
endif
|
|
enddo
|
|
|
|
end subroutine smooth_heapify
|
|
subroutine smooth_extract(array,mask,offset)
|
|
real(prec), intent(inout) :: array(0:)
|
|
integer :: i,j,last,ch(2)
|
|
integer(i64),value :: mask
|
|
integer, value :: offset
|
|
|
|
last=size(array)
|
|
|
|
do i=last-1,2,-1
|
|
if (offset.lt.2) then
|
|
do
|
|
mask=ishft(mask,-1)
|
|
offset=offset+1
|
|
if (iand(mask,1_i64).eq.1) exit
|
|
enddo
|
|
else
|
|
ch(2)=i-1
|
|
ch(1)=ch(2)-L(offset-2)
|
|
mask=iand(mask,not(1_i64))
|
|
do j=1,2
|
|
mask=ior(ishft(mask,1),1_i64)
|
|
offset=offset-1
|
|
call interheap_sift(array,ch(j),mask,offset)
|
|
end do
|
|
end if
|
|
end do
|
|
|
|
|
|
end subroutine smooth_extract
|
|
|
|
|
|
subroutine sift_in(array,root,sz)
|
|
real(prec), intent(inout) :: array(0:)
|
|
real(prec) :: tmp
|
|
integer :: left,right,next
|
|
integer,value :: root,sz
|
|
|
|
if (sz.lt.2) return
|
|
|
|
tmp=array(root)
|
|
|
|
do
|
|
right=root-1
|
|
left=right-L(sz-2)
|
|
|
|
if (array(right).lt.array(left)) then
|
|
next=left
|
|
sz=sz-1
|
|
else
|
|
next=right
|
|
sz=sz-2
|
|
end if
|
|
|
|
if (array(next).le.tmp) exit
|
|
|
|
array(root)=array(next)
|
|
|
|
root=next
|
|
if (sz.le.1) exit
|
|
end do
|
|
|
|
array(root)=tmp
|
|
|
|
end subroutine sift_in
|
|
|
|
subroutine interheap_sift(array,root,mask,offset)
|
|
real(prec), intent(inout) :: array(0:)
|
|
real(prec) :: tmp,mx
|
|
integer :: left,right,next
|
|
integer(i64),value :: mask
|
|
integer, value :: root,offset
|
|
|
|
tmp=array(root)
|
|
|
|
do while (mask.ne.1)
|
|
mx=tmp
|
|
|
|
if (offset.gt.1) then
|
|
right=root-1
|
|
left=right-L(offset-2)
|
|
mx=max(mx,array(left),array(right))
|
|
end if
|
|
|
|
next=root-L(offset)
|
|
|
|
if (array(next).le.mx) exit
|
|
|
|
array(root)=array(next)
|
|
root=next
|
|
|
|
do
|
|
mask=ishft(mask,-1)
|
|
offset=offset+1
|
|
if (iand(mask,1_i64).eq.1) exit
|
|
end do
|
|
end do
|
|
|
|
|
|
array(root)=tmp
|
|
call sift_in(array,root,offset)
|
|
end subroutine interheap_sift
|
|
|
|
|
|
end module sorts
|