m_mrgrnk.f90

! Author: Michel Olagnon
! Code downloaded from: http://fortran-2000.com/rank/
! License: not specified, but public it seems
module m_mrgrnk
Integer, Parameter :: kdp = selected_real_kind(15)
public :: mrgrnk
private :: kdp
private :: r_mrgrnk, i_mrgrnk, d_mrgrnk
interface mrgrnk
  module procedure d_mrgrnk, r_mrgrnk, i_mrgrnk
end interface mrgrnk
contains
 
Subroutine d_mrgrnk (XDONT, IRNGT)
! __________________________________________________________
!   MRGRNK = Merge-sort ranking of an array
!   For performance reasons, the first 2 passes are taken
!   out of the standard loop, and use dedicated coding.
! __________________________________________________________
! __________________________________________________________
      Real (kind=kdp), Dimension (:), Intent (In) :: xdont
      Integer, Dimension (:), Intent (Out) :: IRNGT
! __________________________________________________________
      Real (kind=kdp) :: xvala, xvalb
!
      Integer, Dimension (SIZE(IRNGT)) :: JWRKT
      Integer :: LMTNA, LMTNC, IRNG1, IRNG2
      Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
!
      nval = min(SIZE(xdont), SIZE(irngt))
      Select Case (nval)
      Case (:0)
         Return
      Case (1)
         irngt(1) = 1
         Return
      Case Default
         Continue
      End Select
!
!  Fill-in the index array, creating ordered couples
!
      Do iind = 2, nval, 2
         If (xdont(iind-1) <= xdont(iind)) Then
            irngt(iind-1) = iind - 1
            irngt(iind) = iind
         Else
            irngt(iind-1) = iind
            irngt(iind) = iind - 1
         End If
      End Do
      If (modulo(nval, 2) /= 0) Then
         irngt(nval) = nval
      End If
!
!  We will now have ordered subsets A - B - A - B - ...
!  and merge A and B couples into     C   -   C   - ...
!
      lmtna = 2
      lmtnc = 4
!
!  First iteration. The length of the ordered subsets goes from 2 to 4
!
      Do
         If (nval <= 2) Exit
!
!   Loop on merges of A and B into C
!
         Do iwrkd = 0, nval - 1, 4
            If ((iwrkd+4) > nval) Then
               If ((iwrkd+2) >= nval) Exit
!
!   1 2 3
!
               If (xdont(irngt(iwrkd+2)) <= xdont(irngt(iwrkd+3))) Exit
!
!   1 3 2
!
               If (xdont(irngt(iwrkd+1)) <= xdont(irngt(iwrkd+3))) Then
                  irng2 = irngt(iwrkd+2)
                  irngt(iwrkd+2) = irngt(iwrkd+3)
                  irngt(iwrkd+3) = irng2
!
!   3 1 2
!
               Else
                  irng1 = irngt(iwrkd+1)
                  irngt(iwrkd+1) = irngt(iwrkd+3)
                  irngt(iwrkd+3) = irngt(iwrkd+2)
                  irngt(iwrkd+2) = irng1
               End If
               Exit
            End If
!
!   1 2 3 4
!
            If (xdont(irngt(iwrkd+2)) <= xdont(irngt(iwrkd+3))) cycle
!
!   1 3 x x
!
            If (xdont(irngt(iwrkd+1)) <= xdont(irngt(iwrkd+3))) Then
               irng2 = irngt(iwrkd+2)
               irngt(iwrkd+2) = irngt(iwrkd+3)
               If (xdont(irng2) <= xdont(irngt(iwrkd+4))) Then
!   1 3 2 4
                  irngt(iwrkd+3) = irng2
               Else
!   1 3 4 2
                  irngt(iwrkd+3) = irngt(iwrkd+4)
                  irngt(iwrkd+4) = irng2
               End If
!
!   3 x x x
!
            Else
               irng1 = irngt(iwrkd+1)
               irng2 = irngt(iwrkd+2)
               irngt(iwrkd+1) = irngt(iwrkd+3)
               If (xdont(irng1) <= xdont(irngt(iwrkd+4))) Then
                  irngt(iwrkd+2) = irng1
                  If (xdont(irng2) <= xdont(irngt(iwrkd+4))) Then
!   3 1 2 4
                     irngt(iwrkd+3) = irng2
                  Else
!   3 1 4 2
                     irngt(iwrkd+3) = irngt(iwrkd+4)
                     irngt(iwrkd+4) = irng2
                  End If
               Else
!   3 4 1 2
                  irngt(iwrkd+2) = irngt(iwrkd+4)
                  irngt(iwrkd+3) = irng1
                  irngt(iwrkd+4) = irng2
               End If
            End If
         End Do
!
!  The Cs become As and Bs
!
         lmtna = 4
         Exit
      End Do
