sort.f90 2.86 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_sort
Matthias Redies's avatar
Matthias Redies committed
7
   USE m_judft
8 9
CONTAINS

Matthias Redies's avatar
Matthias Redies committed
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
   SUBROUTINE sort(ind,lv,lv1)
      !********************************************************************
      !     heapsort routine
      !     input:   lv    = array of objects to be sorted
      !              lv1   = second array to use as secondary sort key
      !     output:  ind(i) = index of i'th smallest object
      !********************************************************************
      IMPLICIT NONE
      !
      REAL,    INTENT (IN) :: lv(:)
      INTEGER, INTENT (OUT) :: ind(:)
      REAL,INTENT(IN),OPTIONAL :: lv1(:)
      !     ..
      !     .. Local Scalars ..
      REAL eps,q,q1
      INTEGER i,idx,ir,j,l,n
      REAL,ALLOCATABLE :: llv(:)
      !     ..
      !     .. Data statements ..
      DATA eps/1.e-10/
      !     ..
      !
      n=SIZE(ind)
      IF (n>SIZE(lv)) CALL judft_error("BUG: incosistent dimensions")
      ALLOCATE(llv(n))
      IF (PRESENT(lv1)) THEN
         IF (n>SIZE(lv1)) CALL judft_error("BUG: incosistent dimensions")
         llv=lv1
      ELSE
         llv=(/(1.*i,i=1,n)/)
      END IF
      IF (n == 0) RETURN ! Nothing to do
      IF (n == 1) THEN   ! Not much to do
         ind(1) = 1
         RETURN
      END IF
46

Matthias Redies's avatar
Matthias Redies committed
47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77
      DO i = 1,n
         ind(i) = i
      ENDDO
      !
      l = n/2 + 1
      ir = n
      DO
         IF (l.GT.1) THEN
            l = l - 1
            idx = ind(l)
            q = lv(idx)
            q1= llv(idx)
         ELSE
            idx = ind(ir)
            q = lv(idx)
            q1= llv(idx)
            ind(ir) = ind(1)
            ir = ir - 1
            IF (ir.EQ.1) THEN
               ind(1) = idx
               RETURN
            END IF
         END IF
         i = l
         j = l + l
         DO WHILE(j.LE.ir)
            IF (j.LT.ir) THEN
               !           if(lv(ind(j)).lt.lv(ind(j+1))) j=j+1
               IF (((lv(ind(j+1))-lv(ind(j))).GE.eps).OR. &!Standard comparison
                   ((ABS((lv(ind(j+1))-lv(ind(j))))<eps).AND.&!Same length, check second key
                    ((llv(ind(j+1))-llv(ind(j))).GE.eps))) &
78
                  j=j+1
Matthias Redies's avatar
Matthias Redies committed
79 80 81 82 83 84 85 86 87 88 89 90 91 92
            END IF
            !        if(q.lt.lv(ind(j))) then
            IF ((lv(ind(j))-q).GE.eps.OR.&
                (ABS((lv(ind(j))-q))<eps.AND.(llv(ind(j))-q1).GE.eps))THEN
               ind(i) = ind(j)
               i = j
               j = j + j
            ELSE
               j = ir + 1
            END IF
         enddo
         ind(i) = idx
      ENDDO
   END SUBROUTINE sort
93 94

END MODULE m_sort