aselimov
/
CACmb
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Replaced sorting method which resulted in broken deletion algorithm

Browse Source
master
Alex Selimov 5 years ago
parent 4f79085956
commit 9ccc5f1caf
3 changed files with 895 additions and 71 deletions
Show Stats Download Patch File Download Diff File

64
src/elements.f90
Unescape View File

@ -5,6 +5,7 @@ module elements
use parameters
use functions
use subroutines
use sorts
use box
implicit none
@ -395,54 +396,69 @@ module elements
return
end subroutine rhombshape
subroutine delete_atoms(num, index)
subroutine delete_atoms(num, in_index)
!This subroutine deletes atoms from the atom arrays
integer, intent(in) :: num
integer, intent(inout), dimension(num) :: index
integer, intent(in), dimension(num) :: in_index
real(kind=dp), dimension(num) :: for_sort
integer, dimension(num) :: sorted_index
integer :: i
call heapsort(index)
for_sort = in_index
call dpquicksort(for_sort)
do i = 1, num
sorted_index(i) = nint(for_sort(i))
end do
!We go from largest index to smallest index just to make sure that we don't miss
!accidentally overwrite values which need to be deleted
do i = num, 1, -1
if(index(i) == atom_num) then
r_atom(:,index(i)) = 0.0_dp
type_atom(index(i)) = 0
if(sorted_index(i) == atom_num) then
r_atom(:,sorted_index(i)) = 0.0_dp
type_atom(sorted_index(i)) = 0
else
r_atom(:,index(i)) = r_atom(:, atom_num)
type_atom(index(i)) = type_atom(atom_num)
r_atom(:,sorted_index(i)) = r_atom(:, atom_num)
type_atom(sorted_index(i)) = type_atom(atom_num)
end if
atom_num = atom_num - 1
end do
end subroutine
subroutine delete_elements(num, index)
subroutine delete_elements(num, in_index)
!This subroutine deletes elements from the element array
integer, intent(in) :: num
integer, intent(inout), dimension(num) :: index
integer, intent(in), dimension(num) :: in_index
integer :: i
real(kind=dp), dimension(num) :: for_sort
integer, dimension(num) :: sorted_index
call heapsort(index)
for_sort = in_index
call dpquicksort(for_sort)
do i = 1, num
sorted_index(i) = nint(for_sort(i))
end do
!We go from largest index to smallest index just to make sure that we don't miss
!accidentally overwrite values which need to be deleted
do i = num, 1, -1
if(index(i) == ele_num) then
r_node(:,:,:,index(i)) = 0.0_dp
type_ele(index(i)) =''
size_ele(index(i)) = 0
lat_ele(index(i)) = 0
sbox_ele(index(i)) = 0
if(sorted_index(i) == ele_num) then
r_node(:,:,:,sorted_index(i)) = 0.0_dp
type_ele(sorted_index(i)) =''
size_ele(sorted_index(i)) = 0
lat_ele(sorted_index(i)) = 0
sbox_ele(sorted_index(i)) = 0
else
node_num = node_num - ng_node(lat_ele(index(i)))
r_node(:,:,:,index(i)) = r_node(:,:,:,ele_num)
type_ele(index(i)) = type_ele(ele_num)
size_ele(index(i)) = size_ele(ele_num)
lat_ele(index(i)) = lat_ele(ele_num)
sbox_ele(index(i)) = sbox_ele(ele_num)
node_num = node_num - ng_node(lat_ele(sorted_index(i)))
r_node(:,:,:,sorted_index(i)) = r_node(:,:,:,ele_num)
type_ele(sorted_index(i)) = type_ele(ele_num)
size_ele(sorted_index(i)) = size_ele(ele_num)
lat_ele(sorted_index(i)) = lat_ele(ele_num)
sbox_ele(sorted_index(i)) = sbox_ele(ele_num)
end if
ele_num = ele_num - 1
end do
@ -632,4 +648,4 @@ module elements
end select
end subroutine
end module elements
end module elements

854
src/sorts.f90
Unescape View File

@ -0,0 +1,854 @@
! 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

48
src/subroutines.f90
Unescape View File

@ -173,52 +173,6 @@ module subroutines
end subroutine parse_ori_vec
subroutine heapsort(a)
integer, intent(in out) :: a(0:)
integer :: start, n, bottom
real :: temp
n = size(a)
do start = (n - 2) / 2, 0, -1
call siftdown(a, start, n);
end do
do bottom = n - 1, 1, -1
temp = a(0)
a(0) = a(bottom)
a(bottom) = temp;
call siftdown(a, 0, bottom)
end do
end subroutine heapsort
subroutine siftdown(a, start, bottom)
integer, intent(in out) :: a(0:)
integer, intent(in) :: start, bottom
integer :: child, root
real :: temp
root = start
do while(root*2 + 1 < bottom)
child = root * 2 + 1
if (child + 1 < bottom) then
if (a(child) < a(child+1)) child = child + 1
end if
if (a(root) < a(child)) then
temp = a(child)
a(child) = a (root)
a(root) = temp
root = child
else
return
end if
end do
end subroutine siftdown
subroutine apply_periodic(r)
@ -344,4 +298,4 @@ module subroutines
return
end subroutine check_right_ortho
end module subroutines
end module subroutines

Write Preview
Loading…
Cancel
Save