! Copyright (C) 2010-2016 Samuel Ponce', Roxana Margine, Carla Verdi, Feliciano Giustino

! Copyright (C) 2007-2009 Jesse Noffsinger, Brad Malone, Feliciano Giustino

!

! This file is distributed under the terms of the GNU General Public

! License. See the file `LICENSE' in the root directory of the

! present distribution, or http://www.gnu.org/copyleft.gpl.txt .

!

! Adapted from flib/hpsort_eps

!---------------------------------------------------------------------

subroutine hpsort_eps_epw (n, ra, ind, eps)

!---------------------------------------------------------------------

! sort an array ra(1:n) into ascending order using heapsort algorithm,

! and considering two elements being equal if their values differ

! for less than "eps".

! n is input, ra is replaced on output by its sorted rearrangement.

! create an index table (ind) by making an exchange in the index array

! whenever an exchange is made on the sorted data array (ra).

! in case of equal values in the data array (ra) the values in the

! index array (ind) are used to order the entries.

! if on input ind(1) = 0 then indices are initialized in the routine,

! if on input ind(1) != 0 then indices are assumed to have been

! initialized before entering the routine and these

! indices are carried around during the sorting process

!

! no work space needed !

! free us from machine-dependent sorting-routines !

!

! adapted from Numerical Recipes pg. 329 (new edition)

!

!use kinds, ONLY : DP

implicit none

integer, parameter :: DP=4 ! RJS !Mon Apr 13 20:10:53 EDT 2020

!-input/output variables

integer, intent(in) :: n

real(DP), intent(in) :: eps

integer :: ind (n)

real(DP) :: ra (n)

!-local variables

integer :: i, ir, j, l, iind

real(DP) :: rra

! initialize index array

IF (ind (1) .eq.0) then

DO i = 1, n

ind (i) = i

ENDDO

ENDIF

! nothing to order

IF (n.lt.2) return

! initialize indices for hiring and retirement-promotion phase

l = n / 2 + 1

ir = n

sorting: do

! still in hiring phase

IF ( l .gt. 1 ) then

l = l - 1

rra = ra (l)

iind = ind (l)

! in retirement-promotion phase.

ELSE

! clear a space at the end of the array

rra = ra (ir)

!

iind = ind (ir)

! retire the top of the heap into it

ra (ir) = ra (1)

!

ind (ir) = ind (1)

! decrease the size of the corporation

ir = ir - 1

! done with the last promotion

IF ( ir .eq. 1 ) then

! the least competent worker at all !

ra (1) = rra

!

ind (1) = iind

exit sorting

ENDIF

ENDIF

! wheter in hiring or promotion phase, we

i = l

! set up to place rra in its proper level

j = l + l

!

DO while ( j .le. ir )

IF ( j .lt. ir ) then

! compare to better underling

IF ( hslt( ra (j), ra (j + 1) ) ) then

j = j + 1

!else if ( .not. hslt( ra (j+1), ra (j) ) ) then

! this means ra(j) == ra(j+1) within tolerance

! if (ind (j) .lt.ind (j + 1) ) j = j + 1

ENDIF

ENDIF

! demote rra

IF ( hslt( rra, ra (j) ) ) then

ra (i) = ra (j)

ind (i) = ind (j)

i = j

j = j + j

!else if ( .not. hslt ( ra(j) , rra ) ) then

!this means rra == ra(j) within tolerance

! demote rra

! if (iind.lt.ind (j) ) then

! ra (i) = ra (j)

! ind (i) = ind (j)

! i = j

! j = j + j

! else

! set j to terminate do-while loop

! j = ir + 1

! endif

! this is the right place for rra

ELSE

! set j to terminate do-while loop

j = ir + 1

ENDIF

ENDDO

ra (i) = rra

ind (i) = iind

END DO sorting

contains

! internal function

! compare two real number and return the result

logical function hslt( a, b )

REAL(DP) :: a, b

IF( abs(a-b) < eps ) then

hslt = .false.

ELSE

hslt = ( a < b )

end if

end function hslt

!

end subroutine hpsort_eps_epw