!
!  Iteration loop. Each time, the length of the ordered subsets
!  is doubled.
!
      Do
         If (lmtna >= nval) Exit
         iwrkf = 0
         lmtnc = 2 * lmtnc
!
!   Loop on merges of A and B into C
!
         Do
            iwrk = iwrkf
            iwrkd = iwrkf + 1
            jinda = iwrkf + lmtna
            iwrkf = iwrkf + lmtnc
            If (iwrkf >= nval) Then
               If (jinda >= nval) Exit
               iwrkf = nval
            End If
            iinda = 1
            iindb = jinda + 1
!
!   Shortcut for the case when the max of A is smaller
!   than the min of B. This line may be activated when the
!   initial set is already close to sorted.
!
!          IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE
!
!  One steps in the C subset, that we build in the final rank array
!
!  Make a copy of the rank array for the merge iteration
!
            jwrkt(1:lmtna) = irngt(iwrkd:jinda)
!
            xvala = xdont(jwrkt(iinda))
            xvalb = xdont(irngt(iindb))
!
            Do
               iwrk = iwrk + 1
!
!  We still have unprocessed values in both A and B
!
               If (xvala > xvalb) Then
                  irngt(iwrk) = irngt(iindb)
                  iindb = iindb + 1
                  If (iindb > iwrkf) Then
!  Only A still with unprocessed values
                     irngt(iwrk+1:iwrkf) = jwrkt(iinda:lmtna)
                     Exit
                  End If
                  xvalb = xdont(irngt(iindb))
               Else
                  irngt(iwrk) = jwrkt(iinda)
                  iinda = iinda + 1
                  If (iinda > lmtna) exit! Only B still with unprocessed values
                  xvala = xdont(jwrkt(iinda))
               End If
!
            End Do
         End Do
!
!  The Cs become As and Bs
!
         lmtna = 2 * lmtna
      End Do
!
      Return
!
End Subroutine d_mrgrnk
 
Subroutine r_mrgrnk (XDONT, IRNGT)
! __________________________________________________________
!   MRGRNK = Merge-sort ranking of an array
!   For performance reasons, the first 2 passes are taken
!   out of the standard loop, and use dedicated coding.
! __________________________________________________________
! _________________________________________________________
      Real, Dimension (:), Intent (In) :: XDONT
      Integer, Dimension (:), Intent (Out) :: IRNGT
! __________________________________________________________
      Real :: XVALA, XVALB
!
      Integer, Dimension (SIZE(IRNGT)) :: JWRKT
      Integer :: LMTNA, LMTNC, IRNG1, IRNG2
      Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
!
      nval = min(SIZE(xdont), SIZE(irngt))
      Select Case (nval)
      Case (:0)
         Return
      Case (1)
         irngt(1) = 1
         Return
      Case Default
         Continue
      End Select
!
!  Fill-in the index array, creating ordered couples
!
      Do iind = 2, nval, 2
         If (xdont(iind-1) <= xdont(iind)) Then
            irngt(iind-1) = iind - 1
            irngt(iind) = iind
         Else
            irngt(iind-1) = iind
            irngt(iind) = iind - 1
         End If
      End Do
      If (modulo(nval, 2) /= 0) Then
         irngt(nval) = nval
      End If
!
!  We will now have ordered subsets A - B - A - B - ...
!  and merge A and B couples into     C   -   C   - ...
!
      lmtna = 2
      lmtnc = 4
!
!  First iteration. The length of the ordered subsets goes from 2 to 4
!
      Do
         If (nval <= 2) Exit
!
!   Loop on merges of A and B into C
!
         Do iwrkd = 0, nval - 1, 4
            If ((iwrkd+4) > nval) Then
               If ((iwrkd+2) >= nval) Exit
!
!   1 2 3
!
               If (xdont(irngt(iwrkd+2)) <= xdont(irngt(iwrkd+3))) Exit
!
!   1 3 2
!
               If (xdont(irngt(iwrkd+1)) <= xdont(irngt(iwrkd+3))) Then
                  irng2 = irngt(iwrkd+2)
                  irngt(iwrkd+2) = irngt(iwrkd+3)
                  irngt(iwrkd+3) = irng2
!
!   3 1 2
!
               Else
                  irng1 = irngt(iwrkd+1)
                  irngt(iwrkd+1) = irngt(iwrkd+3)
                  irngt(iwrkd+3) = irngt(iwrkd+2)
                  irngt(iwrkd+2) = irng1
               End If
               Exit
            End If
!
!   1 2 3 4
!
            If (xdont(irngt(iwrkd+2)) <= xdont(irngt(iwrkd+3))) cycle
