Line data    Source code

       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : MODULE m_sort
       7             :    USE m_judft
       8             : CONTAINS
       9             : 
      10       17468 :    SUBROUTINE sort(ind,lv,lv1)
      11             :       !********************************************************************
      12             :       !     heapsort routine
      13             :       !     input:   lv    = array of objects to be sorted
      14             :       !              lv1   = second array to use as secondary sort key
      15             :       !     output:  ind(i) = index of i'th smallest object
      16             :       !********************************************************************
      17             :       IMPLICIT NONE
      18             :       !
      19             :       REAL,    INTENT (IN) :: lv(:)
      20             :       INTEGER, INTENT (OUT) :: ind(:)
      21             :       REAL,INTENT(IN),OPTIONAL :: lv1(:)
      22             :       !     ..
      23             :       !     .. Local Scalars ..
      24             :       REAL eps,q,q1
      25             :       INTEGER i,idx,ir,j,l,n
      26       17468 :       REAL,ALLOCATABLE :: llv(:)
      27             :       !     ..
      28             :       !     .. Data statements ..
      29             :       DATA eps/1.e-10/
      30             :       !     ..
      31             :       !
      32       17468 :       n=SIZE(ind)
      33       17468 :       IF (n>SIZE(lv)) CALL judft_error("BUG in sort: inconsistent dimensions")
      34       52404 :       ALLOCATE(llv(n))
      35       17468 :       IF (PRESENT(lv1)) THEN
      36       16959 :          IF (n>SIZE(lv1)) CALL judft_error("BUG in sort: inconsistent dimensions")
      37     6373574 :          llv=lv1
      38             :       ELSE
      39      553655 :          llv=(/(1.*i,i=1,n)/)
      40             :       END IF
      41       34936 :       IF (n == 0) RETURN ! Nothing to do
      42       17446 :       IF (n == 1) THEN   ! Not much to do
      43           0 :          ind(1) = 1
      44           0 :          RETURN
      45             :       END IF
      46             : 
      47     6541152 :       DO i = 1,n
      48     6541152 :          ind(i) = i
      49             :       ENDDO
      50             :       !
      51       17446 :       l = n/2 + 1
      52       17446 :       ir = n
      53     9747681 :       DO
      54     9765127 :          IF (l.GT.1) THEN
      55     3258867 :             l = l - 1
      56     3258867 :             idx = ind(l)
      57     3258867 :             q = lv(idx)
      58     3258867 :             q1= llv(idx)
      59             :          ELSE
      60     6506260 :             idx = ind(ir)
      61     6506260 :             q = lv(idx)
      62     6506260 :             q1= llv(idx)
      63     6506260 :             ind(ir) = ind(1)
      64     6506260 :             ir = ir - 1
      65     6506260 :             IF (ir.EQ.1) THEN
      66       17446 :                ind(1) = idx
      67       17446 :                RETURN
      68             :             END IF
      69             :          END IF
      70     9747681 :          i = l
      71     9747681 :          j = l + l
      72    89460295 :          DO WHILE(j.LE.ir)
      73    79712614 :             IF (j.LT.ir) THEN
      74             :                !           if(lv(ind(j)).lt.lv(ind(j+1))) j=j+1
      75    79582347 :                IF (((lv(ind(j+1))-lv(ind(j))).GE.eps).OR. &!Standard comparison
      76             :                    ((ABS((lv(ind(j+1))-lv(ind(j))))<eps).AND.&!Same length, check second key
      77             :                     ((llv(ind(j+1))-llv(ind(j))).GE.eps))) &
      78    79712614 :                   j=j+1
      79             :             END IF
      80             :             !        if(q.lt.lv(ind(j))) then
      81    79712614 :             IF ((lv(ind(j))-q).GE.eps.OR.&
      82     9747681 :                 (ABS((lv(ind(j))-q))<eps.AND.(llv(ind(j))-q1).GE.eps))THEN
      83    77800121 :                ind(i) = ind(j)
      84    77800121 :                i = j
      85    77800121 :                j = j + j
      86             :             ELSE
      87     1912493 :                j = ir + 1
      88             :             END IF
      89             :          enddo
      90     9747681 :          ind(i) = idx
      91             :       ENDDO
      92       34936 :    END SUBROUTINE sort
      93             : 
      94             : END MODULE m_sort