From c01b21ab7aa060139fb6f9d197b4501088d60c9a Mon Sep 17 00:00:00 2001
From: David A Minton <daminton@purdue.edu>
Date: Fri, 22 Oct 2021 22:37:17 -0400
Subject: [PATCH] Added sort with better documentation

---
 src/util/util_sort.f90 | 31 ++++++++++++++++++++++---------
 1 file changed, 22 insertions(+), 9 deletions(-)

diff --git a/src/util/util_sort.f90 b/src/util/util_sort.f90
index a662eb52c..1b09b84f7 100644
--- a/src/util/util_sort.f90
+++ b/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