!******************************************************************************
!
! reference
!
! http://www.ics.kagoshima-u.ac.jp/~fuchida/edu/algorithm/sort-algorithm/
!
! Press et al., Numerical Recipes in C (Japanese version), 1993
!
!******************************************************************************

  module sort_module

    use vtype_module

    implicit none

    private


    public :: sort_quick

    interface
       subroutine sort_quick( im, arr, arr1 )
         use vtype_module
         implicit none
         integer(i4b), intent(in   )           :: im
         real(dp)    , intent(inout)           :: arr ( im )
         real(dp)    , intent(inout), optional :: arr1( im )
       end subroutine sort_quick
    end interface


  end module sort_module

  !****************************************************************************

    subroutine sort_heap( im, arr )

      use vtype_module

      implicit none

      integer(i4b), intent(in   ) :: im
      real(dp)    , intent(inout) :: arr( im )


      !
      ! local variables
      !
      integer(i4b) :: iend


      !
      ! "hiring" phase
      !
      iend = im
      call mkheap( im, arr, iend/2, iend )
      !
      ! "retirement and promotion" phase
      !
      do iend = im, 2, -1
         ! "retirement" phase
         ! The smallest number is "retired" (moved to the end of array).
         call swap( im, arr, 1, iend )
         ! "promotion" phase
         call mkheap( im, arr, 1, iend-1 )
      end do


      !************************************************************************

    contains

      !************************************************************************

      subroutine mkheap( im, arr, istart, iend )

        integer(i4b), intent(in   ) :: im
        real(dp)    , intent(inout) :: arr( im )
        integer(i4b), intent(in   ) :: istart, iend


        !
        ! local variables
        !
        integer(i4b) :: i, j, k


        if( iend .le. 1 ) return


!!$        i = iend / 2
        i = istart
        j = i
        do
           k = 2 * j

           if( k .le. iend ) then

              if( ( k + 1 ) .le. iend ) then
                 if( arr( k + 1 ) .gt. arr( k ) ) then
                    k = k + 1
                 end if
              end if

              if( arr( k ) .gt. arr( j ) ) then
                 call swap( im, arr, j, k )
                 j = k
              else
                 i = i - 1
                 j = i
              end if

           else
              i = i - 1
              j = i
           end if

           if( i .eq. 0 ) exit

        end do

      end subroutine mkheap

      !************************************************************************

      subroutine swap( im, arr, i1, i2 )

        integer(i4b), intent(in   ) :: im
        real(dp)    , intent(inout) :: arr( im )
        integer(i4b), intent(in   ) :: i1, i2


        !
        ! local variables
        !
        real(dp) :: rarr


        rarr      = arr( i1 )
        arr( i1 ) = arr( i2 )
        arr( i2 ) = rarr


      end subroutine swap

      !************************************************************************

    end subroutine sort_heap

    !**************************************************************************

    subroutine sort_quick( im, arr, arr1 )

      use vtype_module

      implicit none

      integer(i4b), intent(in   )           :: im
      real(dp)    , intent(inout)           :: arr ( im )
      real(dp)    , intent(inout), optional :: arr1( im )


      !
      ! local varialbes
      !
      integer(i4b) :: is, ie

      is = 1
      ie = im
      call sort_quick0( im, is, ie, arr, arr1 )


      !************************************************************************

    contains

      !************************************************************************

      recursive subroutine sort_quick0( im, is, ie, arr, arr1 )

        integer(i4b), intent(in   )           :: im
        integer(i4b), intent(in   )           :: is, ie
        real(dp)    , intent(inout)           :: arr ( im )
        real(dp)    , intent(inout), optional :: arr1( im )


        !
        ! local varialbes
        !
        real(dp)     :: pivot
        integer(i4b) :: i1, i2, is1, ie1, is2, ie2

!!$        logical      :: sw_find1, sw_find2


        if( is .eq. ie ) return


        !
        ! select pivot
        !
        do i1 = is+1, ie
           if( arr( is ) .ne. arr( i1 ) ) exit
        end do
        if( i1 .gt. ie ) return
        if( arr( is ) .ge. arr( i1 ) ) then
           pivot = arr( is )
        else
           pivot = arr( i1 )
        end if



