quicksort.F90 Source File

#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
        ! to guess a size of the stack to start with
        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
        ! to guess a size of the stack to start with
        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
        ! to guess a size of the stack to start with
        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
        ! to guess a size of the stack to start with
        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
        ! to guess a size of the stack to start with
        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
        ! to guess a size of the stack to start with
        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
        ! to guess a size of the stack to start with
        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
        ! to guess a size of the stack to start with
        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
        ! to guess a size of the stack to start with
        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
        ! to guess a size of the stack to start with
        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
        ! to guess a size of the stack to start with
        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
        ! to guess a size of the stack to start with
        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
        ! to guess a size of the stack to start with
        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)) =