Module m_uniinv Integer, Parameter :: kdp = selected_real_kind(15) public :: uniinv private :: kdp private :: R_uniinv, I_uniinv, D_uniinv private :: R_nearless, I_nearless, D_nearless, nearless interface uniinv module procedure d_uniinv, r_uniinv, i_uniinv end interface uniinv interface nearless module procedure D_nearless, R_nearless, I_nearless end interface nearless contains Subroutine D_uniinv (XDONT, IGOEST) ! __________________________________________________________ ! UNIINV = Merge-sort inverse ranking of an array, with removal of ! duplicate entries. ! The routine is similar to pure merge-sort ranking, but on ! the last pass, it sets indices in IGOEST to the rank ! of the value in the ordered set with duplicates removed. ! 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) :: IGOEST ! __________________________________________________________ Real (kind=kdp) :: XTST, XDONA, XDONB ! ! __________________________________________________________ Integer, Dimension (SIZE(IGOEST)) :: JWRKT, IRNGT Integer :: LMTNA, LMTNC, IRNG, IRNG1, IRNG2, NUNI Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB ! NVAL = Min (SIZE(XDONT), SIZE(IGOEST)) ! Select Case (NVAL) Case (:0) Return Case (1) IGOEST (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 <= 4) 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 (2*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 ! ! One steps in the C subset, that we create in the final rank array ! ! Make a copy of the rank array for the iteration ! JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA) XDONA = XDONT (JWRKT(IINDA)) XDONB = XDONT (IRNGT(IINDB)) ! Do IWRK = IWRK + 1 ! ! We still have unprocessed values in both A and B ! If (XDONA > XDONB) 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 XDONB = XDONT (IRNGT(IINDB)) Else IRNGT (IWRK) = JWRKT (IINDA) IINDA = IINDA + 1 If (IINDA > LMTNA) Exit! Only B still with unprocessed values XDONA = XDONT (JWRKT(IINDA)) End If ! End Do End Do ! ! The Cs become As and Bs ! LMTNA = 2 * LMTNA End Do ! ! Last merge of A and B into C, with removal of duplicates. ! IINDA = 1 IINDB = LMTNA + 1 NUNI = 0 ! ! One steps in the C subset, that we create in the final rank array ! JWRKT (1:LMTNA) = IRNGT (1:LMTNA) If (IINDB <= NVAL) Then XTST = NEARLESS (Min(XDONT(JWRKT(1)), XDONT(IRNGT(IINDB)))) Else XTST = NEARLESS (XDONT(JWRKT(1))) Endif Do IWRK = 1, NVAL ! ! We still have unprocessed values in both A and B ! If (IINDA <= LMTNA) Then If (IINDB <= NVAL) Then If (XDONT(JWRKT(IINDA)) > XDONT(IRNGT(IINDB))) Then IRNG = IRNGT (IINDB) IINDB = IINDB + 1 Else IRNG = JWRKT (IINDA) IINDA = IINDA + 1 End If Else ! ! Only A still with unprocessed values ! IRNG = JWRKT (IINDA) IINDA = IINDA + 1 End If Else ! ! Only B still with unprocessed values ! IRNG = IRNGT (IWRK) End If If (XDONT(IRNG) > XTST) Then XTST = XDONT (IRNG) NUNI = NUNI + 1 End If IGOEST (IRNG) = NUNI ! End Do ! Return ! End Subroutine D_uniinv Subroutine R_uniinv (XDONT, IGOEST) ! __________________________________________________________ ! UNIINV = Merge-sort inverse ranking of an array, with removal of ! duplicate entries. ! The routine is similar to pure merge-sort ranking, but on ! the last pass, it sets indices in IGOEST to the rank ! of the value in the ordered set with duplicates removed. ! 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) :: IGOEST ! __________________________________________________________ Real :: XTST, XDONA, XDONB ! ! __________________________________________________________ Integer, Dimension (SIZE(IGOEST)) :: JWRKT, IRNGT Integer :: LMTNA, LMTNC, IRNG, IRNG1, IRNG2, NUNI Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB ! NVAL = Min (SIZE(XDONT), SIZE(IGOEST)) ! Select Case (NVAL) Case (:0) Return Case (1) IGOEST (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 <= 4) 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 (2*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 ! ! One steps in the C subset, that we create in the final rank array ! ! Make a copy of the rank array for the iteration ! JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA) XDONA = XDONT (JWRKT(IINDA)) XDONB = XDONT (IRNGT(IINDB)) ! Do IWRK = IWRK + 1 ! ! We still have unprocessed values in both A and B ! If (XDONA > XDONB) 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 XDONB = XDONT (IRNGT(IINDB)) Else IRNGT (IWRK) = JWRKT (IINDA) IINDA = IINDA + 1 If (IINDA > LMTNA) Exit! Only B still with unprocessed values XDONA = XDONT (JWRKT(IINDA)) End If ! End Do End Do ! ! The Cs become As and Bs ! LMTNA = 2 * LMTNA End Do ! ! Last merge of A and B into C, with removal of duplicates. ! IINDA = 1 IINDB = LMTNA + 1 NUNI = 0 ! ! One steps in the C subset, that we create in the final rank array ! JWRKT (1:LMTNA) = IRNGT (1:LMTNA) If (IINDB <= NVAL) Then XTST = NEARLESS (Min(XDONT(JWRKT(1)), XDONT(IRNGT(IINDB)))) Else XTST = NEARLESS (XDONT(JWRKT(1))) Endif Do IWRK = 1, NVAL ! ! We still have unprocessed values in both A and B ! If (IINDA <= LMTNA) Then If (IINDB <= NVAL) Then If (XDONT(JWRKT(IINDA)) > XDONT(IRNGT(IINDB))) Then IRNG = IRNGT (IINDB) IINDB = IINDB + 1 Else IRNG = JWRKT (IINDA) IINDA = IINDA + 1 End If Else ! ! Only A still with unprocessed values ! IRNG = JWRKT (IINDA) IINDA = IINDA + 1 End If Else ! ! Only B still with unprocessed values ! IRNG = IRNGT (IWRK) End If If (XDONT(IRNG) > XTST) Then XTST = XDONT (IRNG) NUNI = NUNI + 1 End If IGOEST (IRNG) = NUNI ! End Do ! Return ! End Subroutine R_uniinv Subroutine I_uniinv (XDONT, IGOEST) ! __________________________________________________________ ! UNIINV = Merge-sort inverse ranking of an array, with removal of ! duplicate entries. ! The routine is similar to pure merge-sort ranking, but on ! the last pass, it sets indices in IGOEST to the rank ! of the value in the ordered set with duplicates removed. ! 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) :: IGOEST ! __________________________________________________________ Integer :: XTST, XDONA, XDONB ! ! __________________________________________________________ Integer, Dimension (SIZE(IGOEST)) :: JWRKT, IRNGT Integer :: LMTNA, LMTNC, IRNG, IRNG1, IRNG2, NUNI Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB ! NVAL = Min (SIZE(XDONT), SIZE(IGOEST)) ! Select Case (NVAL) Case (:0) Return Case (1) IGOEST (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 <= 4) 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 (2*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 ! ! One steps in the C subset, that we create in the final rank array ! ! Make a copy of the rank array for the iteration ! JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA) XDONA = XDONT (JWRKT(IINDA)) XDONB = XDONT (IRNGT(IINDB)) ! Do IWRK = IWRK + 1 ! ! We still have unprocessed values in both A and B ! If (XDONA > XDONB) 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 XDONB = XDONT (IRNGT(IINDB)) Else IRNGT (IWRK) = JWRKT (IINDA) IINDA = IINDA + 1 If (IINDA > LMTNA) Exit! Only B still with unprocessed values XDONA = XDONT (JWRKT(IINDA)) End If ! End Do End Do ! ! The Cs become As and Bs ! LMTNA = 2 * LMTNA End Do ! ! Last merge of A and B into C, with removal of duplicates. ! IINDA = 1 IINDB = LMTNA + 1 NUNI = 0 ! ! One steps in the C subset, that we create in the final rank array ! JWRKT (1:LMTNA) = IRNGT (1:LMTNA) If (IINDB <= NVAL) Then XTST = NEARLESS (Min(XDONT(JWRKT(1)), XDONT(IRNGT(IINDB)))) Else XTST = NEARLESS (XDONT(JWRKT(1))) Endif Do IWRK = 1, NVAL ! ! We still have unprocessed values in both A and B ! If (IINDA <= LMTNA) Then If (IINDB <= NVAL) Then If (XDONT(JWRKT(IINDA)) > XDONT(IRNGT(IINDB))) Then IRNG = IRNGT (IINDB) IINDB = IINDB + 1 Else IRNG = JWRKT (IINDA) IINDA = IINDA + 1 End If Else ! ! Only A still with unprocessed values ! IRNG = JWRKT (IINDA) IINDA = IINDA + 1 End If Else ! ! Only B still with unprocessed values ! IRNG = IRNGT (IWRK) End If If (XDONT(IRNG) > XTST) Then XTST = XDONT (IRNG) NUNI = NUNI + 1 End If IGOEST (IRNG) = NUNI ! End Do ! Return ! End Subroutine I_uniinv Function D_nearless (XVAL) result (D_nl) ! Nearest value less than given value ! __________________________________________________________ Real (kind=kdp), Intent (In) :: XVAL Real (kind=kdp) :: D_nl ! __________________________________________________________ D_nl = nearest (XVAL, -1.0_kdp) return ! End Function D_nearless Function R_nearless (XVAL) result (R_nl) ! Nearest value less than given value ! __________________________________________________________ Real, Intent (In) :: XVAL Real :: R_nl ! __________________________________________________________ R_nl = nearest (XVAL, -1.0) return ! End Function R_nearless Function I_nearless (XVAL) result (I_nl) ! Nearest value less than given value ! __________________________________________________________ Integer, Intent (In) :: XVAL Integer :: I_nl ! __________________________________________________________ I_nl = XVAL - 1 return ! End Function I_nearless end module m_uniinv