From 81ecfaf1a151158cdcb9d64151776634555ece29 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 27 Jul 2021 17:56:13 -0400 Subject: [PATCH] Replaced all the old sorting routines with new ones. util_index is now part of the generic util_sort --- .../9pl_18tp_encounters/pl.swiftest.in | 4 +- src/modules/swiftest_classes.f90 | 26 +++-- src/util/util_index.f90 | 103 ------------------ src/util/util_sort.f90 | 87 ++++++++++++++- 4 files changed, 105 insertions(+), 115 deletions(-) delete mode 100644 src/util/util_index.f90 diff --git a/examples/symba_swifter_comparison/9pl_18tp_encounters/pl.swiftest.in b/examples/symba_swifter_comparison/9pl_18tp_encounters/pl.swiftest.in index cd3c71538..8bc636061 100644 --- a/examples/symba_swifter_comparison/9pl_18tp_encounters/pl.swiftest.in +++ b/examples/symba_swifter_comparison/9pl_18tp_encounters/pl.swiftest.in @@ -1,5 +1,5 @@ 8 -1 4.9125474498983623693e-11 0.0014751237493860230134 +201 4.9125474498983623693e-11 0.0014751237493860230134 1.6306381826061645943e-05 -0.032433320146471017464 0.30732647407569840814 0.0280888997405028297 -0.033622812072399158034 -0.0019305604712619159943 0.0029264451427202888868 @@ -7,7 +7,7 @@ 4.0453784346544178454e-05 -0.6608991468450423623 -0.28805695486041710263 0.034183953683804932377 0.007943018642097033136 -0.018635382188272479886 -0.00071410720992500279457 -3 8.9970113821660187435e-10 0.010044863223462002622 +1003 8.9970113821660187435e-10 0.010044863223462002622 4.25875607065040958e-05 0.5665449483756358484 -0.84285201543201082597 3.8152874628327130158e-05 0.0139986033055793102076 0.009533392738922031109 -5.008237574040859916e-07 diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 313893e5d..b16df1020 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -20,7 +20,7 @@ module swiftest_classes public :: tides_kick_getacch_pl, tides_step_spin_system public :: user_kick_getacch_body public :: util_coord_b2h_pl, util_coord_b2h_tp, util_coord_h2b_pl, util_coord_h2b_tp, util_exit, util_fill_body, util_fill_pl, util_fill_tp, & - util_index, util_peri_tp, util_reverse_status, util_set_beg_end_pl, util_set_ir3h, util_set_msys, util_set_mu_pl, & + util_peri_tp, util_reverse_status, util_set_beg_end_pl, util_set_ir3h, util_set_msys, util_set_mu_pl, & util_set_mu_tp, util_set_rhill, util_set_rhill_approximate, util_sort, util_spill_body, util_spill_pl, util_spill_tp, util_valid, util_version !******************************************************************************************************************************** @@ -789,12 +789,6 @@ module subroutine util_fill_tp(self, inserts, lfill_list) logical, dimension(:), intent(in) :: lfill_list !! Logical array of bodies to merge into the keeps end subroutine util_fill_tp - module subroutine util_index(arr, index) - implicit none - integer(I4B), dimension(:), intent(out) :: index - real(DP), dimension(:), intent(in) :: arr - end subroutine util_index - module subroutine util_peri_tp(self, system, param) implicit none class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object @@ -856,15 +850,33 @@ module subroutine util_sort_i4b(arr) integer(I4B), dimension(:), intent(inout) :: arr end subroutine util_sort_i4b + module subroutine util_sort_index_i4b(arr,ind) + implicit none + integer(I4B), dimension(:), intent(in) :: arr + integer(I4B), dimension(:), intent(out) :: ind + end subroutine util_sort_index_i4b + module subroutine util_sort_sp(arr) implicit none real(SP), dimension(:), intent(inout) :: arr end subroutine util_sort_sp + module subroutine util_sort_index_sp(arr,ind) + implicit none + real(SP), dimension(:), intent(in) :: arr + integer(I4B), dimension(:), intent(out) :: ind + end subroutine util_sort_index_sp + module subroutine util_sort_dp(arr) implicit none real(DP), dimension(:), intent(inout) :: arr end subroutine util_sort_dp + + module subroutine util_sort_index_dp(arr,ind) + implicit none + real(DP), dimension(:), intent(in) :: arr + integer(I4B), dimension(:), intent(out) :: ind + end subroutine util_sort_index_dp end interface interface diff --git a/src/util/util_index.f90 b/src/util/util_index.f90 deleted file mode 100644 index fcece8809..000000000 --- a/src/util/util_index.f90 +++ /dev/null @@ -1,103 +0,0 @@ -submodule (swiftest_classes) s_util_index - use swiftest -contains - module subroutine util_index(arr, index) - !! author: David A. Minton - !! - !! Index input real array into ascending numerical order using Quicksort algorithm - !! - !! Adapted from David E. Kaufmann's Swifter routine: util_index.f90 - !! Adapted from Numerical Recipes in Fortran 90: The Art of Parallel Scientific Computing, by Press, Teukolsky, - !! Vetterling, and Flannery, 2nd ed., pp. 1173-4 - implicit none - ! Arguments - integer(I4B), dimension(:), intent(out) :: index - real(DP), dimension(:), intent(in) :: arr - ! Internals - integer(I4B), parameter :: nn = 15, nstack = 50 - integer(I4B) :: n, k, i, j, indext, jstack, l, r, dum - integer(I4B), dimension(nstack) :: istack - real(DP) :: a - - n = size(arr) - if (n /= size(index)) then - write(*, *) "Swiftest Error:" - write(*, *) " array size mismatch in util_index" - call util_exit(FAILURE) - end if - index = arth(1, 1, n) - jstack = 0 - ! l is the counter ie 'the one we are at' - l = 1 - ! r is the length of the array ie 'the total number of particles' - r = n - do - if ((r - l) < nn) then - do j = l + 1, r - indext = index(j) - a = arr(indext) - do i = j - 1, l, -1 - if (arr(index(i)) <= a) exit - index(i+1) = index(i) - end do - index(i+1) = indext - end do - if (jstack == 0) return - r = istack(jstack) - l = istack(jstack-1) - jstack = jstack - 2 - else - k = (l + r)/2 - dum = index(k); index(k) = index(l+1); index(l+1) = dum - ! if the mass of the particle we are at in our counting is greater than the mass of the last particle then put the particle we are at above the last one - if (arr(index(l)) > arr(index(r))) then - dum = index(l); index(l) = index(r); index(r) = dum - end if - ! if the mass of the particle above the one we are at in our counting is greater than the last particle then put that particle above the last one - if (arr(index(l+1)) > arr(index(r))) then - dum = index(l+1); index(l+1) = index(r); index(r) = dum - end if - ! if the mass of teh particle we are at in our counting is greater than the one above it, then put it above the one above it - if (arr(index(l)) > arr(index(l+1))) then - dum = index(l); index(l) = index(l+1); index(l+1) = dum - end if - i = l + 1 - j = r - indext = index(l+1) - a = arr(indext) - do - do - i = i + 1 - if (arr(index(i)) >= a) exit - end do - do - j = j - 1 - if (arr(index(j)) <= a) exit - end do - if (j < i) exit - dum = index(i); index(i) = index(j); index(j) = dum - end do - index(l+1) = index(j) - index(j) = indext - jstack = jstack + 2 - if (jstack > nstack) then - write(*, *) "Swiftest Error:" - write(*, *) " nstack too small in util_sort" - call util_exit(FAILURE) - end if - if ((r - i + 1) >= (j - l)) then - istack(jstack) = r - istack(jstack-1) = i - r = j - 1 - else - istack(jstack) = j - 1 - istack(jstack-1) = l - l = i - end if - end if - end do - - return - - end subroutine util_index -end submodule s_util_index diff --git a/src/util/util_sort.f90 b/src/util/util_sort.f90 index 108e7941d..73e1a4978 100644 --- a/src/util/util_sort.f90 +++ b/src/util/util_sort.f90 @@ -4,7 +4,7 @@ module subroutine util_sort_dp(arr) !! author: David A. Minton !! - !! Sort input double precision array into ascending numerical order using insertion sort. + !! Sort input double 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 @@ -26,10 +26,37 @@ module subroutine util_sort_dp(arr) return end subroutine util_sort_dp + module 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. + !! This algorithm works well for partially sorted arrays (which is usually the case here) + !! + implicit none + ! Arguments + real(DP), dimension(:), intent(in) :: arr + integer(I4B), dimension(:), intent(out) :: ind + ! Internals + real(DP) :: tmp + integer(I4B) :: n, i, j + + n = size(arr) + ind = [(i, i=1, n)] + do i = 2, n + tmp = arr(ind(i)) + do j = i - 1, 1, -1 + if (arr(ind(j)) <= tmp) exit + ind(j + 1) = ind(j) + end do + ind(j + 1) = i + end do + return + end subroutine util_sort_index_dp + module subroutine util_sort_i4b(arr) !! author: David A. Minton !! - !! Sort input integer array into ascending numerical order using insertion sort. + !! Sort input integer 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 @@ -51,10 +78,37 @@ module subroutine util_sort_i4b(arr) return end subroutine util_sort_i4b + module 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) + !! + implicit none + ! Arguments + integer(I4B), dimension(:), intent(in) :: arr + integer(I4B), dimension(:), intent(out) :: ind + ! Internals + integer(I4B) :: tmp + integer(I4B) :: n, i, j + + n = size(arr) + ind = [(i, i=1, n)] + do i = 2, n + tmp = arr(ind(i)) + do j = i - 1, 1, -1 + if (arr(ind(j)) <= tmp) exit + ind(j + 1) = ind(j) + end do + ind(j + 1) = i + end do + return + end subroutine util_sort_index_i4b + module subroutine util_sort_sp(arr) !! author: David A. Minton !! - !! Sort input single precision array into ascending numerical order using insertion sort. + !! 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) ! implicit none @@ -75,4 +129,31 @@ module subroutine util_sort_sp(arr) end do return end subroutine util_sort_sp + + module 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) + !! + implicit none + ! Arguments + real(SP), dimension(:), intent(in) :: arr + integer(I4B), dimension(:), intent(out) :: ind + ! Internals + real(SP) :: tmp + integer(I4B) :: n, i, j + + n = size(arr) + ind = [(i, i=1, n)] + do i = 2, n + tmp = arr(ind(i)) + do j = i - 1, 1, -1 + if (arr(ind(j)) <= tmp) exit + ind(j + 1) = ind(j) + end do + ind(j + 1) = i + end do + return + end subroutine util_sort_index_sp end submodule s_util_sort