Module m_rapknr Integer, Parameter :: kdp = selected_real_kind(15) public :: rapknr private :: kdp private :: R_rapknr, I_rapknr, D_rapknr interface rapknr module procedure d_rapknr, r_rapknr, i_rapknr end interface rapknr contains Subroutine D_rapknr (XDONT, IRNGT, NORD) ! Ranks partially XDONT by IRNGT, up to order NORD, in decreasing order. ! rapknr = (rnkpar backwards) ! __________________________________________________________ ! This routine uses a pivoting strategy such as the one of ! finding the median based on the quicksort algorithm, but ! we skew the pivot choice to try to bring it to NORD as ! fast as possible. It uses 2 temporary arrays, where it ! stores the indices of the values larger than the pivot ! (IHIGT), and the indices of values smaller than the pivot ! that we might still need later on (ILOWT). It iterates ! until it can bring the number of values in IHIGT to ! exactly NORD, and then uses an insertion sort to rank ! this set, since it is supposedly small. ! Michel Olagnon - Feb. 2011 ! __________________________________________________________ ! __________________________________________________________ Real (kind=kdp), Dimension (:), Intent (In) :: XDONT Integer, Dimension (:), Intent (Out) :: IRNGT Integer, Intent (In) :: NORD ! __________________________________________________________ Real (kind=kdp) :: XPIV, XPIV0, XWRK, XWRK1, XMIN, XMAX ! Integer, Dimension (SIZE(XDONT)) :: ILOWT, IHIGT Integer :: NDON, JHIG, JLOW, IHIG, IWRK, IWRK1, IWRK2, IWRK3 Integer :: IDEB, JDEB, IMIL, IFIN, NWRK, ICRS, IDCR, ILOW Integer :: JLM2, JLM1, JHM2, JHM1 ! NDON = SIZE (XDONT) ! ! First loop is used to fill-in ILOWT, IHIGT at the same time ! If (NDON < 2) Then If (NORD >= 1) IRNGT (1) = 1 Return End If ! ! One chooses a pivot, best estimate possible to put fractile near ! mid-point of the set of high values. ! If (XDONT(2) < XDONT(1)) Then ILOWT (1) = 2 IHIGT (1) = 1 Else ILOWT (1) = 1 IHIGT (1) = 2 End If ! If (NDON < 3) Then If (NORD >= 1) IRNGT (1) = IHIGT (1) If (NORD >= 2) IRNGT (2) = ILOWT (1) Return End If ! --- If (XDONT(3) > XDONT(ILOWT(1))) Then ILOWT (2) = ILOWT (1) If (XDONT(3) > XDONT(IHIGT(1))) Then ILOWT (1) = IHIGT (1) IHIGT (1) = 3 Else ILOWT (1) = 3 End If Else ILOWT (2) = 3 End If ! --- If (NDON < 4) Then If (NORD >= 1) IRNGT (1) = IHIGT (1) If (NORD >= 2) IRNGT (2) = ILOWT (1) If (NORD >= 3) IRNGT (3) = ILOWT (2) Return End If ! If (XDONT(NDON) > XDONT(ILOWT(1))) Then ILOWT (3) = ILOWT (2) ILOWT (2) = ILOWT (1) If (XDONT(NDON) > XDONT(IHIGT(1))) Then ILOWT (1) = IHIGT (1) IHIGT (1) = NDON Else ILOWT (1) = NDON End If Else if (XDONT (NDON) > XDONT (ILOWT(2))) Then ILOWT (3) = ILOWT (2) ILOWT (2) = NDON else ILOWT (3) = NDON endif End If ! If (NDON < 5) Then If (NORD >= 1) IRNGT (1) = IHIGT (1) If (NORD >= 2) IRNGT (2) = ILOWT (1) If (NORD >= 3) IRNGT (3) = ILOWT (2) If (NORD >= 4) IRNGT (4) = ILOWT (3) Return End If ! --- JDEB = 0 IDEB = JDEB + 1 JHIG = IDEB JLOW = 3 XPIV = XDONT (IHIGT(IDEB)) + REAL(2*NORD)/REAL(NDON+NORD) * & (XDONT(ILOWT(3))-XDONT(IHIGT(IDEB))) If (XPIV >= XDONT(ILOWT(1))) Then XPIV = XDONT (IHIGT(IDEB)) + REAL(2*NORD)/REAL(NDON+NORD) * & (XDONT(ILOWT(2))-XDONT(IHIGT(IDEB))) If (XPIV >= XDONT(ILOWT(1))) & XPIV = XDONT (IHIGT(IDEB)) + REAL (2*NORD) / REAL (NDON+NORD) * & (XDONT(ILOWT(1))-XDONT(IHIGT(IDEB))) End If XPIV0 = XPIV ! --- ! One puts values < pivot in the end and those >= pivot ! at the beginning. This is split in 2 cases, so that ! we can skip the loop test a number of times. ! As we are also filling in the work arrays at the same time ! we stop filling in the ILOWT array as soon as we have more ! than enough values in IHIGT. ! ! If (XDONT(NDON) < XPIV) Then ICRS = 3 Do ICRS = ICRS + 1 If (XDONT(ICRS) < XPIV) Then If (ICRS >= NDON) Exit JLOW = JLOW + 1 ILOWT (JLOW) = ICRS Else JHIG = JHIG + 1 IHIGT (JHIG) = ICRS If (JHIG >= NORD) Exit End If End Do ! ! One restricts further processing because it is no use ! to store more low values ! If (ICRS < NDON-1) Then Do ICRS = ICRS + 1 If (XDONT(ICRS) >= XPIV) Then JHIG = JHIG + 1 IHIGT (JHIG) = ICRS Else If (ICRS >= NDON) Then Exit End If End Do End If ! ! --- Else ! ! Same as above, but this is not as easy to optimize, so the ! DO-loop is kept ! Do ICRS = 4, NDON - 1 If (XDONT(ICRS) < XPIV) Then JLOW = JLOW + 1 ILOWT (JLOW) = ICRS Else JHIG = JHIG + 1 IHIGT (JHIG) = ICRS If (JHIG >= NORD) Exit End If End Do ! If (ICRS < NDON-1) Then Do ICRS = ICRS + 1 If (XDONT(ICRS) >= XPIV) Then If (ICRS >= NDON) Exit JHIG = JHIG + 1 IHIGT (JHIG) = ICRS End If End Do End If End If ! --- JLM2 = 0 JLM1 = 0 JHM2 = 0 JHM1 = 0 Do if (JHIG == NORD) Exit If (JHM2 == JHIG .And. JLM2 == JLOW) Then ! ! We are oscillating. Perturbate by bringing JHIG closer by one ! to NORD ! If (NORD > JHIG) Then XMAX = XDONT (ILOWT(1)) ILOW = 1 Do ICRS = 2, JLOW If (XDONT(ILOWT(ICRS)) > XMAX) Then XMAX = XDONT (ILOWT(ICRS)) ILOW = ICRS End If End Do ! JHIG = JHIG + 1 IHIGT (JHIG) = ILOWT (ILOW) ILOWT (ILOW) = ILOWT (JLOW) JLOW = JLOW - 1 Else IHIG = IHIGT (JHIG) XMIN = XDONT (IHIG) Do ICRS = 1, JHIG If (XDONT(IHIGT(ICRS)) < XMIN) Then IWRK = IHIGT (ICRS) XMIN = XDONT (IWRK) IHIGT (ICRS) = IHIG IHIG = IWRK End If End Do JHIG = JHIG - 1 End If End If JLM2 = JLM1 JLM1 = JLOW JHM2 = JHM1 JHM1 = JHIG ! --- ! We try to bring the number of values in the high values set ! closer to NORD. ! Select Case (NORD-JHIG) Case (2:) ! ! Not enough values in low part, at least 2 are missing ! Select Case (JLOW) !!!!! CASE DEFAULT !!!!! write (*,*) "Assertion failed" !!!!! STOP ! ! We make a special case when we have so few values in ! the low values set that it is bad performance to choose a pivot ! and apply the general algorithm. ! Case (2) If (XDONT(ILOWT(1)) >= XDONT(ILOWT(2))) Then JHIG = JHIG + 1 IHIGT (JHIG) = ILOWT (1) JHIG = JHIG + 1 IHIGT (JHIG) = ILOWT (2) Else JHIG = JHIG + 1 IHIGT (JHIG) = ILOWT (2) JHIG = JHIG + 1 IHIGT (JHIG) = ILOWT (1) End If Exit ! --- Case (3) ! ! IWRK1 = ILOWT (1) IWRK2 = ILOWT (2) IWRK3 = ILOWT (3) If (XDONT(IWRK2) > XDONT(IWRK1)) Then ILOWT (1) = IWRK2 ILOWT (2) = IWRK1 IWRK2 = IWRK1 End If If (XDONT(IWRK2) < XDONT(IWRK3)) Then ILOWT (3) = IWRK2 ILOWT (2) = IWRK3 IWRK2 = IWRK3 If (XDONT(IWRK2) > XDONT(ILOWT(1))) Then ILOWT (2) = ILOWT (1) ILOWT (1) = IWRK2 End If End If JLOW = 0 Do ICRS = JHIG + 1, NORD JLOW = JLOW + 1 IHIGT (ICRS) = ILOWT (JLOW) End Do JHIG = NORD Exit ! --- Case (4:) ! ! XPIV0 = XPIV IFIN = JLOW ! ! One chooses a pivot from the 2 first values and the last one. ! This should ensure sufficient renewal between iterations to ! avoid worst case behavior effects. ! IWRK1 = ILOWT (1) IWRK2 = ILOWT (2) IWRK3 = ILOWT (IFIN) If (XDONT(IWRK2) > XDONT(IWRK1)) Then ILOWT (1) = IWRK2 ILOWT (2) = IWRK1 IWRK2 = IWRK1 End If If (XDONT(IWRK2) < XDONT(IWRK3)) Then ILOWT (IFIN) = IWRK2 ILOWT (2) = IWRK3 IWRK2 = IWRK3 If (XDONT(IWRK2) > XDONT(IHIGT(1))) Then ILOWT (2) = ILOWT (1) ILOWT (1) = IWRK2 End If End If ! JDEB = JHIG NWRK = NORD - JHIG IWRK1 = ILOWT (1) JHIG = JHIG + 1 IHIGT (JHIG) = IWRK1 XPIV = XDONT (IWRK1) + REAL (NWRK) / REAL (NORD+NWRK) * & (XDONT(ILOWT(IFIN))-XDONT(IWRK1)) ! ! One takes values >= pivot to IHIGT ! Again, 2 parts, one where we take care of the remaining ! low values because we might still need them, and the ! other when we know that we will have more than enough ! high values in the end. ! --- JLOW = 0 Do ICRS = 2, IFIN If (XDONT(ILOWT(ICRS)) >= XPIV) Then JHIG = JHIG + 1 IHIGT (JHIG) = ILOWT (ICRS) If (JHIG >= NORD) Exit Else JLOW = JLOW + 1 ILOWT (JLOW) = ILOWT (ICRS) End If End Do ! Do ICRS = ICRS + 1, IFIN If (XDONT(ILOWT(ICRS)) >= XPIV) Then JHIG = JHIG + 1 IHIGT (JHIG) = ILOWT (ICRS) End If End Do End Select ! --- ! Case (1) ! ! Only 1 value is missing in high part ! XMAX = XDONT (ILOWT(1)) ILOW = 1 Do ICRS = 2, JLOW If (XDONT(ILOWT(ICRS)) > XMAX) Then XMAX = XDONT (ILOWT(ICRS)) ILOW = ICRS End If End Do ! JHIG = JHIG + 1 IHIGT (JHIG) = ILOWT (ILOW) Exit ! ! Case (0) ! ! Low part is exactly what we want ! Exit ! --- ! Case (-5:-1) ! ! Only few values too many in high part ! IRNGT (1) = IHIGT (1) Do ICRS = 2, NORD IWRK = IHIGT (ICRS) XWRK = XDONT (IWRK) Do IDCR = ICRS - 1, 1, - 1 If (XWRK > XDONT(IRNGT(IDCR))) Then IRNGT (IDCR+1) = IRNGT (IDCR) Else Exit End If End Do IRNGT (IDCR+1) = IWRK End Do ! XWRK1 = XDONT (IRNGT(NORD)) Do ICRS = NORD + 1, JHIG If (XDONT(IHIGT (ICRS)) > XWRK1) Then XWRK = XDONT (IHIGT (ICRS)) Do IDCR = NORD - 1, 1, - 1 If (XWRK <= XDONT(IRNGT(IDCR))) Exit IRNGT (IDCR+1) = IRNGT (IDCR) End Do IRNGT (IDCR+1) = IHIGT (ICRS) XWRK1 = XDONT (IRNGT(NORD)) End If End Do ! Return ! ! Case (:-6) ! ! last case: too many values in high part ! --- IDEB = JDEB + 1 IMIL = (JHIG+IDEB) / 2 IFIN = JHIG ! --- ! One chooses a pivot from 1st, last, and middle values ! If (XDONT(IHIGT(IMIL)) > XDONT(IHIGT(IDEB))) Then IWRK = IHIGT (IDEB) IHIGT (IDEB) = IHIGT (IMIL) IHIGT (IMIL) = IWRK End If If (XDONT(IHIGT(IMIL)) < XDONT(IHIGT(IFIN))) Then IWRK = IHIGT (IFIN) IHIGT (IFIN) = IHIGT (IMIL) IHIGT (IMIL) = IWRK If (XDONT(IHIGT(IMIL)) > XDONT(IHIGT(IDEB))) Then IWRK = IHIGT (IDEB) IHIGT (IDEB) = IHIGT (IMIL) IHIGT (IMIL) = IWRK End If End If If (IFIN <= 3) Exit ! --- XPIV = XDONT (IHIGT(1)) + REAL(NORD)/REAL(JHIG+NORD) * & (XDONT(IHIGT(IFIN))-XDONT(IHIGT(1))) If (JDEB > 0) Then If (XPIV <= XPIV0) & XPIV = XPIV0 + REAL(2*NORD-JDEB)/REAL (JHIG+NORD) * & (XDONT(IHIGT(IFIN))-XPIV0) Else IDEB = 1 End If ! ! One takes values < XPIV to ILOWT ! However, we do not process the first values if we have been ! through the case when we did not have enough high values ! --- JLOW = 0 JHIG = JDEB ! --- If (XDONT(IHIGT(IFIN)) < XPIV) Then ICRS = JDEB Do ICRS = ICRS + 1 If (XDONT(IHIGT(ICRS)) < XPIV) Then JLOW = JLOW + 1 ILOWT (JLOW) = IHIGT (ICRS) If (ICRS >= IFIN) Exit Else JHIG = JHIG + 1 IHIGT (JHIG) = IHIGT (ICRS) If (JHIG >= NORD) Exit End If End Do ! --- If (ICRS < IFIN) Then Do ICRS = ICRS + 1 If (XDONT(IHIGT(ICRS)) >= XPIV) Then JHIG = JHIG + 1 IHIGT (JHIG) = IHIGT (ICRS) Else If (ICRS >= IFIN) Exit End If End Do End If Else Do ICRS = IDEB, IFIN If (XDONT(IHIGT(ICRS)) < XPIV) Then JLOW = JLOW + 1 ILOWT (JLOW) = IHIGT (ICRS) Else JHIG = JHIG + 1 IHIGT (JHIG) = IHIGT (ICRS) If (JHIG >= NORD) Exit End If End Do ! Do ICRS = ICRS + 1, IFIN If (XDONT(IHIGT(ICRS)) >= XPIV) Then JHIG = JHIG + 1 IHIGT (JHIG) = IHIGT (ICRS) End If End Do End If ! End Select ! End Do ! --- ! Now, we only need to complete ranking of the 1:NORD set ! Assuming NORD is small, we use a simple insertion sort ! IRNGT (1) = IHIGT (1) Do ICRS = 2, NORD IWRK = IHIGT (ICRS) XWRK = XDONT (IWRK) Do IDCR = ICRS - 1, 1, - 1 If (XWRK > XDONT(IRNGT(IDCR))) Then IRNGT (IDCR+1) = IRNGT (IDCR) Else Exit End If End Do IRNGT (IDCR+1) = IWRK End Do Return ! ! End Subroutine D_rapknr Subroutine R_rapknr (XDONT, IRNGT, NORD) ! Ranks partially XDONT by IRNGT, up to order NORD, in decreasing order. ! rapknr = (rnkpar backwards) ! __________________________________________________________ ! This routine uses a pivoting strategy such as the one of ! finding the median based on the quicksort algorithm, but ! we skew the pivot choice to try to bring it to NORD as ! fast as possible. It uses 2 temporary arrays, where it ! stores the indices of the values larger than the pivot ! (IHIGT), and the indices of values smaller than the pivot ! that we might still need later on (ILOWT). It iterates ! until it can bring the number of values in IHIGT to ! exactly NORD, and then uses an insertion sort to rank ! this set, since it is supposedly small. ! Michel Olagnon - Feb. 2011 ! __________________________________________________________ ! __________________________________________________________ ! _________________________________________________________ Real, Dimension (:), Intent (In) :: XDONT Integer, Dimension (:), Intent (Out) :: IRNGT Integer, Intent (In) :: NORD ! __________________________________________________________ Real :: XPIV, XPIV0, XWRK, XWRK1, XMIN, XMAX ! Integer, Dimension (SIZE(XDONT)) :: ILOWT, IHIGT Integer :: NDON, JHIG, JLOW, IHIG, IWRK, IWRK1, IWRK2, IWRK3 Integer :: IDEB, JDEB, IMIL, IFIN, NWRK, ICRS, IDCR, ILOW Integer :: JLM2, JLM1, JHM2, JHM1 ! NDON = SIZE (XDONT) ! ! First loop is used to fill-in ILOWT, IHIGT at the same time ! If (NDON < 2) Then If (NORD >= 1) IRNGT (1) = 1 Return End If ! ! One chooses a pivot, best estimate possible to put fractile near ! mid-point of the set of high values. ! If (XDONT(2) < XDONT(1)) Then ILOWT (1) = 2 IHIGT (1) = 1 Else ILOWT (1) = 1 IHIGT (1) = 2 End If ! If (NDON < 3) Then If (NORD >= 1) IRNGT (1) = IHIGT (1) If (NORD >= 2) IRNGT (2) = ILOWT (1) Return End If ! --- If (XDONT(3) > XDONT(ILOWT(1))) Then ILOWT (2) = ILOWT (1) If (XDONT(3) > XDONT(IHIGT(1))) Then ILOWT (1) = IHIGT (1) IHIGT (1) = 3 Else ILOWT (1) = 3 End If Else ILOWT (2) = 3 End If ! --- If (NDON < 4) Then If (NORD >= 1) IRNGT (1) = IHIGT (1) If (NORD >= 2) IRNGT (2) = ILOWT (1) If (NORD >= 3) IRNGT (3) = ILOWT (2) Return End If ! If (XDONT(NDON) > XDONT(ILOWT(1))) Then ILOWT (3) = ILOWT (2) ILOWT (2) = ILOWT (1) If (XDONT(NDON) > XDONT(IHIGT(1))) Then ILOWT (1) = IHIGT (1) IHIGT (1) = NDON Else ILOWT (1) = NDON End If Else if (XDONT (NDON) > XDONT (ILOWT(2))) Then ILOWT (3) = ILOWT (2) ILOWT (2) = NDON else ILOWT (3) = NDON endif End If ! If (NDON < 5) Then If (NORD >= 1) IRNGT (1) = IHIGT (1) If (NORD >= 2) IRNGT (2) = ILOWT (1) If (NORD >= 3) IRNGT (3) = ILOWT (2) If (NORD >= 4) IRNGT (4) = ILOWT (3) Return End If ! --- JDEB = 0 IDEB = JDEB + 1 JHIG = IDEB JLOW = 3 XPIV = XDONT (IHIGT(IDEB)) + REAL(2*NORD)/REAL(NDON+NORD) * & (XDONT(ILOWT(3))-XDONT(IHIGT(IDEB))) If (XPIV >= XDONT(ILOWT(1))) Then XPIV = XDONT (IHIGT(IDEB)) + REAL(2*NORD)/REAL(NDON+NORD) * & (XDONT(ILOWT(2))-XDONT(IHIGT(IDEB))) If (XPIV >= XDONT(ILOWT(1))) & XPIV = XDONT (IHIGT(IDEB)) + REAL (2*NORD) / REAL (NDON+NORD) * & (XDONT(ILOWT(1))-XDONT(IHIGT(IDEB))) End If XPIV0 = XPIV ! --- ! One puts values < pivot in the end and those >= pivot ! at the beginning. This is split in 2 cases, so that ! we can skip the loop test a number of times. ! As we are also filling in the work arrays at the same time ! we stop filling in the ILOWT array as soon as we have more ! than enough values in IHIGT. ! ! If (XDONT(NDON) < XPIV) Then ICRS = 3 Do ICRS = ICRS + 1 If (XDONT(ICRS) < XPIV) Then If (ICRS >= NDON) Exit JLOW = JLOW + 1 ILOWT (JLOW) = ICRS Else JHIG = JHIG + 1 IHIGT (JHIG) = ICRS If (JHIG >= NORD) Exit End If End Do ! ! One restricts further processing because it is no use ! to store more low values ! If (ICRS < NDON-1) Then Do ICRS = ICRS + 1 If (XDONT(ICRS) >= XPIV) Then JHIG = JHIG + 1 IHIGT (JHIG) = ICRS Else If (ICRS >= NDON) Then Exit End If End Do End If ! ! --- Else ! ! Same as above, but this is not as easy to optimize, so the ! DO-loop is kept ! Do ICRS = 4, NDON - 1 If (XDONT(ICRS) < XPIV) Then JLOW = JLOW + 1 ILOWT (JLOW) = ICRS Else JHIG = JHIG + 1 IHIGT (JHIG) = ICRS If (JHIG >= NORD) Exit End If End Do ! If (ICRS < NDON-1) Then Do ICRS = ICRS + 1 If (XDONT(ICRS) >= XPIV) Then If (ICRS >= NDON) Exit JHIG = JHIG + 1 IHIGT (JHIG) = ICRS End If End Do End If End If ! --- JLM2 = 0 JLM1 = 0 JHM2 = 0 JHM1 = 0 Do if (JHIG == NORD) Exit If (JHM2 == JHIG .And. JLM2 == JLOW) Then ! ! We are oscillating. Perturbate by bringing JHIG closer by one ! to NORD ! If (NORD > JHIG) Then XMAX = XDONT (ILOWT(1)) ILOW = 1 Do ICRS = 2, JLOW If (XDONT(ILOWT(ICRS)) > XMAX) Then XMAX = XDONT (ILOWT(ICRS)) ILOW = ICRS End If End Do ! JHIG = JHIG + 1 IHIGT (JHIG) = ILOWT (ILOW) ILOWT (ILOW) = ILOWT (JLOW) JLOW = JLOW - 1 Else IHIG = IHIGT (JHIG) XMIN = XDONT (IHIG) Do ICRS = 1, JHIG If (XDONT(IHIGT(ICRS)) < XMIN) Then IWRK = IHIGT (ICRS) XMIN = XDONT (IWRK) IHIGT (ICRS) = IHIG IHIG = IWRK End If End Do JHIG = JHIG - 1 End If End If JLM2 = JLM1 JLM1 = JLOW JHM2 = JHM1 JHM1 = JHIG ! --- ! We try to bring the number of values in the high values set ! closer to NORD. ! Select Case (NORD-JHIG) Case (2:) ! ! Not enough values in low part, at least 2 are missing ! Select Case (JLOW) !!!!! CASE DEFAULT !!!!! write (*,*) "Assertion failed" !!!!! STOP ! ! We make a special case when we have so few values in ! the low values set that it is bad performance to choose a pivot ! and apply the general algorithm. ! Case (2) If (XDONT(ILOWT(1)) >= XDONT(ILOWT(2))) Then JHIG = JHIG + 1 IHIGT (JHIG) = ILOWT (1) JHIG = JHIG + 1 IHIGT (JHIG) = ILOWT (2) Else JHIG = JHIG + 1 IHIGT (JHIG) = ILOWT (2) JHIG = JHIG + 1 IHIGT (JHIG) = ILOWT (1) End If Exit ! --- Case (3) ! ! IWRK1 = ILOWT (1) IWRK2 = ILOWT (2) IWRK3 = ILOWT (3) If (XDONT(IWRK2) > XDONT(IWRK1)) Then ILOWT (1) = IWRK2 ILOWT (2) = IWRK1 IWRK2 = IWRK1 End If If (XDONT(IWRK2) < XDONT(IWRK3)) Then ILOWT (3) = IWRK2 ILOWT (2) = IWRK3 IWRK2 = IWRK3 If (XDONT(IWRK2) > XDONT(ILOWT(1))) Then ILOWT (2) = ILOWT (1) ILOWT (1) = IWRK2 End If End If JLOW = 0 Do ICRS = JHIG + 1, NORD JLOW = JLOW + 1 IHIGT (ICRS) = ILOWT (JLOW) End Do JHIG = NORD Exit ! --- Case (4:) ! ! XPIV0 = XPIV IFIN = JLOW ! ! One chooses a pivot from the 2 first values and the last one. ! This should ensure sufficient renewal between iterations to ! avoid worst case behavior effects. ! IWRK1 = ILOWT (1) IWRK2 = ILOWT (2) IWRK3 = ILOWT (IFIN) If (XDONT(IWRK2) > XDONT(IWRK1)) Then ILOWT (1) = IWRK2 ILOWT (2) = IWRK1 IWRK2 = IWRK1 End If If (XDONT(IWRK2) < XDONT(IWRK3)) Then ILOWT (IFIN) = IWRK2 ILOWT (2) = IWRK3 IWRK2 = IWRK3 If (XDONT(IWRK2) > XDONT(IHIGT(1))) Then ILOWT (2) = ILOWT (1) ILOWT (1) = IWRK2 End If End If ! JDEB = JHIG NWRK = NORD - JHIG IWRK1 = ILOWT (1) JHIG = JHIG + 1 IHIGT (JHIG) = IWRK1 XPIV = XDONT (IWRK1) + REAL (NWRK) / REAL (NORD+NWRK) * & (XDONT(ILOWT(IFIN))-XDONT(IWRK1)) ! ! One takes values >= pivot to IHIGT ! Again, 2 parts, one where we take care of the remaining ! low values because we might still need them, and the ! other when we know that we will have more than enough ! high values in the end. ! --- JLOW = 0 Do ICRS = 2, IFIN If (XDONT(ILOWT(ICRS)) >= XPIV) Then JHIG = JHIG + 1 IHIGT (JHIG) = ILOWT (ICRS) If (JHIG >= NORD) Exit Else JLOW = JLOW + 1 ILOWT (JLOW) = ILOWT (ICRS) End If End Do ! Do ICRS = ICRS + 1, IFIN If (XDONT(ILOWT(ICRS)) >= XPIV) Then JHIG = JHIG + 1 IHIGT (JHIG) = ILOWT (ICRS) End If End Do End Select ! --- ! Case (1) ! ! Only 1 value is missing in high part ! XMAX = XDONT (ILOWT(1)) ILOW = 1 Do ICRS = 2, JLOW If (XDONT(ILOWT(ICRS)) > XMAX) Then XMAX = XDONT (ILOWT(ICRS)) ILOW = ICRS End If End Do ! JHIG = JHIG + 1 IHIGT (JHIG) = ILOWT (ILOW) Exit ! ! Case (0) ! ! Low part is exactly what we want ! Exit ! --- ! Case (-5:-1) ! ! Only few values too many in high part ! IRNGT (1) = IHIGT (1) Do ICRS = 2, NORD IWRK = IHIGT (ICRS) XWRK = XDONT (IWRK) Do IDCR = ICRS - 1, 1, - 1 If (XWRK > XDONT(IRNGT(IDCR))) Then IRNGT (IDCR+1) = IRNGT (IDCR) Else Exit End If End Do IRNGT (IDCR+1) = IWRK End Do ! XWRK1 = XDONT (IRNGT(NORD)) Do ICRS = NORD + 1, JHIG If (XDONT(IHIGT (ICRS)) > XWRK1) Then XWRK = XDONT (IHIGT (ICRS)) Do IDCR = NORD - 1, 1, - 1 If (XWRK <= XDONT(IRNGT(IDCR))) Exit IRNGT (IDCR+1) = IRNGT (IDCR) End Do IRNGT (IDCR+1) = IHIGT (ICRS) XWRK1 = XDONT (IRNGT(NORD)) End If End Do ! Return ! ! Case (:-6) ! ! last case: too many values in high part ! --- IDEB = JDEB + 1 IMIL = (JHIG+IDEB) / 2 IFIN = JHIG ! --- ! One chooses a pivot from 1st, last, and middle values ! If (XDONT(IHIGT(IMIL)) > XDONT(IHIGT(IDEB))) Then IWRK = IHIGT (IDEB) IHIGT (IDEB) = IHIGT (IMIL) IHIGT (IMIL) = IWRK End If If (XDONT(IHIGT(IMIL)) < XDONT(IHIGT(IFIN))) Then IWRK = IHIGT (IFIN) IHIGT (IFIN) = IHIGT (IMIL) IHIGT (IMIL) = IWRK If (XDONT(IHIGT(IMIL)) > XDONT(IHIGT(IDEB))) Then IWRK = IHIGT (IDEB) IHIGT (IDEB) = IHIGT (IMIL) IHIGT (IMIL) = IWRK End If End If If (IFIN <= 3) Exit ! --- XPIV = XDONT (IHIGT(1)) + REAL(NORD)/REAL(JHIG+NORD) * & (XDONT(IHIGT(IFIN))-XDONT(IHIGT(1))) If (JDEB > 0) Then If (XPIV <= XPIV0) & XPIV = XPIV0 + REAL(2*NORD-JDEB)/REAL (JHIG+NORD) * & (XDONT(IHIGT(IFIN))-XPIV0) Else IDEB = 1 End If ! ! One takes values < XPIV to ILOWT ! However, we do not process the first values if we have been ! through the case when we did not have enough high values ! --- JLOW = 0 JHIG = JDEB ! --- If (XDONT(IHIGT(IFIN)) < XPIV) Then ICRS = JDEB Do ICRS = ICRS + 1 If (XDONT(IHIGT(ICRS)) < XPIV) Then JLOW = JLOW + 1 ILOWT (JLOW) = IHIGT (ICRS) If (ICRS >= IFIN) Exit Else JHIG = JHIG + 1 IHIGT (JHIG) = IHIGT (ICRS) If (JHIG >= NORD) Exit End If End Do ! --- If (ICRS < IFIN) Then Do ICRS = ICRS + 1 If (XDONT(IHIGT(ICRS)) >= XPIV) Then JHIG = JHIG + 1 IHIGT (JHIG) = IHIGT (ICRS) Else If (ICRS >= IFIN) Exit End If End Do End If Else Do ICRS = IDEB, IFIN If (XDONT(IHIGT(ICRS)) < XPIV) Then JLOW = JLOW + 1 ILOWT (JLOW) = IHIGT (ICRS) Else JHIG = JHIG + 1 IHIGT (JHIG) = IHIGT (ICRS) If (JHIG >= NORD) Exit End If End Do ! Do ICRS = ICRS + 1, IFIN If (XDONT(IHIGT(ICRS)) >= XPIV) Then JHIG = JHIG + 1 IHIGT (JHIG) = IHIGT (ICRS) End If End Do End If ! End Select ! End Do ! --- ! Now, we only need to complete ranking of the 1:NORD set ! Assuming NORD is small, we use a simple insertion sort ! IRNGT (1) = IHIGT (1) Do ICRS = 2, NORD IWRK = IHIGT (ICRS) XWRK = XDONT (IWRK) Do IDCR = ICRS - 1, 1, - 1 If (XWRK > XDONT(IRNGT(IDCR))) Then IRNGT (IDCR+1) = IRNGT (IDCR) Else Exit End If End Do IRNGT (IDCR+1) = IWRK End Do Return ! ! End Subroutine R_rapknr Subroutine I_rapknr (XDONT, IRNGT, NORD) ! Ranks partially XDONT by IRNGT, up to order NORD, in decreasing order. ! rapknr = (rnkpar backwards) ! __________________________________________________________ ! This routine uses a pivoting strategy such as the one of ! finding the median based on the quicksort algorithm, but ! we skew the pivot choice to try to bring it to NORD as ! fast as possible. It uses 2 temporary arrays, where it ! stores the indices of the values larger than the pivot ! (IHIGT), and the indices of values smaller than the pivot ! that we might still need later on (ILOWT). It iterates ! until it can bring the number of values in IHIGT to ! exactly NORD, and then uses an insertion sort to rank ! this set, since it is supposedly small. ! Michel Olagnon - Feb. 2011 ! __________________________________________________________ ! __________________________________________________________ ! __________________________________________________________ Integer, Dimension (:), Intent (In) :: XDONT Integer, Dimension (:), Intent (Out) :: IRNGT Integer, Intent (In) :: NORD ! __________________________________________________________ Integer :: XPIV, XPIV0, XWRK, XWRK1, XMIN, XMAX ! Integer, Dimension (SIZE(XDONT)) :: ILOWT, IHIGT Integer :: NDON, JHIG, JLOW, IHIG, IWRK, IWRK1, IWRK2, IWRK3 Integer :: IDEB, JDEB, IMIL, IFIN, NWRK, ICRS, IDCR, ILOW Integer :: JLM2, JLM1, JHM2, JHM1 ! NDON = SIZE (XDONT) ! ! First loop is used to fill-in ILOWT, IHIGT at the same time ! If (NDON < 2) Then If (NORD >= 1) IRNGT (1) = 1 Return End If ! ! One chooses a pivot, best estimate possible to put fractile near ! mid-point of the set of high values. ! If (XDONT(2) < XDONT(1)) Then ILOWT (1) = 2 IHIGT (1) = 1 Else ILOWT (1) = 1 IHIGT (1) = 2 End If ! If (NDON < 3) Then If (NORD >= 1) IRNGT (1) = IHIGT (1) If (NORD >= 2) IRNGT (2) = ILOWT (1) Return End If ! --- If (XDONT(3) > XDONT(ILOWT(1))) Then ILOWT (2) = ILOWT (1) If (XDONT(3) > XDONT(IHIGT(1))) Then ILOWT (1) = IHIGT (1) IHIGT (1) = 3 Else ILOWT (1) = 3 End If Else ILOWT (2) = 3 End If ! --- If (NDON < 4) Then If (NORD >= 1) IRNGT (1) = IHIGT (1) If (NORD >= 2) IRNGT (2) = ILOWT (1) If (NORD >= 3) IRNGT (3) = ILOWT (2) Return End If ! If (XDONT(NDON) > XDONT(ILOWT(1))) Then ILOWT (3) = ILOWT (2) ILOWT (2) = ILOWT (1) If (XDONT(NDON) > XDONT(IHIGT(1))) Then ILOWT (1) = IHIGT (1) IHIGT (1) = NDON Else ILOWT (1) = NDON End If Else if (XDONT (NDON) > XDONT (ILOWT(2))) Then ILOWT (3) = ILOWT (2) ILOWT (2) = NDON else ILOWT (3) = NDON endif End If ! If (NDON < 5) Then If (NORD >= 1) IRNGT (1) = IHIGT (1) If (NORD >= 2) IRNGT (2) = ILOWT (1) If (NORD >= 3) IRNGT (3) = ILOWT (2) If (NORD >= 4) IRNGT (4) = ILOWT (3) Return End If ! --- JDEB = 0 IDEB = JDEB + 1 JHIG = IDEB JLOW = 3 XPIV = XDONT (IHIGT(IDEB)) + REAL(2*NORD)/REAL(NDON+NORD) * & (XDONT(ILOWT(3))-XDONT(IHIGT(IDEB))) If (XPIV >= XDONT(ILOWT(1))) Then XPIV = XDONT (IHIGT(IDEB)) + REAL(2*NORD)/REAL(NDON+NORD) * & (XDONT(ILOWT(2))-XDONT(IHIGT(IDEB))) If (XPIV >= XDONT(ILOWT(1))) & XPIV = XDONT (IHIGT(IDEB)) + REAL (2*NORD) / REAL (NDON+NORD) * & (XDONT(ILOWT(1))-XDONT(IHIGT(IDEB))) End If XPIV0 = XPIV ! --- ! One puts values < pivot in the end and those >= pivot ! at the beginning. This is split in 2 cases, so that ! we can skip the loop test a number of times. ! As we are also filling in the work arrays at the same time ! we stop filling in the ILOWT array as soon as we have more ! than enough values in IHIGT. ! ! If (XDONT(NDON) < XPIV) Then ICRS = 3 Do ICRS = ICRS + 1 If (XDONT(ICRS) < XPIV) Then If (ICRS >= NDON) Exit JLOW = JLOW + 1 ILOWT (JLOW) = ICRS Else JHIG = JHIG + 1 IHIGT (JHIG) = ICRS If (JHIG >= NORD) Exit End If End Do ! ! One restricts further processing because it is no use ! to store more low values ! If (ICRS < NDON-1) Then Do ICRS = ICRS + 1 If (XDONT(ICRS) >= XPIV) Then JHIG = JHIG + 1 IHIGT (JHIG) = ICRS Else If (ICRS >= NDON) Then Exit End If End Do End If ! ! --- Else ! ! Same as above, but this is not as easy to optimize, so the ! DO-loop is kept ! Do ICRS = 4, NDON - 1 If (XDONT(ICRS) < XPIV) Then JLOW = JLOW + 1 ILOWT (JLOW) = ICRS Else JHIG = JHIG + 1 IHIGT (JHIG) = ICRS If (JHIG >= NORD) Exit End If End Do ! If (ICRS < NDON-1) Then Do ICRS = ICRS + 1 If (XDONT(ICRS) >= XPIV) Then If (ICRS >= NDON) Exit JHIG = JHIG + 1 IHIGT (JHIG) = ICRS End If End Do End If End If ! --- JLM2 = 0 JLM1 = 0 JHM2 = 0 JHM1 = 0 Do if (JHIG == NORD) Exit If (JHM2 == JHIG .And. JLM2 == JLOW) Then ! ! We are oscillating. Perturbate by bringing JHIG closer by one ! to NORD ! If (NORD > JHIG) Then XMAX = XDONT (ILOWT(1)) ILOW = 1 Do ICRS = 2, JLOW If (XDONT(ILOWT(ICRS)) > XMAX) Then XMAX = XDONT (ILOWT(ICRS)) ILOW = ICRS End If End Do ! JHIG = JHIG + 1 IHIGT (JHIG) = ILOWT (ILOW) ILOWT (ILOW) = ILOWT (JLOW) JLOW = JLOW - 1 Else IHIG = IHIGT (JHIG) XMIN = XDONT (IHIG) Do ICRS = 1, JHIG If (XDONT(IHIGT(ICRS)) < XMIN) Then IWRK = IHIGT (ICRS) XMIN = XDONT (IWRK) IHIGT (ICRS) = IHIG IHIG = IWRK End If End Do JHIG = JHIG - 1 End If End If JLM2 = JLM1 JLM1 = JLOW JHM2 = JHM1 JHM1 = JHIG ! --- ! We try to bring the number of values in the high values set ! closer to NORD. ! Select Case (NORD-JHIG) Case (2:) ! ! Not enough values in low part, at least 2 are missing ! Select Case (JLOW) !!!!! CASE DEFAULT !!!!! write (*,*) "Assertion failed" !!!!! STOP ! ! We make a special case when we have so few values in ! the low values set that it is bad performance to choose a pivot ! and apply the general algorithm. ! Case (2) If (XDONT(ILOWT(1)) >= XDONT(ILOWT(2))) Then JHIG = JHIG + 1 IHIGT (JHIG) = ILOWT (1) JHIG = JHIG + 1 IHIGT (JHIG) = ILOWT (2) Else JHIG = JHIG + 1 IHIGT (JHIG) = ILOWT (2) JHIG = JHIG + 1 IHIGT (JHIG) = ILOWT (1) End If Exit ! --- Case (3) ! ! IWRK1 = ILOWT (1) IWRK2 = ILOWT (2) IWRK3 = ILOWT (3) If (XDONT(IWRK2) > XDONT(IWRK1)) Then ILOWT (1) = IWRK2 ILOWT (2) = IWRK1 IWRK2 = IWRK1 End If If (XDONT(IWRK2) < XDONT(IWRK3)) Then ILOWT (3) = IWRK2 ILOWT (2) = IWRK3 IWRK2 = IWRK3 If (XDONT(IWRK2) > XDONT(ILOWT(1))) Then ILOWT (2) = ILOWT (1) ILOWT (1) = IWRK2 End If End If JLOW = 0 Do ICRS = JHIG + 1, NORD JLOW = JLOW + 1 IHIGT (ICRS) = ILOWT (JLOW) End Do JHIG = NORD Exit ! --- Case (4:) ! ! XPIV0 = XPIV IFIN = JLOW ! ! One chooses a pivot from the 2 first values and the last one. ! This should ensure sufficient renewal between iterations to ! avoid worst case behavior effects. ! IWRK1 = ILOWT (1) IWRK2 = ILOWT (2) IWRK3 = ILOWT (IFIN) If (XDONT(IWRK2) > XDONT(IWRK1)) Then ILOWT (1) = IWRK2 ILOWT (2) = IWRK1 IWRK2 = IWRK1 End If If (XDONT(IWRK2) < XDONT(IWRK3)) Then ILOWT (IFIN) = IWRK2 ILOWT (2) = IWRK3 IWRK2 = IWRK3 If (XDONT(IWRK2) > XDONT(IHIGT(1))) Then ILOWT (2) = ILOWT (1) ILOWT (1) = IWRK2 End If End If ! JDEB = JHIG NWRK = NORD - JHIG IWRK1 = ILOWT (1) JHIG = JHIG + 1 IHIGT (JHIG) = IWRK1 XPIV = XDONT (IWRK1) + REAL (NWRK) / REAL (NORD+NWRK) * & (XDONT(ILOWT(IFIN))-XDONT(IWRK1)) ! ! One takes values >= pivot to IHIGT ! Again, 2 parts, one where we take care of the remaining ! low values because we might still need them, and the ! other when we know that we will have more than enough ! high values in the end. ! --- JLOW = 0 Do ICRS = 2, IFIN If (XDONT(ILOWT(ICRS)) >= XPIV) Then JHIG = JHIG + 1 IHIGT (JHIG) = ILOWT (ICRS) If (JHIG >= NORD) Exit Else JLOW = JLOW + 1 ILOWT (JLOW) = ILOWT (ICRS) End If End Do ! Do ICRS = ICRS + 1, IFIN If (XDONT(ILOWT(ICRS)) >= XPIV) Then JHIG = JHIG + 1 IHIGT (JHIG) = ILOWT (ICRS) End If End Do End Select ! --- ! Case (1) ! ! Only 1 value is missing in high part ! XMAX = XDONT (ILOWT(1)) ILOW = 1 Do ICRS = 2, JLOW If (XDONT(ILOWT(ICRS)) > XMAX) Then XMAX = XDONT (ILOWT(ICRS)) ILOW = ICRS End If End Do ! JHIG = JHIG + 1 IHIGT (JHIG) = ILOWT (ILOW) Exit ! ! Case (0) ! ! Low part is exactly what we want ! Exit ! --- ! Case (-5:-1) ! ! Only few values too many in high part ! IRNGT (1) = IHIGT (1) Do ICRS = 2, NORD IWRK = IHIGT (ICRS) XWRK = XDONT (IWRK) Do IDCR = ICRS - 1, 1, - 1 If (XWRK > XDONT(IRNGT(IDCR))) Then IRNGT (IDCR+1) = IRNGT (IDCR) Else Exit End If End Do IRNGT (IDCR+1) = IWRK End Do ! XWRK1 = XDONT (IRNGT(NORD)) Do ICRS = NORD + 1, JHIG If (XDONT(IHIGT (ICRS)) > XWRK1) Then XWRK = XDONT (IHIGT (ICRS)) Do IDCR = NORD - 1, 1, - 1 If (XWRK <= XDONT(IRNGT(IDCR))) Exit IRNGT (IDCR+1) = IRNGT (IDCR) End Do IRNGT (IDCR+1) = IHIGT (ICRS) XWRK1 = XDONT (IRNGT(NORD)) End If End Do ! Return ! ! Case (:-6) ! ! last case: too many values in high part ! --- IDEB = JDEB + 1 IMIL = (JHIG+IDEB) / 2 IFIN = JHIG ! --- ! One chooses a pivot from 1st, last, and middle values ! If (XDONT(IHIGT(IMIL)) > XDONT(IHIGT(IDEB))) Then IWRK = IHIGT (IDEB) IHIGT (IDEB) = IHIGT (IMIL) IHIGT (IMIL) = IWRK End If If (XDONT(IHIGT(IMIL)) < XDONT(IHIGT(IFIN))) Then IWRK = IHIGT (IFIN) IHIGT (IFIN) = IHIGT (IMIL) IHIGT (IMIL) = IWRK If (XDONT(IHIGT(IMIL)) > XDONT(IHIGT(IDEB))) Then IWRK = IHIGT (IDEB) IHIGT (IDEB) = IHIGT (IMIL) IHIGT (IMIL) = IWRK End If End If If (IFIN <= 3) Exit ! --- XPIV = XDONT (IHIGT(1)) + REAL(NORD)/REAL(JHIG+NORD) * & (XDONT(IHIGT(IFIN))-XDONT(IHIGT(1))) If (JDEB > 0) Then If (XPIV <= XPIV0) & XPIV = XPIV0 + REAL(2*NORD-JDEB)/REAL (JHIG+NORD) * & (XDONT(IHIGT(IFIN))-XPIV0) Else IDEB = 1 End If ! ! One takes values < XPIV to ILOWT ! However, we do not process the first values if we have been ! through the case when we did not have enough high values ! --- JLOW = 0 JHIG = JDEB ! --- If (XDONT(IHIGT(IFIN)) < XPIV) Then ICRS = JDEB Do ICRS = ICRS + 1 If (XDONT(IHIGT(ICRS)) < XPIV) Then JLOW = JLOW + 1 ILOWT (JLOW) = IHIGT (ICRS) If (ICRS >= IFIN) Exit Else JHIG = JHIG + 1 IHIGT (JHIG) = IHIGT (ICRS) If (JHIG >= NORD) Exit End If End Do ! --- If (ICRS < IFIN) Then Do ICRS = ICRS + 1 If (XDONT(IHIGT(ICRS)) >= XPIV) Then JHIG = JHIG + 1 IHIGT (JHIG) = IHIGT (ICRS) Else If (ICRS >= IFIN) Exit End If End Do End If Else Do ICRS = IDEB, IFIN If (XDONT(IHIGT(ICRS)) < XPIV) Then JLOW = JLOW + 1 ILOWT (JLOW) = IHIGT (ICRS) Else JHIG = JHIG + 1 IHIGT (JHIG) = IHIGT (ICRS) If (JHIG >= NORD) Exit End If End Do ! Do ICRS = ICRS + 1, IFIN If (XDONT(IHIGT(ICRS)) >= XPIV) Then JHIG = JHIG + 1 IHIGT (JHIG) = IHIGT (ICRS) End If End Do End If ! End Select ! End Do ! --- ! Now, we only need to complete ranking of the 1:NORD set ! Assuming NORD is small, we use a simple insertion sort ! IRNGT (1) = IHIGT (1) Do ICRS = 2, NORD IWRK = IHIGT (ICRS) XWRK = XDONT (IWRK) Do IDCR = ICRS - 1, 1, - 1 If (XWRK > XDONT(IRNGT(IDCR))) Then IRNGT (IDCR+1) = IRNGT (IDCR) Else Exit End If End Do IRNGT (IDCR+1) = IWRK End Do Return ! ! End Subroutine I_rapknr end module m_rapknr