psmile_quicksort_index.F90

Go to the documentation of this file.

!-----------------------------------------------------------------------
! Copyright 2006-2010, NEC Europe Ltd., London, UK.
! All rights reserved. Use is subject to OASIS4 license terms.
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: PSMILe_Quicksort_index
!
! !INTERFACE:

      subroutine psmile_quicksort_index(a, n, index)
!
! !USES:
!
        use PSMILe, dummy_interface => PSMILe_Quicksort_index

        implicit none
!
! !INPUT PARAMETERS:
!
        integer, intent(In)    :: n
!
! !INPUT/OUTPUT PARAMETERS:
!
        integer, intent(InOut) :: a(n)
        integer, intent(InOut) :: index(n)
!
! !LOCAL VARIABLES
!
        integer, Parameter     :: pointer_inc = 64
        integer                :: pointer_size
        integer, pointer       :: stackl(:), stackr(:)
        integer, pointer       :: new_stackl(:), new_stackr(:)

        integer                :: i, j, k, l, r, s
        integer                :: ww
        integer                :: w, x
!
! !DESCRIPTION:
!
! Non-recursive stack version of Quicksort from N. Wirth's Pascal Book,
! 'Algorithms + Data Structures = Programms'.
!
! taken from:
!    http://www.nag.com/nagware/examples.asp
!    http://www.nag.com/nagware/Examples/nur.f90
!
! see also:
!    http://en.wikipedia.org/wiki/Quicksort
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 19.07.95    Alan Miller  created
! 13.02.08    R. Redler    revised
!
!----------------------------------------------------------------------
!
!  Initialization
!
#ifdef VERBOSE
        print 9990, trim(ch_id)

        call psmile_flushstd
#endif /* VERBOSE */

        pointer_size = pointer_inc

        allocate(new_stackl(pointer_inc), new_stackr(pointer_inc))

        stackl => new_stackl
        stackr => new_stackr

        s = 1
        stackl(1) = 1
        stackr(1) = n
        !
        !  Start sorting
        !        
        ! ... keep taking the top request from the stack until s = 0.

10      continue
        l = stackl(s)
        r = stackr(s)
        s = s - 1

        ! ... keep splitting a(l), ... ,a(r) until l>= r.

20      continue
        i = l
        j = r
        k = (l+r) / 2
        x = a(k)

        ! ... repeat until i > j.

        do
           do
              if (a(i) < x) then ! Search from lower end
                 i = i + 1
                 cycle
              else
                 exit
              end if
           end do

           do
              if (x < a(j)) then ! Search from upper end
                 j = j - 1
                 cycle
              else
                 exit
              end if
           end do

           if (i <= j) then ! Swap positions i & j
              w = a(i)
              ww = index(i)
              a(i) = a(j)
              index(i) = index(j)
              a(j) = w
              index(j) = ww
              i = i + 1
              j = j - 1
              if (i.gt.j) exit
           else
              exit
           end if
        end do

        if (j-l >= r-i) then

           if (l < j) then

              if ( s+1 > pointer_size ) then
                 pointer_size = pointer_size + pointer_inc
                 allocate(new_stackl(pointer_size), new_stackr(pointer_size))
                 new_stackl(1:pointer_size-pointer_inc) = &
                     stackl(1:pointer_size-pointer_inc)
                 new_stackr(1:pointer_size-pointer_inc) = &
                     stackr(1:pointer_size-pointer_inc)
                 deallocate(stackl, stackr)
                 stackl => new_stackl
                 stackr => new_stackr
              endif

              s = s + 1
              stackl(s) = l
              stackr(s) = j
           end if
           l = i

        else

           if (i < r) then

              if ( s+1 > pointer_size ) then
                 pointer_size = pointer_size + pointer_inc
                 allocate(new_stackl(pointer_size), new_stackr(pointer_size))
                 new_stackl(1:pointer_size-pointer_inc) = &
                     stackl(1:pointer_size-pointer_inc)
                 new_stackr(1:pointer_size-pointer_inc) = &
                     stackr(1:pointer_size-pointer_inc)
                 deallocate(stackl, stackr)
                 stackl => new_stackl
                 stackr => new_stackr
              endif

              s = s + 1
              stackl(s) = i
              stackr(s) = r
           end if
           r = j

        end if

        if (l <  r) GO TO 20
        if (s /= 0) GO TO 10

#ifdef VERBOSE
        print 9980, trim(ch_id), pointer_size

        call psmile_flushstd
#endif /* VERBOSE */
        !
        !  Formats:
        !
9990    format (1x, a, ': psmile_quicksort_index')
9980    format (1x, a, ': psmile_quicksort_index: eof, stack size was ', i8)

      end subroutine PSMILe_Quicksort_index