Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Replaced NR quicksort with simple insertion sort, as it performs bett…
Browse files Browse the repository at this point in the history 
…er for partially sorted lists, and nearly every time we need to sort in SyMBA it's a nearly sorted list
  • Loading branch information
@daminton
daminton committed Jul 27, 2021
1 parent f38d19e commit 49b18b0
Show file tree
Hide file tree
Showing 4 changed files with 50 additions and 236 deletions.
4 changes: 2 additions & 2 deletions examples/symba_swifter_comparison/9pl_18tp_encounters/swiftest_symba_vs_swifter_symba.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -622,7 +622,7 @@
" 1.25413096e+00])\n",
"Coordinates:\n",
" id int64 2\n",
" * time (d) (time (d)) float64 0.0 11.0 22.0 ... 3.63e+03 3.641e+03 3.652e+03</pre><div class='xr-wrap' hidden><div class='xr-header'><div class='xr-obj-type'>xarray.DataArray</div><div class='xr-array-name'>'rmag'</div><ul class='xr-dim-list'><li><span class='xr-has-index'>time (d)</span>: 333</li></ul></div><ul class='xr-sections'><li class='xr-section-item'><div class='xr-array-wrap'><input id='section-292226df-2122-47b0-8b1c-f02a3035728d' class='xr-array-in' type='checkbox' checked><label for='section-292226df-2122-47b0-8b1c-f02a3035728d' title='Show/hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-array-preview xr-preview'><span>0.0 2.817e-10 0.1966 0.341 0.3666 ... 1.253 1.253 1.253 1.253 1.254</span></div><div class='xr-array-data'><pre>array([0.00000000e+00, 2.81661760e-10, 1.96588167e-01, 3.41033868e-01,\n",
" * time (d) (time (d)) float64 0.0 11.0 22.0 ... 3.63e+03 3.641e+03 3.652e+03</pre><div class='xr-wrap' hidden><div class='xr-header'><div class='xr-obj-type'>xarray.DataArray</div><div class='xr-array-name'>'rmag'</div><ul class='xr-dim-list'><li><span class='xr-has-index'>time (d)</span>: 333</li></ul></div><ul class='xr-sections'><li class='xr-section-item'><div class='xr-array-wrap'><input id='section-934c4bab-67f7-4643-9182-363295c7239b' class='xr-array-in' type='checkbox' checked><label for='section-934c4bab-67f7-4643-9182-363295c7239b' title='Show/hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-array-preview xr-preview'><span>0.0 2.817e-10 0.1966 0.341 0.3666 ... 1.253 1.253 1.253 1.253 1.254</span></div><div class='xr-array-data'><pre>array([0.00000000e+00, 2.81661760e-10, 1.96588167e-01, 3.41033868e-01,\n",
" 3.66560841e-01, 3.66301948e-01, 3.66249335e-01, 3.66396298e-01,\n",
" 3.66718394e-01, 3.67175496e-01, 3.67715021e-01, 3.68275152e-01,\n",
" 3.68790522e-01, 3.69199134e-01, 3.69448202e-01, 3.69499736e-01,\n",
Expand Down Expand Up @@ -662,7 +662,7 @@
" 1.26066111e+00, 1.26026226e+00, 1.25953643e+00, 1.25855297e+00,\n",
" 1.25740517e+00, 1.25620218e+00, 1.25505935e+00, 1.25408756e+00,\n",
" 1.25338232e+00, 1.25301396e+00, 1.25302049e+00, 1.25340406e+00,\n",
" 1.25413096e+00])</pre></div></div></li><li class='xr-section-item'><input id='section-6e81a036-0504-45c3-97e0-e261861e9005' class='xr-section-summary-in' type='checkbox' checked><label for='section-6e81a036-0504-45c3-97e0-e261861e9005' class='xr-section-summary' >Coordinates: <span>(2)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span>id</span></div><div class='xr-var-dims'>()</div><div class='xr-var-dtype'>int64</div><div class='xr-var-preview xr-preview'>2</div><input id='attrs-9e79e1bd-7853-4607-8b6b-c61cdac53b49' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-9e79e1bd-7853-4607-8b6b-c61cdac53b49' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-67d0383c-a38c-4c85-8c7d-b50aaca9b0dc' class='xr-var-data-in' type='checkbox'><label for='data-67d0383c-a38c-4c85-8c7d-b50aaca9b0dc' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array(2)</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>time (d)</span></div><div class='xr-var-dims'>(time (d))</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>0.0 11.0 ... 3.641e+03 3.652e+03</div><input id='attrs-7a219eba-d644-43b9-88fc-39fee45c01e3' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-7a219eba-d644-43b9-88fc-39fee45c01e3' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-63b1d9bd-4611-430c-b768-8a30cc3c5b83' class='xr-var-data-in' type='checkbox'><label for='data-63b1d9bd-4611-430c-b768-8a30cc3c5b83' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([ 0., 11., 22., ..., 3630., 3641., 3652.])</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-4965cf8c-f6f8-45fe-80e6-f0e14fb44ccb' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-4965cf8c-f6f8-45fe-80e6-f0e14fb44ccb' class='xr-section-summary' title='Expand/collapse section'>Attributes: <span>(0)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><dl class='xr-attrs'></dl></div></li></ul></div></div>"
" 1.25413096e+00])</pre></div></div></li><li class='xr-section-item'><input id='section-8c0d58c8-a907-42a2-908d-53c85ec339a2' class='xr-section-summary-in' type='checkbox' checked><label for='section-8c0d58c8-a907-42a2-908d-53c85ec339a2' class='xr-section-summary' >Coordinates: <span>(2)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span>id</span></div><div class='xr-var-dims'>()</div><div class='xr-var-dtype'>int64</div><div class='xr-var-preview xr-preview'>2</div><input id='attrs-82f0be55-5674-40bd-86cc-ee0f43675535' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-82f0be55-5674-40bd-86cc-ee0f43675535' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-b221e6ba-0310-46c9-88b0-d34cc6482ca6' class='xr-var-data-in' type='checkbox'><label for='data-b221e6ba-0310-46c9-88b0-d34cc6482ca6' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array(2)</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>time (d)</span></div><div class='xr-var-dims'>(time (d))</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>0.0 11.0 ... 3.641e+03 3.652e+03</div><input id='attrs-58db77b4-04a6-4a79-9adf-42a7047607ed' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-58db77b4-04a6-4a79-9adf-42a7047607ed' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-23b94260-97d0-463c-aa88-bb24dd635e26' class='xr-var-data-in' type='checkbox'><label for='data-23b94260-97d0-463c-aa88-bb24dd635e26' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([ 0., 11., 22., ..., 3630., 3641., 3652.])</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-ec271324-fb01-4492-a273-5faa3bacbe74' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-ec271324-fb01-4492-a273-5faa3bacbe74' class='xr-section-summary' title='Expand/collapse section'>Attributes: <span>(0)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><dl class='xr-attrs'></dl></div></li></ul></div></div>"
],
"text/plain": [
"<xarray.DataArray 'rmag' (time (d): 333)>\n",
Expand Down
2 changes: 1 addition & 1 deletion src/setup/setup.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,6 @@ module subroutine setup_construct_system(system, param)
return
end subroutine setup_construct_system


module subroutine setup_initialize_system(self, param)
!! author: David A. Minton
!!
Expand All @@ -83,6 +82,7 @@ module subroutine setup_initialize_system(self, param)
call self%cb%initialize(param)
call self%pl%initialize(param)
call self%tp%initialize(param)
call util_valid(self%pl, self%tp)
call self%set_msys()
call self%pl%set_mu(self%cb)
call self%tp%set_mu(self%cb)
Expand Down
279 changes: 47 additions & 232 deletions src/util/util_sort.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,260 +4,75 @@
module subroutine util_sort_dp(arr)
!! author: David A. Minton
!!
!! Sort input double precision array into ascending numerical order using Quicksort algorithm
!! Sort input double precision array into ascending numerical order using insertion sort.
!! This algorithm works well for partially sorted arrays (which is usually the case here)
!!
!! Adapted from David E. Kaufmann's Swifter routine: util_sort_dp.f90
!! Adapted from Numerical Recipes in Fortran 90: The Art of Parallel Scientific Computing, by Press, Teukolsky,
!! Vetterling, and Flannery, 2nd ed., pp. 1169-70
implicit none
! Arguments
real(DP), dimension(:), intent(inout) :: arr
! Internals
integer(I4B), parameter :: NN = 15, NSTACK = 50
real(DP) :: a, dum
integer(I4B) :: n, k, i, j, jstack, l, r
integer(I4B), dimension(NSTACK) :: istack
real(DP) :: tmp
integer(I4B) :: n, i, j

! executable code
n = size(arr)
jstack = 0
l = 1
r = n
do
if ((r - l) < NN) then
do j = l + 1, r
a = arr(j)
do i = j - 1, l, -1
if (arr(i) <= a) exit
arr(i+1) = arr(i)
end do
arr(i+1) = a
end do
if (jstack == 0) return
r = istack(jstack)
l = istack(jstack-1)
jstack = jstack - 2
else
k = (l + r)/2
dum = arr(k); arr(k) = arr(l+1); arr(l+1) = dum
if (arr(l) > arr(r)) then
dum = arr(l); arr(l) = arr(r); arr(r) = dum
end if
if (arr(l+1) > arr(r)) then
dum = arr(l+1); arr(l+1) = arr(r); arr(r) = dum
end if
if (arr(l) > arr(l+1)) then
dum = arr(l); arr(l) = arr(l+1); arr(l+1) = dum
end if
i = l + 1
j = r
a = arr(l+1)
do
do
i = i + 1
if (arr(i) >= a) exit
end do
do
j = j - 1
if (arr(j) <= a) exit
end do
if (j < i) exit
dum = arr(i); arr(i) = arr(j); arr(j) = dum
end do
arr(l+1) = arr(j)
arr(j) = a
jstack = jstack + 2
if (jstack > NSTACK) then
write(*, *) "Swiftest Error:"
write(*, *) " NSTACK too small in util_sort_I4B"
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
do i = 2, n
tmp = arr(i)
do j = i - 1, 1, -1
if (arr(j) <= tmp) exit
arr(j + 1) = arr(j)
end do
arr(j + 1) = tmp
end do

return

end subroutine util_sort_dp

module subroutine util_sort_i4b(arr)
!! author: David A. Minton
!!
!! Sort input double precision array into ascending numerical order using Quicksort algorithm
!! Sort input integer array into ascending numerical order using insertion sort.
!! This algorithm works well for partially sorted arrays (which is usually the case here)
!!
!! Adapted from David E. Kaufmann's Swifter routine: util_sort_i4b.f90
!! Adapted from Numerical Recipes in Fortran 90: The Art of Parallel Scientific Computing, by Press, Teukolsky,
!! Vetterling, and Flannery, 2nd ed., pp. 1169-70
implicit none
! Arguments
integer(I4B), dimension(:), intent(inout) :: arr
! Internals
integer(I4B), parameter :: NN = 15, NSTACK = 50
integer(I4B) :: a, n, k, i, j, jstack, l, r, dum
integer(I4B), dimension(NSTACK) :: istack

! executable code
integer(I4B) :: tmp
integer(I4B) :: n, i, j

n = size(arr)
jstack = 0
l = 1
r = n
do
if ((r - l) < NN) then
do j = l + 1, r
a = arr(j)
do i = j - 1, l, -1
if (arr(i) <= a) exit
arr(i+1) = arr(i)
end do
arr(i+1) = a
end do
if (jstack == 0) return
r = istack(jstack)
l = istack(jstack-1)
jstack = jstack - 2
else
k = (l + r)/2
dum = arr(k); arr(k) = arr(l+1); arr(l+1) = dum
if (arr(l) > arr(r)) then
dum = arr(l); arr(l) = arr(r); arr(r) = dum
end if
if (arr(l+1) > arr(r)) then
dum = arr(l+1); arr(l+1) = arr(r); arr(r) = dum
end if
if (arr(l) > arr(l+1)) then
dum = arr(l); arr(l) = arr(l+1); arr(l+1) = dum
end if
i = l + 1
j = r
a = arr(l+1)
do
do
i = i + 1
if (arr(i) >= a) exit
end do
do
j = j - 1
if (arr(j) <= a) exit
end do
if (j < i) exit
dum = arr(i); arr(i) = arr(j); arr(j) = dum
end do
arr(l+1) = arr(j)
arr(j) = a
jstack = jstack + 2
if (jstack > NSTACK) then
write(*, *) "Swiftest Error:"
write(*, *) " NSTACK too small in util_sort_i4b"
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
do i = 2, n
tmp = arr(i)
do j = i - 1, 1, -1
if (arr(j) <= tmp) exit
arr(j + 1) = arr(j)
end do
arr(j + 1) = tmp
end do

return

end subroutine util_sort_i4b
end subroutine util_sort_i4b

module subroutine util_sort_sp(arr)
!! author: David A. Minton
!!
!! Sort input single precision array into 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(inout) :: arr
! Internals
real(SP) :: tmp
integer(I4B) :: n, i, j

module subroutine util_sort_sp(arr)
!! author: David A. Minton
!!
!! Sort input single precision array into ascending numerical order using Quicksort algorithm
!!
!! Adapted from David E. Kaufmann's Swifter routine: util_sort_DP.f90
!! Adapted from Numerical Recipes in Fortran 90: The Art of Parallel Scientific Computing, by Press, Teukolsky,
!! Vetterling, and Flannery, 2nd ed., pp. 1169-70
implicit none
! Arguments
real(SP), dimension(:), intent(inout) :: arr
! Internals
integer(I4B), parameter :: NN = 15, NSTACK = 50
real(SP) :: a, dum
integer(I4B) :: n, k, i, j, jstack, l, r
integer(I4B), dimension(NSTACK) :: istack

! executable code
n = size(arr)
jstack = 0
l = 1
r = n
do
if ((r - l) < NN) then
do j = l + 1, r
a = arr(j)
do i = j - 1, l, -1
if (arr(i) <= a) exit
arr(i+1) = arr(i)
end do
arr(i+1) = a
end do
if (jstack == 0) return
r = istack(jstack)
l = istack(jstack-1)
jstack = jstack - 2
else
k = (l + r)/2
dum = arr(k); arr(k) = arr(l+1); arr(l+1) = dum
if (arr(l) > arr(r)) then
dum = arr(l); arr(l) = arr(r); arr(r) = dum
end if
if (arr(l+1) > arr(r)) then
dum = arr(l+1); arr(l+1) = arr(r); arr(r) = dum
end if
if (arr(l) > arr(l+1)) then
dum = arr(l); arr(l) = arr(l+1); arr(l+1) = dum
end if
i = l + 1
j = r
a = arr(l+1)
do
do
i = i + 1
if (arr(i) >= a) exit
end do
do
j = j - 1
if (arr(j) <= a) exit
end do
if (j < i) exit
dum = arr(i); arr(i) = arr(j); arr(j) = dum
end do
arr(l+1) = arr(j)
arr(j) = a
jstack = jstack + 2
if (jstack > NSTACK) then
write(*, *) "Swiftest Error:"
write(*, *) " NSTACK too small in util_sort_I4B"
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
n = size(arr)
do i = 2, n
tmp = arr(i)
do j = i - 1, 1, -1
if (arr(j) <= tmp) exit
arr(j + 1) = arr(j)
end do

return

end subroutine util_sort_sp
arr(j + 1) = tmp
end do
return
end subroutine util_sort_sp
end submodule s_util_sort
1 change: 0 additions & 1 deletion src/util/util_valid.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,6 @@ module subroutine util_valid(pl, tp)
call util_exit(FAILURE)
end if
end do
deallocate(idarr)
end associate

return
Expand Down

0 comments on commit 49b18b0

Please sign in to comment.