!
!   1 3 x x
!
            If (xdont(irngt(iwrkd+1)) <= xdont(irngt(iwrkd+3))) Then
               irng2 = irngt(iwrkd+2)
               irngt(iwrkd+2) = irngt(iwrkd+3)
               If (xdont(irng2) <= xdont(irngt(iwrkd+4))) Then
!   1 3 2 4
                  irngt(iwrkd+3) = irng2
               Else
!   1 3 4 2
                  irngt(iwrkd+3) = irngt(iwrkd+4)
                  irngt(iwrkd+4) = irng2
               End If
!
!   3 x x x
!
            Else
               irng1 = irngt(iwrkd+1)
               irng2 = irngt(iwrkd+2)
               irngt(iwrkd+1) = irngt(iwrkd+3)
               If (xdont(irng1) <= xdont(irngt(iwrkd+4))) Then
                  irngt(iwrkd+2) = irng1
                  If (xdont(irng2) <= xdont(irngt(iwrkd+4))) Then
!   3 1 2 4
                     irngt(iwrkd+3) = irng2
                  Else
!   3 1 4 2
                     irngt(iwrkd+3) = irngt(iwrkd+4)
                     irngt(iwrkd+4) = irng2
                  End If
               Else
!   3 4 1 2
                  irngt(iwrkd+2) = irngt(iwrkd+4)
                  irngt(iwrkd+3) = irng1
                  irngt(iwrkd+4) = irng2
               End If
            End If
         End Do
!
!  The Cs become As and Bs
!
         lmtna = 4
         Exit
      End Do
!
!  Iteration loop. Each time, the length of the ordered subsets
!  is doubled.
!
      Do
         If (lmtna >= nval) Exit
         iwrkf = 0
         lmtnc = 2 * lmtnc
!
!   Loop on merges of A and B into C
!
         Do
            iwrk = iwrkf
            iwrkd = iwrkf + 1
            jinda = iwrkf + lmtna
            iwrkf = iwrkf + lmtnc
            If (iwrkf >= nval) Then
               If (jinda >= nval) Exit
               iwrkf = nval
            End If
            iinda = 1
            iindb = jinda + 1
!
!   Shortcut for the case when the max of A is smaller
!   than the min of B. This line may be activated when the
!   initial set is already close to sorted.
!
!          IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE
!
!  One steps in the C subset, that we build in the final rank array
!
!  Make a copy of the rank array for the merge iteration
!
            jwrkt(1:lmtna) = irngt(iwrkd:jinda)
!
            xvala = xdont(jwrkt(iinda))
            xvalb = xdont(irngt(iindb))
!
            Do
               iwrk = iwrk + 1
!
!  We still have unprocessed values in both A and B
!
               If (xvala > xvalb) Then
                  irngt(iwrk) = irngt(iindb)
                  iindb = iindb + 1
                  If (iindb > iwrkf) Then
!  Only A still with unprocessed values
                     irngt(iwrk+1:iwrkf) = jwrkt(iinda:lmtna)
                     Exit
                  End If
                  xvalb = xdont(irngt(iindb))
               Else
                  irngt(iwrk) = jwrkt(iinda)
                  iinda = iinda + 1
                  If (iinda > lmtna) exit! Only B still with unprocessed values
                  xvala = xdont(jwrkt(iinda))
               End If
!
            End Do
         End Do
!
!  The Cs become As and Bs
!
         lmtna = 2 * lmtna
      End Do
!
      Return
!
End Subroutine r_mrgrnk
Subroutine i_mrgrnk (XDONT, IRNGT)
! __________________________________________________________
!   MRGRNK = Merge-sort ranking of an array
!   For performance reasons, the first 2 passes are taken
!   out of the standard loop, and use dedicated coding.
! __________________________________________________________
! __________________________________________________________
      Integer, Dimension (:), Intent (In)  :: XDONT
      Integer, Dimension (:), Intent (Out) :: IRNGT
! __________________________________________________________
      Integer :: XVALA, XVALB
!
      Integer, Dimension (SIZE(IRNGT)) :: JWRKT
      Integer :: LMTNA, LMTNC, IRNG1, IRNG2
      Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
!
      nval = min(SIZE(xdont), SIZE(irngt))
      Select Case (nval)
      Case (:0)
         Return
      Case (1)
         irngt(1) = 1
         Return
      Case Default
         Continue
      End Select
!
!  Fill-in the index array, creating ordered couples
!
      Do iind = 2, nval, 2
         If (xdont(iind-1) <= xdont(iind)) Then
            irngt(iind-1) = iind - 1
            irngt(iind) = iind
         Else
            irngt(iind-1) = iind
            irngt(iind) = iind - 1
         End If
      End Do
      If (modulo(nval, 2) /= 0) Then
         irngt(nval) = nval
      End If
!
!  We will now have ordered subsets A - B - A - B - ...
!  and merge A and B couples into     C   -   C   - ...
!
      lmtna = 2
      lmtnc = 4
