From ccc88d7f42c5c43f0737ec7586d29e3afc71e035 Mon Sep 17 00:00:00 2001 From: David Minton Date: Fri, 22 Oct 2021 16:39:36 -0400 Subject: [PATCH 1/4] Updated this branch with the changes made on the profiling branch --- src/encounter/encounter_check.f90 | 149 +++++++++++++++++------------- 1 file changed, 84 insertions(+), 65 deletions(-) diff --git a/src/encounter/encounter_check.f90 b/src/encounter/encounter_check.f90 index c64dbd68a..a77025a73 100644 --- a/src/encounter/encounter_check.f90 +++ b/src/encounter/encounter_check.f90 @@ -447,12 +447,12 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt ntot = nplm + nplt n = 2 * ntot if (ntot /= ntot_last) then - call boundingbox%setup(ntot, ntot_last) - ntot_last = ntot end if + call timer%reset() + call timer%start() !$omp parallel do default(private) schedule(static) & !$omp shared(xplm, xplt, vplm, vplt, rencm, renct, boundingbox) & !$omp firstprivate(dt, nplm, nplt, ntot) @@ -929,16 +929,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) @@ -966,16 +972,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) @@ -990,34 +1003,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 @@ -1029,33 +1042,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 @@ -1067,60 +1080,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 From b09010f9e4c23cc9cbd9247cbc798bbe3badcdf5 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 22 Oct 2021 16:43:30 -0400 Subject: [PATCH 2/4] Added quicksort to this branch --- src/util/util_sort.f90 | 85 +++++++++++++++++++++++++++++++++++------- 1 file changed, 72 insertions(+), 13 deletions(-) diff --git a/src/util/util_sort.f90 b/src/util/util_sort.f90 index 8ea7a3026..a662eb52c 100644 --- a/src/util/util_sort.f90 +++ b/src/util/util_sort.f90 @@ -60,7 +60,7 @@ end subroutine util_sort_body module pure subroutine util_sort_dp(arr) !! author: David A. Minton !! - !! Sort input double precision array in place into ascending numerical order using insertion sort. + !! 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) !! implicit none @@ -89,7 +89,7 @@ end subroutine util_sort_dp module pure subroutine util_sort_index_dp(arr, ind) !! author: David A. Minton !! - !! Sort input double precision array by index in ascending numerical order using insertion sort. + !! Sort input DP 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). !! 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. @@ -99,27 +99,86 @@ module pure subroutine util_sort_index_dp(arr, ind) real(DP), dimension(:), intent(in) :: arr integer(I4B), dimension(:), allocatable, intent(inout) :: ind ! Internals - integer(I4B) :: n, i, j, jp, itmp + integer(I4B) :: n, i, j, itmp + real(DP), 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_index_DP(tmparr, ind) + + return + end subroutine util_sort_index_dp + + + 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 + real(DP), dimension(:), intent(inout) :: arr + integer(I4B),dimension(:),intent(out) :: ind + 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:)) + end if + + return + end subroutine qsort_index_DP + + + pure subroutine partition_DP(arr, marker, ind) + !! author: David A. Minton + !! + !! Partition function for quicksort on DP type + real(DP), intent(inout), dimension(:) :: arr + integer(I4B), intent(inout), dimension(:), optional :: ind + integer(I4B), intent(out) :: marker + integer(I4B) :: i, j, itmp + real(DP) :: temp + real(DP) :: x ! pivot point + + x = arr(1) + i = 0 + j = size(arr) + 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 + itmp = ind(i) + ind(i) = ind(j) + ind(j) = itmp + else if (i == j) then + marker = i + 1 + return + else + marker = i + return + endif end do - + return - end subroutine util_sort_index_dp - + end subroutine Partition_DP + module pure subroutine util_sort_i4b(arr) !! author: David A. Minton From 396da9a71dd72d7cd6810aef4be78b1467d480fc Mon Sep 17 00:00:00 2001 From: David Minton Date: Fri, 22 Oct 2021 22:32:17 -0400 Subject: [PATCH 3/4] Added in improvements from the profiling branch --- src/encounter/encounter_check.f90 | 62 +++++++++++++++++-------------- src/modules/encounter_classes.f90 | 2 +- 2 files changed, 36 insertions(+), 28 deletions(-) diff --git a/src/encounter/encounter_check.f90 b/src/encounter/encounter_check.f90 index a77025a73..e8115be32 100644 --- a/src/encounter/encounter_check.f90 +++ b/src/encounter/encounter_check.f90 @@ -1020,10 +1020,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 @@ -1058,10 +1057,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 @@ -1080,7 +1079,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) @@ -1095,20 +1094,23 @@ 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 + 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) + 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 - lenc(dim) = (iend(dim,j) > ibegi(dim)) .and. (ibeg(dim,j) < iendi(dim)) - end do - lencounteri(j) = all(lenc(:)) + ! Check the other dimensions + lencounteri(j) = all((iend(2:SWEEPDIM,j) > ibegi(2:SWEEPDIM)) .and. (ibeg(2:SWEEPDIM,j) < iendi(2:SWEEPDIM))) end do return @@ -1127,19 +1129,25 @@ 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 + 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)) - end do - lencounteri(j) = all(lenc(:)) + lencounteri(j) = all((iend(2:SWEEPDIM,j) > ibegi(2:SWEEPDIM)) .and. (ibeg(2:SWEEPDIM,j) < iendi(2:SWEEPDIM))) end do return diff --git a/src/modules/encounter_classes.f90 b/src/modules/encounter_classes.f90 index 4965f14f1..1c182619b 100644 --- a/src/modules/encounter_classes.f90 +++ b/src/modules/encounter_classes.f90 @@ -7,7 +7,7 @@ module encounter_classes implicit none public - integer(I4B), parameter :: SWEEPDIM = 2 + integer(I4B), parameter :: SWEEPDIM = 3 type :: encounter_list integer(I4B) :: nenc = 0 !! Total number of encounters From a802bccfa3b798bfd83f55d38f1a0ef7d853bde8 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 22 Oct 2021 22:35:01 -0400 Subject: [PATCH 4/4] Added pure back to the sweep subroutine --- src/encounter/encounter_check.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/encounter/encounter_check.f90 b/src/encounter/encounter_check.f90 index e8115be32..298bf7d4f 100644 --- a/src/encounter/encounter_check.f90 +++ b/src/encounter/encounter_check.f90 @@ -1079,7 +1079,7 @@ subroutine encounter_check_sweep_aabb_all_single_list(n, ext_ind, ibeg, iend, in end subroutine encounter_check_sweep_aabb_all_single_list - subroutine encounter_check_sweep_aabb_one_double_list(i, n1, n2, ntot, ext_ind, ibegi, iendi, ibeg, iend, 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)