# quicksort.F90 Source File

## Source Code

```#if !defined(SX)

#define srt_ind(i) (((i)*nskip_)+1)

module sort_mod_int
use util_mod, only: stop_all
use util_mod_comparisons, only: operator(.arrgt.), operator(.arrlt.)

use constants
use SystemData, only: Symmetry, operator(.eq.), operator(.ne.), &
operator(.gt.), operator(.lt.), assignment(=)
use symdata, only: SymPairProd, operator(.eq.), operator(.ne.), &
operator(.gt.), operator(.lt.), assignment(=)
implicit none

private
public :: sort

! Private operator for sorting complex numbers by their real values
interface operator(.gt.)
module procedure cmplx_gt_int
end interface

interface operator(.lt.)
module procedure cmplx_lt_int
end interface

interface sort
module procedure sort_int
end interface

interface cmplx_gt
module procedure cmplx_gt_int
end interface

interface cmplx_lt
module procedure cmplx_lt_int
end interface

contains
pure subroutine sort_int (arr, nskip, par)

! Perform a quicksort on an array of integers, arr(n). Uses the
! sample code in NumericalRecipies as a base.
! Optionally sort arr2 in parallel (in the routines it is enabled)

! Custom comparisons. Use the operator .locallt. & .localgt.
! interface operator(.locallt.)
!     pure function custom_lt (b, c) result (ret)
!         use constants
!         integer(kind=int32), intent(in) :: b(:)
!         integer(kind=int32), intent(in) :: c(:)
!         logical :: ret
!     end function
! end interface

! interface operator(.localgt.)
!     pure function custom_gt (b, c) result (ret)
!         use constants
!         integer(kind=int32), intent(in) :: b(:)
!         integer(kind=int32), intent(in) :: c(:)
!         logical :: ret
!     end function
! end interface

! sort needs auxiliary storage of length 2*log_2(n).
integer, parameter :: nStackMax = 50
integer, parameter :: nInsertionSort = 7
integer, intent(in), optional :: nskip
integer, intent(out), optional :: par

integer(kind=int32), intent(inout) :: arr(:)
!, intent(inout) :: arr2(:)
!, intent(inout) :: arr3(:)
!, intent(inout) :: arr4(:)
!, intent(inout) :: arr5(:)
!, intent(inout) :: arr6(:)

! Oh, how lovely it would be to be able to use push/pop and not need
integer :: stack(nStackMax), nstack, nskip_
integer :: pivot, lo, hi, n, i, j, par_
! n.b. This size statement is removed if type1 is scalar ...
integer(kind=int32) :: a
! :: a2(size(arr2(1)))
! :: a3(size(arr3(1)))
! :: a4(size(arr4(1)))
! :: a5(size(arr5(1)))
! :: a6(size(arr6(1)))

! Temporary variables for swapping
integer(kind=int32) :: tmp1
! :: tmp2(size(arr2(1)))
! :: tmp3(size(arr3(1)))
! :: tmp4(size(arr4(1)))
! :: tmp5(size(arr5(1)))
! :: tmp6(size(arr6(1)))
character(*), parameter :: this_routine = 'sort'

! Initialise temporary variables if required
!(a)
!(tmp1)
!(a2)
!(tmp2)
!(a3)
!(tmp3)
!(a4)
!(tmp4)
!(a5)
!(tmp5)
!(a6)
!(tmp6)

! *** HACK ***
! Workaround for gfortran compiler bug
! n.b. This will produce spurious warnings if run in valgrind, as
!      a is not initialised. Not a problem in optimised build.
i = 0
#ifdef DEBUG_
! it IS a problem in debug build, however
! *** HACK MOD ***
!to prevent crashes in debug mode
if((ubound(arr,1) - lbound(arr,1) + 1) &
> 0 ) a = arr(1)
#endif
! if (custom_lt(a, a)) i = i
! if (custom_gt(a, a)) i = i

if (present(nskip)) then
nskip_ = nskip
else
nskip_ = 1
endif

! Initialise parity calculation
par_ = 1

! The size of the array to sort.
! N.B. Use zero based indices. To obtain ! the entry into the actual
! array, multiply indices by nskip and add ! 1 (hence zero based)
! **** See local macro srt_ind() ******
lo = lbound(arr, 1) - 1
n = ((ubound(arr, 1) - lo -1)/nskip_) + 1
hi = lo + n - 1

nstack = 0
do while (.true.)
! If the section/partition we are looking at is smaller than
! nInsertSort then perform an insertion sort
if (hi - lo < nInsertionSort) then
do j = lo + 1, hi
a = arr(srt_ind(j))
! a2 = arr2(srt_ind(j))
! a3 = arr3(srt_ind(j))
! a4 = arr4(srt_ind(j))
! a5 = arr5(srt_ind(j))
! a6 = arr6(srt_ind(j))
do i = j - 1, 0, -1
if (arr(srt_ind(i)) < a) exit
arr(srt_ind(i+1)) = arr(srt_ind(i))
! arr2(srt_ind(i+1)) = arr2(srt_ind(i))
! arr3(srt_ind(i+1)) = arr3(srt_ind(i))
! arr4(srt_ind(i+1)) = arr4(srt_ind(i))
! arr5(srt_ind(i+1)) = arr5(srt_ind(i))
! arr6(srt_ind(i+1)) = arr6(srt_ind(i))
par_ = -par_
enddo
arr(srt_ind(i+1)) = a
! arr2(srt_ind(i+1)) = a2
! arr3(srt_ind(i+1)) = a3
! arr4(srt_ind(i+1)) = a4
! arr5(srt_ind(i+1)) = a5
! arr6(srt_ind(i+1)) = a6
enddo

if (nstack == 0) exit
hi = stack(nstack)
lo = stack(nstack-1)
nstack = nstack - 2

! Otherwise start partitioning with quicksort.
else
! Pick a partitioning element, and arrange such that
! arr(lo) <= arr(lo+1) <= arr(hi)
pivot = (lo + hi) / 2
tmp1 = arr(srt_ind(pivot))
arr(srt_ind(pivot)) = arr(srt_ind(lo + 1))
arr(srt_ind(lo + 1)) = tmp1
! tmp2 = arr2(srt_ind(pivot))
! arr2(srt_ind(pivot)) = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = tmp2
! tmp3 = arr3(srt_ind(pivot))
! arr3(srt_ind(pivot)) = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = tmp3
! tmp4 = arr4(srt_ind(pivot))
! arr4(srt_ind(pivot)) = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = tmp4
! tmp5 = arr5(srt_ind(pivot))
! arr5(srt_ind(pivot)) = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = tmp5
! tmp6 = arr6(srt_ind(pivot))
! arr6(srt_ind(pivot)) = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = tmp6
par_ = -par_

if (arr(srt_ind(lo)) > arr(srt_ind(hi))) then
tmp1 = arr(srt_ind(lo))
arr(srt_ind(lo)) = arr(srt_ind(hi))
arr(srt_ind(hi)) = tmp1
! tmp2 = arr2(srt_ind(lo))
! arr2(srt_ind(lo)) = arr2(srt_ind(hi))
! arr2(srt_ind(hi)) = tmp2
! tmp3 = arr3(srt_ind(lo))
! arr3(srt_ind(lo)) = arr3(srt_ind(hi))
! arr3(srt_ind(hi)) = tmp3
! tmp4 = arr4(srt_ind(lo))
! arr4(srt_ind(lo)) = arr4(srt_ind(hi))
! arr4(srt_ind(hi)) = tmp4
! tmp5 = arr5(srt_ind(lo))
! arr5(srt_ind(lo)) = arr5(srt_ind(hi))
! arr5(srt_ind(hi)) = tmp5
! tmp6 = arr6(srt_ind(lo))
! arr6(srt_ind(lo)) = arr6(srt_ind(hi))
! arr6(srt_ind(hi)) = tmp6
par_ = -par_
endif
if (arr(srt_ind(lo+1)) > arr(srt_ind(hi))) then
tmp1 = arr(srt_ind(lo+1))
arr(srt_ind(lo+1)) = arr(srt_ind(hi))
arr(srt_ind(hi)) = tmp1
! tmp2 = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = arr2(srt_ind(hi))
! arr2(srt_ind(hi)) = tmp2
! tmp3 = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = arr3(srt_ind(hi))
! arr3(srt_ind(hi)) = tmp3
! tmp4 = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = arr4(srt_ind(hi))
! arr4(srt_ind(hi)) = tmp4
! tmp5 = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = arr5(srt_ind(hi))
! arr5(srt_ind(hi)) = tmp5
! tmp6 = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = arr6(srt_ind(hi))
! arr6(srt_ind(hi)) = tmp6
par_ = -par_
endif
if (arr(srt_ind(lo)) > arr(srt_ind(lo+1))) then
tmp1 = arr(srt_ind(lo))
arr(srt_ind(lo)) = arr(srt_ind(lo+1))
arr(srt_ind(lo+1)) = tmp1
! tmp2 = arr2(srt_ind(lo))
! arr2(srt_ind(lo)) = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = tmp2
! tmp3 = arr3(srt_ind(lo))
! arr3(srt_ind(lo)) = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = tmp3
! tmp4 = arr4(srt_ind(lo))
! arr4(srt_ind(lo)) = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = tmp4
! tmp5 = arr5(srt_ind(lo))
! arr5(srt_ind(lo)) = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = tmp5
! tmp6 = arr6(srt_ind(lo))
! arr6(srt_ind(lo)) = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = tmp6
par_ = -par_
endif

i = lo + 1
j = hi
a = arr(srt_ind(lo + 1)) !! a is the pivot value
! a2 = arr2(srt_ind(lo + 1))
! a3 = arr3(srt_ind(lo + 1))
! a4 = arr4(srt_ind(lo + 1))
! a5 = arr5(srt_ind(lo + 1))
! a6 = arr6(srt_ind(lo + 1))
do while (.true.)
! Scand down list to find element > a
i = i + 1
do while (arr(srt_ind(i)) < a)
i = i + 1
enddo

! Scan down list to find element < a
j = j - 1
do while (arr(srt_ind(j)) > a)
j = j - 1
enddo

! When the pointers crossed, partitioning is complete.
if (j < i) exit

! Swap the elements, so that all elements < a end up
! in lower indexed variables.
tmp1 = arr(srt_ind(i))
arr(srt_ind(i)) = arr(srt_ind(j))
arr(srt_ind(j)) = tmp1
! tmp2 = arr2(srt_ind(i))
! arr2(srt_ind(i)) = arr2(srt_ind(j))
! arr2(srt_ind(j)) = tmp2
! tmp3 = arr3(srt_ind(i))
! arr3(srt_ind(i)) = arr3(srt_ind(j))
! arr3(srt_ind(j)) = tmp3
! tmp4 = arr4(srt_ind(i))
! arr4(srt_ind(i)) = arr4(srt_ind(j))
! arr4(srt_ind(j)) = tmp4
! tmp5 = arr5(srt_ind(i))
! arr5(srt_ind(i)) = arr5(srt_ind(j))
! arr5(srt_ind(j)) = tmp5
! tmp6 = arr6(srt_ind(i))
! arr6(srt_ind(i)) = arr6(srt_ind(j))
! arr6(srt_ind(j)) = tmp6
par_ = -par_
enddo;

! Insert partitioning element
arr(srt_ind(lo + 1)) = arr(srt_ind(j))
arr(srt_ind(j)) = a
! arr2(srt_ind(lo + 1)) = arr2(srt_ind(j))
! arr3(srt_ind(lo + 1)) = arr3(srt_ind(j))
! arr4(srt_ind(lo + 1)) = arr4(srt_ind(j))
! arr5(srt_ind(lo + 1)) = arr5(srt_ind(j))
! arr6(srt_ind(lo + 1)) = arr6(srt_ind(j))
! arr2(srt_ind(j)) = a2
! arr3(srt_ind(j)) = a3
! arr4(srt_ind(j)) = a4
! arr5(srt_ind(j)) = a5
! arr6(srt_ind(j)) = a6
par_ = -par_

! Push the larger of the partitioned sections onto the stack
! of sections to look at later.
! --> need fewest stack elements.
nstack = nstack + 2
if (nstack > nStackMax) then
call stop_all (this_routine, &
"parameter nStackMax too small")
endif
if (hi - i + 1 >= j - lo) then
stack(nstack) = hi
stack(nstack-1) = i
hi = j - 1
else
stack(nstack) = j - 1
stack(nstack-1) = lo
lo = i
endif
endif
enddo

! Destroy temporary variables if required
!(a)
!(tmp1)
!(a2)
!(tmp2)
!(a3)
!(tmp3)
!(a4)
!(tmp4)
!(a5)
!(tmp5)
!(a6)
!(tmp6)

! Return the parity if required
if (present(par)) &
par = par_

end subroutine

! A private comparison. This should not be available outside of the
! module!
elemental function cmplx_gt_int (a, b) result (bGt)
complex(dp), intent(in) :: a, b
logical :: bGt

bGt = real(a, dp) > real(b, dp)
end function

elemental function cmplx_lt_int (a, b) result (bLt)
complex(dp), intent(in) :: a, b
logical :: bLt

bLt = real(a, dp) < real(b, dp)
end function

end module

#endif

#define srt_ind(i) (((i)*nskip_)+1)

module sort_mod_int64
use util_mod, only: stop_all
use util_mod_comparisons, only: operator(.arrgt.), operator(.arrlt.)

use constants
use SystemData, only: Symmetry, operator(.eq.), operator(.ne.), &
operator(.gt.), operator(.lt.), assignment(=)
use symdata, only: SymPairProd, operator(.eq.), operator(.ne.), &
operator(.gt.), operator(.lt.), assignment(=)
implicit none

private
public :: sort

! Private operator for sorting complex numbers by their real values
interface operator(.gt.)
module procedure cmplx_gt_int64
end interface

interface operator(.lt.)
module procedure cmplx_lt_int64
end interface

interface sort
module procedure sort_int64
end interface

interface cmplx_gt
module procedure cmplx_gt_int64
end interface

interface cmplx_lt
module procedure cmplx_lt_int64
end interface

contains
pure subroutine sort_int64 (arr, nskip, par)

! Perform a quicksort on an array of integers, arr(n). Uses the
! sample code in NumericalRecipies as a base.
! Optionally sort arr2 in parallel (in the routines it is enabled)

! Custom comparisons. Use the operator .locallt. & .localgt.
! interface operator(.locallt.)
!     pure function custom_lt (b, c) result (ret)
!         use constants
!         integer(kind=int64), intent(in) :: b(:)
!         integer(kind=int64), intent(in) :: c(:)
!         logical :: ret
!     end function
! end interface

! interface operator(.localgt.)
!     pure function custom_gt (b, c) result (ret)
!         use constants
!         integer(kind=int64), intent(in) :: b(:)
!         integer(kind=int64), intent(in) :: c(:)
!         logical :: ret
!     end function
! end interface

! sort needs auxiliary storage of length 2*log_2(n).
integer, parameter :: nStackMax = 50
integer, parameter :: nInsertionSort = 7
integer, intent(in), optional :: nskip
integer, intent(out), optional :: par

integer(kind=int64), intent(inout) :: arr(:)
!, intent(inout) :: arr2(:)
!, intent(inout) :: arr3(:)
!, intent(inout) :: arr4(:)
!, intent(inout) :: arr5(:)
!, intent(inout) :: arr6(:)

! Oh, how lovely it would be to be able to use push/pop and not need
integer :: stack(nStackMax), nstack, nskip_
integer :: pivot, lo, hi, n, i, j, par_
! n.b. This size statement is removed if type1 is scalar ...
integer(kind=int64) :: a
! :: a2(size(arr2(1)))
! :: a3(size(arr3(1)))
! :: a4(size(arr4(1)))
! :: a5(size(arr5(1)))
! :: a6(size(arr6(1)))

! Temporary variables for swapping
integer(kind=int64) :: tmp1
! :: tmp2(size(arr2(1)))
! :: tmp3(size(arr3(1)))
! :: tmp4(size(arr4(1)))
! :: tmp5(size(arr5(1)))
! :: tmp6(size(arr6(1)))
character(*), parameter :: this_routine = 'sort'

! Initialise temporary variables if required
!(a)
!(tmp1)
!(a2)
!(tmp2)
!(a3)
!(tmp3)
!(a4)
!(tmp4)
!(a5)
!(tmp5)
!(a6)
!(tmp6)

! *** HACK ***
! Workaround for gfortran compiler bug
! n.b. This will produce spurious warnings if run in valgrind, as
!      a is not initialised. Not a problem in optimised build.
i = 0
#ifdef DEBUG_
! it IS a problem in debug build, however
! *** HACK MOD ***
!to prevent crashes in debug mode
if((ubound(arr,1) - lbound(arr,1) + 1) &
> 0 ) a = arr(1)
#endif
! if (custom_lt(a, a)) i = i
! if (custom_gt(a, a)) i = i

if (present(nskip)) then
nskip_ = nskip
else
nskip_ = 1
endif

! Initialise parity calculation
par_ = 1

! The size of the array to sort.
! N.B. Use zero based indices. To obtain ! the entry into the actual
! array, multiply indices by nskip and add ! 1 (hence zero based)
! **** See local macro srt_ind() ******
lo = lbound(arr, 1) - 1
n = ((ubound(arr, 1) - lo -1)/nskip_) + 1
hi = lo + n - 1

nstack = 0
do while (.true.)
! If the section/partition we are looking at is smaller than
! nInsertSort then perform an insertion sort
if (hi - lo < nInsertionSort) then
do j = lo + 1, hi
a = arr(srt_ind(j))
! a2 = arr2(srt_ind(j))
! a3 = arr3(srt_ind(j))
! a4 = arr4(srt_ind(j))
! a5 = arr5(srt_ind(j))
! a6 = arr6(srt_ind(j))
do i = j - 1, 0, -1
if (arr(srt_ind(i)) < a) exit
arr(srt_ind(i+1)) = arr(srt_ind(i))
! arr2(srt_ind(i+1)) = arr2(srt_ind(i))
! arr3(srt_ind(i+1)) = arr3(srt_ind(i))
! arr4(srt_ind(i+1)) = arr4(srt_ind(i))
! arr5(srt_ind(i+1)) = arr5(srt_ind(i))
! arr6(srt_ind(i+1)) = arr6(srt_ind(i))
par_ = -par_
enddo
arr(srt_ind(i+1)) = a
! arr2(srt_ind(i+1)) = a2
! arr3(srt_ind(i+1)) = a3
! arr4(srt_ind(i+1)) = a4
! arr5(srt_ind(i+1)) = a5
! arr6(srt_ind(i+1)) = a6
enddo

if (nstack == 0) exit
hi = stack(nstack)
lo = stack(nstack-1)
nstack = nstack - 2

! Otherwise start partitioning with quicksort.
else
! Pick a partitioning element, and arrange such that
! arr(lo) <= arr(lo+1) <= arr(hi)
pivot = (lo + hi) / 2
tmp1 = arr(srt_ind(pivot))
arr(srt_ind(pivot)) = arr(srt_ind(lo + 1))
arr(srt_ind(lo + 1)) = tmp1
! tmp2 = arr2(srt_ind(pivot))
! arr2(srt_ind(pivot)) = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = tmp2
! tmp3 = arr3(srt_ind(pivot))
! arr3(srt_ind(pivot)) = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = tmp3
! tmp4 = arr4(srt_ind(pivot))
! arr4(srt_ind(pivot)) = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = tmp4
! tmp5 = arr5(srt_ind(pivot))
! arr5(srt_ind(pivot)) = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = tmp5
! tmp6 = arr6(srt_ind(pivot))
! arr6(srt_ind(pivot)) = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = tmp6
par_ = -par_

if (arr(srt_ind(lo)) > arr(srt_ind(hi))) then
tmp1 = arr(srt_ind(lo))
arr(srt_ind(lo)) = arr(srt_ind(hi))
arr(srt_ind(hi)) = tmp1
! tmp2 = arr2(srt_ind(lo))
! arr2(srt_ind(lo)) = arr2(srt_ind(hi))
! arr2(srt_ind(hi)) = tmp2
! tmp3 = arr3(srt_ind(lo))
! arr3(srt_ind(lo)) = arr3(srt_ind(hi))
! arr3(srt_ind(hi)) = tmp3
! tmp4 = arr4(srt_ind(lo))
! arr4(srt_ind(lo)) = arr4(srt_ind(hi))
! arr4(srt_ind(hi)) = tmp4
! tmp5 = arr5(srt_ind(lo))
! arr5(srt_ind(lo)) = arr5(srt_ind(hi))
! arr5(srt_ind(hi)) = tmp5
! tmp6 = arr6(srt_ind(lo))
! arr6(srt_ind(lo)) = arr6(srt_ind(hi))
! arr6(srt_ind(hi)) = tmp6
par_ = -par_
endif
if (arr(srt_ind(lo+1)) > arr(srt_ind(hi))) then
tmp1 = arr(srt_ind(lo+1))
arr(srt_ind(lo+1)) = arr(srt_ind(hi))
arr(srt_ind(hi)) = tmp1
! tmp2 = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = arr2(srt_ind(hi))
! arr2(srt_ind(hi)) = tmp2
! tmp3 = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = arr3(srt_ind(hi))
! arr3(srt_ind(hi)) = tmp3
! tmp4 = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = arr4(srt_ind(hi))
! arr4(srt_ind(hi)) = tmp4
! tmp5 = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = arr5(srt_ind(hi))
! arr5(srt_ind(hi)) = tmp5
! tmp6 = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = arr6(srt_ind(hi))
! arr6(srt_ind(hi)) = tmp6
par_ = -par_
endif
if (arr(srt_ind(lo)) > arr(srt_ind(lo+1))) then
tmp1 = arr(srt_ind(lo))
arr(srt_ind(lo)) = arr(srt_ind(lo+1))
arr(srt_ind(lo+1)) = tmp1
! tmp2 = arr2(srt_ind(lo))
! arr2(srt_ind(lo)) = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = tmp2
! tmp3 = arr3(srt_ind(lo))
! arr3(srt_ind(lo)) = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = tmp3
! tmp4 = arr4(srt_ind(lo))
! arr4(srt_ind(lo)) = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = tmp4
! tmp5 = arr5(srt_ind(lo))
! arr5(srt_ind(lo)) = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = tmp5
! tmp6 = arr6(srt_ind(lo))
! arr6(srt_ind(lo)) = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = tmp6
par_ = -par_
endif

i = lo + 1
j = hi
a = arr(srt_ind(lo + 1)) !! a is the pivot value
! a2 = arr2(srt_ind(lo + 1))
! a3 = arr3(srt_ind(lo + 1))
! a4 = arr4(srt_ind(lo + 1))
! a5 = arr5(srt_ind(lo + 1))
! a6 = arr6(srt_ind(lo + 1))
do while (.true.)
! Scand down list to find element > a
i = i + 1
do while (arr(srt_ind(i)) < a)
i = i + 1
enddo

! Scan down list to find element < a
j = j - 1
do while (arr(srt_ind(j)) > a)
j = j - 1
enddo

! When the pointers crossed, partitioning is complete.
if (j < i) exit

! Swap the elements, so that all elements < a end up
! in lower indexed variables.
tmp1 = arr(srt_ind(i))
arr(srt_ind(i)) = arr(srt_ind(j))
arr(srt_ind(j)) = tmp1
! tmp2 = arr2(srt_ind(i))
! arr2(srt_ind(i)) = arr2(srt_ind(j))
! arr2(srt_ind(j)) = tmp2
! tmp3 = arr3(srt_ind(i))
! arr3(srt_ind(i)) = arr3(srt_ind(j))
! arr3(srt_ind(j)) = tmp3
! tmp4 = arr4(srt_ind(i))
! arr4(srt_ind(i)) = arr4(srt_ind(j))
! arr4(srt_ind(j)) = tmp4
! tmp5 = arr5(srt_ind(i))
! arr5(srt_ind(i)) = arr5(srt_ind(j))
! arr5(srt_ind(j)) = tmp5
! tmp6 = arr6(srt_ind(i))
! arr6(srt_ind(i)) = arr6(srt_ind(j))
! arr6(srt_ind(j)) = tmp6
par_ = -par_
enddo;

! Insert partitioning element
arr(srt_ind(lo + 1)) = arr(srt_ind(j))
arr(srt_ind(j)) = a
! arr2(srt_ind(lo + 1)) = arr2(srt_ind(j))
! arr3(srt_ind(lo + 1)) = arr3(srt_ind(j))
! arr4(srt_ind(lo + 1)) = arr4(srt_ind(j))
! arr5(srt_ind(lo + 1)) = arr5(srt_ind(j))
! arr6(srt_ind(lo + 1)) = arr6(srt_ind(j))
! arr2(srt_ind(j)) = a2
! arr3(srt_ind(j)) = a3
! arr4(srt_ind(j)) = a4
! arr5(srt_ind(j)) = a5
! arr6(srt_ind(j)) = a6
par_ = -par_

! Push the larger of the partitioned sections onto the stack
! of sections to look at later.
! --> need fewest stack elements.
nstack = nstack + 2
if (nstack > nStackMax) then
call stop_all (this_routine, &
"parameter nStackMax too small")
endif
if (hi - i + 1 >= j - lo) then
stack(nstack) = hi
stack(nstack-1) = i
hi = j - 1
else
stack(nstack) = j - 1
stack(nstack-1) = lo
lo = i
endif
endif
enddo

! Destroy temporary variables if required
!(a)
!(tmp1)
!(a2)
!(tmp2)
!(a3)
!(tmp3)
!(a4)
!(tmp4)
!(a5)
!(tmp5)
!(a6)
!(tmp6)

! Return the parity if required
if (present(par)) &
par = par_

end subroutine

! A private comparison. This should not be available outside of the
! module!
elemental function cmplx_gt_int64 (a, b) result (bGt)
complex(dp), intent(in) :: a, b
logical :: bGt

bGt = real(a, dp) > real(b, dp)
end function

elemental function cmplx_lt_int64 (a, b) result (bLt)
complex(dp), intent(in) :: a, b
logical :: bLt

bLt = real(a, dp) < real(b, dp)
end function

end module

#define srt_ind(i) (((i)*nskip_)+1)

module sort_mod_doub
use util_mod, only: stop_all
use util_mod_comparisons, only: operator(.arrgt.), operator(.arrlt.)

use constants
use SystemData, only: Symmetry, operator(.eq.), operator(.ne.), &
operator(.gt.), operator(.lt.), assignment(=)
use symdata, only: SymPairProd, operator(.eq.), operator(.ne.), &
operator(.gt.), operator(.lt.), assignment(=)
implicit none

private
public :: sort

! Private operator for sorting complex numbers by their real values
interface operator(.gt.)
module procedure cmplx_gt_doub
end interface

interface operator(.lt.)
module procedure cmplx_lt_doub
end interface

interface sort
module procedure sort_doub
end interface

interface cmplx_gt
module procedure cmplx_gt_doub
end interface

interface cmplx_lt
module procedure cmplx_lt_doub
end interface

contains
pure subroutine sort_doub (arr, nskip, par)

! Perform a quicksort on an array of integers, arr(n). Uses the
! sample code in NumericalRecipies as a base.
! Optionally sort arr2 in parallel (in the routines it is enabled)

! Custom comparisons. Use the operator .locallt. & .localgt.
! interface operator(.locallt.)
!     pure function custom_lt (b, c) result (ret)
!         use constants
!         real(dp), intent(in) :: b(:)
!         real(dp), intent(in) :: c(:)
!         logical :: ret
!     end function
! end interface

! interface operator(.localgt.)
!     pure function custom_gt (b, c) result (ret)
!         use constants
!         real(dp), intent(in) :: b(:)
!         real(dp), intent(in) :: c(:)
!         logical :: ret
!     end function
! end interface

! sort needs auxiliary storage of length 2*log_2(n).
integer, parameter :: nStackMax = 50
integer, parameter :: nInsertionSort = 7
integer, intent(in), optional :: nskip
integer, intent(out), optional :: par

real(dp), intent(inout) :: arr(:)
!, intent(inout) :: arr2(:)
!, intent(inout) :: arr3(:)
!, intent(inout) :: arr4(:)
!, intent(inout) :: arr5(:)
!, intent(inout) :: arr6(:)

! Oh, how lovely it would be to be able to use push/pop and not need
integer :: stack(nStackMax), nstack, nskip_
integer :: pivot, lo, hi, n, i, j, par_
! n.b. This size statement is removed if type1 is scalar ...
real(dp) :: a
! :: a2(size(arr2(1)))
! :: a3(size(arr3(1)))
! :: a4(size(arr4(1)))
! :: a5(size(arr5(1)))
! :: a6(size(arr6(1)))

! Temporary variables for swapping
real(dp) :: tmp1
! :: tmp2(size(arr2(1)))
! :: tmp3(size(arr3(1)))
! :: tmp4(size(arr4(1)))
! :: tmp5(size(arr5(1)))
! :: tmp6(size(arr6(1)))
character(*), parameter :: this_routine = 'sort'

! Initialise temporary variables if required
!(a)
!(tmp1)
!(a2)
!(tmp2)
!(a3)
!(tmp3)
!(a4)
!(tmp4)
!(a5)
!(tmp5)
!(a6)
!(tmp6)

! *** HACK ***
! Workaround for gfortran compiler bug
! n.b. This will produce spurious warnings if run in valgrind, as
!      a is not initialised. Not a problem in optimised build.
i = 0
#ifdef DEBUG_
! it IS a problem in debug build, however
! *** HACK MOD ***
!to prevent crashes in debug mode
if((ubound(arr,1) - lbound(arr,1) + 1) &
> 0 ) a = arr(1)
#endif
! if (custom_lt(a, a)) i = i
! if (custom_gt(a, a)) i = i

if (present(nskip)) then
nskip_ = nskip
else
nskip_ = 1
endif

! Initialise parity calculation
par_ = 1

! The size of the array to sort.
! N.B. Use zero based indices. To obtain ! the entry into the actual
! array, multiply indices by nskip and add ! 1 (hence zero based)
! **** See local macro srt_ind() ******
lo = lbound(arr, 1) - 1
n = ((ubound(arr, 1) - lo -1)/nskip_) + 1
hi = lo + n - 1

nstack = 0
do while (.true.)
! If the section/partition we are looking at is smaller than
! nInsertSort then perform an insertion sort
if (hi - lo < nInsertionSort) then
do j = lo + 1, hi
a = arr(srt_ind(j))
! a2 = arr2(srt_ind(j))
! a3 = arr3(srt_ind(j))
! a4 = arr4(srt_ind(j))
! a5 = arr5(srt_ind(j))
! a6 = arr6(srt_ind(j))
do i = j - 1, 0, -1
if (arr(srt_ind(i)) < a) exit
arr(srt_ind(i+1)) = arr(srt_ind(i))
! arr2(srt_ind(i+1)) = arr2(srt_ind(i))
! arr3(srt_ind(i+1)) = arr3(srt_ind(i))
! arr4(srt_ind(i+1)) = arr4(srt_ind(i))
! arr5(srt_ind(i+1)) = arr5(srt_ind(i))
! arr6(srt_ind(i+1)) = arr6(srt_ind(i))
par_ = -par_
enddo
arr(srt_ind(i+1)) = a
! arr2(srt_ind(i+1)) = a2
! arr3(srt_ind(i+1)) = a3
! arr4(srt_ind(i+1)) = a4
! arr5(srt_ind(i+1)) = a5
! arr6(srt_ind(i+1)) = a6
enddo

if (nstack == 0) exit
hi = stack(nstack)
lo = stack(nstack-1)
nstack = nstack - 2

! Otherwise start partitioning with quicksort.
else
! Pick a partitioning element, and arrange such that
! arr(lo) <= arr(lo+1) <= arr(hi)
pivot = (lo + hi) / 2
tmp1 = arr(srt_ind(pivot))
arr(srt_ind(pivot)) = arr(srt_ind(lo + 1))
arr(srt_ind(lo + 1)) = tmp1
! tmp2 = arr2(srt_ind(pivot))
! arr2(srt_ind(pivot)) = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = tmp2
! tmp3 = arr3(srt_ind(pivot))
! arr3(srt_ind(pivot)) = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = tmp3
! tmp4 = arr4(srt_ind(pivot))
! arr4(srt_ind(pivot)) = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = tmp4
! tmp5 = arr5(srt_ind(pivot))
! arr5(srt_ind(pivot)) = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = tmp5
! tmp6 = arr6(srt_ind(pivot))
! arr6(srt_ind(pivot)) = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = tmp6
par_ = -par_

if (arr(srt_ind(lo)) > arr(srt_ind(hi))) then
tmp1 = arr(srt_ind(lo))
arr(srt_ind(lo)) = arr(srt_ind(hi))
arr(srt_ind(hi)) = tmp1
! tmp2 = arr2(srt_ind(lo))
! arr2(srt_ind(lo)) = arr2(srt_ind(hi))
! arr2(srt_ind(hi)) = tmp2
! tmp3 = arr3(srt_ind(lo))
! arr3(srt_ind(lo)) = arr3(srt_ind(hi))
! arr3(srt_ind(hi)) = tmp3
! tmp4 = arr4(srt_ind(lo))
! arr4(srt_ind(lo)) = arr4(srt_ind(hi))
! arr4(srt_ind(hi)) = tmp4
! tmp5 = arr5(srt_ind(lo))
! arr5(srt_ind(lo)) = arr5(srt_ind(hi))
! arr5(srt_ind(hi)) = tmp5
! tmp6 = arr6(srt_ind(lo))
! arr6(srt_ind(lo)) = arr6(srt_ind(hi))
! arr6(srt_ind(hi)) = tmp6
par_ = -par_
endif
if (arr(srt_ind(lo+1)) > arr(srt_ind(hi))) then
tmp1 = arr(srt_ind(lo+1))
arr(srt_ind(lo+1)) = arr(srt_ind(hi))
arr(srt_ind(hi)) = tmp1
! tmp2 = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = arr2(srt_ind(hi))
! arr2(srt_ind(hi)) = tmp2
! tmp3 = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = arr3(srt_ind(hi))
! arr3(srt_ind(hi)) = tmp3
! tmp4 = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = arr4(srt_ind(hi))
! arr4(srt_ind(hi)) = tmp4
! tmp5 = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = arr5(srt_ind(hi))
! arr5(srt_ind(hi)) = tmp5
! tmp6 = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = arr6(srt_ind(hi))
! arr6(srt_ind(hi)) = tmp6
par_ = -par_
endif
if (arr(srt_ind(lo)) > arr(srt_ind(lo+1))) then
tmp1 = arr(srt_ind(lo))
arr(srt_ind(lo)) = arr(srt_ind(lo+1))
arr(srt_ind(lo+1)) = tmp1
! tmp2 = arr2(srt_ind(lo))
! arr2(srt_ind(lo)) = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = tmp2
! tmp3 = arr3(srt_ind(lo))
! arr3(srt_ind(lo)) = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = tmp3
! tmp4 = arr4(srt_ind(lo))
! arr4(srt_ind(lo)) = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = tmp4
! tmp5 = arr5(srt_ind(lo))
! arr5(srt_ind(lo)) = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = tmp5
! tmp6 = arr6(srt_ind(lo))
! arr6(srt_ind(lo)) = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = tmp6
par_ = -par_
endif

i = lo + 1
j = hi
a = arr(srt_ind(lo + 1)) !! a is the pivot value
! a2 = arr2(srt_ind(lo + 1))
! a3 = arr3(srt_ind(lo + 1))
! a4 = arr4(srt_ind(lo + 1))
! a5 = arr5(srt_ind(lo + 1))
! a6 = arr6(srt_ind(lo + 1))
do while (.true.)
! Scand down list to find element > a
i = i + 1
do while (arr(srt_ind(i)) < a)
i = i + 1
enddo

! Scan down list to find element < a
j = j - 1
do while (arr(srt_ind(j)) > a)
j = j - 1
enddo

! When the pointers crossed, partitioning is complete.
if (j < i) exit

! Swap the elements, so that all elements < a end up
! in lower indexed variables.
tmp1 = arr(srt_ind(i))
arr(srt_ind(i)) = arr(srt_ind(j))
arr(srt_ind(j)) = tmp1
! tmp2 = arr2(srt_ind(i))
! arr2(srt_ind(i)) = arr2(srt_ind(j))
! arr2(srt_ind(j)) = tmp2
! tmp3 = arr3(srt_ind(i))
! arr3(srt_ind(i)) = arr3(srt_ind(j))
! arr3(srt_ind(j)) = tmp3
! tmp4 = arr4(srt_ind(i))
! arr4(srt_ind(i)) = arr4(srt_ind(j))
! arr4(srt_ind(j)) = tmp4
! tmp5 = arr5(srt_ind(i))
! arr5(srt_ind(i)) = arr5(srt_ind(j))
! arr5(srt_ind(j)) = tmp5
! tmp6 = arr6(srt_ind(i))
! arr6(srt_ind(i)) = arr6(srt_ind(j))
! arr6(srt_ind(j)) = tmp6
par_ = -par_
enddo;

! Insert partitioning element
arr(srt_ind(lo + 1)) = arr(srt_ind(j))
arr(srt_ind(j)) = a
! arr2(srt_ind(lo + 1)) = arr2(srt_ind(j))
! arr3(srt_ind(lo + 1)) = arr3(srt_ind(j))
! arr4(srt_ind(lo + 1)) = arr4(srt_ind(j))
! arr5(srt_ind(lo + 1)) = arr5(srt_ind(j))
! arr6(srt_ind(lo + 1)) = arr6(srt_ind(j))
! arr2(srt_ind(j)) = a2
! arr3(srt_ind(j)) = a3
! arr4(srt_ind(j)) = a4
! arr5(srt_ind(j)) = a5
! arr6(srt_ind(j)) = a6
par_ = -par_

! Push the larger of the partitioned sections onto the stack
! of sections to look at later.
! --> need fewest stack elements.
nstack = nstack + 2
if (nstack > nStackMax) then
call stop_all (this_routine, &
"parameter nStackMax too small")
endif
if (hi - i + 1 >= j - lo) then
stack(nstack) = hi
stack(nstack-1) = i
hi = j - 1
else
stack(nstack) = j - 1
stack(nstack-1) = lo
lo = i
endif
endif
enddo

! Destroy temporary variables if required
!(a)
!(tmp1)
!(a2)
!(tmp2)
!(a3)
!(tmp3)
!(a4)
!(tmp4)
!(a5)
!(tmp5)
!(a6)
!(tmp6)

! Return the parity if required
if (present(par)) &
par = par_

end subroutine

! A private comparison. This should not be available outside of the
! module!
elemental function cmplx_gt_doub (a, b) result (bGt)
complex(dp), intent(in) :: a, b
logical :: bGt

bGt = real(a, dp) > real(b, dp)
end function

elemental function cmplx_lt_doub (a, b) result (bLt)
complex(dp), intent(in) :: a, b
logical :: bLt

bLt = real(a, dp) < real(b, dp)
end function

end module

#define srt_ind(i) (((i)*nskip_)+1)

module sort_mod_sym
use util_mod, only: stop_all
use util_mod_comparisons, only: operator(.arrgt.), operator(.arrlt.)

use constants
use SystemData, only: Symmetry, operator(.eq.), operator(.ne.), &
operator(.gt.), operator(.lt.), assignment(=)
use symdata, only: SymPairProd, operator(.eq.), operator(.ne.), &
operator(.gt.), operator(.lt.), assignment(=)
implicit none

private
public :: sort

! Private operator for sorting complex numbers by their real values
interface operator(.gt.)
module procedure cmplx_gt_sym
end interface

interface operator(.lt.)
module procedure cmplx_lt_sym
end interface

interface sort
module procedure sort_sym
end interface

interface cmplx_gt
module procedure cmplx_gt_sym
end interface

interface cmplx_lt
module procedure cmplx_lt_sym
end interface

contains
pure subroutine sort_sym (arr, nskip, par)

! Perform a quicksort on an array of integers, arr(n). Uses the
! sample code in NumericalRecipies as a base.
! Optionally sort arr2 in parallel (in the routines it is enabled)

! Custom comparisons. Use the operator .locallt. & .localgt.
! interface operator(.locallt.)
!     pure function custom_lt (b, c) result (ret)
!         use constants
!         type(Symmetry), intent(in) :: b(:)
!         type(Symmetry), intent(in) :: c(:)
!         logical :: ret
!     end function
! end interface

! interface operator(.localgt.)
!     pure function custom_gt (b, c) result (ret)
!         use constants
!         type(Symmetry), intent(in) :: b(:)
!         type(Symmetry), intent(in) :: c(:)
!         logical :: ret
!     end function
! end interface

! sort needs auxiliary storage of length 2*log_2(n).
integer, parameter :: nStackMax = 50
integer, parameter :: nInsertionSort = 7
integer, intent(in), optional :: nskip
integer, intent(out), optional :: par

type(Symmetry), intent(inout) :: arr(:)
!, intent(inout) :: arr2(:)
!, intent(inout) :: arr3(:)
!, intent(inout) :: arr4(:)
!, intent(inout) :: arr5(:)
!, intent(inout) :: arr6(:)

! Oh, how lovely it would be to be able to use push/pop and not need
integer :: stack(nStackMax), nstack, nskip_
integer :: pivot, lo, hi, n, i, j, par_
! n.b. This size statement is removed if type1 is scalar ...
type(Symmetry) :: a
! :: a2(size(arr2(1)))
! :: a3(size(arr3(1)))
! :: a4(size(arr4(1)))
! :: a5(size(arr5(1)))
! :: a6(size(arr6(1)))

! Temporary variables for swapping
type(Symmetry) :: tmp1
! :: tmp2(size(arr2(1)))
! :: tmp3(size(arr3(1)))
! :: tmp4(size(arr4(1)))
! :: tmp5(size(arr5(1)))
! :: tmp6(size(arr6(1)))
character(*), parameter :: this_routine = 'sort'

! Initialise temporary variables if required
!(a)
!(tmp1)
!(a2)
!(tmp2)
!(a3)
!(tmp3)
!(a4)
!(tmp4)
!(a5)
!(tmp5)
!(a6)
!(tmp6)

! *** HACK ***
! Workaround for gfortran compiler bug
! n.b. This will produce spurious warnings if run in valgrind, as
!      a is not initialised. Not a problem in optimised build.
i = 0
#ifdef DEBUG_
! it IS a problem in debug build, however
! *** HACK MOD ***
!to prevent crashes in debug mode
if((ubound(arr,1) - lbound(arr,1) + 1) &
> 0 ) a = arr(1)
#endif
! if (custom_lt(a, a)) i = i
! if (custom_gt(a, a)) i = i

if (present(nskip)) then
nskip_ = nskip
else
nskip_ = 1
endif

! Initialise parity calculation
par_ = 1

! The size of the array to sort.
! N.B. Use zero based indices. To obtain ! the entry into the actual
! array, multiply indices by nskip and add ! 1 (hence zero based)
! **** See local macro srt_ind() ******
lo = lbound(arr, 1) - 1
n = ((ubound(arr, 1) - lo -1)/nskip_) + 1
hi = lo + n - 1

nstack = 0
do while (.true.)
! If the section/partition we are looking at is smaller than
! nInsertSort then perform an insertion sort
if (hi - lo < nInsertionSort) then
do j = lo + 1, hi
a = arr(srt_ind(j))
! a2 = arr2(srt_ind(j))
! a3 = arr3(srt_ind(j))
! a4 = arr4(srt_ind(j))
! a5 = arr5(srt_ind(j))
! a6 = arr6(srt_ind(j))
do i = j - 1, 0, -1
if (arr(srt_ind(i)) < a) exit
arr(srt_ind(i+1)) = arr(srt_ind(i))
! arr2(srt_ind(i+1)) = arr2(srt_ind(i))
! arr3(srt_ind(i+1)) = arr3(srt_ind(i))
! arr4(srt_ind(i+1)) = arr4(srt_ind(i))
! arr5(srt_ind(i+1)) = arr5(srt_ind(i))
! arr6(srt_ind(i+1)) = arr6(srt_ind(i))
par_ = -par_
enddo
arr(srt_ind(i+1)) = a
! arr2(srt_ind(i+1)) = a2
! arr3(srt_ind(i+1)) = a3
! arr4(srt_ind(i+1)) = a4
! arr5(srt_ind(i+1)) = a5
! arr6(srt_ind(i+1)) = a6
enddo

if (nstack == 0) exit
hi = stack(nstack)
lo = stack(nstack-1)
nstack = nstack - 2

! Otherwise start partitioning with quicksort.
else
! Pick a partitioning element, and arrange such that
! arr(lo) <= arr(lo+1) <= arr(hi)
pivot = (lo + hi) / 2
tmp1 = arr(srt_ind(pivot))
arr(srt_ind(pivot)) = arr(srt_ind(lo + 1))
arr(srt_ind(lo + 1)) = tmp1
! tmp2 = arr2(srt_ind(pivot))
! arr2(srt_ind(pivot)) = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = tmp2
! tmp3 = arr3(srt_ind(pivot))
! arr3(srt_ind(pivot)) = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = tmp3
! tmp4 = arr4(srt_ind(pivot))
! arr4(srt_ind(pivot)) = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = tmp4
! tmp5 = arr5(srt_ind(pivot))
! arr5(srt_ind(pivot)) = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = tmp5
! tmp6 = arr6(srt_ind(pivot))
! arr6(srt_ind(pivot)) = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = tmp6
par_ = -par_

if (arr(srt_ind(lo)) > arr(srt_ind(hi))) then
tmp1 = arr(srt_ind(lo))
arr(srt_ind(lo)) = arr(srt_ind(hi))
arr(srt_ind(hi)) = tmp1
! tmp2 = arr2(srt_ind(lo))
! arr2(srt_ind(lo)) = arr2(srt_ind(hi))
! arr2(srt_ind(hi)) = tmp2
! tmp3 = arr3(srt_ind(lo))
! arr3(srt_ind(lo)) = arr3(srt_ind(hi))
! arr3(srt_ind(hi)) = tmp3
! tmp4 = arr4(srt_ind(lo))
! arr4(srt_ind(lo)) = arr4(srt_ind(hi))
! arr4(srt_ind(hi)) = tmp4
! tmp5 = arr5(srt_ind(lo))
! arr5(srt_ind(lo)) = arr5(srt_ind(hi))
! arr5(srt_ind(hi)) = tmp5
! tmp6 = arr6(srt_ind(lo))
! arr6(srt_ind(lo)) = arr6(srt_ind(hi))
! arr6(srt_ind(hi)) = tmp6
par_ = -par_
endif
if (arr(srt_ind(lo+1)) > arr(srt_ind(hi))) then
tmp1 = arr(srt_ind(lo+1))
arr(srt_ind(lo+1)) = arr(srt_ind(hi))
arr(srt_ind(hi)) = tmp1
! tmp2 = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = arr2(srt_ind(hi))
! arr2(srt_ind(hi)) = tmp2
! tmp3 = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = arr3(srt_ind(hi))
! arr3(srt_ind(hi)) = tmp3
! tmp4 = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = arr4(srt_ind(hi))
! arr4(srt_ind(hi)) = tmp4
! tmp5 = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = arr5(srt_ind(hi))
! arr5(srt_ind(hi)) = tmp5
! tmp6 = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = arr6(srt_ind(hi))
! arr6(srt_ind(hi)) = tmp6
par_ = -par_
endif
if (arr(srt_ind(lo)) > arr(srt_ind(lo+1))) then
tmp1 = arr(srt_ind(lo))
arr(srt_ind(lo)) = arr(srt_ind(lo+1))
arr(srt_ind(lo+1)) = tmp1
! tmp2 = arr2(srt_ind(lo))
! arr2(srt_ind(lo)) = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = tmp2
! tmp3 = arr3(srt_ind(lo))
! arr3(srt_ind(lo)) = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = tmp3
! tmp4 = arr4(srt_ind(lo))
! arr4(srt_ind(lo)) = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = tmp4
! tmp5 = arr5(srt_ind(lo))
! arr5(srt_ind(lo)) = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = tmp5
! tmp6 = arr6(srt_ind(lo))
! arr6(srt_ind(lo)) = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = tmp6
par_ = -par_
endif

i = lo + 1
j = hi
a = arr(srt_ind(lo + 1)) !! a is the pivot value
! a2 = arr2(srt_ind(lo + 1))
! a3 = arr3(srt_ind(lo + 1))
! a4 = arr4(srt_ind(lo + 1))
! a5 = arr5(srt_ind(lo + 1))
! a6 = arr6(srt_ind(lo + 1))
do while (.true.)
! Scand down list to find element > a
i = i + 1
do while (arr(srt_ind(i)) < a)
i = i + 1
enddo

! Scan down list to find element < a
j = j - 1
do while (arr(srt_ind(j)) > a)
j = j - 1
enddo

! When the pointers crossed, partitioning is complete.
if (j < i) exit

! Swap the elements, so that all elements < a end up
! in lower indexed variables.
tmp1 = arr(srt_ind(i))
arr(srt_ind(i)) = arr(srt_ind(j))
arr(srt_ind(j)) = tmp1
! tmp2 = arr2(srt_ind(i))
! arr2(srt_ind(i)) = arr2(srt_ind(j))
! arr2(srt_ind(j)) = tmp2
! tmp3 = arr3(srt_ind(i))
! arr3(srt_ind(i)) = arr3(srt_ind(j))
! arr3(srt_ind(j)) = tmp3
! tmp4 = arr4(srt_ind(i))
! arr4(srt_ind(i)) = arr4(srt_ind(j))
! arr4(srt_ind(j)) = tmp4
! tmp5 = arr5(srt_ind(i))
! arr5(srt_ind(i)) = arr5(srt_ind(j))
! arr5(srt_ind(j)) = tmp5
! tmp6 = arr6(srt_ind(i))
! arr6(srt_ind(i)) = arr6(srt_ind(j))
! arr6(srt_ind(j)) = tmp6
par_ = -par_
enddo;

! Insert partitioning element
arr(srt_ind(lo + 1)) = arr(srt_ind(j))
arr(srt_ind(j)) = a
! arr2(srt_ind(lo + 1)) = arr2(srt_ind(j))
! arr3(srt_ind(lo + 1)) = arr3(srt_ind(j))
! arr4(srt_ind(lo + 1)) = arr4(srt_ind(j))
! arr5(srt_ind(lo + 1)) = arr5(srt_ind(j))
! arr6(srt_ind(lo + 1)) = arr6(srt_ind(j))
! arr2(srt_ind(j)) = a2
! arr3(srt_ind(j)) = a3
! arr4(srt_ind(j)) = a4
! arr5(srt_ind(j)) = a5
! arr6(srt_ind(j)) = a6
par_ = -par_

! Push the larger of the partitioned sections onto the stack
! of sections to look at later.
! --> need fewest stack elements.
nstack = nstack + 2
if (nstack > nStackMax) then
call stop_all (this_routine, &
"parameter nStackMax too small")
endif
if (hi - i + 1 >= j - lo) then
stack(nstack) = hi
stack(nstack-1) = i
hi = j - 1
else
stack(nstack) = j - 1
stack(nstack-1) = lo
lo = i
endif
endif
enddo

! Destroy temporary variables if required
!(a)
!(tmp1)
!(a2)
!(tmp2)
!(a3)
!(tmp3)
!(a4)
!(tmp4)
!(a5)
!(tmp5)
!(a6)
!(tmp6)

! Return the parity if required
if (present(par)) &
par = par_

end subroutine

! A private comparison. This should not be available outside of the
! module!
elemental function cmplx_gt_sym (a, b) result (bGt)
complex(dp), intent(in) :: a, b
logical :: bGt

bGt = real(a, dp) > real(b, dp)
end function

elemental function cmplx_lt_sym (a, b) result (bLt)
complex(dp), intent(in) :: a, b
logical :: bLt

bLt = real(a, dp) < real(b, dp)
end function

end module

#define srt_ind(i) (((i)*nskip_)+1)

module sort_mod_sympairprod
use util_mod, only: stop_all
use util_mod_comparisons, only: operator(.arrgt.), operator(.arrlt.)

use constants
use SystemData, only: Symmetry, operator(.eq.), operator(.ne.), &
operator(.gt.), operator(.lt.), assignment(=)
use symdata, only: SymPairProd, operator(.eq.), operator(.ne.), &
operator(.gt.), operator(.lt.), assignment(=)
implicit none

private
public :: sort

! Private operator for sorting complex numbers by their real values
interface operator(.gt.)
module procedure cmplx_gt_sympairprod
end interface

interface operator(.lt.)
module procedure cmplx_lt_sympairprod
end interface

interface sort
module procedure sort_sympairprod
end interface

interface cmplx_gt
module procedure cmplx_gt_sympairprod
end interface

interface cmplx_lt
module procedure cmplx_lt_sympairprod
end interface

contains
pure subroutine sort_sympairprod (arr, nskip, par)

! Perform a quicksort on an array of integers, arr(n). Uses the
! sample code in NumericalRecipies as a base.
! Optionally sort arr2 in parallel (in the routines it is enabled)

! Custom comparisons. Use the operator .locallt. & .localgt.
! interface operator(.locallt.)
!     pure function custom_lt (b, c) result (ret)
!         use constants
!         type(SymPairProd), intent(in) :: b(:)
!         type(SymPairProd), intent(in) :: c(:)
!         logical :: ret
!     end function
! end interface

! interface operator(.localgt.)
!     pure function custom_gt (b, c) result (ret)
!         use constants
!         type(SymPairProd), intent(in) :: b(:)
!         type(SymPairProd), intent(in) :: c(:)
!         logical :: ret
!     end function
! end interface

! sort needs auxiliary storage of length 2*log_2(n).
integer, parameter :: nStackMax = 50
integer, parameter :: nInsertionSort = 7
integer, intent(in), optional :: nskip
integer, intent(out), optional :: par

type(SymPairProd), intent(inout) :: arr(:)
!, intent(inout) :: arr2(:)
!, intent(inout) :: arr3(:)
!, intent(inout) :: arr4(:)
!, intent(inout) :: arr5(:)
!, intent(inout) :: arr6(:)

! Oh, how lovely it would be to be able to use push/pop and not need
integer :: stack(nStackMax), nstack, nskip_
integer :: pivot, lo, hi, n, i, j, par_
! n.b. This size statement is removed if type1 is scalar ...
type(SymPairProd) :: a
! :: a2(size(arr2(1)))
! :: a3(size(arr3(1)))
! :: a4(size(arr4(1)))
! :: a5(size(arr5(1)))
! :: a6(size(arr6(1)))

! Temporary variables for swapping
type(SymPairProd) :: tmp1
! :: tmp2(size(arr2(1)))
! :: tmp3(size(arr3(1)))
! :: tmp4(size(arr4(1)))
! :: tmp5(size(arr5(1)))
! :: tmp6(size(arr6(1)))
character(*), parameter :: this_routine = 'sort'

! Initialise temporary variables if required
!(a)
!(tmp1)
!(a2)
!(tmp2)
!(a3)
!(tmp3)
!(a4)
!(tmp4)
!(a5)
!(tmp5)
!(a6)
!(tmp6)

! *** HACK ***
! Workaround for gfortran compiler bug
! n.b. This will produce spurious warnings if run in valgrind, as
!      a is not initialised. Not a problem in optimised build.
i = 0
#ifdef DEBUG_
! it IS a problem in debug build, however
! *** HACK MOD ***
!to prevent crashes in debug mode
if((ubound(arr,1) - lbound(arr,1) + 1) &
> 0 ) a = arr(1)
#endif
! if (custom_lt(a, a)) i = i
! if (custom_gt(a, a)) i = i

if (present(nskip)) then
nskip_ = nskip
else
nskip_ = 1
endif

! Initialise parity calculation
par_ = 1

! The size of the array to sort.
! N.B. Use zero based indices. To obtain ! the entry into the actual
! array, multiply indices by nskip and add ! 1 (hence zero based)
! **** See local macro srt_ind() ******
lo = lbound(arr, 1) - 1
n = ((ubound(arr, 1) - lo -1)/nskip_) + 1
hi = lo + n - 1

nstack = 0
do while (.true.)
! If the section/partition we are looking at is smaller than
! nInsertSort then perform an insertion sort
if (hi - lo < nInsertionSort) then
do j = lo + 1, hi
a = arr(srt_ind(j))
! a2 = arr2(srt_ind(j))
! a3 = arr3(srt_ind(j))
! a4 = arr4(srt_ind(j))
! a5 = arr5(srt_ind(j))
! a6 = arr6(srt_ind(j))
do i = j - 1, 0, -1
if (arr(srt_ind(i)) < a) exit
arr(srt_ind(i+1)) = arr(srt_ind(i))
! arr2(srt_ind(i+1)) = arr2(srt_ind(i))
! arr3(srt_ind(i+1)) = arr3(srt_ind(i))
! arr4(srt_ind(i+1)) = arr4(srt_ind(i))
! arr5(srt_ind(i+1)) = arr5(srt_ind(i))
! arr6(srt_ind(i+1)) = arr6(srt_ind(i))
par_ = -par_
enddo
arr(srt_ind(i+1)) = a
! arr2(srt_ind(i+1)) = a2
! arr3(srt_ind(i+1)) = a3
! arr4(srt_ind(i+1)) = a4
! arr5(srt_ind(i+1)) = a5
! arr6(srt_ind(i+1)) = a6
enddo

if (nstack == 0) exit
hi = stack(nstack)
lo = stack(nstack-1)
nstack = nstack - 2

! Otherwise start partitioning with quicksort.
else
! Pick a partitioning element, and arrange such that
! arr(lo) <= arr(lo+1) <= arr(hi)
pivot = (lo + hi) / 2
tmp1 = arr(srt_ind(pivot))
arr(srt_ind(pivot)) = arr(srt_ind(lo + 1))
arr(srt_ind(lo + 1)) = tmp1
! tmp2 = arr2(srt_ind(pivot))
! arr2(srt_ind(pivot)) = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = tmp2
! tmp3 = arr3(srt_ind(pivot))
! arr3(srt_ind(pivot)) = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = tmp3
! tmp4 = arr4(srt_ind(pivot))
! arr4(srt_ind(pivot)) = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = tmp4
! tmp5 = arr5(srt_ind(pivot))
! arr5(srt_ind(pivot)) = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = tmp5
! tmp6 = arr6(srt_ind(pivot))
! arr6(srt_ind(pivot)) = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = tmp6
par_ = -par_

if (arr(srt_ind(lo)) > arr(srt_ind(hi))) then
tmp1 = arr(srt_ind(lo))
arr(srt_ind(lo)) = arr(srt_ind(hi))
arr(srt_ind(hi)) = tmp1
! tmp2 = arr2(srt_ind(lo))
! arr2(srt_ind(lo)) = arr2(srt_ind(hi))
! arr2(srt_ind(hi)) = tmp2
! tmp3 = arr3(srt_ind(lo))
! arr3(srt_ind(lo)) = arr3(srt_ind(hi))
! arr3(srt_ind(hi)) = tmp3
! tmp4 = arr4(srt_ind(lo))
! arr4(srt_ind(lo)) = arr4(srt_ind(hi))
! arr4(srt_ind(hi)) = tmp4
! tmp5 = arr5(srt_ind(lo))
! arr5(srt_ind(lo)) = arr5(srt_ind(hi))
! arr5(srt_ind(hi)) = tmp5
! tmp6 = arr6(srt_ind(lo))
! arr6(srt_ind(lo)) = arr6(srt_ind(hi))
! arr6(srt_ind(hi)) = tmp6
par_ = -par_
endif
if (arr(srt_ind(lo+1)) > arr(srt_ind(hi))) then
tmp1 = arr(srt_ind(lo+1))
arr(srt_ind(lo+1)) = arr(srt_ind(hi))
arr(srt_ind(hi)) = tmp1
! tmp2 = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = arr2(srt_ind(hi))
! arr2(srt_ind(hi)) = tmp2
! tmp3 = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = arr3(srt_ind(hi))
! arr3(srt_ind(hi)) = tmp3
! tmp4 = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = arr4(srt_ind(hi))
! arr4(srt_ind(hi)) = tmp4
! tmp5 = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = arr5(srt_ind(hi))
! arr5(srt_ind(hi)) = tmp5
! tmp6 = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = arr6(srt_ind(hi))
! arr6(srt_ind(hi)) = tmp6
par_ = -par_
endif
if (arr(srt_ind(lo)) > arr(srt_ind(lo+1))) then
tmp1 = arr(srt_ind(lo))
arr(srt_ind(lo)) = arr(srt_ind(lo+1))
arr(srt_ind(lo+1)) = tmp1
! tmp2 = arr2(srt_ind(lo))
! arr2(srt_ind(lo)) = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = tmp2
! tmp3 = arr3(srt_ind(lo))
! arr3(srt_ind(lo)) = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = tmp3
! tmp4 = arr4(srt_ind(lo))
! arr4(srt_ind(lo)) = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = tmp4
! tmp5 = arr5(srt_ind(lo))
! arr5(srt_ind(lo)) = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = tmp5
! tmp6 = arr6(srt_ind(lo))
! arr6(srt_ind(lo)) = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = tmp6
par_ = -par_
endif

i = lo + 1
j = hi
a = arr(srt_ind(lo + 1)) !! a is the pivot value
! a2 = arr2(srt_ind(lo + 1))
! a3 = arr3(srt_ind(lo + 1))
! a4 = arr4(srt_ind(lo + 1))
! a5 = arr5(srt_ind(lo + 1))
! a6 = arr6(srt_ind(lo + 1))
do while (.true.)
! Scand down list to find element > a
i = i + 1
do while (arr(srt_ind(i)) < a)
i = i + 1
enddo

! Scan down list to find element < a
j = j - 1
do while (arr(srt_ind(j)) > a)
j = j - 1
enddo

! When the pointers crossed, partitioning is complete.
if (j < i) exit

! Swap the elements, so that all elements < a end up
! in lower indexed variables.
tmp1 = arr(srt_ind(i))
arr(srt_ind(i)) = arr(srt_ind(j))
arr(srt_ind(j)) = tmp1
! tmp2 = arr2(srt_ind(i))
! arr2(srt_ind(i)) = arr2(srt_ind(j))
! arr2(srt_ind(j)) = tmp2
! tmp3 = arr3(srt_ind(i))
! arr3(srt_ind(i)) = arr3(srt_ind(j))
! arr3(srt_ind(j)) = tmp3
! tmp4 = arr4(srt_ind(i))
! arr4(srt_ind(i)) = arr4(srt_ind(j))
! arr4(srt_ind(j)) = tmp4
! tmp5 = arr5(srt_ind(i))
! arr5(srt_ind(i)) = arr5(srt_ind(j))
! arr5(srt_ind(j)) = tmp5
! tmp6 = arr6(srt_ind(i))
! arr6(srt_ind(i)) = arr6(srt_ind(j))
! arr6(srt_ind(j)) = tmp6
par_ = -par_
enddo;

! Insert partitioning element
arr(srt_ind(lo + 1)) = arr(srt_ind(j))
arr(srt_ind(j)) = a
! arr2(srt_ind(lo + 1)) = arr2(srt_ind(j))
! arr3(srt_ind(lo + 1)) = arr3(srt_ind(j))
! arr4(srt_ind(lo + 1)) = arr4(srt_ind(j))
! arr5(srt_ind(lo + 1)) = arr5(srt_ind(j))
! arr6(srt_ind(lo + 1)) = arr6(srt_ind(j))
! arr2(srt_ind(j)) = a2
! arr3(srt_ind(j)) = a3
! arr4(srt_ind(j)) = a4
! arr5(srt_ind(j)) = a5
! arr6(srt_ind(j)) = a6
par_ = -par_

! Push the larger of the partitioned sections onto the stack
! of sections to look at later.
! --> need fewest stack elements.
nstack = nstack + 2
if (nstack > nStackMax) then
call stop_all (this_routine, &
"parameter nStackMax too small")
endif
if (hi - i + 1 >= j - lo) then
stack(nstack) = hi
stack(nstack-1) = i
hi = j - 1
else
stack(nstack) = j - 1
stack(nstack-1) = lo
lo = i
endif
endif
enddo

! Destroy temporary variables if required
!(a)
!(tmp1)
!(a2)
!(tmp2)
!(a3)
!(tmp3)
!(a4)
!(tmp4)
!(a5)
!(tmp5)
!(a6)
!(tmp6)

! Return the parity if required
if (present(par)) &
par = par_

end subroutine

! A private comparison. This should not be available outside of the
! module!
elemental function cmplx_gt_sympairprod (a, b) result (bGt)
complex(dp), intent(in) :: a, b
logical :: bGt

bGt = real(a, dp) > real(b, dp)
end function

elemental function cmplx_lt_sympairprod (a, b) result (bLt)
complex(dp), intent(in) :: a, b
logical :: bLt

bLt = real(a, dp) < real(b, dp)
end function

end module

#define srt_ind(i) (((i)*nskip_)+1)

module sort_mod_cmplx
use util_mod, only: stop_all
use util_mod_comparisons, only: operator(.arrgt.), operator(.arrlt.)

use constants
use SystemData, only: Symmetry, operator(.eq.), operator(.ne.), &
operator(.gt.), operator(.lt.), assignment(=)
use symdata, only: SymPairProd, operator(.eq.), operator(.ne.), &
operator(.gt.), operator(.lt.), assignment(=)
implicit none

private
public :: sort

! Private operator for sorting complex numbers by their real values
interface operator(.gt.)
module procedure cmplx_gt_cmplx
end interface

interface operator(.lt.)
module procedure cmplx_lt_cmplx
end interface

interface sort
module procedure sort_cmplx
end interface

interface cmplx_gt
module procedure cmplx_gt_cmplx
end interface

interface cmplx_lt
module procedure cmplx_lt_cmplx
end interface

contains
pure subroutine sort_cmplx (arr, nskip, par)

! Perform a quicksort on an array of integers, arr(n). Uses the
! sample code in NumericalRecipies as a base.
! Optionally sort arr2 in parallel (in the routines it is enabled)

! Custom comparisons. Use the operator .locallt. & .localgt.
! interface operator(.locallt.)
!     pure function custom_lt (b, c) result (ret)
!         use constants
!         complex(dp), intent(in) :: b(:)
!         complex(dp), intent(in) :: c(:)
!         logical :: ret
!     end function
! end interface

! interface operator(.localgt.)
!     pure function custom_gt (b, c) result (ret)
!         use constants
!         complex(dp), intent(in) :: b(:)
!         complex(dp), intent(in) :: c(:)
!         logical :: ret
!     end function
! end interface

! sort needs auxiliary storage of length 2*log_2(n).
integer, parameter :: nStackMax = 50
integer, parameter :: nInsertionSort = 7
integer, intent(in), optional :: nskip
integer, intent(out), optional :: par

complex(dp), intent(inout) :: arr(:)
!, intent(inout) :: arr2(:)
!, intent(inout) :: arr3(:)
!, intent(inout) :: arr4(:)
!, intent(inout) :: arr5(:)
!, intent(inout) :: arr6(:)

! Oh, how lovely it would be to be able to use push/pop and not need
integer :: stack(nStackMax), nstack, nskip_
integer :: pivot, lo, hi, n, i, j, par_
! n.b. This size statement is removed if type1 is scalar ...
complex(dp) :: a
! :: a2(size(arr2(1)))
! :: a3(size(arr3(1)))
! :: a4(size(arr4(1)))
! :: a5(size(arr5(1)))
! :: a6(size(arr6(1)))

! Temporary variables for swapping
complex(dp) :: tmp1
! :: tmp2(size(arr2(1)))
! :: tmp3(size(arr3(1)))
! :: tmp4(size(arr4(1)))
! :: tmp5(size(arr5(1)))
! :: tmp6(size(arr6(1)))
character(*), parameter :: this_routine = 'sort'

! Initialise temporary variables if required
!(a)
!(tmp1)
!(a2)
!(tmp2)
!(a3)
!(tmp3)
!(a4)
!(tmp4)
!(a5)
!(tmp5)
!(a6)
!(tmp6)

! *** HACK ***
! Workaround for gfortran compiler bug
! n.b. This will produce spurious warnings if run in valgrind, as
!      a is not initialised. Not a problem in optimised build.
i = 0
#ifdef DEBUG_
! it IS a problem in debug build, however
! *** HACK MOD ***
!to prevent crashes in debug mode
if((ubound(arr,1) - lbound(arr,1) + 1) &
> 0 ) a = arr(1)
#endif
! if (custom_lt(a, a)) i = i
! if (custom_gt(a, a)) i = i

if (present(nskip)) then
nskip_ = nskip
else
nskip_ = 1
endif

! Initialise parity calculation
par_ = 1

! The size of the array to sort.
! N.B. Use zero based indices. To obtain ! the entry into the actual
! array, multiply indices by nskip and add ! 1 (hence zero based)
! **** See local macro srt_ind() ******
lo = lbound(arr, 1) - 1
n = ((ubound(arr, 1) - lo -1)/nskip_) + 1
hi = lo + n - 1

nstack = 0
do while (.true.)
! If the section/partition we are looking at is smaller than
! nInsertSort then perform an insertion sort
if (hi - lo < nInsertionSort) then
do j = lo + 1, hi
a = arr(srt_ind(j))
! a2 = arr2(srt_ind(j))
! a3 = arr3(srt_ind(j))
! a4 = arr4(srt_ind(j))
! a5 = arr5(srt_ind(j))
! a6 = arr6(srt_ind(j))
do i = j - 1, 0, -1
if (arr(srt_ind(i)) < a) exit
arr(srt_ind(i+1)) = arr(srt_ind(i))
! arr2(srt_ind(i+1)) = arr2(srt_ind(i))
! arr3(srt_ind(i+1)) = arr3(srt_ind(i))
! arr4(srt_ind(i+1)) = arr4(srt_ind(i))
! arr5(srt_ind(i+1)) = arr5(srt_ind(i))
! arr6(srt_ind(i+1)) = arr6(srt_ind(i))
par_ = -par_
enddo
arr(srt_ind(i+1)) = a
! arr2(srt_ind(i+1)) = a2
! arr3(srt_ind(i+1)) = a3
! arr4(srt_ind(i+1)) = a4
! arr5(srt_ind(i+1)) = a5
! arr6(srt_ind(i+1)) = a6
enddo

if (nstack == 0) exit
hi = stack(nstack)
lo = stack(nstack-1)
nstack = nstack - 2

! Otherwise start partitioning with quicksort.
else
! Pick a partitioning element, and arrange such that
! arr(lo) <= arr(lo+1) <= arr(hi)
pivot = (lo + hi) / 2
tmp1 = arr(srt_ind(pivot))
arr(srt_ind(pivot)) = arr(srt_ind(lo + 1))
arr(srt_ind(lo + 1)) = tmp1
! tmp2 = arr2(srt_ind(pivot))
! arr2(srt_ind(pivot)) = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = tmp2
! tmp3 = arr3(srt_ind(pivot))
! arr3(srt_ind(pivot)) = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = tmp3
! tmp4 = arr4(srt_ind(pivot))
! arr4(srt_ind(pivot)) = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = tmp4
! tmp5 = arr5(srt_ind(pivot))
! arr5(srt_ind(pivot)) = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = tmp5
! tmp6 = arr6(srt_ind(pivot))
! arr6(srt_ind(pivot)) = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = tmp6
par_ = -par_

if (arr(srt_ind(lo)) > arr(srt_ind(hi))) then
tmp1 = arr(srt_ind(lo))
arr(srt_ind(lo)) = arr(srt_ind(hi))
arr(srt_ind(hi)) = tmp1
! tmp2 = arr2(srt_ind(lo))
! arr2(srt_ind(lo)) = arr2(srt_ind(hi))
! arr2(srt_ind(hi)) = tmp2
! tmp3 = arr3(srt_ind(lo))
! arr3(srt_ind(lo)) = arr3(srt_ind(hi))
! arr3(srt_ind(hi)) = tmp3
! tmp4 = arr4(srt_ind(lo))
! arr4(srt_ind(lo)) = arr4(srt_ind(hi))
! arr4(srt_ind(hi)) = tmp4
! tmp5 = arr5(srt_ind(lo))
! arr5(srt_ind(lo)) = arr5(srt_ind(hi))
! arr5(srt_ind(hi)) = tmp5
! tmp6 = arr6(srt_ind(lo))
! arr6(srt_ind(lo)) = arr6(srt_ind(hi))
! arr6(srt_ind(hi)) = tmp6
par_ = -par_
endif
if (arr(srt_ind(lo+1)) > arr(srt_ind(hi))) then
tmp1 = arr(srt_ind(lo+1))
arr(srt_ind(lo+1)) = arr(srt_ind(hi))
arr(srt_ind(hi)) = tmp1
! tmp2 = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = arr2(srt_ind(hi))
! arr2(srt_ind(hi)) = tmp2
! tmp3 = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = arr3(srt_ind(hi))
! arr3(srt_ind(hi)) = tmp3
! tmp4 = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = arr4(srt_ind(hi))
! arr4(srt_ind(hi)) = tmp4
! tmp5 = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = arr5(srt_ind(hi))
! arr5(srt_ind(hi)) = tmp5
! tmp6 = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = arr6(srt_ind(hi))
! arr6(srt_ind(hi)) = tmp6
par_ = -par_
endif
if (arr(srt_ind(lo)) > arr(srt_ind(lo+1))) then
tmp1 = arr(srt_ind(lo))
arr(srt_ind(lo)) = arr(srt_ind(lo+1))
arr(srt_ind(lo+1)) = tmp1
! tmp2 = arr2(srt_ind(lo))
! arr2(srt_ind(lo)) = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = tmp2
! tmp3 = arr3(srt_ind(lo))
! arr3(srt_ind(lo)) = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = tmp3
! tmp4 = arr4(srt_ind(lo))
! arr4(srt_ind(lo)) = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = tmp4
! tmp5 = arr5(srt_ind(lo))
! arr5(srt_ind(lo)) = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = tmp5
! tmp6 = arr6(srt_ind(lo))
! arr6(srt_ind(lo)) = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = tmp6
par_ = -par_
endif

i = lo + 1
j = hi
a = arr(srt_ind(lo + 1)) !! a is the pivot value
! a2 = arr2(srt_ind(lo + 1))
! a3 = arr3(srt_ind(lo + 1))
! a4 = arr4(srt_ind(lo + 1))
! a5 = arr5(srt_ind(lo + 1))
! a6 = arr6(srt_ind(lo + 1))
do while (.true.)
! Scand down list to find element > a
i = i + 1
do while (arr(srt_ind(i)) < a)
i = i + 1
enddo

! Scan down list to find element < a
j = j - 1
do while (arr(srt_ind(j)) > a)
j = j - 1
enddo

! When the pointers crossed, partitioning is complete.
if (j < i) exit

! Swap the elements, so that all elements < a end up
! in lower indexed variables.
tmp1 = arr(srt_ind(i))
arr(srt_ind(i)) = arr(srt_ind(j))
arr(srt_ind(j)) = tmp1
! tmp2 = arr2(srt_ind(i))
! arr2(srt_ind(i)) = arr2(srt_ind(j))
! arr2(srt_ind(j)) = tmp2
! tmp3 = arr3(srt_ind(i))
! arr3(srt_ind(i)) = arr3(srt_ind(j))
! arr3(srt_ind(j)) = tmp3
! tmp4 = arr4(srt_ind(i))
! arr4(srt_ind(i)) = arr4(srt_ind(j))
! arr4(srt_ind(j)) = tmp4
! tmp5 = arr5(srt_ind(i))
! arr5(srt_ind(i)) = arr5(srt_ind(j))
! arr5(srt_ind(j)) = tmp5
! tmp6 = arr6(srt_ind(i))
! arr6(srt_ind(i)) = arr6(srt_ind(j))
! arr6(srt_ind(j)) = tmp6
par_ = -par_
enddo;

! Insert partitioning element
arr(srt_ind(lo + 1)) = arr(srt_ind(j))
arr(srt_ind(j)) = a
! arr2(srt_ind(lo + 1)) = arr2(srt_ind(j))
! arr3(srt_ind(lo + 1)) = arr3(srt_ind(j))
! arr4(srt_ind(lo + 1)) = arr4(srt_ind(j))
! arr5(srt_ind(lo + 1)) = arr5(srt_ind(j))
! arr6(srt_ind(lo + 1)) = arr6(srt_ind(j))
! arr2(srt_ind(j)) = a2
! arr3(srt_ind(j)) = a3
! arr4(srt_ind(j)) = a4
! arr5(srt_ind(j)) = a5
! arr6(srt_ind(j)) = a6
par_ = -par_

! Push the larger of the partitioned sections onto the stack
! of sections to look at later.
! --> need fewest stack elements.
nstack = nstack + 2
if (nstack > nStackMax) then
call stop_all (this_routine, &
"parameter nStackMax too small")
endif
if (hi - i + 1 >= j - lo) then
stack(nstack) = hi
stack(nstack-1) = i
hi = j - 1
else
stack(nstack) = j - 1
stack(nstack-1) = lo
lo = i
endif
endif
enddo

! Destroy temporary variables if required
!(a)
!(tmp1)
!(a2)
!(tmp2)
!(a3)
!(tmp3)
!(a4)
!(tmp4)
!(a5)
!(tmp5)
!(a6)
!(tmp6)

! Return the parity if required
if (present(par)) &
par = par_

end subroutine

! A private comparison. This should not be available outside of the
! module!
elemental function cmplx_gt_cmplx (a, b) result (bGt)
complex(dp), intent(in) :: a, b
logical :: bGt

bGt = real(a, dp) > real(b, dp)
end function

elemental function cmplx_lt_cmplx (a, b) result (bLt)
complex(dp), intent(in) :: a, b
logical :: bLt

bLt = real(a, dp) < real(b, dp)
end function

end module

#if !defined(SX)

#define srt_ind(i) (((i)*nskip_)+1)

module sort_mod_a_i
use util_mod, only: stop_all
use util_mod_comparisons, only: operator(.arrgt.), operator(.arrlt.)

use constants
use SystemData, only: Symmetry, operator(.eq.), operator(.ne.), &
operator(.gt.), operator(.lt.), assignment(=)
use symdata, only: SymPairProd, operator(.eq.), operator(.ne.), &
operator(.gt.), operator(.lt.), assignment(=)
implicit none

private
public :: sort

! Private operator for sorting complex numbers by their real values
interface operator(.gt.)
module procedure cmplx_gt_a_i
end interface

interface operator(.lt.)
module procedure cmplx_lt_a_i
end interface

interface sort
module procedure sort_a_i
end interface

interface cmplx_gt
module procedure cmplx_gt_a_i
end interface

interface cmplx_lt
module procedure cmplx_lt_a_i
end interface

contains
pure subroutine sort_a_i (arr, nskip, par)

! Perform a quicksort on an array of integers, arr(n). Uses the
! sample code in NumericalRecipies as a base.
! Optionally sort arr2 in parallel (in the routines it is enabled)

! Custom comparisons. Use the operator .locallt. & .localgt.
! interface operator(.locallt.)
!     pure function custom_lt (b, c) result (ret)
!         use constants
!         integer(kind=int32), intent(in) :: b(:)
!         integer(kind=int32), intent(in) :: c(:)
!         logical :: ret
!     end function
! end interface
! interface operator(.localgt.)
!     pure function custom_gt (b, c) result (ret)
!         use constants
!         integer(kind=int32), intent(in) :: b(:)
!         integer(kind=int32), intent(in) :: c(:)
!         logical :: ret
!     end function
! end interface

! sort needs auxiliary storage of length 2*log_2(n).
integer, parameter :: nStackMax = 50
integer, parameter :: nInsertionSort = 7
integer, intent(in), optional :: nskip
integer, intent(out), optional :: par

integer(kind=int32), intent(inout) :: arr(:,:)
!, intent(inout) :: arr2(:)
!, intent(inout) :: arr3(:)
!, intent(inout) :: arr4(:)
!, intent(inout) :: arr5(:)
!, intent(inout) :: arr6(:)

! Oh, how lovely it would be to be able to use push/pop and not need
integer :: stack(nStackMax), nstack, nskip_
integer :: pivot, lo, hi, n, i, j, par_
! n.b. This size statement is removed if type1 is scalar ...
integer(kind=int32) :: a(size(arr(:,1)))
! :: a2(size(arr2(1)))
! :: a3(size(arr3(1)))
! :: a4(size(arr4(1)))
! :: a5(size(arr5(1)))
! :: a6(size(arr6(1)))

! Temporary variables for swapping
integer(kind=int32) :: tmp1(size(arr(:,1)))
! :: tmp2(size(arr2(1)))
! :: tmp3(size(arr3(1)))
! :: tmp4(size(arr4(1)))
! :: tmp5(size(arr5(1)))
! :: tmp6(size(arr6(1)))
character(*), parameter :: this_routine = 'sort'

! Initialise temporary variables if required
!(a)
!(tmp1)
!(a2)
!(tmp2)
!(a3)
!(tmp3)
!(a4)
!(tmp4)
!(a5)
!(tmp5)
!(a6)
!(tmp6)

! *** HACK ***
! Workaround for gfortran compiler bug
! n.b. This will produce spurious warnings if run in valgrind, as
!      a is not initialised. Not a problem in optimised build.
i = 0
#ifdef DEBUG_
! it IS a problem in debug build, however
! *** HACK MOD ***
!to prevent crashes in debug mode
if((ubound(arr,2) - lbound(arr,2) + 1) &
> 0 ) a = arr(:,1)
#endif
! if (custom_lt(a, a)) i = i
! if (custom_gt(a, a)) i = i

if (present(nskip)) then
nskip_ = nskip
else
nskip_ = 1
endif

! Initialise parity calculation
par_ = 1

! The size of the array to sort.
! N.B. Use zero based indices. To obtain ! the entry into the actual
! array, multiply indices by nskip and add ! 1 (hence zero based)
! **** See local macro srt_ind() ******
lo = lbound(arr, 2) - 1
n = ((ubound(arr, 2) - lo -1)/nskip_) + 1
hi = lo + n - 1

nstack = 0
do while (.true.)
! If the section/partition we are looking at is smaller than
! nInsertSort then perform an insertion sort
if (hi - lo < nInsertionSort) then
do j = lo + 1, hi
a = arr(:,srt_ind(j))
! a2 = arr2(srt_ind(j))
! a3 = arr3(srt_ind(j))
! a4 = arr4(srt_ind(j))
! a5 = arr5(srt_ind(j))
! a6 = arr6(srt_ind(j))
do i = j - 1, 0, -1
if (arr(:,srt_ind(i)) .arrlt. a) exit
arr(:,srt_ind(i+1)) = arr(:,srt_ind(i))
! arr2(srt_ind(i+1)) = arr2(srt_ind(i))
! arr3(srt_ind(i+1)) = arr3(srt_ind(i))
! arr4(srt_ind(i+1)) = arr4(srt_ind(i))
! arr5(srt_ind(i+1)) = arr5(srt_ind(i))
! arr6(srt_ind(i+1)) = arr6(srt_ind(i))
par_ = -par_
enddo
arr(:,srt_ind(i+1)) = a
! arr2(srt_ind(i+1)) = a2
! arr3(srt_ind(i+1)) = a3
! arr4(srt_ind(i+1)) = a4
! arr5(srt_ind(i+1)) = a5
! arr6(srt_ind(i+1)) = a6
enddo

if (nstack == 0) exit
hi = stack(nstack)
lo = stack(nstack-1)
nstack = nstack - 2

! Otherwise start partitioning with quicksort.
else
! Pick a partitioning element, and arrange such that
! arr(:,lo) <= arr(:,lo+1) <= arr(:,hi)
pivot = (lo + hi) / 2
tmp1 = arr(:,srt_ind(pivot))
arr(:,srt_ind(pivot)) = arr(:,srt_ind(lo + 1))
arr(:,srt_ind(lo + 1)) = tmp1
! tmp2 = arr2(srt_ind(pivot))
! arr2(srt_ind(pivot)) = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = tmp2
! tmp3 = arr3(srt_ind(pivot))
! arr3(srt_ind(pivot)) = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = tmp3
! tmp4 = arr4(srt_ind(pivot))
! arr4(srt_ind(pivot)) = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = tmp4
! tmp5 = arr5(srt_ind(pivot))
! arr5(srt_ind(pivot)) = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = tmp5
! tmp6 = arr6(srt_ind(pivot))
! arr6(srt_ind(pivot)) = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = tmp6
par_ = -par_

if (arr(:,srt_ind(lo)) .arrgt. arr(:,srt_ind(hi))) then
tmp1 = arr(:,srt_ind(lo))
arr(:,srt_ind(lo)) = arr(:,srt_ind(hi))
arr(:,srt_ind(hi)) = tmp1
! tmp2 = arr2(srt_ind(lo))
! arr2(srt_ind(lo)) = arr2(srt_ind(hi))
! arr2(srt_ind(hi)) = tmp2
! tmp3 = arr3(srt_ind(lo))
! arr3(srt_ind(lo)) = arr3(srt_ind(hi))
! arr3(srt_ind(hi)) = tmp3
! tmp4 = arr4(srt_ind(lo))
! arr4(srt_ind(lo)) = arr4(srt_ind(hi))
! arr4(srt_ind(hi)) = tmp4
! tmp5 = arr5(srt_ind(lo))
! arr5(srt_ind(lo)) = arr5(srt_ind(hi))
! arr5(srt_ind(hi)) = tmp5
! tmp6 = arr6(srt_ind(lo))
! arr6(srt_ind(lo)) = arr6(srt_ind(hi))
! arr6(srt_ind(hi)) = tmp6
par_ = -par_
endif
if (arr(:,srt_ind(lo+1)) .arrgt. arr(:,srt_ind(hi))) then
tmp1 = arr(:,srt_ind(lo+1))
arr(:,srt_ind(lo+1)) = arr(:,srt_ind(hi))
arr(:,srt_ind(hi)) = tmp1
! tmp2 = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = arr2(srt_ind(hi))
! arr2(srt_ind(hi)) = tmp2
! tmp3 = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = arr3(srt_ind(hi))
! arr3(srt_ind(hi)) = tmp3
! tmp4 = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = arr4(srt_ind(hi))
! arr4(srt_ind(hi)) = tmp4
! tmp5 = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = arr5(srt_ind(hi))
! arr5(srt_ind(hi)) = tmp5
! tmp6 = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = arr6(srt_ind(hi))
! arr6(srt_ind(hi)) = tmp6
par_ = -par_
endif
if (arr(:,srt_ind(lo)) .arrgt. arr(:,srt_ind(lo+1))) then
tmp1 = arr(:,srt_ind(lo))
arr(:,srt_ind(lo)) = arr(:,srt_ind(lo+1))
arr(:,srt_ind(lo+1)) = tmp1
! tmp2 = arr2(srt_ind(lo))
! arr2(srt_ind(lo)) = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = tmp2
! tmp3 = arr3(srt_ind(lo))
! arr3(srt_ind(lo)) = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = tmp3
! tmp4 = arr4(srt_ind(lo))
! arr4(srt_ind(lo)) = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = tmp4
! tmp5 = arr5(srt_ind(lo))
! arr5(srt_ind(lo)) = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = tmp5
! tmp6 = arr6(srt_ind(lo))
! arr6(srt_ind(lo)) = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = tmp6
par_ = -par_
endif

i = lo + 1
j = hi
a = arr(:,srt_ind(lo + 1)) !! a is the pivot value
! a2 = arr2(srt_ind(lo + 1))
! a3 = arr3(srt_ind(lo + 1))
! a4 = arr4(srt_ind(lo + 1))
! a5 = arr5(srt_ind(lo + 1))
! a6 = arr6(srt_ind(lo + 1))
do while (.true.)
! Scand down list to find element > a
i = i + 1
do while (arr(:,srt_ind(i)) .arrlt. a)
i = i + 1
enddo

! Scan down list to find element < a
j = j - 1
do while (arr(:,srt_ind(j)) .arrgt. a)
j = j - 1
enddo

! When the pointers crossed, partitioning is complete.
if (j < i) exit

! Swap the elements, so that all elements < a end up
! in lower indexed variables.
tmp1 = arr(:,srt_ind(i))
arr(:,srt_ind(i)) = arr(:,srt_ind(j))
arr(:,srt_ind(j)) = tmp1
! tmp2 = arr2(srt_ind(i))
! arr2(srt_ind(i)) = arr2(srt_ind(j))
! arr2(srt_ind(j)) = tmp2
! tmp3 = arr3(srt_ind(i))
! arr3(srt_ind(i)) = arr3(srt_ind(j))
! arr3(srt_ind(j)) = tmp3
! tmp4 = arr4(srt_ind(i))
! arr4(srt_ind(i)) = arr4(srt_ind(j))
! arr4(srt_ind(j)) = tmp4
! tmp5 = arr5(srt_ind(i))
! arr5(srt_ind(i)) = arr5(srt_ind(j))
! arr5(srt_ind(j)) = tmp5
! tmp6 = arr6(srt_ind(i))
! arr6(srt_ind(i)) = arr6(srt_ind(j))
! arr6(srt_ind(j)) = tmp6
par_ = -par_
enddo;

! Insert partitioning element
arr(:,srt_ind(lo + 1)) = arr(:,srt_ind(j))
arr(:,srt_ind(j)) = a
! arr2(srt_ind(lo + 1)) = arr2(srt_ind(j))
! arr3(srt_ind(lo + 1)) = arr3(srt_ind(j))
! arr4(srt_ind(lo + 1)) = arr4(srt_ind(j))
! arr5(srt_ind(lo + 1)) = arr5(srt_ind(j))
! arr6(srt_ind(lo + 1)) = arr6(srt_ind(j))
! arr2(srt_ind(j)) = a2
! arr3(srt_ind(j)) = a3
! arr4(srt_ind(j)) = a4
! arr5(srt_ind(j)) = a5
! arr6(srt_ind(j)) = a6
par_ = -par_

! Push the larger of the partitioned sections onto the stack
! of sections to look at later.
! --> need fewest stack elements.
nstack = nstack + 2
if (nstack > nStackMax) then
call stop_all (this_routine, &
"parameter nStackMax too small")
endif
if (hi - i + 1 >= j - lo) then
stack(nstack) = hi
stack(nstack-1) = i
hi = j - 1
else
stack(nstack) = j - 1
stack(nstack-1) = lo
lo = i
endif
endif
enddo

! Destroy temporary variables if required
!(a)
!(tmp1)
!(a2)
!(tmp2)
!(a3)
!(tmp3)
!(a4)
!(tmp4)
!(a5)
!(tmp5)
!(a6)
!(tmp6)

! Return the parity if required
if (present(par)) &
par = par_

end subroutine

! A private comparison. This should not be available outside of the
! module!
elemental function cmplx_gt_a_i (a, b) result (bGt)
complex(dp), intent(in) :: a, b
logical :: bGt

bGt = real(a, dp) > real(b, dp)
end function

elemental function cmplx_lt_a_i (a, b) result (bLt)
complex(dp), intent(in) :: a, b
logical :: bLt

bLt = real(a, dp) < real(b, dp)
end function

end module

#endif

#define srt_ind(i) (((i)*nskip_)+1)

module sort_mod_a_i64
use util_mod, only: stop_all
use util_mod_comparisons, only: operator(.arrgt.), operator(.arrlt.)

use constants
use SystemData, only: Symmetry, operator(.eq.), operator(.ne.), &
operator(.gt.), operator(.lt.), assignment(=)
use symdata, only: SymPairProd, operator(.eq.), operator(.ne.), &
operator(.gt.), operator(.lt.), assignment(=)
implicit none

private
public :: sort

! Private operator for sorting complex numbers by their real values
interface operator(.gt.)
module procedure cmplx_gt_a_i64
end interface

interface operator(.lt.)
module procedure cmplx_lt_a_i64
end interface

interface sort
module procedure sort_a_i64
end interface

interface cmplx_gt
module procedure cmplx_gt_a_i64
end interface

interface cmplx_lt
module procedure cmplx_lt_a_i64
end interface

contains
pure subroutine sort_a_i64 (arr, nskip, par)

! Perform a quicksort on an array of integers, arr(n). Uses the
! sample code in NumericalRecipies as a base.
! Optionally sort arr2 in parallel (in the routines it is enabled)

! Custom comparisons. Use the operator .locallt. & .localgt.
! interface operator(.locallt.)
!     pure function custom_lt (b, c) result (ret)
!         use constants
!         integer(kind=int64), intent(in) :: b(:)
!         integer(kind=int64), intent(in) :: c(:)
!         logical :: ret
!     end function
! end interface
! interface operator(.localgt.)
!     pure function custom_gt (b, c) result (ret)
!         use constants
!         integer(kind=int64), intent(in) :: b(:)
!         integer(kind=int64), intent(in) :: c(:)
!         logical :: ret
!     end function
! end interface

! sort needs auxiliary storage of length 2*log_2(n).
integer, parameter :: nStackMax = 50
integer, parameter :: nInsertionSort = 7
integer, intent(in), optional :: nskip
integer, intent(out), optional :: par

integer(kind=int64), intent(inout) :: arr(:,:)
!, intent(inout) :: arr2(:)
!, intent(inout) :: arr3(:)
!, intent(inout) :: arr4(:)
!, intent(inout) :: arr5(:)
!, intent(inout) :: arr6(:)

! Oh, how lovely it would be to be able to use push/pop and not need
integer :: stack(nStackMax), nstack, nskip_
integer :: pivot, lo, hi, n, i, j, par_
! n.b. This size statement is removed if type1 is scalar ...
integer(kind=int64) :: a(size(arr(:,1)))
! :: a2(size(arr2(1)))
! :: a3(size(arr3(1)))
! :: a4(size(arr4(1)))
! :: a5(size(arr5(1)))
! :: a6(size(arr6(1)))

! Temporary variables for swapping
integer(kind=int64) :: tmp1(size(arr(:,1)))
! :: tmp2(size(arr2(1)))
! :: tmp3(size(arr3(1)))
! :: tmp4(size(arr4(1)))
! :: tmp5(size(arr5(1)))
! :: tmp6(size(arr6(1)))
character(*), parameter :: this_routine = 'sort'

! Initialise temporary variables if required
!(a)
!(tmp1)
!(a2)
!(tmp2)
!(a3)
!(tmp3)
!(a4)
!(tmp4)
!(a5)
!(tmp5)
!(a6)
!(tmp6)

! *** HACK ***
! Workaround for gfortran compiler bug
! n.b. This will produce spurious warnings if run in valgrind, as
!      a is not initialised. Not a problem in optimised build.
i = 0
#ifdef DEBUG_
! it IS a problem in debug build, however
! *** HACK MOD ***
!to prevent crashes in debug mode
if((ubound(arr,2) - lbound(arr,2) + 1) &
> 0 ) a = arr(:,1)
#endif
! if (custom_lt(a, a)) i = i
! if (custom_gt(a, a)) i = i

if (present(nskip)) then
nskip_ = nskip
else
nskip_ = 1
endif

! Initialise parity calculation
par_ = 1

! The size of the array to sort.
! N.B. Use zero based indices. To obtain ! the entry into the actual
! array, multiply indices by nskip and add ! 1 (hence zero based)
! **** See local macro srt_ind() ******
lo = lbound(arr, 2) - 1
n = ((ubound(arr, 2) - lo -1)/nskip_) + 1
hi = lo + n - 1

nstack = 0
do while (.true.)
! If the section/partition we are looking at is smaller than
! nInsertSort then perform an insertion sort
if (hi - lo < nInsertionSort) then
do j = lo + 1, hi
a = arr(:,srt_ind(j))
! a2 = arr2(srt_ind(j))
! a3 = arr3(srt_ind(j))
! a4 = arr4(srt_ind(j))
! a5 = arr5(srt_ind(j))
! a6 = arr6(srt_ind(j))
do i = j - 1, 0, -1
if (arr(:,srt_ind(i)) .arrlt. a) exit
arr(:,srt_ind(i+1)) = arr(:,srt_ind(i))
! arr2(srt_ind(i+1)) = arr2(srt_ind(i))
! arr3(srt_ind(i+1)) = arr3(srt_ind(i))
! arr4(srt_ind(i+1)) = arr4(srt_ind(i))
! arr5(srt_ind(i+1)) = arr5(srt_ind(i))
! arr6(srt_ind(i+1)) = arr6(srt_ind(i))
par_ = -par_
enddo
arr(:,srt_ind(i+1)) = a
! arr2(srt_ind(i+1)) = a2
! arr3(srt_ind(i+1)) = a3
! arr4(srt_ind(i+1)) = a4
! arr5(srt_ind(i+1)) = a5
! arr6(srt_ind(i+1)) = a6
enddo

if (nstack == 0) exit
hi = stack(nstack)
lo = stack(nstack-1)
nstack = nstack - 2

! Otherwise start partitioning with quicksort.
else
! Pick a partitioning element, and arrange such that
! arr(:,lo) <= arr(:,lo+1) <= arr(:,hi)
pivot = (lo + hi) / 2
tmp1 = arr(:,srt_ind(pivot))
arr(:,srt_ind(pivot)) = arr(:,srt_ind(lo + 1))
arr(:,srt_ind(lo + 1)) = tmp1
! tmp2 = arr2(srt_ind(pivot))
! arr2(srt_ind(pivot)) = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = tmp2
! tmp3 = arr3(srt_ind(pivot))
! arr3(srt_ind(pivot)) = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = tmp3
! tmp4 = arr4(srt_ind(pivot))
! arr4(srt_ind(pivot)) = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = tmp4
! tmp5 = arr5(srt_ind(pivot))
! arr5(srt_ind(pivot)) = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = tmp5
! tmp6 = arr6(srt_ind(pivot))
! arr6(srt_ind(pivot)) = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = tmp6
par_ = -par_

if (arr(:,srt_ind(lo)) .arrgt. arr(:,srt_ind(hi))) then
tmp1 = arr(:,srt_ind(lo))
arr(:,srt_ind(lo)) = arr(:,srt_ind(hi))
arr(:,srt_ind(hi)) = tmp1
! tmp2 = arr2(srt_ind(lo))
! arr2(srt_ind(lo)) = arr2(srt_ind(hi))
! arr2(srt_ind(hi)) = tmp2
! tmp3 = arr3(srt_ind(lo))
! arr3(srt_ind(lo)) = arr3(srt_ind(hi))
! arr3(srt_ind(hi)) = tmp3
! tmp4 = arr4(srt_ind(lo))
! arr4(srt_ind(lo)) = arr4(srt_ind(hi))
! arr4(srt_ind(hi)) = tmp4
! tmp5 = arr5(srt_ind(lo))
! arr5(srt_ind(lo)) = arr5(srt_ind(hi))
! arr5(srt_ind(hi)) = tmp5
! tmp6 = arr6(srt_ind(lo))
! arr6(srt_ind(lo)) = arr6(srt_ind(hi))
! arr6(srt_ind(hi)) = tmp6
par_ = -par_
endif
if (arr(:,srt_ind(lo+1)) .arrgt. arr(:,srt_ind(hi))) then
tmp1 = arr(:,srt_ind(lo+1))
arr(:,srt_ind(lo+1)) = arr(:,srt_ind(hi))
arr(:,srt_ind(hi)) = tmp1
! tmp2 = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = arr2(srt_ind(hi))
! arr2(srt_ind(hi)) = tmp2
! tmp3 = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = arr3(srt_ind(hi))
! arr3(srt_ind(hi)) = tmp3
! tmp4 = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = arr4(srt_ind(hi))
! arr4(srt_ind(hi)) = tmp4
! tmp5 = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = arr5(srt_ind(hi))
! arr5(srt_ind(hi)) = tmp5
! tmp6 = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = arr6(srt_ind(hi))
! arr6(srt_ind(hi)) = tmp6
par_ = -par_
endif
if (arr(:,srt_ind(lo)) .arrgt. arr(:,srt_ind(lo+1))) then
tmp1 = arr(:,srt_ind(lo))
arr(:,srt_ind(lo)) = arr(:,srt_ind(lo+1))
arr(:,srt_ind(lo+1)) = tmp1
! tmp2 = arr2(srt_ind(lo))
! arr2(srt_ind(lo)) = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = tmp2
! tmp3 = arr3(srt_ind(lo))
! arr3(srt_ind(lo)) = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = tmp3
! tmp4 = arr4(srt_ind(lo))
! arr4(srt_ind(lo)) = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = tmp4
! tmp5 = arr5(srt_ind(lo))
! arr5(srt_ind(lo)) = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = tmp5
! tmp6 = arr6(srt_ind(lo))
! arr6(srt_ind(lo)) = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = tmp6
par_ = -par_
endif

i = lo + 1
j = hi
a = arr(:,srt_ind(lo + 1)) !! a is the pivot value
! a2 = arr2(srt_ind(lo + 1))
! a3 = arr3(srt_ind(lo + 1))
! a4 = arr4(srt_ind(lo + 1))
! a5 = arr5(srt_ind(lo + 1))
! a6 = arr6(srt_ind(lo + 1))
do while (.true.)
! Scand down list to find element > a
i = i + 1
do while (arr(:,srt_ind(i)) .arrlt. a)
i = i + 1
enddo

! Scan down list to find element < a
j = j - 1
do while (arr(:,srt_ind(j)) .arrgt. a)
j = j - 1
enddo

! When the pointers crossed, partitioning is complete.
if (j < i) exit

! Swap the elements, so that all elements < a end up
! in lower indexed variables.
tmp1 = arr(:,srt_ind(i))
arr(:,srt_ind(i)) = arr(:,srt_ind(j))
arr(:,srt_ind(j)) = tmp1
! tmp2 = arr2(srt_ind(i))
! arr2(srt_ind(i)) = arr2(srt_ind(j))
! arr2(srt_ind(j)) = tmp2
! tmp3 = arr3(srt_ind(i))
! arr3(srt_ind(i)) = arr3(srt_ind(j))
! arr3(srt_ind(j)) = tmp3
! tmp4 = arr4(srt_ind(i))
! arr4(srt_ind(i)) = arr4(srt_ind(j))
! arr4(srt_ind(j)) = tmp4
! tmp5 = arr5(srt_ind(i))
! arr5(srt_ind(i)) = arr5(srt_ind(j))
! arr5(srt_ind(j)) = tmp5
! tmp6 = arr6(srt_ind(i))
! arr6(srt_ind(i)) = arr6(srt_ind(j))
! arr6(srt_ind(j)) = tmp6
par_ = -par_
enddo;

! Insert partitioning element
arr(:,srt_ind(lo + 1)) = arr(:,srt_ind(j))
arr(:,srt_ind(j)) = a
! arr2(srt_ind(lo + 1)) = arr2(srt_ind(j))
! arr3(srt_ind(lo + 1)) = arr3(srt_ind(j))
! arr4(srt_ind(lo + 1)) = arr4(srt_ind(j))
! arr5(srt_ind(lo + 1)) = arr5(srt_ind(j))
! arr6(srt_ind(lo + 1)) = arr6(srt_ind(j))
! arr2(srt_ind(j)) = a2
! arr3(srt_ind(j)) = a3
! arr4(srt_ind(j)) = a4
! arr5(srt_ind(j)) = a5
! arr6(srt_ind(j)) = a6
par_ = -par_

! Push the larger of the partitioned sections onto the stack
! of sections to look at later.
! --> need fewest stack elements.
nstack = nstack + 2
if (nstack > nStackMax) then
call stop_all (this_routine, &
"parameter nStackMax too small")
endif
if (hi - i + 1 >= j - lo) then
stack(nstack) = hi
stack(nstack-1) = i
hi = j - 1
else
stack(nstack) = j - 1
stack(nstack-1) = lo
lo = i
endif
endif
enddo

! Destroy temporary variables if required
!(a)
!(tmp1)
!(a2)
!(tmp2)
!(a3)
!(tmp3)
!(a4)
!(tmp4)
!(a5)
!(tmp5)
!(a6)
!(tmp6)

! Return the parity if required
if (present(par)) &
par = par_

end subroutine

! A private comparison. This should not be available outside of the
! module!
elemental function cmplx_gt_a_i64 (a, b) result (bGt)
complex(dp), intent(in) :: a, b
logical :: bGt

bGt = real(a, dp) > real(b, dp)
end function

elemental function cmplx_lt_a_i64 (a, b) result (bLt)
complex(dp), intent(in) :: a, b
logical :: bLt

bLt = real(a, dp) < real(b, dp)
end function

end module

#define srt_ind(i) (((i)*nskip_)+1)

module sort_mod_a_i64_custom
use util_mod, only: stop_all
use util_mod_comparisons, only: operator(.arrgt.), operator(.arrlt.)

use constants
use SystemData, only: Symmetry, operator(.eq.), operator(.ne.), &
operator(.gt.), operator(.lt.), assignment(=)
use symdata, only: SymPairProd, operator(.eq.), operator(.ne.), &
operator(.gt.), operator(.lt.), assignment(=)
implicit none

private
public :: sort

! Private operator for sorting complex numbers by their real values
interface operator(.gt.)
module procedure cmplx_gt_a_i64_custom
end interface

interface operator(.lt.)
module procedure cmplx_lt_a_i64_custom
end interface

interface sort
module procedure sort_a_i64_custom
end interface

interface cmplx_gt
module procedure cmplx_gt_a_i64_custom
end interface

interface cmplx_lt
module procedure cmplx_lt_a_i64_custom
end interface

contains
pure subroutine sort_a_i64_custom (arr, custom_lt, custom_gt,nskip, par)

! Perform a quicksort on an array of integers, arr(n). Uses the
! sample code in NumericalRecipies as a base.
! Optionally sort arr2 in parallel (in the routines it is enabled)

! Custom comparisons. Use the operator .locallt. & .localgt.
interface operator(.locallt.)
pure function custom_lt (b, c) result (ret)
use constants
integer(kind=int64), intent(in) :: b(:)
integer(kind=int64), intent(in) :: c(:)
logical :: ret
end function
end interface
interface operator(.localgt.)
pure function custom_gt (b, c) result (ret)
use constants
integer(kind=int64), intent(in) :: b(:)
integer(kind=int64), intent(in) :: c(:)
logical :: ret
end function
end interface

! sort needs auxiliary storage of length 2*log_2(n).
integer, parameter :: nStackMax = 50
integer, parameter :: nInsertionSort = 7
integer, intent(in), optional :: nskip
integer, intent(out), optional :: par

integer(kind=int64), intent(inout) :: arr(:,:)
!, intent(inout) :: arr2(:)
!, intent(inout) :: arr3(:)
!, intent(inout) :: arr4(:)
!, intent(inout) :: arr5(:)
!, intent(inout) :: arr6(:)

! Oh, how lovely it would be to be able to use push/pop and not need
integer :: stack(nStackMax), nstack, nskip_
integer :: pivot, lo, hi, n, i, j, par_
! n.b. This size statement is removed if type1 is scalar ...
integer(kind=int64) :: a(size(arr(:,1)))
! :: a2(size(arr2(1)))
! :: a3(size(arr3(1)))
! :: a4(size(arr4(1)))
! :: a5(size(arr5(1)))
! :: a6(size(arr6(1)))

! Temporary variables for swapping
integer(kind=int64) :: tmp1(size(arr(:,1)))
! :: tmp2(size(arr2(1)))
! :: tmp3(size(arr3(1)))
! :: tmp4(size(arr4(1)))
! :: tmp5(size(arr5(1)))
! :: tmp6(size(arr6(1)))
character(*), parameter :: this_routine = 'sort'

! Initialise temporary variables if required
!(a)
!(tmp1)
!(a2)
!(tmp2)
!(a3)
!(tmp3)
!(a4)
!(tmp4)
!(a5)
!(tmp5)
!(a6)
!(tmp6)

! *** HACK ***
! Workaround for gfortran compiler bug
! n.b. This will produce spurious warnings if run in valgrind, as
!      a is not initialised. Not a problem in optimised build.
i = 0
#ifdef DEBUG_
! it IS a problem in debug build, however
! *** HACK MOD ***
!to prevent crashes in debug mode
if((ubound(arr,2) - lbound(arr,2) + 1) &
> 0 ) a = arr(:,1)
#endif
if (custom_lt(a, a)) i = i
if (custom_gt(a, a)) i = i

if (present(nskip)) then
nskip_ = nskip
else
nskip_ = 1
endif

! Initialise parity calculation
par_ = 1

! The size of the array to sort.
! N.B. Use zero based indices. To obtain ! the entry into the actual
! array, multiply indices by nskip and add ! 1 (hence zero based)
! **** See local macro srt_ind() ******
lo = lbound(arr, 2) - 1
n = ((ubound(arr, 2) - lo -1)/nskip_) + 1
hi = lo + n - 1

nstack = 0
do while (.true.)
! If the section/partition we are looking at is smaller than
! nInsertSort then perform an insertion sort
if (hi - lo < nInsertionSort) then
do j = lo + 1, hi
a = arr(:,srt_ind(j))
! a2 = arr2(srt_ind(j))
! a3 = arr3(srt_ind(j))
! a4 = arr4(srt_ind(j))
! a5 = arr5(srt_ind(j))
! a6 = arr6(srt_ind(j))
do i = j - 1, 0, -1
if (arr(:,srt_ind(i)) .locallt. a) exit
arr(:,srt_ind(i+1)) = arr(:,srt_ind(i))
! arr2(srt_ind(i+1)) = arr2(srt_ind(i))
! arr3(srt_ind(i+1)) = arr3(srt_ind(i))
! arr4(srt_ind(i+1)) = arr4(srt_ind(i))
! arr5(srt_ind(i+1)) = arr5(srt_ind(i))
! arr6(srt_ind(i+1)) = arr6(srt_ind(i))
par_ = -par_
enddo
arr(:,srt_ind(i+1)) = a
! arr2(srt_ind(i+1)) = a2
! arr3(srt_ind(i+1)) = a3
! arr4(srt_ind(i+1)) = a4
! arr5(srt_ind(i+1)) = a5
! arr6(srt_ind(i+1)) = a6
enddo

if (nstack == 0) exit
hi = stack(nstack)
lo = stack(nstack-1)
nstack = nstack - 2

! Otherwise start partitioning with quicksort.
else
! Pick a partitioning element, and arrange such that
! arr(:,lo) <= arr(:,lo+1) <= arr(:,hi)
pivot = (lo + hi) / 2
tmp1 = arr(:,srt_ind(pivot))
arr(:,srt_ind(pivot)) = arr(:,srt_ind(lo + 1))
arr(:,srt_ind(lo + 1)) = tmp1
! tmp2 = arr2(srt_ind(pivot))
! arr2(srt_ind(pivot)) = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = tmp2
! tmp3 = arr3(srt_ind(pivot))
! arr3(srt_ind(pivot)) = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = tmp3
! tmp4 = arr4(srt_ind(pivot))
! arr4(srt_ind(pivot)) = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = tmp4
! tmp5 = arr5(srt_ind(pivot))
! arr5(srt_ind(pivot)) = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = tmp5
! tmp6 = arr6(srt_ind(pivot))
! arr6(srt_ind(pivot)) = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = tmp6
par_ = -par_

if (arr(:,srt_ind(lo)) .localgt. arr(:,srt_ind(hi))) then
tmp1 = arr(:,srt_ind(lo))
arr(:,srt_ind(lo)) = arr(:,srt_ind(hi))
arr(:,srt_ind(hi)) = tmp1
! tmp2 = arr2(srt_ind(lo))
! arr2(srt_ind(lo)) = arr2(srt_ind(hi))
! arr2(srt_ind(hi)) = tmp2
! tmp3 = arr3(srt_ind(lo))
! arr3(srt_ind(lo)) = arr3(srt_ind(hi))
! arr3(srt_ind(hi)) = tmp3
! tmp4 = arr4(srt_ind(lo))
! arr4(srt_ind(lo)) = arr4(srt_ind(hi))
! arr4(srt_ind(hi)) = tmp4
! tmp5 = arr5(srt_ind(lo))
! arr5(srt_ind(lo)) = arr5(srt_ind(hi))
! arr5(srt_ind(hi)) = tmp5
! tmp6 = arr6(srt_ind(lo))
! arr6(srt_ind(lo)) = arr6(srt_ind(hi))
! arr6(srt_ind(hi)) = tmp6
par_ = -par_
endif
if (arr(:,srt_ind(lo+1)) .localgt. arr(:,srt_ind(hi))) then
tmp1 = arr(:,srt_ind(lo+1))
arr(:,srt_ind(lo+1)) = arr(:,srt_ind(hi))
arr(:,srt_ind(hi)) = tmp1
! tmp2 = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = arr2(srt_ind(hi))
! arr2(srt_ind(hi)) = tmp2
! tmp3 = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = arr3(srt_ind(hi))
! arr3(srt_ind(hi)) = tmp3
! tmp4 = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = arr4(srt_ind(hi))
! arr4(srt_ind(hi)) = tmp4
! tmp5 = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = arr5(srt_ind(hi))
! arr5(srt_ind(hi)) = tmp5
! tmp6 = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = arr6(srt_ind(hi))
! arr6(srt_ind(hi)) = tmp6
par_ = -par_
endif
if (arr(:,srt_ind(lo)) .localgt. arr(:,srt_ind(lo+1))) then
tmp1 = arr(:,srt_ind(lo))
arr(:,srt_ind(lo)) = arr(:,srt_ind(lo+1))
arr(:,srt_ind(lo+1)) = tmp1
! tmp2 = arr2(srt_ind(lo))
! arr2(srt_ind(lo)) = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = tmp2
! tmp3 = arr3(srt_ind(lo))
! arr3(srt_ind(lo)) = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = tmp3
! tmp4 = arr4(srt_ind(lo))
! arr4(srt_ind(lo)) = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = tmp4
! tmp5 = arr5(srt_ind(lo))
! arr5(srt_ind(lo)) = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = tmp5
! tmp6 = arr6(srt_ind(lo))
! arr6(srt_ind(lo)) = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = tmp6
par_ = -par_
endif

i = lo + 1
j = hi
a = arr(:,srt_ind(lo + 1)) !! a is the pivot value
! a2 = arr2(srt_ind(lo + 1))
! a3 = arr3(srt_ind(lo + 1))
! a4 = arr4(srt_ind(lo + 1))
! a5 = arr5(srt_ind(lo + 1))
! a6 = arr6(srt_ind(lo + 1))
do while (.true.)
! Scand down list to find element > a
i = i + 1
do while (arr(:,srt_ind(i)) .locallt. a)
i = i + 1
enddo

! Scan down list to find element < a
j = j - 1
do while (arr(:,srt_ind(j)) .localgt. a)
j = j - 1
enddo

! When the pointers crossed, partitioning is complete.
if (j < i) exit

! Swap the elements, so that all elements < a end up
! in lower indexed variables.
tmp1 = arr(:,srt_ind(i))
arr(:,srt_ind(i)) = arr(:,srt_ind(j))
arr(:,srt_ind(j)) = tmp1
! tmp2 = arr2(srt_ind(i))
! arr2(srt_ind(i)) = arr2(srt_ind(j))
! arr2(srt_ind(j)) = tmp2
! tmp3 = arr3(srt_ind(i))
! arr3(srt_ind(i)) = arr3(srt_ind(j))
! arr3(srt_ind(j)) = tmp3
! tmp4 = arr4(srt_ind(i))
! arr4(srt_ind(i)) = arr4(srt_ind(j))
! arr4(srt_ind(j)) = tmp4
! tmp5 = arr5(srt_ind(i))
! arr5(srt_ind(i)) = arr5(srt_ind(j))
! arr5(srt_ind(j)) = tmp5
! tmp6 = arr6(srt_ind(i))
! arr6(srt_ind(i)) = arr6(srt_ind(j))
! arr6(srt_ind(j)) = tmp6
par_ = -par_
enddo;

! Insert partitioning element
arr(:,srt_ind(lo + 1)) = arr(:,srt_ind(j))
arr(:,srt_ind(j)) = a
! arr2(srt_ind(lo + 1)) = arr2(srt_ind(j))
! arr3(srt_ind(lo + 1)) = arr3(srt_ind(j))
! arr4(srt_ind(lo + 1)) = arr4(srt_ind(j))
! arr5(srt_ind(lo + 1)) = arr5(srt_ind(j))
! arr6(srt_ind(lo + 1)) = arr6(srt_ind(j))
! arr2(srt_ind(j)) = a2
! arr3(srt_ind(j)) = a3
! arr4(srt_ind(j)) = a4
! arr5(srt_ind(j)) = a5
! arr6(srt_ind(j)) = a6
par_ = -par_

! Push the larger of the partitioned sections onto the stack
! of sections to look at later.
! --> need fewest stack elements.
nstack = nstack + 2
if (nstack > nStackMax) then
call stop_all (this_routine, &
"parameter nStackMax too small")
endif
if (hi - i + 1 >= j - lo) then
stack(nstack) = hi
stack(nstack-1) = i
hi = j - 1
else
stack(nstack) = j - 1
stack(nstack-1) = lo
lo = i
endif
endif
enddo

! Destroy temporary variables if required
!(a)
!(tmp1)
!(a2)
!(tmp2)
!(a3)
!(tmp3)
!(a4)
!(tmp4)
!(a5)
!(tmp5)
!(a6)
!(tmp6)

! Return the parity if required
if (present(par)) &
par = par_

end subroutine

! A private comparison. This should not be available outside of the
! module!
elemental function cmplx_gt_a_i64_custom (a, b) result (bGt)
complex(dp), intent(in) :: a, b
logical :: bGt

bGt = real(a, dp) > real(b, dp)
end function

elemental function cmplx_lt_a_i64_custom (a, b) result (bLt)
complex(dp), intent(in) :: a, b
logical :: bLt

bLt = real(a, dp) < real(b, dp)
end function

end module

#define srt_ind(i) (((i)*nskip_)+1)

module sort_mod_a_i_custom
use util_mod, only: stop_all
use util_mod_comparisons, only: operator(.arrgt.), operator(.arrlt.)

use constants
use SystemData, only: Symmetry, operator(.eq.), operator(.ne.), &
operator(.gt.), operator(.lt.), assignment(=)
use symdata, only: SymPairProd, operator(.eq.), operator(.ne.), &
operator(.gt.), operator(.lt.), assignment(=)
implicit none

private
public :: sort

! Private operator for sorting complex numbers by their real values
interface operator(.gt.)
module procedure cmplx_gt_a_i_custom
end interface

interface operator(.lt.)
module procedure cmplx_lt_a_i_custom
end interface

interface sort
module procedure sort_a_i_custom
end interface

interface cmplx_gt
module procedure cmplx_gt_a_i_custom
end interface

interface cmplx_lt
module procedure cmplx_lt_a_i_custom
end interface

contains
pure subroutine sort_a_i_custom (arr, custom_lt, custom_gt,nskip, par)

! Perform a quicksort on an array of integers, arr(n). Uses the
! sample code in NumericalRecipies as a base.
! Optionally sort arr2 in parallel (in the routines it is enabled)

! Custom comparisons. Use the operator .locallt. & .localgt.
interface operator(.locallt.)
pure function custom_lt (b, c) result (ret)
use constants
integer(kind=int32), intent(in) :: b(:)
integer(kind=int32), intent(in) :: c(:)
logical :: ret
end function
end interface
interface operator(.localgt.)
pure function custom_gt (b, c) result (ret)
use constants
integer(kind=int32), intent(in) :: b(:)
integer(kind=int32), intent(in) :: c(:)
logical :: ret
end function
end interface

! sort needs auxiliary storage of length 2*log_2(n).
integer, parameter :: nStackMax = 50
integer, parameter :: nInsertionSort = 7
integer, intent(in), optional :: nskip
integer, intent(out), optional :: par

integer(kind=int32), intent(inout) :: arr(:,:)
!, intent(inout) :: arr2(:)
!, intent(inout) :: arr3(:)
!, intent(inout) :: arr4(:)
!, intent(inout) :: arr5(:)
!, intent(inout) :: arr6(:)

! Oh, how lovely it would be to be able to use push/pop and not need
integer :: stack(nStackMax), nstack, nskip_
integer :: pivot, lo, hi, n, i, j, par_
! n.b. This size statement is removed if type1 is scalar ...
integer(kind=int32) :: a(size(arr(:,1)))
! :: a2(size(arr2(1)))
! :: a3(size(arr3(1)))
! :: a4(size(arr4(1)))
! :: a5(size(arr5(1)))
! :: a6(size(arr6(1)))

! Temporary variables for swapping
integer(kind=int32) :: tmp1(size(arr(:,1)))
! :: tmp2(size(arr2(1)))
! :: tmp3(size(arr3(1)))
! :: tmp4(size(arr4(1)))
! :: tmp5(size(arr5(1)))
! :: tmp6(size(arr6(1)))
character(*), parameter :: this_routine = 'sort'

! Initialise temporary variables if required
!(a)
!(tmp1)
!(a2)
!(tmp2)
!(a3)
!(tmp3)
!(a4)
!(tmp4)
!(a5)
!(tmp5)
!(a6)
!(tmp6)

! *** HACK ***
! Workaround for gfortran compiler bug
! n.b. This will produce spurious warnings if run in valgrind, as
!      a is not initialised. Not a problem in optimised build.
i = 0
#ifdef DEBUG_
! it IS a problem in debug build, however
! *** HACK MOD ***
!to prevent crashes in debug mode
if((ubound(arr,2) - lbound(arr,2) + 1) &
> 0 ) a = arr(:,1)
#endif
if (custom_lt(a, a)) i = i
if (custom_gt(a, a)) i = i

if (present(nskip)) then
nskip_ = nskip
else
nskip_ = 1
endif

! Initialise parity calculation
par_ = 1

! The size of the array to sort.
! N.B. Use zero based indices. To obtain ! the entry into the actual
! array, multiply indices by nskip and add ! 1 (hence zero based)
! **** See local macro srt_ind() ******
lo = lbound(arr, 2) - 1
n = ((ubound(arr, 2) - lo -1)/nskip_) + 1
hi = lo + n - 1

nstack = 0
do while (.true.)
! If the section/partition we are looking at is smaller than
! nInsertSort then perform an insertion sort
if (hi - lo < nInsertionSort) then
do j = lo + 1, hi
a = arr(:,srt_ind(j))
! a2 = arr2(srt_ind(j))
! a3 = arr3(srt_ind(j))
! a4 = arr4(srt_ind(j))
! a5 = arr5(srt_ind(j))
! a6 = arr6(srt_ind(j))
do i = j - 1, 0, -1
if (arr(:,srt_ind(i)) .locallt. a) exit
arr(:,srt_ind(i+1)) = arr(:,srt_ind(i))
! arr2(srt_ind(i+1)) = arr2(srt_ind(i))
! arr3(srt_ind(i+1)) = arr3(srt_ind(i))
! arr4(srt_ind(i+1)) = arr4(srt_ind(i))
! arr5(srt_ind(i+1)) = arr5(srt_ind(i))
! arr6(srt_ind(i+1)) = arr6(srt_ind(i))
par_ = -par_
enddo
arr(:,srt_ind(i+1)) = a
! arr2(srt_ind(i+1)) = a2
! arr3(srt_ind(i+1)) = a3
! arr4(srt_ind(i+1)) = a4
! arr5(srt_ind(i+1)) = a5
! arr6(srt_ind(i+1)) = a6
enddo

if (nstack == 0) exit
hi = stack(nstack)
lo = stack(nstack-1)
nstack = nstack - 2

! Otherwise start partitioning with quicksort.
else
! Pick a partitioning element, and arrange such that
! arr(:,lo) <= arr(:,lo+1) <= arr(:,hi)
pivot = (lo + hi) / 2
tmp1 = arr(:,srt_ind(pivot))
arr(:,srt_ind(pivot)) = arr(:,srt_ind(lo + 1))
arr(:,srt_ind(lo + 1)) = tmp1
! tmp2 = arr2(srt_ind(pivot))
! arr2(srt_ind(pivot)) = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = tmp2
! tmp3 = arr3(srt_ind(pivot))
! arr3(srt_ind(pivot)) = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = tmp3
! tmp4 = arr4(srt_ind(pivot))
! arr4(srt_ind(pivot)) = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = tmp4
! tmp5 = arr5(srt_ind(pivot))
! arr5(srt_ind(pivot)) = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = tmp5
! tmp6 = arr6(srt_ind(pivot))
! arr6(srt_ind(pivot)) = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = tmp6
par_ = -par_

if (arr(:,srt_ind(lo)) .localgt. arr(:,srt_ind(hi))) then
tmp1 = arr(:,srt_ind(lo))
arr(:,srt_ind(lo)) = arr(:,srt_ind(hi))
arr(:,srt_ind(hi)) = tmp1
! tmp2 = arr2(srt_ind(lo))
! arr2(srt_ind(lo)) = arr2(srt_ind(hi))
! arr2(srt_ind(hi)) = tmp2
! tmp3 = arr3(srt_ind(lo))
! arr3(srt_ind(lo)) = arr3(srt_ind(hi))
! arr3(srt_ind(hi)) = tmp3
! tmp4 = arr4(srt_ind(lo))
! arr4(srt_ind(lo)) = arr4(srt_ind(hi))
! arr4(srt_ind(hi)) = tmp4
! tmp5 = arr5(srt_ind(lo))
! arr5(srt_ind(lo)) = arr5(srt_ind(hi))
! arr5(srt_ind(hi)) = tmp5
! tmp6 = arr6(srt_ind(lo))
! arr6(srt_ind(lo)) = arr6(srt_ind(hi))
! arr6(srt_ind(hi)) = tmp6
par_ = -par_
endif
if (arr(:,srt_ind(lo+1)) .localgt. arr(:,srt_ind(hi))) then
tmp1 = arr(:,srt_ind(lo+1))
arr(:,srt_ind(lo+1)) = arr(:,srt_ind(hi))
arr(:,srt_ind(hi)) = tmp1
! tmp2 = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = arr2(srt_ind(hi))
! arr2(srt_ind(hi)) = tmp2
! tmp3 = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = arr3(srt_ind(hi))
! arr3(srt_ind(hi)) = tmp3
! tmp4 = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = arr4(srt_ind(hi))
! arr4(srt_ind(hi)) = tmp4
! tmp5 = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = arr5(srt_ind(hi))
! arr5(srt_ind(hi)) = tmp5
! tmp6 = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = arr6(srt_ind(hi))
! arr6(srt_ind(hi)) = tmp6
par_ = -par_
endif
if (arr(:,srt_ind(lo)) .localgt. arr(:,srt_ind(lo+1))) then
tmp1 = arr(:,srt_ind(lo))
arr(:,srt_ind(lo)) = arr(:,srt_ind(lo+1))
arr(:,srt_ind(lo+1)) = tmp1
! tmp2 = arr2(srt_ind(lo))
! arr2(srt_ind(lo)) = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = tmp2
! tmp3 = arr3(srt_ind(lo))
! arr3(srt_ind(lo)) = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = tmp3
! tmp4 = arr4(srt_ind(lo))
! arr4(srt_ind(lo)) = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = tmp4
! tmp5 = arr5(srt_ind(lo))
! arr5(srt_ind(lo)) = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = tmp5
! tmp6 = arr6(srt_ind(lo))
! arr6(srt_ind(lo)) = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = tmp6
par_ = -par_
endif

i = lo + 1
j = hi
a = arr(:,srt_ind(lo + 1)) !! a is the pivot value
! a2 = arr2(srt_ind(lo + 1))
! a3 = arr3(srt_ind(lo + 1))
! a4 = arr4(srt_ind(lo + 1))
! a5 = arr5(srt_ind(lo + 1))
! a6 = arr6(srt_ind(lo + 1))
do while (.true.)
! Scand down list to find element > a
i = i + 1
do while (arr(:,srt_ind(i)) .locallt. a)
i = i + 1
enddo

! Scan down list to find element < a
j = j - 1
do while (arr(:,srt_ind(j)) .localgt. a)
j = j - 1
enddo

! When the pointers crossed, partitioning is complete.
if (j < i) exit

! Swap the elements, so that all elements < a end up
! in lower indexed variables.
tmp1 = arr(:,srt_ind(i))
arr(:,srt_ind(i)) = arr(:,srt_ind(j))
arr(:,srt_ind(j)) = tmp1
! tmp2 = arr2(srt_ind(i))
! arr2(srt_ind(i)) = arr2(srt_ind(j))
! arr2(srt_ind(j)) = tmp2
! tmp3 = arr3(srt_ind(i))
! arr3(srt_ind(i)) = arr3(srt_ind(j))
! arr3(srt_ind(j)) = tmp3
! tmp4 = arr4(srt_ind(i))
! arr4(srt_ind(i)) = arr4(srt_ind(j))
! arr4(srt_ind(j)) = tmp4
! tmp5 = arr5(srt_ind(i))
! arr5(srt_ind(i)) = arr5(srt_ind(j))
! arr5(srt_ind(j)) = tmp5
! tmp6 = arr6(srt_ind(i))
! arr6(srt_ind(i)) = arr6(srt_ind(j))
! arr6(srt_ind(j)) = tmp6
par_ = -par_
enddo;

! Insert partitioning element
arr(:,srt_ind(lo + 1)) = arr(:,srt_ind(j))
arr(:,srt_ind(j)) = a
! arr2(srt_ind(lo + 1)) = arr2(srt_ind(j))
! arr3(srt_ind(lo + 1)) = arr3(srt_ind(j))
! arr4(srt_ind(lo + 1)) = arr4(srt_ind(j))
! arr5(srt_ind(lo + 1)) = arr5(srt_ind(j))
! arr6(srt_ind(lo + 1)) = arr6(srt_ind(j))
! arr2(srt_ind(j)) = a2
! arr3(srt_ind(j)) = a3
! arr4(srt_ind(j)) = a4
! arr5(srt_ind(j)) = a5
! arr6(srt_ind(j)) = a6
par_ = -par_

! Push the larger of the partitioned sections onto the stack
! of sections to look at later.
! --> need fewest stack elements.
nstack = nstack + 2
if (nstack > nStackMax) then
call stop_all (this_routine, &
"parameter nStackMax too small")
endif
if (hi - i + 1 >= j - lo) then
stack(nstack) = hi
stack(nstack-1) = i
hi = j - 1
else
stack(nstack) = j - 1
stack(nstack-1) = lo
lo = i
endif
endif
enddo

! Destroy temporary variables if required
!(a)
!(tmp1)
!(a2)
!(tmp2)
!(a3)
!(tmp3)
!(a4)
!(tmp4)
!(a5)
!(tmp5)
!(a6)
!(tmp6)

! Return the parity if required
if (present(par)) &
par = par_

end subroutine

! A private comparison. This should not be available outside of the
! module!
elemental function cmplx_gt_a_i_custom (a, b) result (bGt)
complex(dp), intent(in) :: a, b
logical :: bGt

bGt = real(a, dp) > real(b, dp)
end function

elemental function cmplx_lt_a_i_custom (a, b) result (bLt)
complex(dp), intent(in) :: a, b
logical :: bLt

bLt = real(a, dp) < real(b, dp)
end function

end module

#define srt_ind(i) (((i)*nskip_)+1)

module sort_mod_a_d
use util_mod, only: stop_all
use util_mod_comparisons, only: operator(.arrgt.), operator(.arrlt.)

use constants
use SystemData, only: Symmetry, operator(.eq.), operator(.ne.), &
operator(.gt.), operator(.lt.), assignment(=)
use symdata, only: SymPairProd, operator(.eq.), operator(.ne.), &
operator(.gt.), operator(.lt.), assignment(=)
implicit none

private
public :: sort

! Private operator for sorting complex numbers by their real values
interface operator(.gt.)
module procedure cmplx_gt_a_d
end interface

interface operator(.lt.)
module procedure cmplx_lt_a_d
end interface

interface sort
module procedure sort_a_d
end interface

interface cmplx_gt
module procedure cmplx_gt_a_d
end interface

interface cmplx_lt
module procedure cmplx_lt_a_d
end interface

contains
pure subroutine sort_a_d (arr, nskip, par)

! Perform a quicksort on an array of integers, arr(n). Uses the
! sample code in NumericalRecipies as a base.
! Optionally sort arr2 in parallel (in the routines it is enabled)

! Custom comparisons. Use the operator .locallt. & .localgt.
! interface operator(.locallt.)
!     pure function custom_lt (b, c) result (ret)
!         use constants
!         real(dp), intent(in) :: b(:)
!         real(dp), intent(in) :: c(:)
!         logical :: ret
!     end function
! end interface
! interface operator(.localgt.)
!     pure function custom_gt (b, c) result (ret)
!         use constants
!         real(dp), intent(in) :: b(:)
!         real(dp), intent(in) :: c(:)
!         logical :: ret
!     end function
! end interface

! sort needs auxiliary storage of length 2*log_2(n).
integer, parameter :: nStackMax = 50
integer, parameter :: nInsertionSort = 7
integer, intent(in), optional :: nskip
integer, intent(out), optional :: par

real(dp), intent(inout) :: arr(:,:)
!, intent(inout) :: arr2(:)
!, intent(inout) :: arr3(:)
!, intent(inout) :: arr4(:)
!, intent(inout) :: arr5(:)
!, intent(inout) :: arr6(:)

! Oh, how lovely it would be to be able to use push/pop and not need
integer :: stack(nStackMax), nstack, nskip_
integer :: pivot, lo, hi, n, i, j, par_
! n.b. This size statement is removed if type1 is scalar ...
real(dp) :: a(size(arr(:,1)))
! :: a2(size(arr2(1)))
! :: a3(size(arr3(1)))
! :: a4(size(arr4(1)))
! :: a5(size(arr5(1)))
! :: a6(size(arr6(1)))

! Temporary variables for swapping
real(dp) :: tmp1(size(arr(:,1)))
! :: tmp2(size(arr2(1)))
! :: tmp3(size(arr3(1)))
! :: tmp4(size(arr4(1)))
! :: tmp5(size(arr5(1)))
! :: tmp6(size(arr6(1)))
character(*), parameter :: this_routine = 'sort'

! Initialise temporary variables if required
!(a)
!(tmp1)
!(a2)
!(tmp2)
!(a3)
!(tmp3)
!(a4)
!(tmp4)
!(a5)
!(tmp5)
!(a6)
!(tmp6)

! *** HACK ***
! Workaround for gfortran compiler bug
! n.b. This will produce spurious warnings if run in valgrind, as
!      a is not initialised. Not a problem in optimised build.
i = 0
#ifdef DEBUG_
! it IS a problem in debug build, however
! *** HACK MOD ***
!to prevent crashes in debug mode
if((ubound(arr,2) - lbound(arr,2) + 1) &
> 0 ) a = arr(:,1)
#endif
! if (custom_lt(a, a)) i = i
! if (custom_gt(a, a)) i = i

if (present(nskip)) then
nskip_ = nskip
else
nskip_ = 1
endif

! Initialise parity calculation
par_ = 1

! The size of the array to sort.
! N.B. Use zero based indices. To obtain ! the entry into the actual
! array, multiply indices by nskip and add ! 1 (hence zero based)
! **** See local macro srt_ind() ******
lo = lbound(arr, 2) - 1
n = ((ubound(arr, 2) - lo -1)/nskip_) + 1
hi = lo + n - 1

nstack = 0
do while (.true.)
! If the section/partition we are looking at is smaller than
! nInsertSort then perform an insertion sort
if (hi - lo < nInsertionSort) then
do j = lo + 1, hi
a = arr(:,srt_ind(j))
! a2 = arr2(srt_ind(j))
! a3 = arr3(srt_ind(j))
! a4 = arr4(srt_ind(j))
! a5 = arr5(srt_ind(j))
! a6 = arr6(srt_ind(j))
do i = j - 1, 0, -1
if (arr(:,srt_ind(i)) .arrlt. a) exit
arr(:,srt_ind(i+1)) = arr(:,srt_ind(i))
! arr2(srt_ind(i+1)) = arr2(srt_ind(i))
! arr3(srt_ind(i+1)) = arr3(srt_ind(i))
! arr4(srt_ind(i+1)) = arr4(srt_ind(i))
! arr5(srt_ind(i+1)) = arr5(srt_ind(i))
! arr6(srt_ind(i+1)) = arr6(srt_ind(i))
par_ = -par_
enddo
arr(:,srt_ind(i+1)) = a
! arr2(srt_ind(i+1)) = a2
! arr3(srt_ind(i+1)) = a3
! arr4(srt_ind(i+1)) = a4
! arr5(srt_ind(i+1)) = a5
! arr6(srt_ind(i+1)) = a6
enddo

if (nstack == 0) exit
hi = stack(nstack)
lo = stack(nstack-1)
nstack = nstack - 2

! Otherwise start partitioning with quicksort.
else
! Pick a partitioning element, and arrange such that
! arr(:,lo) <= arr(:,lo+1) <= arr(:,hi)
pivot = (lo + hi) / 2
tmp1 = arr(:,srt_ind(pivot))
arr(:,srt_ind(pivot)) = arr(:,srt_ind(lo + 1))
arr(:,srt_ind(lo + 1)) = tmp1
! tmp2 = arr2(srt_ind(pivot))
! arr2(srt_ind(pivot)) = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = tmp2
! tmp3 = arr3(srt_ind(pivot))
! arr3(srt_ind(pivot)) = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = tmp3
! tmp4 = arr4(srt_ind(pivot))
! arr4(srt_ind(pivot)) = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = tmp4
! tmp5 = arr5(srt_ind(pivot))
! arr5(srt_ind(pivot)) = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = tmp5
! tmp6 = arr6(srt_ind(pivot))
! arr6(srt_ind(pivot)) = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = tmp6
par_ = -par_

if (arr(:,srt_ind(lo)) .arrgt. arr(:,srt_ind(hi))) then
tmp1 = arr(:,srt_ind(lo))
arr(:,srt_ind(lo)) = arr(:,srt_ind(hi))
arr(:,srt_ind(hi)) = tmp1
! tmp2 = arr2(srt_ind(lo))
! arr2(srt_ind(lo)) = arr2(srt_ind(hi))
! arr2(srt_ind(hi)) = tmp2
! tmp3 = arr3(srt_ind(lo))
! arr3(srt_ind(lo)) = arr3(srt_ind(hi))
! arr3(srt_ind(hi)) = tmp3
! tmp4 = arr4(srt_ind(lo))
! arr4(srt_ind(lo)) = arr4(srt_ind(hi))
! arr4(srt_ind(hi)) = tmp4
! tmp5 = arr5(srt_ind(lo))
! arr5(srt_ind(lo)) = arr5(srt_ind(hi))
! arr5(srt_ind(hi)) = tmp5
! tmp6 = arr6(srt_ind(lo))
! arr6(srt_ind(lo)) = arr6(srt_ind(hi))
! arr6(srt_ind(hi)) = tmp6
par_ = -par_
endif
if (arr(:,srt_ind(lo+1)) .arrgt. arr(:,srt_ind(hi))) then
tmp1 = arr(:,srt_ind(lo+1))
arr(:,srt_ind(lo+1)) = arr(:,srt_ind(hi))
arr(:,srt_ind(hi)) = tmp1
! tmp2 = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = arr2(srt_ind(hi))
! arr2(srt_ind(hi)) = tmp2
! tmp3 = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = arr3(srt_ind(hi))
! arr3(srt_ind(hi)) = tmp3
! tmp4 = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = arr4(srt_ind(hi))
! arr4(srt_ind(hi)) = tmp4
! tmp5 = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = arr5(srt_ind(hi))
! arr5(srt_ind(hi)) = tmp5
! tmp6 = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = arr6(srt_ind(hi))
! arr6(srt_ind(hi)) = tmp6
par_ = -par_
endif
if (arr(:,srt_ind(lo)) .arrgt. arr(:,srt_ind(lo+1))) then
tmp1 = arr(:,srt_ind(lo))
arr(:,srt_ind(lo)) = arr(:,srt_ind(lo+1))
arr(:,srt_ind(lo+1)) = tmp1
! tmp2 = arr2(srt_ind(lo))
! arr2(srt_ind(lo)) = arr2(srt_ind(lo+1))
! arr2(srt_ind(lo+1)) = tmp2
! tmp3 = arr3(srt_ind(lo))
! arr3(srt_ind(lo)) = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = tmp3
! tmp4 = arr4(srt_ind(lo))
! arr4(srt_ind(lo)) = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = tmp4
! tmp5 = arr5(srt_ind(lo))
! arr5(srt_ind(lo)) = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = tmp5
! tmp6 = arr6(srt_ind(lo))
! arr6(srt_ind(lo)) = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = tmp6
par_ = -par_
endif

i = lo + 1
j = hi
a = arr(:,srt_ind(lo + 1)) !! a is the pivot value
! a2 = arr2(srt_ind(lo + 1))
! a3 = arr3(srt_ind(lo + 1))
! a4 = arr4(srt_ind(lo + 1))
! a5 = arr5(srt_ind(lo + 1))
! a6 = arr6(srt_ind(lo + 1))
do while (.true.)
! Scand down list to find element > a
i = i + 1
do while (arr(:,srt_ind(i)) .arrlt. a)
i = i + 1
enddo

! Scan down list to find element < a
j = j - 1
do while (arr(:,srt_ind(j)) .arrgt. a)
j = j - 1
enddo

! When the pointers crossed, partitioning is complete.
if (j < i) exit

! Swap the elements, so that all elements < a end up
! in lower indexed variables.
tmp1 = arr(:,srt_ind(i))
arr(:,srt_ind(i)) = arr(:,srt_ind(j))
arr(:,srt_ind(j)) = tmp1
! tmp2 = arr2(srt_ind(i))
! arr2(srt_ind(i)) = arr2(srt_ind(j))
! arr2(srt_ind(j)) = tmp2
! tmp3 = arr3(srt_ind(i))
! arr3(srt_ind(i)) = arr3(srt_ind(j))
! arr3(srt_ind(j)) = tmp3
! tmp4 = arr4(srt_ind(i))
! arr4(srt_ind(i)) = arr4(srt_ind(j))
! arr4(srt_ind(j)) = tmp4
! tmp5 = arr5(srt_ind(i))
! arr5(srt_ind(i)) = arr5(srt_ind(j))
! arr5(srt_ind(j)) = tmp5
! tmp6 = arr6(srt_ind(i))
! arr6(srt_ind(i)) = arr6(srt_ind(j))
! arr6(srt_ind(j)) = tmp6
par_ = -par_
enddo;

! Insert partitioning element
arr(:,srt_ind(lo + 1)) = arr(:,srt_ind(j))
arr(:,srt_ind(j)) = a
! arr2(srt_ind(lo + 1)) = arr2(srt_ind(j))
! arr3(srt_ind(lo + 1)) = arr3(srt_ind(j))
! arr4(srt_ind(lo + 1)) = arr4(srt_ind(j))
! arr5(srt_ind(lo + 1)) = arr5(srt_ind(j))
! arr6(srt_ind(lo + 1)) = arr6(srt_ind(j))
! arr2(srt_ind(j)) = a2
! arr3(srt_ind(j)) = a3
! arr4(srt_ind(j)) = a4
! arr5(srt_ind(j)) = a5
! arr6(srt_ind(j)) = a6
par_ = -par_

! Push the larger of the partitioned sections onto the stack
! of sections to look at later.
! --> need fewest stack elements.
nstack = nstack + 2
if (nstack > nStackMax) then
call stop_all (this_routine, &
"parameter nStackMax too small")
endif
if (hi - i + 1 >= j - lo) then
stack(nstack) = hi
stack(nstack-1) = i
hi = j - 1
else
stack(nstack) = j - 1
stack(nstack-1) = lo
lo = i
endif
endif
enddo

! Destroy temporary variables if required
!(a)
!(tmp1)
!(a2)
!(tmp2)
!(a3)
!(tmp3)
!(a4)
!(tmp4)
!(a5)
!(tmp5)
!(a6)
!(tmp6)

! Return the parity if required
if (present(par)) &
par = par_

end subroutine

! A private comparison. This should not be available outside of the
! module!
elemental function cmplx_gt_a_d (a, b) result (bGt)
complex(dp), intent(in) :: a, b
logical :: bGt

bGt = real(a, dp) > real(b, dp)
end function

elemental function cmplx_lt_a_d (a, b) result (bLt)
complex(dp), intent(in) :: a, b
logical :: bLt

bLt = real(a, dp) < real(b, dp)
end function

end module

#if !defined(SX)

#define srt_ind(i) (((i)*nskip_)+1)

module sort_mod_i_i
use util_mod, only: stop_all
use util_mod_comparisons, only: operator(.arrgt.), operator(.arrlt.)

use constants
use SystemData, only: Symmetry, operator(.eq.), operator(.ne.), &
operator(.gt.), operator(.lt.), assignment(=)
use symdata, only: SymPairProd, operator(.eq.), operator(.ne.), &
operator(.gt.), operator(.lt.), assignment(=)
implicit none

private
public :: sort

! Private operator for sorting complex numbers by their real values
interface operator(.gt.)
module procedure cmplx_gt_i_i
end interface

interface operator(.lt.)
module procedure cmplx_lt_i_i
end interface

interface sort
module procedure sort_i_i
end interface

interface cmplx_gt
module procedure cmplx_gt_i_i
end interface

interface cmplx_lt
module procedure cmplx_lt_i_i
end interface

contains
pure subroutine sort_i_i (arr, arr2,nskip, par)

! Perform a quicksort on an array of integers, arr(n). Uses the
! sample code in NumericalRecipies as a base.
! Optionally sort arr2 in parallel (in the routines it is enabled)

! Custom comparisons. Use the operator .locallt. & .localgt.
! interface operator(.locallt.)
!     pure function custom_lt (b, c) result (ret)
!         use constants
!         integer(kind=int32), intent(in) :: b(:)
!         integer(kind=int32), intent(in) :: c(:)
!         logical :: ret
!     end function
! end interface

! interface operator(.localgt.)
!     pure function custom_gt (b, c) result (ret)
!         use constants
!         integer(kind=int32), intent(in) :: b(:)
!         integer(kind=int32), intent(in) :: c(:)
!         logical :: ret
!     end function
! end interface

! sort needs auxiliary storage of length 2*log_2(n).
integer, parameter :: nStackMax = 50
integer, parameter :: nInsertionSort = 7
integer, intent(in), optional :: nskip
integer, intent(out), optional :: par

integer(kind=int32), intent(inout) :: arr(:)
integer(kind=int32), intent(inout) :: arr2(:)
!, intent(inout) :: arr3(:)
!, intent(inout) :: arr4(:)
!, intent(inout) :: arr5(:)
!, intent(inout) :: arr6(:)

! Oh, how lovely it would be to be able to use push/pop and not need
integer :: stack(nStackMax), nstack, nskip_
integer :: pivot, lo, hi, n, i, j, par_
! n.b. This size statement is removed if type1 is scalar ...
integer(kind=int32) :: a
integer(kind=int32) :: a2
! :: a3(size(arr3(1)))
! :: a4(size(arr4(1)))
! :: a5(size(arr5(1)))
! :: a6(size(arr6(1)))

! Temporary variables for swapping
integer(kind=int32) :: tmp1
integer(kind=int32) :: tmp2
! :: tmp3(size(arr3(1)))
! :: tmp4(size(arr4(1)))
! :: tmp5(size(arr5(1)))
! :: tmp6(size(arr6(1)))
character(*), parameter :: this_routine = 'sort'

! Initialise temporary variables if required
!(a)
!(tmp1)
!(a2)
!(tmp2)
!(a3)
!(tmp3)
!(a4)
!(tmp4)
!(a5)
!(tmp5)
!(a6)
!(tmp6)

! *** HACK ***
! Workaround for gfortran compiler bug
! n.b. This will produce spurious warnings if run in valgrind, as
!      a is not initialised. Not a problem in optimised build.
i = 0
#ifdef DEBUG_
! it IS a problem in debug build, however
! *** HACK MOD ***
!to prevent crashes in debug mode
if((ubound(arr,1) - lbound(arr,1) + 1) &
> 0 ) a = arr(1)
#endif
! if (custom_lt(a, a)) i = i
! if (custom_gt(a, a)) i = i

if (present(nskip)) then
nskip_ = nskip
else
nskip_ = 1
endif

! Initialise parity calculation
par_ = 1

! The size of the array to sort.
! N.B. Use zero based indices. To obtain ! the entry into the actual
! array, multiply indices by nskip and add ! 1 (hence zero based)
! **** See local macro srt_ind() ******
lo = lbound(arr, 1) - 1
n = ((ubound(arr, 1) - lo -1)/nskip_) + 1
hi = lo + n - 1

nstack = 0
do while (.true.)
! If the section/partition we are looking at is smaller than
! nInsertSort then perform an insertion sort
if (hi - lo < nInsertionSort) then
do j = lo + 1, hi
a = arr(srt_ind(j))
a2 = arr2(srt_ind(j))
! a3 = arr3(srt_ind(j))
! a4 = arr4(srt_ind(j))
! a5 = arr5(srt_ind(j))
! a6 = arr6(srt_ind(j))
do i = j - 1, 0, -1
if (arr(srt_ind(i)) < a) exit
arr(srt_ind(i+1)) = arr(srt_ind(i))
arr2(srt_ind(i+1)) = arr2(srt_ind(i))
! arr3(srt_ind(i+1)) = arr3(srt_ind(i))
! arr4(srt_ind(i+1)) = arr4(srt_ind(i))
! arr5(srt_ind(i+1)) = arr5(srt_ind(i))
! arr6(srt_ind(i+1)) = arr6(srt_ind(i))
par_ = -par_
enddo
arr(srt_ind(i+1)) = a
arr2(srt_ind(i+1)) = a2
! arr3(srt_ind(i+1)) = a3
! arr4(srt_ind(i+1)) = a4
! arr5(srt_ind(i+1)) = a5
! arr6(srt_ind(i+1)) = a6
enddo

if (nstack == 0) exit
hi = stack(nstack)
lo = stack(nstack-1)
nstack = nstack - 2

! Otherwise start partitioning with quicksort.
else
! Pick a partitioning element, and arrange such that
! arr(lo) <= arr(lo+1) <= arr(hi)
pivot = (lo + hi) / 2
tmp1 = arr(srt_ind(pivot))
arr(srt_ind(pivot)) = arr(srt_ind(lo + 1))
arr(srt_ind(lo + 1)) = tmp1
tmp2 = arr2(srt_ind(pivot))
arr2(srt_ind(pivot)) = arr2(srt_ind(lo+1))
arr2(srt_ind(lo+1)) = tmp2
! tmp3 = arr3(srt_ind(pivot))
! arr3(srt_ind(pivot)) = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = tmp3
! tmp4 = arr4(srt_ind(pivot))
! arr4(srt_ind(pivot)) = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = tmp4
! tmp5 = arr5(srt_ind(pivot))
! arr5(srt_ind(pivot)) = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = tmp5
! tmp6 = arr6(srt_ind(pivot))
! arr6(srt_ind(pivot)) = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = tmp6
par_ = -par_

if (arr(srt_ind(lo)) > arr(srt_ind(hi))) then
tmp1 = arr(srt_ind(lo))
arr(srt_ind(lo)) = arr(srt_ind(hi))
arr(srt_ind(hi)) = tmp1
tmp2 = arr2(srt_ind(lo))
arr2(srt_ind(lo)) = arr2(srt_ind(hi))
arr2(srt_ind(hi)) = tmp2
! tmp3 = arr3(srt_ind(lo))
! arr3(srt_ind(lo)) = arr3(srt_ind(hi))
! arr3(srt_ind(hi)) = tmp3
! tmp4 = arr4(srt_ind(lo))
! arr4(srt_ind(lo)) = arr4(srt_ind(hi))
! arr4(srt_ind(hi)) = tmp4
! tmp5 = arr5(srt_ind(lo))
! arr5(srt_ind(lo)) = arr5(srt_ind(hi))
! arr5(srt_ind(hi)) = tmp5
! tmp6 = arr6(srt_ind(lo))
! arr6(srt_ind(lo)) = arr6(srt_ind(hi))
! arr6(srt_ind(hi)) = tmp6
par_ = -par_
endif
if (arr(srt_ind(lo+1)) > arr(srt_ind(hi))) then
tmp1 = arr(srt_ind(lo+1))
arr(srt_ind(lo+1)) = arr(srt_ind(hi))
arr(srt_ind(hi)) = tmp1
tmp2 = arr2(srt_ind(lo+1))
arr2(srt_ind(lo+1)) = arr2(srt_ind(hi))
arr2(srt_ind(hi)) = tmp2
! tmp3 = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = arr3(srt_ind(hi))
! arr3(srt_ind(hi)) = tmp3
! tmp4 = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = arr4(srt_ind(hi))
! arr4(srt_ind(hi)) = tmp4
! tmp5 = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = arr5(srt_ind(hi))
! arr5(srt_ind(hi)) = tmp5
! tmp6 = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = arr6(srt_ind(hi))
! arr6(srt_ind(hi)) = tmp6
par_ = -par_
endif
if (arr(srt_ind(lo)) > arr(srt_ind(lo+1))) then
tmp1 = arr(srt_ind(lo))
arr(srt_ind(lo)) = arr(srt_ind(lo+1))
arr(srt_ind(lo+1)) = tmp1
tmp2 = arr2(srt_ind(lo))
arr2(srt_ind(lo)) = arr2(srt_ind(lo+1))
arr2(srt_ind(lo+1)) = tmp2
! tmp3 = arr3(srt_ind(lo))
! arr3(srt_ind(lo)) = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = tmp3
! tmp4 = arr4(srt_ind(lo))
! arr4(srt_ind(lo)) = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = tmp4
! tmp5 = arr5(srt_ind(lo))
! arr5(srt_ind(lo)) = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = tmp5
! tmp6 = arr6(srt_ind(lo))
! arr6(srt_ind(lo)) = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = tmp6
par_ = -par_
endif

i = lo + 1
j = hi
a = arr(srt_ind(lo + 1)) !! a is the pivot value
a2 = arr2(srt_ind(lo + 1))
! a3 = arr3(srt_ind(lo + 1))
! a4 = arr4(srt_ind(lo + 1))
! a5 = arr5(srt_ind(lo + 1))
! a6 = arr6(srt_ind(lo + 1))
do while (.true.)
! Scand down list to find element > a
i = i + 1
do while (arr(srt_ind(i)) < a)
i = i + 1
enddo

! Scan down list to find element < a
j = j - 1
do while (arr(srt_ind(j)) > a)
j = j - 1
enddo

! When the pointers crossed, partitioning is complete.
if (j < i) exit

! Swap the elements, so that all elements < a end up
! in lower indexed variables.
tmp1 = arr(srt_ind(i))
arr(srt_ind(i)) = arr(srt_ind(j))
arr(srt_ind(j)) = tmp1
tmp2 = arr2(srt_ind(i))
arr2(srt_ind(i)) = arr2(srt_ind(j))
arr2(srt_ind(j)) = tmp2
! tmp3 = arr3(srt_ind(i))
! arr3(srt_ind(i)) = arr3(srt_ind(j))
! arr3(srt_ind(j)) = tmp3
! tmp4 = arr4(srt_ind(i))
! arr4(srt_ind(i)) = arr4(srt_ind(j))
! arr4(srt_ind(j)) = tmp4
! tmp5 = arr5(srt_ind(i))
! arr5(srt_ind(i)) = arr5(srt_ind(j))
! arr5(srt_ind(j)) = tmp5
! tmp6 = arr6(srt_ind(i))
! arr6(srt_ind(i)) = arr6(srt_ind(j))
! arr6(srt_ind(j)) = tmp6
par_ = -par_
enddo;

! Insert partitioning element
arr(srt_ind(lo + 1)) = arr(srt_ind(j))
arr(srt_ind(j)) = a
arr2(srt_ind(lo + 1)) = arr2(srt_ind(j))
! arr3(srt_ind(lo + 1)) = arr3(srt_ind(j))
! arr4(srt_ind(lo + 1)) = arr4(srt_ind(j))
! arr5(srt_ind(lo + 1)) = arr5(srt_ind(j))
! arr6(srt_ind(lo + 1)) = arr6(srt_ind(j))
arr2(srt_ind(j)) = a2
! arr3(srt_ind(j)) = a3
! arr4(srt_ind(j)) = a4
! arr5(srt_ind(j)) = a5
! arr6(srt_ind(j)) = a6
par_ = -par_

! Push the larger of the partitioned sections onto the stack
! of sections to look at later.
! --> need fewest stack elements.
nstack = nstack + 2
if (nstack > nStackMax) then
call stop_all (this_routine, &
"parameter nStackMax too small")
endif
if (hi - i + 1 >= j - lo) then
stack(nstack) = hi
stack(nstack-1) = i
hi = j - 1
else
stack(nstack) = j - 1
stack(nstack-1) = lo
lo = i
endif
endif
enddo

! Destroy temporary variables if required
!(a)
!(tmp1)
!(a2)
!(tmp2)
!(a3)
!(tmp3)
!(a4)
!(tmp4)
!(a5)
!(tmp5)
!(a6)
!(tmp6)

! Return the parity if required
if (present(par)) &
par = par_

end subroutine

! A private comparison. This should not be available outside of the
! module!
elemental function cmplx_gt_i_i (a, b) result (bGt)
complex(dp), intent(in) :: a, b
logical :: bGt

bGt = real(a, dp) > real(b, dp)
end function

elemental function cmplx_lt_i_i (a, b) result (bLt)
complex(dp), intent(in) :: a, b
logical :: bLt

bLt = real(a, dp) < real(b, dp)
end function

end module

#endif

#if !defined(SX)

#define srt_ind(i) (((i)*nskip_)+1)

module sort_mod_i_d
use util_mod, only: stop_all
use util_mod_comparisons, only: operator(.arrgt.), operator(.arrlt.)

use constants
use SystemData, only: Symmetry, operator(.eq.), operator(.ne.), &
operator(.gt.), operator(.lt.), assignment(=)
use symdata, only: SymPairProd, operator(.eq.), operator(.ne.), &
operator(.gt.), operator(.lt.), assignment(=)
implicit none

private
public :: sort

! Private operator for sorting complex numbers by their real values
interface operator(.gt.)
module procedure cmplx_gt_i_d
end interface

interface operator(.lt.)
module procedure cmplx_lt_i_d
end interface

interface sort
module procedure sort_i_d
end interface

interface cmplx_gt
module procedure cmplx_gt_i_d
end interface

interface cmplx_lt
module procedure cmplx_lt_i_d
end interface

contains
pure subroutine sort_i_d (arr, arr2,nskip, par)

! Perform a quicksort on an array of integers, arr(n). Uses the
! sample code in NumericalRecipies as a base.
! Optionally sort arr2 in parallel (in the routines it is enabled)

! Custom comparisons. Use the operator .locallt. & .localgt.
! interface operator(.locallt.)
!     pure function custom_lt (b, c) result (ret)
!         use constants
!         integer(kind=int32), intent(in) :: b(:)
!         integer(kind=int32), intent(in) :: c(:)
!         logical :: ret
!     end function
! end interface

! interface operator(.localgt.)
!     pure function custom_gt (b, c) result (ret)
!         use constants
!         integer(kind=int32), intent(in) :: b(:)
!         integer(kind=int32), intent(in) :: c(:)
!         logical :: ret
!     end function
! end interface

! sort needs auxiliary storage of length 2*log_2(n).
integer, parameter :: nStackMax = 50
integer, parameter :: nInsertionSort = 7
integer, intent(in), optional :: nskip
integer, intent(out), optional :: par

integer(kind=int32), intent(inout) :: arr(:)
real(dp), intent(inout) :: arr2(:)
!, intent(inout) :: arr3(:)
!, intent(inout) :: arr4(:)
!, intent(inout) :: arr5(:)
!, intent(inout) :: arr6(:)

! Oh, how lovely it would be to be able to use push/pop and not need
integer :: stack(nStackMax), nstack, nskip_
integer :: pivot, lo, hi, n, i, j, par_
! n.b. This size statement is removed if type1 is scalar ...
integer(kind=int32) :: a
real(dp) :: a2
! :: a3(size(arr3(1)))
! :: a4(size(arr4(1)))
! :: a5(size(arr5(1)))
! :: a6(size(arr6(1)))

! Temporary variables for swapping
integer(kind=int32) :: tmp1
real(dp) :: tmp2
! :: tmp3(size(arr3(1)))
! :: tmp4(size(arr4(1)))
! :: tmp5(size(arr5(1)))
! :: tmp6(size(arr6(1)))
character(*), parameter :: this_routine = 'sort'

! Initialise temporary variables if required
!(a)
!(tmp1)
!(a2)
!(tmp2)
!(a3)
!(tmp3)
!(a4)
!(tmp4)
!(a5)
!(tmp5)
!(a6)
!(tmp6)

! *** HACK ***
! Workaround for gfortran compiler bug
! n.b. This will produce spurious warnings if run in valgrind, as
!      a is not initialised. Not a problem in optimised build.
i = 0
#ifdef DEBUG_
! it IS a problem in debug build, however
! *** HACK MOD ***
!to prevent crashes in debug mode
if((ubound(arr,1) - lbound(arr,1) + 1) &
> 0 ) a = arr(1)
#endif
! if (custom_lt(a, a)) i = i
! if (custom_gt(a, a)) i = i

if (present(nskip)) then
nskip_ = nskip
else
nskip_ = 1
endif

! Initialise parity calculation
par_ = 1

! The size of the array to sort.
! N.B. Use zero based indices. To obtain ! the entry into the actual
! array, multiply indices by nskip and add ! 1 (hence zero based)
! **** See local macro srt_ind() ******
lo = lbound(arr, 1) - 1
n = ((ubound(arr, 1) - lo -1)/nskip_) + 1
hi = lo + n - 1

nstack = 0
do while (.true.)
! If the section/partition we are looking at is smaller than
! nInsertSort then perform an insertion sort
if (hi - lo < nInsertionSort) then
do j = lo + 1, hi
a = arr(srt_ind(j))
a2 = arr2(srt_ind(j))
! a3 = arr3(srt_ind(j))
! a4 = arr4(srt_ind(j))
! a5 = arr5(srt_ind(j))
! a6 = arr6(srt_ind(j))
do i = j - 1, 0, -1
if (arr(srt_ind(i)) < a) exit
arr(srt_ind(i+1)) = arr(srt_ind(i))
arr2(srt_ind(i+1)) = arr2(srt_ind(i))
! arr3(srt_ind(i+1)) = arr3(srt_ind(i))
! arr4(srt_ind(i+1)) = arr4(srt_ind(i))
! arr5(srt_ind(i+1)) = arr5(srt_ind(i))
! arr6(srt_ind(i+1)) = arr6(srt_ind(i))
par_ = -par_
enddo
arr(srt_ind(i+1)) = a
arr2(srt_ind(i+1)) = a2
! arr3(srt_ind(i+1)) = a3
! arr4(srt_ind(i+1)) = a4
! arr5(srt_ind(i+1)) = a5
! arr6(srt_ind(i+1)) = a6
enddo

if (nstack == 0) exit
hi = stack(nstack)
lo = stack(nstack-1)
nstack = nstack - 2

! Otherwise start partitioning with quicksort.
else
! Pick a partitioning element, and arrange such that
! arr(lo) <= arr(lo+1) <= arr(hi)
pivot = (lo + hi) / 2
tmp1 = arr(srt_ind(pivot))
arr(srt_ind(pivot)) = arr(srt_ind(lo + 1))
arr(srt_ind(lo + 1)) = tmp1
tmp2 = arr2(srt_ind(pivot))
arr2(srt_ind(pivot)) = arr2(srt_ind(lo+1))
arr2(srt_ind(lo+1)) = tmp2
! tmp3 = arr3(srt_ind(pivot))
! arr3(srt_ind(pivot)) = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = tmp3
! tmp4 = arr4(srt_ind(pivot))
! arr4(srt_ind(pivot)) = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = tmp4
! tmp5 = arr5(srt_ind(pivot))
! arr5(srt_ind(pivot)) = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = tmp5
! tmp6 = arr6(srt_ind(pivot))
! arr6(srt_ind(pivot)) = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = tmp6
par_ = -par_

if (arr(srt_ind(lo)) > arr(srt_ind(hi))) then
tmp1 = arr(srt_ind(lo))
arr(srt_ind(lo)) = arr(srt_ind(hi))
arr(srt_ind(hi)) = tmp1
tmp2 = arr2(srt_ind(lo))
arr2(srt_ind(lo)) = arr2(srt_ind(hi))
arr2(srt_ind(hi)) = tmp2
! tmp3 = arr3(srt_ind(lo))
! arr3(srt_ind(lo)) = arr3(srt_ind(hi))
! arr3(srt_ind(hi)) = tmp3
! tmp4 = arr4(srt_ind(lo))
! arr4(srt_ind(lo)) = arr4(srt_ind(hi))
! arr4(srt_ind(hi)) = tmp4
! tmp5 = arr5(srt_ind(lo))
! arr5(srt_ind(lo)) = arr5(srt_ind(hi))
! arr5(srt_ind(hi)) = tmp5
! tmp6 = arr6(srt_ind(lo))
! arr6(srt_ind(lo)) = arr6(srt_ind(hi))
! arr6(srt_ind(hi)) = tmp6
par_ = -par_
endif
if (arr(srt_ind(lo+1)) > arr(srt_ind(hi))) then
tmp1 = arr(srt_ind(lo+1))
arr(srt_ind(lo+1)) = arr(srt_ind(hi))
arr(srt_ind(hi)) = tmp1
tmp2 = arr2(srt_ind(lo+1))
arr2(srt_ind(lo+1)) = arr2(srt_ind(hi))
arr2(srt_ind(hi)) = tmp2
! tmp3 = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = arr3(srt_ind(hi))
! arr3(srt_ind(hi)) = tmp3
! tmp4 = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = arr4(srt_ind(hi))
! arr4(srt_ind(hi)) = tmp4
! tmp5 = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = arr5(srt_ind(hi))
! arr5(srt_ind(hi)) = tmp5
! tmp6 = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = arr6(srt_ind(hi))
! arr6(srt_ind(hi)) = tmp6
par_ = -par_
endif
if (arr(srt_ind(lo)) > arr(srt_ind(lo+1))) then
tmp1 = arr(srt_ind(lo))
arr(srt_ind(lo)) = arr(srt_ind(lo+1))
arr(srt_ind(lo+1)) = tmp1
tmp2 = arr2(srt_ind(lo))
arr2(srt_ind(lo)) = arr2(srt_ind(lo+1))
arr2(srt_ind(lo+1)) = tmp2
! tmp3 = arr3(srt_ind(lo))
! arr3(srt_ind(lo)) = arr3(srt_ind(lo+1))
! arr3(srt_ind(lo+1)) = tmp3
! tmp4 = arr4(srt_ind(lo))
! arr4(srt_ind(lo)) = arr4(srt_ind(lo+1))
! arr4(srt_ind(lo+1)) = tmp4
! tmp5 = arr5(srt_ind(lo))
! arr5(srt_ind(lo)) = arr5(srt_ind(lo+1))
! arr5(srt_ind(lo+1)) = tmp5
! tmp6 = arr6(srt_ind(lo))
! arr6(srt_ind(lo)) = arr6(srt_ind(lo+1))
! arr6(srt_ind(lo+1)) = tmp6
par_ = -par_
endif

i = lo + 1
j = hi
a = arr(srt_ind(lo + 1)) !! a is the pivot value
a2 = arr2(srt_ind(lo + 1))
! a3 = arr3(srt_ind(lo + 1))
! a4 = arr4(srt_ind(lo + 1))
! a5 = arr5(srt_ind(lo + 1))
! a6 = arr6(srt_ind(lo + 1))
do while (.true.)
! Scand down list to find element > a
i = i + 1
do while (arr(srt_ind(i)) < a)
i = i + 1
enddo

! Scan down list to find element < a
j = j - 1
do while (arr(srt_ind(j)) > a)
j = j - 1
enddo

! When the pointers crossed, partitioning is complete.
if (j < i) exit

! Swap the elements, so that all elements < a end up
! in lower indexed variables.
tmp1 = arr(srt_ind(i))
arr(srt_ind(i)) = arr(srt_ind(j))
arr(srt_ind(j)) = tmp1
tmp2 = arr2(srt_ind(i))
arr2(srt_ind(i)) = arr2(srt_ind(j))
arr2(srt_ind(j)) = tmp2
! tmp3 = arr3(srt_ind(i))
! arr3(srt_ind(i)) = arr3(srt_ind(j))
! arr3(srt_ind(j)) = tmp3
! tmp4 = arr4(srt_ind(i))
! arr4(srt_ind(i)) = arr4(srt_ind(j))
! arr4(srt_ind(j)) = tmp4
! tmp5 = arr5(srt_ind(i))
! arr5(srt_ind(i)) = arr5(srt_ind(j))
! arr5(srt_ind(j)) = tmp5
! tmp6 = arr6(srt_ind(i))
! arr6(srt_ind(i)) = arr6(srt_ind(j))
! arr6(srt_ind(j)) = tmp6
par_ = -par_
enddo;

! Insert partitioning element
arr(srt_ind(lo + 1)) = arr(srt_ind(j))
arr(srt_ind(j)) =```