dqsort Subroutine

RECURSIVE SUBROUTINE dqsort(list, order)
    !!  Sorts a sequence of reals
    !""

IMPLICIT NONE

REAL(RP), DIMENSION (:), INTENT(IN OUT)  :: list
    !!  Sequence of reals to be sorted
INTEGER, DIMENSION (:), INTENT(OUT), OPTIONAL  :: order
    !!  Indices of the sorted sequence

! Local variable
INTEGER :: i

IF (PRESENT(order)) THEN
    DO i = 1, SIZE(list)
      order(i) = i
    END DO
END IF

CALL quick_sort_1(1, SIZE(list))

CONTAINS

    RECURSIVE SUBROUTINE quick_sort_1(left_end, right_end)
    
    INTEGER, INTENT(IN) :: left_end, right_end
    
    !     Local variables
    INTEGER             :: i, j, itemp
    REAL(RP)            :: reference, temp
    INTEGER, PARAMETER  :: max_simple_sort_size = 8
    
    IF (right_end < left_end + max_simple_sort_size) THEN
        ! Use interchange sort for small lists
        CALL interchange_sort(left_end, right_end)
    
    ELSE
        ! Use partition ("quick") sort
        reference = list((left_end + right_end)/2)
        i = left_end - 1; j = right_end + 1
    
        DO
            ! Scan list from left end until element >= reference is found
            DO
                i = i + 1
                IF (list(i) >= reference) EXIT
            END DO
            ! Scan list from right end until element <= reference is found
            DO
                j = j - 1
                IF (list(j) <= reference) EXIT
            END DO
    
            IF (i < j) THEN
                ! Swap two out-of-order elements
                temp = list(i); list(i) = list(j); list(j) = temp
                IF (PRESENT(order)) THEN
                    itemp = order(i); order(i) = order(j); order(j) = itemp
                END IF
            ELSE IF (i == j) THEN
                i = i + 1
                EXIT
            ELSE
                EXIT
            END IF
        END DO
    
        IF (left_end < j) CALL quick_sort_1(left_end, j)
        IF (i < right_end) CALL quick_sort_1(i, right_end)
    END IF
    
    END SUBROUTINE quick_sort_1


    SUBROUTINE interchange_sort(left_end, right_end)
    
    INTEGER, INTENT(IN) :: left_end, right_end
    
    !     Local variables
    INTEGER             :: i, j, itemp
    REAL(RP)            :: temp
    
    DO i = left_end, right_end - 1
        DO j = i+1, right_end
            IF (list(i) > list(j)) THEN
                temp = list(i); list(i) = list(j); list(j) = temp
                IF (PRESENT(order)) THEN
                    itemp = order(i); order(i) = order(j); order(j) = itemp
                END IF
            END IF
        END DO
    END DO
    
    END SUBROUTINE interchange_sort

END SUBROUTINE dqsort