Improved performance of sweep phase

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))

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