diff --git a/src/encounter/encounter_check.f90 b/src/encounter/encounter_check.f90 index 93580eb33..9c5570508 100644 --- a/src/encounter/encounter_check.f90 +++ b/src/encounter/encounter_check.f90 @@ -971,6 +971,7 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, nenc, ind 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 + type(walltimer) :: timer ntot = n1 + n2 call util_index_array(ind_arr, ntot) @@ -982,16 +983,28 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, nenc, ind 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 timer%reset() + call timer%start() 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 timer%stop() + write(*,*) "sweep double: ",timer%count_stop_step - timer%count_start_step + + call timer%reset() + call timer%start() call encounter_check_collapse_ragged_list(lenc, ntot, nenc, index1, index2) + call timer%stop() + write(*,*) "collapse rag: ",timer%count_stop_step - timer%count_start_step + call timer%reset() + call timer%start() ! Reorder the pairs and sort the first index in order to remove any duplicates do concurrent(k = 1:nenc, index2(k) < index1(k)) i = index2(k) index2(k) = index1(k) index1(k) = i end do + call timer%stop() + write(*,*) "reorder arr : ",timer%count_stop_step - timer%count_start_step return end subroutine encounter_check_sweep_aabb_double_list @@ -1058,10 +1071,9 @@ subroutine encounter_check_sweep_aabb_all_double_list(n1, n2, ext_ind, ibeg, ien integer(I4B), dimension(SWEEPDIM) :: ibegi, iendi ntot = n1 + n2 - !$omp parallel do default(private) schedule(guided)& + !$omp parallel do default(private) schedule(dynamic)& !$omp shared(ext_ind, ibeg, iend, ind_arr, lenc) & - !$omp firstprivate(ntot, n1, n2) & - !$omp lastprivate(ibegi, iendi) + !$omp firstprivate(ntot, n1, n2) do i = 1, ntot ibegi(1) = ibeg(1,i) + 1 iendi(1) = iend(1,i) - 1 @@ -1096,10 +1108,10 @@ subroutine encounter_check_sweep_aabb_all_single_list(n, ext_ind, ibeg, iend, in logical, dimension(n) :: lencounteri integer(I4B), dimension(SWEEPDIM) :: ibegi, iendi - !$omp parallel do default(private) schedule(guided)& + !$omp parallel do default(private) schedule(dynamic)& !$omp shared(ext_ind, ibeg, iend, ind_arr, lenc) & !$omp firstprivate(n) & - !$omp lastprivate(ibegi, iendi) + !$omp lastprivate(ibegi, iendi, lencounteri) do i = 1, n ibegi(1) = ibeg(1,i) + 1 iendi(1) = iend(1,i) - 1 @@ -1118,7 +1130,7 @@ subroutine encounter_check_sweep_aabb_all_single_list(n, ext_ind, ibeg, iend, in end subroutine encounter_check_sweep_aabb_all_single_list - pure subroutine encounter_check_sweep_aabb_one_double_list(i, n1, n2, ntot, ext_ind, ibegi, iendi, ibeg, iend, lencounteri) + 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) @@ -1133,14 +1145,20 @@ pure subroutine encounter_check_sweep_aabb_one_double_list(i, n1, n2, ntot, ext_ 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, dim - logical, dimension(SWEEPDIM) :: lenc + integer(I4B) :: j, jbox, dim, jlo, jhi + logical, dimension(2:SWEEPDIM) :: lenc + integer(I4B), dimension(:), allocatable :: box lencounteri(:) = .false. - 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 + jlo = ibegi(1) + jhi = iendi(1) + allocate(box(jlo:jhi)) + box(:) = ext_ind(jlo:jhi) + where(box(:) > ntot) + box(:) = box(:) - ntot + endwhere + do concurrent(jbox = jlo:jhi) ! Sweep forward until the end of the interval + j = box(jbox) 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 do dim = 2, SWEEPDIM @@ -1165,14 +1183,24 @@ pure subroutine encounter_check_sweep_aabb_one_single_list(n, ext_ind, ibegi, ie 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, dim - logical, dimension(SWEEPDIM) :: lenc + integer(I4B) :: j, jbox, dim, jlo, jhi + logical, dimension(2:SWEEPDIM) :: lenc + integer(I4B), dimension(:), allocatable :: box lencounteri(:) = .false. - 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 + + jlo = ibegi(1) + jhi = iendi(1) + + allocate(box(jlo:jhi)) + box(:) = ext_ind(jlo:jhi) + + where(box(:) > n) + box(:) = box(:) - n + endwhere + + do concurrent(jbox = jlo:jhi) ! Sweep forward until the end of the interval + j = box(jbox) ! Check the other dimensions do dim = 2, SWEEPDIM lenc(dim) = (iend(dim,j) > ibegi(dim)) .and. (ibeg(dim,j) < iendi(dim)) diff --git a/src/util/util_sort.f90 b/src/util/util_sort.f90 index a662eb52c..1b09b84f7 100644 --- a/src/util/util_sort.f90 +++ b/src/util/util_sort.f90 @@ -118,9 +118,12 @@ recursive pure subroutine qsort_index_DP(arr, ind) !! author: David A. Minton !! !! Sort input DP precision array by index in ascending numerical order using quicksort sort. + !! implicit none + ! Arguments real(DP), dimension(:), intent(inout) :: arr integer(I4B),dimension(:),intent(out) :: ind + !! Internals integer :: iq if (size(arr) > 1) then @@ -137,16 +140,24 @@ pure subroutine partition_DP(arr, marker, ind) !! author: David A. Minton !! !! Partition function for quicksort on DP type - real(DP), intent(inout), dimension(:) :: arr + !! + implicit none + ! Arguments + real(DP), intent(inout), dimension(:) :: arr integer(I4B), intent(inout), dimension(:), optional :: ind - integer(I4B), intent(out) :: marker - integer(I4B) :: i, j, itmp + integer(I4B), intent(out) :: marker + ! Internals + integer(I4B) :: i, j, itmp, narr, ipiv real(DP) :: temp - real(DP) :: x ! pivot point + real(DP) :: x ! pivot point + + narr = size(arr) - x = arr(1) + ! Get center as pivot, as this is likely partially sorted + ipiv = narr / 2 + x = arr(ipiv) i = 0 - j = size(arr) + 1 + j = narr + 1 do j = j - 1 @@ -164,9 +175,11 @@ pure subroutine partition_DP(arr, marker, ind) temp = arr(i) arr(i) = arr(j) arr(j) = temp - itmp = ind(i) - ind(i) = ind(j) - ind(j) = itmp + if (present(ind)) then + itmp = ind(i) + ind(i) = ind(j) + ind(j) = itmp + end if else if (i == j) then marker = i + 1 return