Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Swapped insertion sort for quicksort for better performance in the so…
Browse files Browse the repository at this point in the history 
…rt/sweep algorithm
  • Loading branch information
@daminton
daminton committed Oct 22, 2021
1 parent a42ac97 commit ed2bdaf
Show file tree
Hide file tree
Showing 2 changed files with 155 additions and 85 deletions.
155 changes: 83 additions & 72 deletions src/encounter/encounter_check.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -497,7 +497,7 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt
!$omp end parallel do
call timer%stop()
write(*,*) "plplm sort : ",timer%count_stop_step - timer%count_start_step

call timer%reset()
call timer%start()
call boundingbox%sweep(nplm, nplt, nenc, index1, index2)
Expand Down Expand Up @@ -702,11 +702,7 @@ subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, lvdotr, inde
call timer%stop()
write(*,*) "plpl triang: ",timer%count_stop_step - timer%count_start_step

call timer%reset()
call timer%start()
call encounter_check_collapse_ragged_list(lenc, npl, nenc, index1, index2, lvdotr)
call timer%stop()
write(*,*) "plpl tricol: ",timer%count_stop_step - timer%count_start_step

return
end subroutine encounter_check_all_triangular_plpl
Expand Down Expand Up @@ -758,11 +754,7 @@ subroutine encounter_check_all_triangular_plplm(nplm, nplt, xplm, vplm, xplt, vp
call timer%stop()
write(*,*) "plplm triang: ",timer%count_stop_step - timer%count_start_step

call timer%reset()
call timer%start()
call encounter_check_collapse_ragged_list(lenc, nplm, nenc, index1, index2, lvdotr)
call timer%stop()
write(*,*) "plplm tricol: ",timer%count_stop_step - timer%count_start_step

return
end subroutine encounter_check_all_triangular_plplm
Expand Down Expand Up @@ -975,16 +967,22 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, nenc, ind
integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter candidate pair
integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter candidate pair
!Internals
integer(I4B) :: i, k, ntot
integer(I4B) :: i, k, ntot, dim
type(encounter_list), dimension(n1+n2) :: lenc !! Array of encounter lists (one encounter list per body)
integer(I4B), dimension(:), allocatable, save :: ind_arr
integer(I4B), dimension(:), allocatable :: ibeg, iend

ntot = n1 + n2
call util_index_array(ind_arr, ntot)
allocate(ibeg(SWEEPDIM * ntot))
allocate(iend(SWEEPDIM * ntot))
do dim = 1, SWEEPDIM
ibeg((dim - 1) * ntot + 1:dim * ntot) = self%aabb(dim)%ibeg(:)
iend((dim - 1) * ntot + 1:dim * ntot) = self%aabb(dim)%iend(:)
end do
! Sweep the intervals for each of the massive bodies along one dimension
! This will build a ragged pair of index lists inside of the lenc data structure

call encounter_check_sweep_aabb_all_double_list(n1, n2, self%aabb(1)%ind(:), self%aabb(1)%ibeg(:), self%aabb(1)%iend(:), self%aabb(2)%ibeg(:), self%aabb(2)%iend(:), ind_arr(:), lenc(:))
call encounter_check_sweep_aabb_all_double_list(n1, n2, self%aabb(1)%ind(:), reshape(ibeg(:), [SWEEPDIM, ntot]), reshape(iend(:), [SWEEPDIM, ntot]), ind_arr(:), lenc(:))

call encounter_check_collapse_ragged_list(lenc, ntot, nenc, index1, index2)

Expand Down Expand Up @@ -1012,16 +1010,23 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, nenc, index1,
integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter candidate pair
integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter candidate pair
!Internals
Integer(I4B) :: i, k
Integer(I4B) :: i, k, dim
type(encounter_list), dimension(n) :: lenc !! Array of encounter lists (one encounter list per body)
integer(I4B), dimension(:), allocatable, save :: ind_arr
type(walltimer) :: timer
integer(I4B), dimension(:), allocatable :: ibeg, iend

call util_index_array(ind_arr, n)
allocate(ibeg(SWEEPDIM * n))
allocate(iend(SWEEPDIM * n))
do dim = 1, SWEEPDIM
ibeg((dim - 1) * n + 1:dim * n) = self%aabb(dim)%ibeg(:)
iend((dim - 1) * n + 1:dim * n) = self%aabb(dim)%iend(:)
end do

! Sweep the intervals for each of the massive bodies along one dimension
! This will build a ragged pair of index lists inside of the lenc data structure
call encounter_check_sweep_aabb_all_single_list(n, self%aabb(1)%ind(:), self%aabb(1)%ibeg(:), self%aabb(1)%iend(:), self%aabb(2)%ibeg(:), self%aabb(2)%iend(:), ind_arr(:), lenc(:))
call encounter_check_sweep_aabb_all_single_list(n, self%aabb(1)%ind(:), reshape(ibeg(:), [SWEEPDIM, n]), reshape(iend(:), [SWEEPDIM, n]), ind_arr(:), lenc(:))

call encounter_check_collapse_ragged_list(lenc, n, nenc, index1, index2)

Expand All @@ -1036,34 +1041,34 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, nenc, index1,
end subroutine encounter_check_sweep_aabb_single_list


subroutine encounter_check_sweep_aabb_all_double_list(n1, n2, ext_ind, ibegx, iendx, ibegy, iendy, ind_arr, lenc)
subroutine encounter_check_sweep_aabb_all_double_list(n1, n2, ext_ind, ibeg, iend, ind_arr, lenc)
!! author: David A. Minton
!!
!! Performs the loop part of the sweep operation. Double list version (e.g. pl-tp or plm-plt)
implicit none
! Arguments
integer(I4B), intent(in) :: n1, n2 !! Number of bodies
integer(I4B), dimension(:), intent(in) :: ext_ind !! Sorted index array of extents
integer(I4B), dimension(:), intent(in) :: ibegx, iendx !! Beginning and ending index lists in the x-dimension
integer(I4B), dimension(:), intent(in) :: ibegy, iendy !! Beginning and ending index lists in the y-dimension
integer(I4B), dimension(:), intent(in) :: ind_arr !! index array for mapping body 2 indexes
type(encounter_list), dimension(:), intent(inout) :: lenc !! Encounter list for the ith body
integer(I4B), intent(in) :: n1, n2 !! Number of bodies
integer(I4B), dimension(:), intent(in) :: ext_ind !! Sorted index array of extents
integer(I4B), dimension(:,:), intent(in) :: ibeg, iend !! Beginning and ending index lists in the n-dimensions
integer(I4B), dimension(:), intent(in) :: ind_arr !! index array for mapping body 2 indexes
type(encounter_list), dimension(:), intent(inout) :: lenc !! Encounter list for the ith body
! Internals
integer(I4B) :: i, ntot
logical, dimension(n1+n2) :: lencounteri
integer(I4B) :: ibegxi, iendxi, ibegyi, iendyi
integer(I4B), dimension(SWEEPDIM) :: ibegi, iendi

ntot = n1 + n2
!$omp parallel do default(private) schedule(guided)&
!$omp shared(ext_ind, ibegx, iendx, ibegy, iendy, ind_arr, lenc) &
!$omp firstprivate(ntot, n1, n2)
!$omp shared(ext_ind, ibeg, iend, ind_arr, lenc) &
!$omp firstprivate(ntot, n1, n2) &
!$omp lastprivate(ibegi, iendi)
do i = 1, ntot
ibegxi = ibegx(i) + 1
iendxi = iendx(i) - 1
if (iendxi >= ibegxi) then
ibegyi = ibegy(i)
iendyi = iendy(i)
call encounter_check_sweep_aabb_one_double_list(i, n1, n2, ntot, ext_ind(:), ibegxi, iendxi, ibegyi, iendyi, ibegx(:), iendx(:), ibegy(:), iendy(:), lencounteri(:))
ibegi(1) = ibeg(1,i) + 1
iendi(1) = iend(1,i) - 1
if (iendi(1) >= ibegi(1)) then
ibegi(2:SWEEPDIM) = ibeg(2:SWEEPDIM,i)
iendi(2:SWEEPDIM) = iend(2:SWEEPDIM,i)
call encounter_check_sweep_aabb_one_double_list(i, n1, n2, ntot, ext_ind(:), ibegi(:), iendi(:), ibeg(:,:), iend(:,:), lencounteri(:))
call encounter_check_make_ragged_list(lencounteri(:), ind_arr(:), lenc(i)%nenc, lenc(i)%index2)
else
lenc(i)%nenc = 0
Expand All @@ -1075,33 +1080,33 @@ subroutine encounter_check_sweep_aabb_all_double_list(n1, n2, ext_ind, ibegx, ie
end subroutine encounter_check_sweep_aabb_all_double_list


subroutine encounter_check_sweep_aabb_all_single_list(n, ext_ind, ibegx, iendx, ibegy, iendy, ind_arr, lenc)
subroutine encounter_check_sweep_aabb_all_single_list(n, ext_ind, ibeg, iend, ind_arr, lenc)
!! author: David A. Minton
!!
!! Performs the loop part of the sweep operation. Single list version (e.g. pl-pl)
implicit none
! Arguments
integer(I4B), intent(in) :: n !! Number of bodies
integer(I4B), dimension(:), intent(in) :: ext_ind !! Sorted index array of extents
integer(I4B), dimension(:), intent(in) :: ibegx, iendx !! Beginning and ending index lists in the x-dimension
integer(I4B), dimension(:), intent(in) :: ibegy, iendy !! Beginning and ending index lists in the y-dimension
integer(I4B), dimension(:), intent(in) :: ind_arr !! index array for mapping body 2 indexes
type(encounter_list), dimension(:), intent(inout) :: lenc !! Encounter list for the ith body
integer(I4B), intent(in) :: n !! Number of bodies
integer(I4B), dimension(:), intent(in) :: ext_ind !! Sorted index array of extents
integer(I4B), dimension(:,:), intent(in) :: ibeg, iend !! Beginning and ending index lists in the n-dimension
integer(I4B), dimension(:), intent(in) :: ind_arr !! index array for mapping body 2 indexes
type(encounter_list), dimension(:), intent(inout) :: lenc !! Encounter list for the ith body
! Internals
integer(I4B) :: i
logical, dimension(n) :: lencounteri
integer(I4B) :: ibegxi, iendxi, ibegyi, iendyi
integer(I4B), dimension(SWEEPDIM) :: ibegi, iendi

!$omp parallel do default(private) schedule(guided)&
!$omp shared(ext_ind, ibegx, iendx, ibegy, iendy, ind_arr, lenc) &
!$omp firstprivate(n)
!$omp shared(ext_ind, ibeg, iend, ind_arr, lenc) &
!$omp firstprivate(n) &
!$omp lastprivate(ibegi, iendi)
do i = 1, n
ibegxi = ibegx(i) + 1
iendxi = iendx(i) - 1
if (iendxi >= ibegxi) then
ibegyi = ibegy(i)
iendyi = iendy(i)
call encounter_check_sweep_aabb_one_single_list(n, ext_ind(:), ibegxi, iendxi, ibegyi, iendyi, ibegx(:), iendx(:), ibegy(:), iendy(:), lencounteri(:))
ibegi(1) = ibeg(1,i) + 1
iendi(1) = iend(1,i) - 1
if (iendi(1) >= ibegi(1)) then
ibegi(2:SWEEPDIM) = ibeg(2:SWEEPDIM,i)
iendi(2:SWEEPDIM) = iend(2:SWEEPDIM,i)
call encounter_check_sweep_aabb_one_single_list(n, ext_ind(:), ibegi(:), iendi(:), ibeg(:,:), iend(:,:), lencounteri(:))
call encounter_check_make_ragged_list(lencounteri(:), ind_arr(:), lenc(i)%nenc, lenc(i)%index2)
else
lenc(i)%nenc = 0
Expand All @@ -1113,60 +1118,66 @@ subroutine encounter_check_sweep_aabb_all_single_list(n, ext_ind, ibegx, iendx,
end subroutine encounter_check_sweep_aabb_all_single_list


pure subroutine encounter_check_sweep_aabb_one_double_list(i, n1, n2, ntot, ext_ind, ibegxi, iendxi, ibegyi, iendyi, ibegx, iendx, ibegy, iendy, lencounteri)
pure subroutine encounter_check_sweep_aabb_one_double_list(i, n1, n2, ntot, ext_ind, ibegi, iendi, ibeg, iend, lencounteri)
!! author: David A. Minton
!!
!! Performs a sweep operation on a single body. Encounters from the same lists not allowed (e.g. pl-tp encounters only)
implicit none
! Arguments
integer(I4B), intent(in) :: i !! The current index of the ith body
integer(I4B), intent(in) :: n1 !! Number of bodies 1
integer(I4B), intent(in) :: n2 !! Number of bodies 2
integer(I4B), intent(in) :: ntot !! n1 + n2
integer(I4B), dimension(:), intent(in) :: ext_ind !! Sorted index array of extents
integer(I4B), intent(in) :: ibegxi, iendxi !! The beginning and ending indices of the ith bounding box in the x-dimension
integer(I4B), intent(in) :: ibegyi, iendyi !! The beginning and ending indices of the ith bounding box in the y-dimension
integer(I4B), dimension(:), intent(in) :: ibegx, iendx !! Beginning and ending index lists in the x-dimension
integer(I4B), dimension(:), intent(in) :: ibegy, iendy !! Beginning and ending index lists in the y-dimensio
logical, dimension(:), intent(out) :: lencounteri !! Encounter list for the ith body
integer(I4B), intent(in) :: i !! The current index of the ith body
integer(I4B), intent(in) :: n1 !! Number of bodies 1
integer(I4B), intent(in) :: n2 !! Number of bodies 2
integer(I4B), intent(in) :: ntot !! n1 + n2
integer(I4B), dimension(:), intent(in) :: ext_ind !! Sorted index array of extents
integer(I4B), dimension(:), intent(in) :: ibegi, iendi !! The beginning and ending indices of the ith bounding box in the n-dimensions
integer(I4B), dimension(:,:), intent(in) :: ibeg, iend !! Beginning and ending index lists in the n-dimensions
logical, dimension(:), intent(out) :: lencounteri !! Encounter list for the ith body
! Internals
integer(I4B) :: j, jbox
integer(I4B) :: j, jbox, dim
logical, dimension(SWEEPDIM) :: lenc

lencounteri(:) = .false.
do concurrent(jbox = ibegxi:iendxi) ! Sweep forward until the end of the interval
lenc(1) = .true.
do concurrent(jbox = ibegi(1):iendi(1)) ! Sweep forward until the end of the interval
j = ext_ind(jbox)
if (j > ntot) j = j - ntot ! If this is an endpoint index, shift it to the correct range
if (((i <= n1) .and. (j <= n1)) .or. ((i > n1) .and. (j > n1))) cycle ! only pairs from the two different lists allowed
! Check the y-dimension
lencounteri(j) = (iendy(j) > ibegyi) .and. (ibegy(j) < iendyi)
do dim = 2, SWEEPDIM
lenc(dim) = (iend(dim,j) > ibegi(dim)) .and. (ibeg(dim,j) < iendi(dim))
end do
lencounteri(j) = all(lenc(:))
end do

return
end subroutine encounter_check_sweep_aabb_one_double_list


pure subroutine encounter_check_sweep_aabb_one_single_list(n, ext_ind, ibegxi, iendxi, ibegyi, iendyi, ibegx, iendx, ibegy, iendy, lencounteri)
pure subroutine encounter_check_sweep_aabb_one_single_list(n, ext_ind, ibegi, iendi, ibeg, iend, lencounteri)
!! author: David A. Minton
!!
!! Performs a sweep operation on a single body. Mutual encounters allowed (e.g. pl-pl)
implicit none
! Arguments
integer(I4B), intent(in) :: n !! Number of bodies
integer(I4B), dimension(:), intent(in) :: ext_ind !! Sorted index array of extents
integer(I4B), intent(in) :: ibegxi, iendxi !! The beginning and ending indices of the ith bounding box in the x-dimension
integer(I4B), intent(in) :: ibegyi, iendyi !! The beginning and ending indices of the ith bounding box in the y-dimension
integer(I4B), dimension(:), intent(in) :: ibegx, iendx !! Beginning and ending index lists in the x-dimension
integer(I4B), dimension(:), intent(in) :: ibegy, iendy !! Beginning and ending index lists in the y-dimension
logical, dimension(:), intent(out) :: lencounteri !! Encounter list for the ith body
integer(I4B), intent(in) :: n !! Number of bodies
integer(I4B), dimension(:), intent(in) :: ext_ind !! Sorted index array of extents
integer(I4B), dimension(:), intent(in) :: ibegi, iendi !! The beginning and ending indices of the ith bounding box in the n-dimensions
integer(I4B), dimension(:,:), intent(in) :: ibeg, iend !! Beginning and ending index lists in the x-dimension
logical, dimension(:), intent(out) :: lencounteri !! Encounter list for the ith body
! Internals
integer(I4B) :: j, jbox
integer(I4B) :: j, jbox, dim
logical, dimension(SWEEPDIM) :: lenc

lencounteri(:) = .false.
do concurrent(jbox = ibegxi:iendxi) ! Sweep forward until the end of the interval
lenc(1) = .true.
do concurrent(jbox = ibegi(1):iendi(1)) ! Sweep forward until the end of the interval
j = ext_ind(jbox)
if (j > n) j = j - n ! If this is an endpoint index, shift it to the correct range
! Check the y-dimension
lencounteri(j) = (iendy(j) > ibegyi) .and. (ibegy(j) < iendyi)
! Check the other dimensions
do dim = 2, SWEEPDIM
lenc(dim) = (iend(dim,j) > ibegi(dim)) .and. (ibeg(dim,j) < iendi(dim))
end do
lencounteri(j) = all(lenc(:))
end do

return
Expand Down
Loading

0 comments on commit ed2bdaf

Please sign in to comment.