!
!  First iteration. The length of the ordered subsets goes from 2 to 4
!
      Do
         If (nval <= 2) Exit
!
!   Loop on merges of A and B into C
!
         Do iwrkd = 0, nval - 1, 4
            If ((iwrkd+4) > nval) Then
               If ((iwrkd+2) >= nval) Exit
!
!   1 2 3
!
               If (xdont(irngt(iwrkd+2)) <= xdont(irngt(iwrkd+3))) Exit
!
!   1 3 2
!
               If (xdont(irngt(iwrkd+1)) <= xdont(irngt(iwrkd+3))) Then
                  irng2 = irngt(iwrkd+2)
                  irngt(iwrkd+2) = irngt(iwrkd+3)
                  irngt(iwrkd+3) = irng2
!
!   3 1 2
!
               Else
                  irng1 = irngt(iwrkd+1)
                  irngt(iwrkd+1) = irngt(iwrkd+3)
                  irngt(iwrkd+3) = irngt(iwrkd+2)
                  irngt(iwrkd+2) = irng1
               End If
               Exit
            End If
!
!   1 2 3 4
!
            If (xdont(irngt(iwrkd+2)) <= xdont(irngt(iwrkd+3))) cycle
!
!   1 3 x x
!
            If (xdont(irngt(iwrkd+1)) <= xdont(irngt(iwrkd+3))) Then
               irng2 = irngt(iwrkd+2)
               irngt(iwrkd+2) = irngt(iwrkd+3)
               If (xdont(irng2) <= xdont(irngt(iwrkd+4))) Then
!   1 3 2 4
                  irngt(iwrkd+3) = irng2
               Else
!   1 3 4 2
                  irngt(iwrkd+3) = irngt(iwrkd+4)
                  irngt(iwrkd+4) = irng2
               End If
!
!   3 x x x
!
            Else
               irng1 = irngt(iwrkd+1)
               irng2 = irngt(iwrkd+2)
               irngt(iwrkd+1) = irngt(iwrkd+3)
               If (xdont(irng1) <= xdont(irngt(iwrkd+4))) Then
                  irngt(iwrkd+2) = irng1
                  If (xdont(irng2) <= xdont(irngt(iwrkd+4))) Then
!   3 1 2 4
                     irngt(iwrkd+3) = irng2
                  Else
!   3 1 4 2
                     irngt(iwrkd+3) = irngt(iwrkd+4)
                     irngt(iwrkd+4) = irng2
                  End If
               Else
!   3 4 1 2
                  irngt(iwrkd+2) = irngt(iwrkd+4)
                  irngt(iwrkd+3) = irng1
                  irngt(iwrkd+4) = irng2
               End If
            End If
         End Do
!
!  The Cs become As and Bs
!
         lmtna = 4
         Exit
      End Do
!
!  Iteration loop. Each time, the length of the ordered subsets
!  is doubled.
!
      Do
         If (lmtna >= nval) Exit
         iwrkf = 0
         lmtnc = 2 * lmtnc
!
!   Loop on merges of A and B into C
!
         Do
            iwrk = iwrkf
            iwrkd = iwrkf + 1
            jinda = iwrkf + lmtna
            iwrkf = iwrkf + lmtnc
            If (iwrkf >= nval) Then
               If (jinda >= nval) Exit
               iwrkf = nval
            End If
            iinda = 1
            iindb = jinda + 1
!
!   Shortcut for the case when the max of A is smaller
!   than the min of B. This line may be activated when the
!   initial set is already close to sorted.
!
!          IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE
!
!  One steps in the C subset, that we build in the final rank array
!
!  Make a copy of the rank array for the merge iteration
!
            jwrkt(1:lmtna) = irngt(iwrkd:jinda)
!
            xvala = xdont(jwrkt(iinda))
            xvalb = xdont(irngt(iindb))
!
            Do
               iwrk = iwrk + 1
!
!  We still have unprocessed values in both A and B
!
               If (xvala > xvalb) Then
                  irngt(iwrk) = irngt(iindb)
                  iindb = iindb + 1
                  If (iindb > iwrkf) Then
!  Only A still with unprocessed values
                     irngt(iwrk+1:iwrkf) = jwrkt(iinda:lmtna)
                     Exit
                  End If
                  xvalb = xdont(irngt(iindb))
               Else
                  irngt(iwrk) = jwrkt(iinda)
                  iinda = iinda + 1
                  If (iinda > lmtna) exit! Only B still with unprocessed values
                  xvala = xdont(jwrkt(iinda))
               End If
!
            End Do
         End Do
!
!  The Cs become As and Bs
!
         lmtna = 2 * lmtna
      End Do
!
      Return
!
End Subroutine i_mrgrnk
end module m_mrgrnk