diff --git a/src/encounter/encounter_check.f90 b/src/encounter/encounter_check.f90 index 99466dc20..ea42073a0 100644 --- a/src/encounter/encounter_check.f90 +++ b/src/encounter/encounter_check.f90 @@ -282,10 +282,10 @@ subroutine encounter_check_reduce_broadphase(n, nenc, index1, index2, lencounter logical, dimension(:), allocatable, intent(inout) :: lencounter !! Logical flag indicating which of the pairs identified in the broad phase was selected in the narrow phase logical, dimension(:), allocatable, intent(inout) :: lvdotr !! Logical flag indicating the sign of v .dot. x ! Internals - integer(I4B) :: i, i0, j, k + integer(I4B) :: i, i0, j, k, nenci, klo, khi integer(I4B), dimension(:), allocatable :: itmp - logical, dimension(n) :: lfresh integer(I4B), dimension(:), allocatable :: ind + integer(I4B), dimension(n) :: ibeg, iend logical, dimension(:), allocatable :: ltmp nenc = count(lencounter(:)) ! Count the true number of encounters @@ -307,25 +307,44 @@ subroutine encounter_check_reduce_broadphase(n, nenc, index1, index2, lencounter lencounter(:) = .true. call util_sort(index1, ind) - lfresh(:) = .true. - - ! Remove duplicates - do k = 1, nenc - i = index1(ind(k)) - j = index2(ind(k)) - if (k == 1) then - lfresh(j) = .false. - else - i0 = index1(ind(k-1)) - if (i /= i0) lfresh(:) = .true. - if (.not.lfresh(j)) lencounter(ind(k)) = .false. - lfresh(j) = .false. + call util_sort_rearrange(index1, ind, nenc) + call util_sort_rearrange(index2, ind, nenc) + call util_sort_rearrange(lvdotr, ind, nenc) + + ! Get the bounds on the bodies in the first index + ibeg(:) = n + iend(:) = 1 + i0 = index1(1) + ibeg(i0) = 1 + do k = 2, nenc + i = index1(k) + if (i /= i0) then + iend(i0) = k - 1 + ibeg(i) = k + i0 = i end if + if (k == nenc) iend(i) = k + end do + + ! Sort on the second index and remove duplicates + if (allocated(itmp)) deallocate(itmp) + allocate(itmp, source=index2) + do concurrent(i = 1:n, iend(i) - ibeg(i) > 0) + klo = ibeg(i) + khi = iend(i) + nenci = khi - klo + 1 + if (allocated(ind)) deallocate(ind) + call util_sort(index2(klo:khi), ind) + index2(klo:khi) = itmp(klo - 1 + ind(:)) + do j = klo + 1, khi + if (index2(j) == index2(j - 1)) lencounter(j) = .false. + end do end do if (count(lencounter(:)) == nenc) return nenc = count(lencounter(:)) ! Count the true number of encounters + if (allocated(itmp)) deallocate(itmp) allocate(itmp(nenc)) itmp(:) = pack(index1(:), lencounter(:)) call move_alloc(itmp, index1) @@ -875,6 +894,7 @@ pure subroutine encounter_check_make_ragged_list(lencounteri, ind_arr, nenc,inde integer(I4B), dimension(:), intent(in) :: ind_arr integer(I4B), intent(out) :: nenc integer(I4B), dimension(:), allocatable, intent(out) :: index2 + ! Internals nenc = count(lencounteri(:)) if (nenc > 0) then diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index cc1b9f66d..7c40516b6 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -35,15 +35,12 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l nplm = pl%nplm nplt = npl - nplm - ! call timer%reset() - ! call timer%start() if (nplt == 0) then call encounter_check_all_plpl(param, npl, pl%xh, pl%vh, pl%renc, dt, lvdotr, index1, index2, nenc) else call encounter_check_all_plplm(param, nplm, nplt, pl%xh(:,1:nplm), pl%vh(:,1:nplm), pl%xh(:,nplm+1:npl), pl%vh(:,nplm+1:npl), pl%renc(1:nplm), pl%renc(nplm+1:npl), dt, lvdotr, index1, index2, nenc) end if - ! call timer%stop() - ! call timer%report(nsubsteps=nthreads, message="Encounter Check Total:") + lany_encounter = nenc > 0 if (lany_encounter) then call plplenc_list%resize(nenc) diff --git a/src/util/util_sort.f90 b/src/util/util_sort.f90 index 1b09b84f7..a3ad9a15c 100644 --- a/src/util/util_sort.f90 +++ b/src/util/util_sort.f90 @@ -60,27 +60,13 @@ end subroutine util_sort_body module pure subroutine util_sort_dp(arr) !! author: David A. Minton !! - !! Sort input DP precision array in place into ascending numerical order using insertion sort. - !! This algorithm works well for partially sorted arrays (which is usually the case here) + !! Sort input DP precision array in place into ascending numerical order using quicksort. !! implicit none ! Arguments real(DP), dimension(:), intent(inout) :: arr - ! Internals - real(DP) :: tmp - integer(I4B) :: n, i, j - n = size(arr) - do i = 2, n - tmp = arr(i) - j = i - 1 - do while (j >= 1) - if (arr(j) <= tmp) exit - arr(j + 1) = arr(j) - j = j - 1 - end do - arr(j + 1) = tmp - end do + call qsort_DP(arr) return end subroutine util_sort_dp @@ -89,17 +75,17 @@ end subroutine util_sort_dp module pure subroutine util_sort_index_dp(arr, ind) !! author: David A. Minton !! - !! Sort input DP precision array by index in ascending numerical order using insertion sort. + !! Sort input DP precision array by index in ascending numerical order using quick sort. !! This algorithm works well for partially sorted arrays (which is usually the case here). !! If ind is supplied already allocated, we assume it is an existing index array (e.g. a previously !! sorted array). If it is not allocated, this subroutine allocates it. !! implicit none ! Arguments - real(DP), dimension(:), intent(in) :: arr + real(DP), dimension(:), intent(in) :: arr integer(I4B), dimension(:), allocatable, intent(inout) :: ind ! Internals - integer(I4B) :: n, i, j, itmp + integer(I4B) :: n, i real(DP), dimension(:), allocatable :: tmparr n = size(arr) @@ -108,32 +94,38 @@ module pure subroutine util_sort_index_dp(arr, ind) ind = [(i, i=1, n)] end if allocate(tmparr, source=arr) - call qsort_index_DP(tmparr, ind) + call qsort_DP(tmparr, ind) return end subroutine util_sort_index_dp - recursive pure subroutine qsort_index_DP(arr, ind) + recursive pure subroutine qsort_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 + real(DP), dimension(:), intent(inout) :: arr + integer(I4B),dimension(:),intent(out), optional :: ind !! Internals integer :: iq if (size(arr) > 1) then - call partition_DP(arr, iq, ind) - call qsort_index_DP(arr(:iq-1),ind(:iq-1)) - call qsort_index_DP(arr(iq:), ind(iq:)) + if (present(ind)) then + call partition_DP(arr, iq, ind) + call qsort_DP(arr(:iq-1),ind(:iq-1)) + call qsort_DP(arr(iq:), ind(iq:)) + else + call partition_DP(arr, iq) + call qsort_DP(arr(:iq-1)) + call qsort_DP(arr(iq:)) + end if end if return - end subroutine qsort_index_DP + end subroutine qsort_DP pure subroutine partition_DP(arr, marker, ind) @@ -190,33 +182,20 @@ pure subroutine partition_DP(arr, marker, ind) end do return - end subroutine Partition_DP + end subroutine partition_DP module pure subroutine util_sort_i4b(arr) !! author: David A. Minton !! - !! Sort input integer array in place into ascending numerical order using insertion sort. + !! Sort input integer array in place into ascending numerical order using quick sort. !! This algorithm works well for partially sorted arrays (which is usually the case here) !! implicit none ! Arguments integer(I4B), dimension(:), intent(inout) :: arr - ! Internals - integer(I4B) :: tmp - integer(I4B) :: n, i, j - n = size(arr) - do i = 2, n - tmp = arr(i) - j = i - 1 - do while (j >= 1) - if (arr(j) <= tmp) exit - arr(j + 1) = arr(j) - j = j - 1 - end do - arr(j + 1) = tmp - end do + call qsort_I4B(arr) return end subroutine util_sort_i4b @@ -225,62 +204,125 @@ end subroutine util_sort_i4b module pure subroutine util_sort_index_i4b(arr, ind) !! author: David A. Minton !! - !! Sort input integer array by index in ascending numerical order using insertion sort. - !! This algorithm works well for partially sorted arrays (which is usually the case here). + !! Sort input integer array by index in ascending numerical order using quicksort. !! If ind is supplied already allocated, we assume it is an existing index array (e.g. a previously !! sorted array). If it is not allocated, this subroutine allocates it. !! implicit none ! Arguments - integer(I4B), dimension(:), intent(in) :: arr + integer(I4B), dimension(:), intent(in) :: arr integer(I4B), dimension(:), allocatable, intent(inout) :: ind ! Internals - integer(I4B) :: n, i, j, itmp + integer(I4B) :: n, i + integer(I4B), dimension(:), allocatable :: tmparr n = size(arr) if (.not.allocated(ind)) then allocate(ind(n)) ind = [(i, i=1, n)] end if - do i = 2, n - itmp = ind(i) - j = i - 1 - do while (j >= 1) - if (arr(ind(j)) <= arr(itmp)) exit - ind(j + 1) = ind(j) - j = j - 1 - end do - ind(j + 1) = itmp - end do + allocate(tmparr, source=arr) + call qsort_I4B(tmparr, ind) return end subroutine util_sort_index_i4b - module pure subroutine util_sort_sp(arr) + recursive pure subroutine qsort_I4B(arr, ind) !! author: David A. Minton !! - !! Sort input single precision array in place into ascending numerical order using insertion sort. - !! This algorithm works well for partially sorted arrays (which is usually the case here) - ! + !! Sort input DP precision array by index in ascending numerical order using quicksort. + !! implicit none ! Arguments - real(SP), dimension(:), intent(inout) :: arr + integer(I4B), dimension(:), intent(inout) :: arr + integer(I4B), dimension(:), intent(out), optional :: ind + !! Internals + integer :: iq + + if (size(arr) > 1) then + if (present(ind)) then + call partition_I4B(arr, iq, ind) + call qsort_I4B(arr(:iq-1),ind(:iq-1)) + call qsort_I4B(arr(iq:), ind(iq:)) + else + call partition_I4B(arr, iq) + call qsort_I4B(arr(:iq-1)) + call qsort_I4B(arr(iq:)) + end if + end if + + return + end subroutine qsort_I4B + + + pure subroutine partition_I4B(arr, marker, ind) + !! author: David A. Minton + !! + !! Partition function for quicksort on I4B type + !! + implicit none + ! Arguments + integer(I4B), intent(inout), dimension(:) :: arr + integer(I4B), intent(inout), dimension(:), optional :: ind + integer(I4B), intent(out) :: marker ! Internals - real(SP) :: tmp - integer(I4B) :: n, i, j + integer(I4B) :: i, j, itmp, narr, ipiv + integer(I4B) :: temp + integer(I4B) :: x ! pivot point - n = size(arr) - do i = 2, n - tmp = arr(i) - j = i - 1 - do while (j >= 1) - if (arr(j) <= tmp) exit - arr(j + 1) = arr(j) + narr = size(arr) + + ! Get center as pivot, as this is likely partially sorted + ipiv = narr / 2 + x = arr(ipiv) + i = 0 + j = narr + 1 + + do + j = j - 1 + do + if (arr(j) <= x) exit j = j - 1 end do - arr(j + 1) = tmp + i = i + 1 + do + if (arr(i) >= x) exit + i = i + 1 + end do + if (i < j) then + ! exchange A(i) and A(j) + temp = arr(i) + arr(i) = arr(j) + arr(j) = temp + 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 + else + marker = i + return + endif end do + + return + end subroutine partition_I4B + + + module pure subroutine util_sort_sp(arr) + !! author: David A. Minton + !! + !! Sort input DP precision array in place into ascending numerical order using quicksort. + !! + implicit none + ! Arguments + real(SP), dimension(:), intent(inout) :: arr + + call qsort_SP(arr) return end subroutine util_sort_sp @@ -289,36 +331,113 @@ end subroutine util_sort_sp module pure subroutine util_sort_index_sp(arr, ind) !! author: David A. Minton !! - !! Sort input single precision array by index in ascending numerical order using insertion sort. - !! This algorithm works well for partially sorted arrays (which is usually the case here). + !! Sort input DP precision array by index in ascending numerical order using quicksort. !! If ind is supplied already allocated, we assume it is an existing index array (e.g. a previously !! sorted array). If it is not allocated, this subroutine allocates it. !! implicit none ! Arguments - real(SP), dimension(:), intent(in) :: arr + real(SP), dimension(:), intent(in) :: arr integer(I4B), dimension(:), allocatable, intent(inout) :: ind ! Internals - integer(I4B) :: n, i, j, itmp + integer(I4B) :: n, i + real(SP), dimension(:), allocatable :: tmparr n = size(arr) if (.not.allocated(ind)) then allocate(ind(n)) ind = [(i, i=1, n)] end if - do i = 2, n - itmp = ind(i) - j = i - 1 - do while (j >= 1) - if (arr(ind(j)) <= arr(itmp)) exit - ind(j + 1) = ind(j) + allocate(tmparr, source=arr) + call qsort_SP(tmparr, ind) + + return + end subroutine util_sort_index_sp + + + recursive pure subroutine qsort_SP(arr, ind) + !! author: David A. Minton + !! + !! Sort input DP precision array by index in ascending numerical order using quicksort. + !! + implicit none + ! Arguments + real(SP), dimension(:), intent(inout) :: arr + integer(I4B),dimension(:),intent(out), optional :: ind + !! Internals + integer :: iq + + if (size(arr) > 1) then + if (present(ind)) then + call partition_SP(arr, iq, ind) + call qsort_SP(arr(:iq-1),ind(:iq-1)) + call qsort_SP(arr(iq:), ind(iq:)) + else + call partition_SP(arr, iq) + call qsort_SP(arr(:iq-1)) + call qsort_SP(arr(iq:)) + end if + end if + + return + end subroutine qsort_SP + + + pure subroutine partition_SP(arr, marker, ind) + !! author: David A. Minton + !! + !! Partition function for quicksort on SP type + !! + implicit none + ! Arguments + real(SP), intent(inout), dimension(:) :: arr + integer(I4B), intent(inout), dimension(:), optional :: ind + integer(I4B), intent(out) :: marker + ! Internals + integer(I4B) :: i, j, itmp, narr, ipiv + real(SP) :: temp + real(SP) :: x ! pivot point + + narr = size(arr) + + ! Get center as pivot, as this is likely partially sorted + ipiv = narr / 2 + x = arr(ipiv) + i = 0 + j = narr + 1 + + do + j = j - 1 + do + if (arr(j) <= x) exit j = j - 1 end do - ind(j + 1) = itmp + i = i + 1 + do + if (arr(i) >= x) exit + i = i + 1 + end do + if (i < j) then + ! exchange A(i) and A(j) + temp = arr(i) + arr(i) = arr(j) + arr(j) = temp + 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 + else + marker = i + return + endif end do - + return - end subroutine util_sort_index_sp + end subroutine partition_SP module subroutine util_sort_pl(self, sortby, ascending)