!!$        sw_find1 = .false. ; sw_find2 = .false.
!!$
!!$        i1 = is
!!$        i2 = ie
!!$
!!$        search_loop : do
!!$           if( arr( i1 ) .ge. pivot ) sw_find1 = .true.
!!$           if( arr( i2 ) .lt. pivot ) sw_find2 = .true.
!!$
!!$           if( ( sw_find1 ) .and. ( sw_find2 ) ) then
!!$              call swap( im, arr, i1, i2 )
!!$              sw_find1 = .false. ; sw_find2 = .false.
!!$           end if
!!$
!!$
!!$!           if( .not. sw_find1 ) i1 = i1 + 1
!!$!           if( i1 .eq. i2 ) then
!!$!              if( arr( i1 ) .ge. pivot ) then
!!$!                 ie1 = i1 - 1
!!$!              else
!!$!                 ie1 = i1
!!$!              end if
!!$!              exit search_loop
!!$!           end if
!!$!
!!$!           if( .not. sw_find2 ) i2 = i2 - 1
!!$!           if( i1 .eq. i2 ) then
!!$!              if( arr( i1 ) .ge. pivot ) then
!!$!                 ie1 = i1 - 1
!!$!              else
!!$!                 ie1 = i1
!!$!              end if
!!$!              exit search_loop
!!$!           end if
!!$
!!$
!!$
!!$           if( .not. sw_find1 ) i1 = i1 + 1
!!$           if( .not. sw_find2 ) i2 = i2 - 1
!!$           if( i1 .ge. i2 ) then
!!$              if( arr( i1 ) .ge. pivot ) then
!!$                 ie1 = i1 - 1
!!$              else
!!$                 ie1 = i1
!!$              end if
!!$              exit search_loop
!!$           end if
!!$
!!$        end do search_loop



        i1 = is
        i2 = ie

        do

           do
              if( arr( i1 ) .ge. pivot ) exit
              i1 = i1 + 1
              ! MEMO:
              ! I never fail to find a value greater than or equal to pivot.
!!$              if( i1 .ge. ie ) exit
           end do
           do
              if( arr( i2 ) .lt. pivot ) exit
              i2 = i2 - 1
              ! MEMO:
              ! I never fail to find a value less than or equal to pivot.
              ! This is due to a way to determine pivot.
              ! A value different from pivot MUST be included in array.
!!$              if( i2 .le. is ) exit
           end do

           ! MEMO:
           ! The values i1 and i2 must not be the same. 
           ! This is due to a way to determine pivot.
           ! A value different from pivot MUST be included in array.
!!$           if( i1 .ge. i2 ) then
           if( i1 .gt. i2 ) then

              ! MEMO:
              ! The value arr( i1 ) must not be greater than or equal to pivot.
              ! This is due to a way to determine pivot.
              ! A value different from pivot MUST be included in array.
!!$              if( arr( i1 ) .ge. pivot ) then
!!$                 ie1 = i1 - 1
!!$              else
!!$                 ie1 = i1
!!$              end if

              ie1 = i1 - 1

              exit
           end if

           call swap( im, arr, i1, i2 )
           if( present( arr1 ) ) call swap( im, arr1, i1, i2 )


        end do


        is1 = is
!!$        ie1 = i1 - 1
        call sort_quick0( im, is1, ie1, arr, arr1 )
        is2 = ie1 + 1
        ie2 = ie
        call sort_quick0( im, is2, ie2, arr, arr1 )


      end subroutine sort_quick0

      !************************************************************************

      subroutine swap( im, arr, i1, i2 )

        integer(i4b), intent(in   ) :: im
        real(dp)    , intent(inout) :: arr( im )
        integer(i4b), intent(in   ) :: i1, i2


        !
        ! local variables
        !
        real(dp) :: rarr


        rarr      = arr( i1 )
        arr( i1 ) = arr( i2 )
        arr( i2 ) = rarr


      end subroutine swap

      !************************************************************************

    end subroutine sort_quick
