PROGRAM Main
    IMPLICIT NONE
    INTEGER, PARAMETER :: dp = KIND(0.0D0), max_bubble = 15,    &
        min_task = 1000
    CHARACTER(LEN=100) :: buffer
    INTEGER :: length, k, n
    REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: vector
    REAL(KIND=dp) :: check, tol

    IF (COMMAND_ARGUMENT_COUNT() /= 1) THEN
        PRINT *, "Wrong number of arguments"
        STOP 1
    END IF
    CALL GET_COMMAND_ARGUMENT(1,buffer)
    READ (buffer,"(I100)",IOSTAT=k) length
    IF (k /= 0 .OR. length < 10) THEN
        PRINT '("Invalid number of elements: ",A)', buffer
        STOP 
    END IF
    ALLOCATE (vector(length))
    check = 0.0_dp
    CALL RANDOM_NUMBER(vector)
    check = SUM(vector)
    PRINT "(I0,' elements with sum ',F0.3)",length, check

!$OMP PARALLEL DEFAULT(NONE), SHARED(vector,k)
!$OMP MASTER
    CALL Quicksort(vector,k)
!$OMP END MASTER
!$OMP END PARALLEL
    IF (k /= 0) THEN
        PRINT "('Quicksort failed: ',I0)", k
        STOP 1
    END IF
    DO n = 2,length
        IF (vector(n) < vector(n-1)) THEN
            PRINT "('Element ',I0,' out of order: ',F0.3,1X,F0.3)",    &
                n, vector(n), vector(n-1)
            STOP 1
        END IF
    END DO
    check = check-SUM(vector)
    tol = 1.0e-12_dp*SQRT(REAL(length,KIND=dp)**3/20)
    IF (ABS(check) > tol) THEN
        PRINT *, "Contents of vector have been mangled"
        STOP 1
    END IF

CONTAINS


    SUBROUTINE Quicksort (array, code)

! This sorts a vector of doubles using Maclaren's variation (CACM ....)
! of Quicksort (CACM 12, p. 185).  Most of the changes from the latter
! are for clarity and speed, at the expense of code.  It is considerably
! better than most implementations of qsort, and has a slightly enhanced
! specification.

! The variation mentioned is to note when a split is very uneven, and mark
! the larger section for special attention on the next pass.  It then
! exchanges the end elements with randomly chosen points before choosing
! the split point.  This adds very little expense, and makes O(N^2) cases
! vanishingly unlikely.  The quality of the random numbers does not matter
! much.

! It returns one of the following values:
!
!     0  -  sort successful
!     1  -  length of vector < 1

        REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: array
        INTEGER, INTENT(OUT) :: code

        code = 0
        if (length < 1) code = 1
        if (length < 2) RETURN
        CALL pass(array,.FALSE.)
    END SUBROUTINE Quicksort



    RECURSIVE SUBROUTINE pass (array, randomise)
        REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: array
        LOGICAL, INTENT(in) :: randomise
        INTEGER :: length, k, m, n
        REAL(KIND=dp) :: x, a

        length = UBOUND(array,1)

! Sort very small arrays with a bubble sort.

        IF (length <= max_bubble) THEN
            DO k = 2,length
                m = k
                IF (array(m-1) > array(m)) THEN
                    x = array(m)
                    DO
                       array(m) = array(m-1)
                       m = m-1                       
                       IF (m < 2) EXIT
                       IF (array(m-1) <= x) EXIT
                    END DO
                    array(m) = x
                END IF
            END DO
            RETURN
        END IF

! Set up for the quicksort.

        n = length
        m = 1

! If the previous pass was bad, exchange the end values with random
! elements.  Note that the generator cannot return 1.0, which avoids
! problems with the result going out of range.

        IF (randomise) THEN
            CALL RANDOM_NUMBER(x)
            k = MAX(INT(n*x)+1,n)
            x = array(m)
            array(m) = array(k)
            array(k) = x
            k = n-k+1
            x = array(n)
            array(n) = array(k)
            array(k) = x
        END IF

! Calculate a median by Singleton's method.

        k = n/2
        IF (array(m) > array(n)) THEN
            x = array(m)
            array(m) = array(n)
            array(n) = x
        END IF
        IF (array(k) < array(m)) THEN
            x = array(k)
            array(k) = array(m)
            array(m) = x
        ELSEIF (array(k) > array(n)) THEN
            x = array(k)
            array(k) = array(n)
            array(n) = x
        END IF
        a = array(k)
 
! Split the vector into two parts; this is where the time is spent.

        DO
            DO
                m = m+1
                IF (array(m) >= a) EXIT
            END DO
            DO
                n = n-1
                IF (array(n) <= a) EXIT
            END DO
            IF (m >= n) EXIT
            x = array(m)
            array(m) = array(n)
            array(n) = x
        END DO

! Sort the subfiles, by recursion.  The OpenMP clauses and use of
! TASKWAIT are utterly ridiculous, but are because OpenMP's model of
! array arguments is incompatible with Fortran's (and C's).  There is
! a way of avoiding the TASKWAIT, but it's even messier than the C
! code.  Contrast this with the C, which is equally disgusting, in a
! different way.

        IF (m-1 >= n+1) THEN
            PRINT *, "Internal inconsistency"
            STOP 1
        END IF
        k = length-n
        IF (m-1 < k) THEN
            IF (m-1 > min_task) THEN
!$OMP TASK DEFAULT(NONE), SHARED(array,m)
                CALL pass(array(:m-1),.FALSE.)
!$OMP END TASK
!$OMP TASK DEFAULT(NONE), SHARED(array,n,length,k,m)
                CALL pass(array(n+1:length),(k > 5*(m+5)))
!$OMP END TASK
!$OMP TASKWAIT
            ELSE
                IF (m-1 > 1) CALL pass(array(:m-1),.FALSE.)
                CALL pass(array(n+1:length),(k > 5*(m+5)))
            END IF 
        ELSE
            IF (k > min_task) THEN
!$OMP TASK DEFAULT(NONE), SHARED(array,n,length)
                CALL pass(array(n+1:length),.FALSE.)
!$OMP END TASK
!$OMP TASK DEFAULT(NONE), SHARED(array,m,k)
                CALL pass(array(:m-1),(m > 5*(k+5)))
!$OMP END TASK
!$OMP TASKWAIT
            ELSE
                IF (k > 1) CALL pass(array(n+1:length),.FALSE.)
                CALL pass(array(:m-1),(m > 5*(k+5)))
            END IF 
        END IF 
    END SUBROUTINE pass

END PROGRAM Main
