Module m_indmed Integer, Parameter :: kdp = selected_real_kind(15) public :: indmed private :: kdp private :: R_indmed, I_indmed, D_indmed private :: r_med, i_med, d_med Integer, Allocatable, Dimension(:), Private, Save :: IDONT interface indmed module procedure d_indmed, r_indmed, i_indmed end interface indmed contains Subroutine D_indmed (XDONT, INDM) ! Returns index of median value of XDONT. ! __________________________________________________________ Real (kind=kdp), Dimension (:), Intent (In) :: XDONT Integer, Intent (Out) :: INDM ! __________________________________________________________ Integer :: IDON ! Allocate (IDONT (SIZE(XDONT))) Do IDON = 1, SIZE(XDONT) IDONT (IDON) = IDON End Do ! Call d_med (XDONT, IDONT, INDM) ! Deallocate (IDONT) End Subroutine D_indmed Recursive Subroutine d_med (XDATT, IDATT, ires_med) ! Finds the index of the median of XDONT using the recursive procedure ! described in Knuth, The Art of Computer Programming, ! vol. 3, 5.3.3 - This procedure is linear in time, and ! does not require to be able to interpolate in the ! set as the one used in INDNTH. It also has better worst ! case behavior than INDNTH, but is about 30% slower in ! average for random uniformly distributed values. ! __________________________________________________________ Real (kind=kdp), Dimension (:), Intent (In) :: XDATT Integer, Dimension (:), Intent (In) :: IDATT Integer, Intent (Out):: ires_med ! __________________________________________________________ ! Real (kind=kdp), Parameter :: XHUGE = HUGE (XDATT) Real (kind=kdp) :: XWRK, XWRK1, XMED7, XMAX, XMIN ! Integer, Dimension (7*(((Size (IDATT)+6)/7+6)/7)) :: ISTRT, IENDT, IMEDT Integer, Dimension (7*((Size(IDATT)+6)/7)) :: IWRKT Integer :: NTRI, NMED, NORD, NEQU, NLEQ, IMED, IDON, IDON1 Integer :: IDEB, ITMP, IDCR, ICRS, ICRS1, ICRS2, IMAX, IMIN Integer :: IWRK, IWRK1, IMED1, IMED7, NDAT ! NDAT = Size (IDATT) NMED = (NDAT+1) / 2 IWRKT = IDATT ! ! If the number of values is small, then use insertion sort ! If (NDAT < 35) Then ! ! Bring minimum to first location to save test in decreasing loop ! IDCR = NDAT If (XDATT (IWRKT (1)) < XDATT (IWRKT (IDCR))) Then IWRK = IWRKT (1) Else IWRK = IWRKT (IDCR) IWRKT (IDCR) = IWRKT (1) Endif XWRK = XDATT (IWRK) Do ITMP = 1, NDAT - 2 IDCR = IDCR - 1 IWRK1 = IWRKT (IDCR) XWRK1 = XDATT (IWRK1) If (XWRK1 < XWRK) Then IWRKT (IDCR) = IWRK XWRK = XWRK1 IWRK = IWRK1 Endif End Do IWRKT (1) = IWRK ! ! Sort the first half, until we have NMED sorted values ! Do ICRS = 3, NMED XWRK = XDATT (IWRKT (ICRS)) IWRK = IWRKT (ICRS) IDCR = ICRS - 1 Do If (XWRK >= XDATT (IWRKT(IDCR))) Exit IWRKT (IDCR+1) = IWRKT (IDCR) IDCR = IDCR - 1 End Do IWRKT (IDCR+1) = IWRK End Do ! ! Insert any value less than the current median in the first half ! XWRK1 = XDATT (IWRKT (NMED)) Do ICRS = NMED+1, NDAT XWRK = XDATT (IWRKT (ICRS)) IWRK = IWRKT (ICRS) If (XWRK < XWRK1) Then IDCR = NMED - 1 Do If (XWRK >= XDATT (IWRKT(IDCR))) Exit IWRKT (IDCR+1) = IWRKT (IDCR) IDCR = IDCR - 1 End Do IWRKT (IDCR+1) = IWRK XWRK1 = XDATT (IWRKT (NMED)) End If End Do ires_med = IWRKT (NMED) Return End If ! ! Make sorted subsets of 7 elements ! This is done by a variant of insertion sort where a first ! pass is used to bring the smallest element to the first position ! decreasing disorder at the same time, so that we may remove ! remove the loop test in the insertion loop. ! IMAX = 1 IMIN = 1 XMAX = XDATT (IWRKT(IMAX)) XMIN = XDATT (IWRKT(IMIN)) DO IDEB = 1, NDAT-6, 7 IDCR = IDEB + 6 If (XDATT (IWRKT(IDEB)) < XDATT (IWRKT(IDCR))) Then IWRK = IWRKT(IDEB) Else IWRK = IWRKT (IDCR) IWRKT (IDCR) = IWRKT(IDEB) Endif XWRK = XDATT (IWRK) Do ITMP = 1, 5 IDCR = IDCR - 1 IWRK1 = IWRKT (IDCR) XWRK1 = XDATT (IWRK1) If (XWRK1 < XWRK) Then IWRKT (IDCR) = IWRK IWRK = IWRK1 XWRK = XWRK1 Endif End Do IWRKT (IDEB) = IWRK If (XWRK < XMIN) Then IMIN = IWRK XMIN = XWRK End If Do ICRS = IDEB+1, IDEB+5 IWRK = IWRKT (ICRS+1) XWRK = XDATT (IWRK) IDON = IWRKT(ICRS) If (XWRK < XDATT(IDON)) Then IWRKT (ICRS+1) = IDON IDCR = ICRS IWRK1 = IWRKT (IDCR-1) XWRK1 = XDATT (IWRK1) Do If (XWRK >= XWRK1) Exit IWRKT (IDCR) = IWRK1 IDCR = IDCR - 1 IWRK1 = IWRKT (IDCR-1) XWRK1 = XDATT (IWRK1) End Do IWRKT (IDCR) = IWRK EndIf End Do If (XWRK > XMAX) Then IMAX = IWRK XMAX = XWRK End If End Do ! ! Add-up alternatively MAX and MIN values to make the number of data ! an exact multiple of 7. ! IDEB = 7 * (NDAT/7) NTRI = NDAT If (IDEB < NDAT) Then ! Do ICRS = IDEB+1, NDAT XWRK1 = XDATT (IWRKT (ICRS)) IF (XWRK1 > XMAX) Then IMAX = IWRKT (ICRS) XMAX = XWRK1 End If IF (XWRK1 < XMIN) Then IMIN = IWRKT (ICRS) XMIN = XWRK1 End If End Do IWRK1 = IMAX Do ICRS = NDAT+1, IDEB+7 IWRKT (ICRS) = IWRK1 If (IWRK1 == IMAX) Then IWRK1 = IMIN Else NMED = NMED + 1 IWRK1 = IMAX End If End Do ! Do ICRS = IDEB+2, IDEB+7 IWRK = IWRKT (ICRS) XWRK = XDATT (IWRK) Do IDCR = ICRS - 1, IDEB+1, - 1 If (XWRK >= XDATT (IWRKT(IDCR))) Exit IWRKT (IDCR+1) = IWRKT (IDCR) End Do IWRKT (IDCR+1) = IWRK End Do ! NTRI = IDEB+7 End If ! ! Make the set of the indices of median values of each sorted subset ! IDON1 = 0 Do IDON = 1, NTRI, 7 IDON1 = IDON1 + 1 IMEDT (IDON1) = IWRKT (IDON + 3) End Do ! ! Find XMED7, the median of the medians ! Call d_med (XDATT, IMEDT(1:IDON1), IMED7) XMED7 = XDATT (IMED7) ! ! Count how many values are not higher than (and how many equal to) XMED7 ! This number is at least 4 * 1/2 * (N/7) : 4 values in each of the ! subsets where the median is lower than the median of medians. For similar ! reasons, we also have at least 2N/7 values not lower than XMED7. At the ! same time, we find in each subset the index of the last value < XMED7, ! and that of the first > XMED7. These indices will be used to restrict the ! search for the median as the Kth element in the subset (> or <) where ! we know it to be. ! IDON1 = 1 NLEQ = 0 NEQU = 0 Do IDON = 1, NTRI, 7 IMED = IDON+3 If (XDATT (IWRKT (IMED)) > XMED7) Then IMED = IMED - 2 If (XDATT (IWRKT (IMED)) > XMED7) Then IMED = IMED - 1 Else If (XDATT (IWRKT (IMED)) < XMED7) Then IMED = IMED + 1 Endif Else If (XDATT (IWRKT (IMED)) < XMED7) Then IMED = IMED + 2 If (XDATT (IWRKT (IMED)) > XMED7) Then IMED = IMED - 1 Else If (XDATT (IWRKT (IMED)) < XMED7) Then IMED = IMED + 1 Endif Endif If (XDATT (IWRKT (IMED)) > XMED7) Then NLEQ = NLEQ + IMED - IDON IENDT (IDON1) = IMED - 1 ISTRT (IDON1) = IMED Else If (XDATT (IWRKT (IMED)) < XMED7) Then NLEQ = NLEQ + IMED - IDON + 1 IENDT (IDON1) = IMED ISTRT (IDON1) = IMED + 1 Else ! If (XDATT (IWRKT (IMED)) == XMED7) NLEQ = NLEQ + IMED - IDON + 1 NEQU = NEQU + 1 IENDT (IDON1) = IMED - 1 Do IMED1 = IMED - 1, IDON, -1 If (XDATT (IWRKT (IMED1)) == XMED7) Then NEQU = NEQU + 1 IENDT (IDON1) = IMED1 - 1 Else Exit End If End Do ISTRT (IDON1) = IMED + 1 Do IMED1 = IMED + 1, IDON + 6 If (XDATT (IWRKT (IMED1)) == XMED7) Then NEQU = NEQU + 1 NLEQ = NLEQ + 1 ISTRT (IDON1) = IMED1 + 1 Else Exit End If End Do Endif IDON1 = IDON1 + 1 End Do ! ! Carry out a partial insertion sort to find the Kth smallest of the ! large values, or the Kth largest of the small values, according to ! what is needed. ! ! If (NLEQ - NEQU + 1 <= NMED) Then If (NLEQ < NMED) Then ! Not enough low values IWRK1 = IMAX XWRK1 = XDATT (IWRK1) NORD = NMED - NLEQ IDON1 = 0 ICRS1 = 1 ICRS2 = 0 IDCR = 0 Do IDON = 1, NTRI, 7 IDON1 = IDON1 + 1 If (ICRS2 < NORD) Then Do ICRS = ISTRT (IDON1), IDON + 6 If (XDATT (IWRKT (ICRS)) < XWRK1) Then IWRK = IWRKT (ICRS) XWRK = XDATT (IWRK) Do IDCR = ICRS1 - 1, 1, - 1 If (XWRK >= XDATT (IWRKT (IDCR))) Exit IWRKT (IDCR+1) = IWRKT (IDCR) End Do IWRKT (IDCR+1) = IWRK IWRK1 = IWRKT (ICRS1) XWRK1 = XDATT (IWRK1) Else If (ICRS2 < NORD) Then IWRKT (ICRS1) = IWRKT (ICRS) IWRK1 = IWRKT (ICRS1) XWRK1 = XDATT (IWRK1) Endif End If ICRS1 = MIN (NORD, ICRS1 + 1) ICRS2 = MIN (NORD, ICRS2 + 1) End Do Else Do ICRS = ISTRT (IDON1), IDON + 6 If (XDATT (IWRKT (ICRS)) >= XWRK1) Exit IWRK = IWRKT (ICRS) XWRK = XDATT (IWRK) Do IDCR = ICRS1 - 1, 1, - 1 If (XWRK >= XDATT (IWRKT (IDCR))) Exit IWRKT (IDCR+1) = IWRKT (IDCR) End Do IWRKT (IDCR+1) = IWRK IWRK1 = IWRKT (ICRS1) XWRK1 = XDATT (IWRK1) End Do End If End Do ires_med = IWRK1 Return Else ires_med = IMED7 Return End If Else ! If (NLEQ > NMED) ! Not enough high values XWRK1 = -XHUGE NORD = NLEQ - NEQU - NMED + 1 IDON1 = 0 ICRS1 = 1 ICRS2 = 0 Do IDON = 1, NTRI, 7 IDON1 = IDON1 + 1 If (ICRS2 < NORD) Then ! Do ICRS = IDON, IENDT (IDON1) If (XDATT(IWRKT (ICRS)) > XWRK1) Then IWRK = IWRKT (ICRS) XWRK = XDATT (IWRK) IDCR = ICRS1 - 1 Do IDCR = ICRS1 - 1, 1, - 1 If (XWRK <= XDATT(IWRKT (IDCR))) Exit IWRKT (IDCR+1) = IWRKT (IDCR) End Do IWRKT (IDCR+1) = IWRK IWRK1 = IWRKT(ICRS1) XWRK1 = XDATT(IWRK1) Else If (ICRS2 < NORD) Then IWRKT (ICRS1) = IWRKT (ICRS) IWRK1 = IWRKT(ICRS1) XWRK1 = XDATT(IWRK1) End If End If ICRS1 = MIN (NORD, ICRS1 + 1) ICRS2 = MIN (NORD, ICRS2 + 1) End Do Else Do ICRS = IENDT (IDON1), IDON, -1 If (XDATT(IWRKT (ICRS)) <= XWRK1) Exit IWRK = IWRKT (ICRS) XWRK = XDATT (IWRK) IDCR = ICRS1 - 1 Do IDCR = ICRS1 - 1, 1, - 1 If (XWRK <= XDATT(IWRKT (IDCR))) Exit IWRKT (IDCR+1) = IWRKT (IDCR) End Do IWRKT (IDCR+1) = IWRK IWRK1 = IWRKT(ICRS1) XWRK1 = XDATT(IWRK1) End Do Endif End Do ! ires_med = IWRK1 Return End If ! END Subroutine d_med ! Subroutine R_indmed (XDONT, INDM) ! Returns index of median value of XDONT. ! __________________________________________________________ Real, Dimension (:), Intent (In) :: XDONT Integer, Intent (Out) :: INDM ! __________________________________________________________ Integer :: IDON ! Allocate (IDONT (SIZE(XDONT))) Do IDON = 1, SIZE(XDONT) IDONT (IDON) = IDON End Do ! Call r_med (XDONT, IDONT, INDM) ! Deallocate (IDONT) End Subroutine R_indmed Recursive Subroutine r_med (XDATT, IDATT, ires_med) ! Finds the index of the median of XDONT using the recursive procedure ! described in Knuth, The Art of Computer Programming, ! vol. 3, 5.3.3 - This procedure is linear in time, and ! does not require to be able to interpolate in the ! set as the one used in INDNTH. It also has better worst ! case behavior than INDNTH, but is about 30% slower in ! average for random uniformly distributed values. ! __________________________________________________________ Real, Dimension (:), Intent (In) :: XDATT Integer, Dimension (:), Intent (In) :: IDATT Integer, Intent (Out) :: ires_med ! __________________________________________________________ ! Real, Parameter :: XHUGE = HUGE (XDATT) Real :: XWRK, XWRK1, XMED7, XMAX, XMIN ! Integer, Dimension (7*(((Size (IDATT)+6)/7+6)/7)) :: ISTRT, IENDT, IMEDT Integer, Dimension (7*((Size(IDATT)+6)/7)) :: IWRKT Integer :: NTRI, NMED, NORD, NEQU, NLEQ, IMED, IDON, IDON1 Integer :: IDEB, ITMP, IDCR, ICRS, ICRS1, ICRS2, IMAX, IMIN Integer :: IWRK, IWRK1, IMED1, IMED7, NDAT ! NDAT = Size (IDATT) NMED = (NDAT+1) / 2 IWRKT = IDATT ! ! If the number of values is small, then use insertion sort ! If (NDAT < 35) Then ! ! Bring minimum to first location to save test in decreasing loop ! IDCR = NDAT If (XDATT (IWRKT (1)) < XDATT (IWRKT (IDCR))) Then IWRK = IWRKT (1) Else IWRK = IWRKT (IDCR) IWRKT (IDCR) = IWRKT (1) Endif XWRK = XDATT (IWRK) Do ITMP = 1, NDAT - 2 IDCR = IDCR - 1 IWRK1 = IWRKT (IDCR) XWRK1 = XDATT (IWRK1) If (XWRK1 < XWRK) Then IWRKT (IDCR) = IWRK XWRK = XWRK1 IWRK = IWRK1 Endif End Do IWRKT (1) = IWRK ! ! Sort the first half, until we have NMED sorted values ! Do ICRS = 3, NMED XWRK = XDATT (IWRKT (ICRS)) IWRK = IWRKT (ICRS) IDCR = ICRS - 1 Do If (XWRK >= XDATT (IWRKT(IDCR))) Exit IWRKT (IDCR+1) = IWRKT (IDCR) IDCR = IDCR - 1 End Do IWRKT (IDCR+1) = IWRK End Do ! ! Insert any value less than the current median in the first half ! XWRK1 = XDATT (IWRKT (NMED)) Do ICRS = NMED+1, NDAT XWRK = XDATT (IWRKT (ICRS)) IWRK = IWRKT (ICRS) If (XWRK < XWRK1) Then IDCR = NMED - 1 Do If (XWRK >= XDATT (IWRKT(IDCR))) Exit IWRKT (IDCR+1) = IWRKT (IDCR) IDCR = IDCR - 1 End Do IWRKT (IDCR+1) = IWRK XWRK1 = XDATT (IWRKT (NMED)) End If End Do ires_med = IWRKT (NMED) Return End If ! ! Make sorted subsets of 7 elements ! This is done by a variant of insertion sort where a first ! pass is used to bring the smallest element to the first position ! decreasing disorder at the same time, so that we may remove ! remove the loop test in the insertion loop. ! IMAX = 1 IMIN = 1 XMAX = XDATT (IWRKT(IMAX)) XMIN = XDATT (IWRKT(IMIN)) DO IDEB = 1, NDAT-6, 7 IDCR = IDEB + 6 If (XDATT (IWRKT(IDEB)) < XDATT (IWRKT(IDCR))) Then IWRK = IWRKT(IDEB) Else IWRK = IWRKT (IDCR) IWRKT (IDCR) = IWRKT(IDEB) Endif XWRK = XDATT (IWRK) Do ITMP = 1, 5 IDCR = IDCR - 1 IWRK1 = IWRKT (IDCR) XWRK1 = XDATT (IWRK1) If (XWRK1 < XWRK) Then IWRKT (IDCR) = IWRK IWRK = IWRK1 XWRK = XWRK1 Endif End Do IWRKT (IDEB) = IWRK If (XWRK < XMIN) Then IMIN = IWRK XMIN = XWRK End If Do ICRS = IDEB+1, IDEB+5 IWRK = IWRKT (ICRS+1) XWRK = XDATT (IWRK) IDON = IWRKT(ICRS) If (XWRK < XDATT(IDON)) Then IWRKT (ICRS+1) = IDON IDCR = ICRS IWRK1 = IWRKT (IDCR-1) XWRK1 = XDATT (IWRK1) Do If (XWRK >= XWRK1) Exit IWRKT (IDCR) = IWRK1 IDCR = IDCR - 1 IWRK1 = IWRKT (IDCR-1) XWRK1 = XDATT (IWRK1) End Do IWRKT (IDCR) = IWRK EndIf End Do If (XWRK > XMAX) Then IMAX = IWRK XMAX = XWRK End If End Do ! ! Add-up alternatively MAX and MIN values to make the number of data ! an exact multiple of 7. ! IDEB = 7 * (NDAT/7) NTRI = NDAT If (IDEB < NDAT) Then ! Do ICRS = IDEB+1, NDAT XWRK1 = XDATT (IWRKT (ICRS)) IF (XWRK1 > XMAX) Then IMAX = IWRKT (ICRS) XMAX = XWRK1 End If IF (XWRK1 < XMIN) Then IMIN = IWRKT (ICRS) XMIN = XWRK1 End If End Do IWRK1 = IMAX Do ICRS = NDAT+1, IDEB+7 IWRKT (ICRS) = IWRK1 If (IWRK1 == IMAX) Then IWRK1 = IMIN Else NMED = NMED + 1 IWRK1 = IMAX End If End Do ! Do ICRS = IDEB+2, IDEB+7 IWRK = IWRKT (ICRS) XWRK = XDATT (IWRK) Do IDCR = ICRS - 1, IDEB+1, - 1 If (XWRK >= XDATT (IWRKT(IDCR))) Exit IWRKT (IDCR+1) = IWRKT (IDCR) End Do IWRKT (IDCR+1) = IWRK End Do ! NTRI = IDEB+7 End If ! ! Make the set of the indices of median values of each sorted subset ! IDON1 = 0 Do IDON = 1, NTRI, 7 IDON1 = IDON1 + 1 IMEDT (IDON1) = IWRKT (IDON + 3) End Do ! ! Find XMED7, the median of the medians ! Call r_med (XDATT, IMEDT(1:IDON1), IMED7) XMED7 = XDATT (IMED7) ! ! Count how many values are not higher than (and how many equal to) XMED7 ! This number is at least 4 * 1/2 * (N/7) : 4 values in each of the ! subsets where the median is lower than the median of medians. For similar ! reasons, we also have at least 2N/7 values not lower than XMED7. At the ! same time, we find in each subset the index of the last value < XMED7, ! and that of the first > XMED7. These indices will be used to restrict the ! search for the median as the Kth element in the subset (> or <) where ! we know it to be. ! IDON1 = 1 NLEQ = 0 NEQU = 0 Do IDON = 1, NTRI, 7 IMED = IDON+3 If (XDATT (IWRKT (IMED)) > XMED7) Then IMED = IMED - 2 If (XDATT (IWRKT (IMED)) > XMED7) Then IMED = IMED - 1 Else If (XDATT (IWRKT (IMED)) < XMED7) Then IMED = IMED + 1 Endif Else If (XDATT (IWRKT (IMED)) < XMED7) Then IMED = IMED + 2 If (XDATT (IWRKT (IMED)) > XMED7) Then IMED = IMED - 1 Else If (XDATT (IWRKT (IMED)) < XMED7) Then IMED = IMED + 1 Endif Endif If (XDATT (IWRKT (IMED)) > XMED7) Then NLEQ = NLEQ + IMED - IDON IENDT (IDON1) = IMED - 1 ISTRT (IDON1) = IMED Else If (XDATT (IWRKT (IMED)) < XMED7) Then NLEQ = NLEQ + IMED - IDON + 1 IENDT (IDON1) = IMED ISTRT (IDON1) = IMED + 1 Else ! If (XDATT (IWRKT (IMED)) == XMED7) NLEQ = NLEQ + IMED - IDON + 1 NEQU = NEQU + 1 IENDT (IDON1) = IMED - 1 Do IMED1 = IMED - 1, IDON, -1 If (XDATT (IWRKT (IMED1)) == XMED7) Then NEQU = NEQU + 1 IENDT (IDON1) = IMED1 - 1 Else Exit End If End Do ISTRT (IDON1) = IMED + 1 Do IMED1 = IMED + 1, IDON + 6 If (XDATT (IWRKT (IMED1)) == XMED7) Then NEQU = NEQU + 1 NLEQ = NLEQ + 1 ISTRT (IDON1) = IMED1 + 1 Else Exit End If End Do Endif IDON1 = IDON1 + 1 End Do ! ! Carry out a partial insertion sort to find the Kth smallest of the ! large values, or the Kth largest of the small values, according to ! what is needed. ! ! If (NLEQ - NEQU + 1 <= NMED) Then If (NLEQ < NMED) Then ! Not enough low values IWRK1 = IMAX XWRK1 = XDATT (IWRK1) NORD = NMED - NLEQ IDON1 = 0 ICRS1 = 1 ICRS2 = 0 IDCR = 0 Do IDON = 1, NTRI, 7 IDON1 = IDON1 + 1 If (ICRS2 < NORD) Then Do ICRS = ISTRT (IDON1), IDON + 6 If (XDATT (IWRKT (ICRS)) < XWRK1) Then IWRK = IWRKT (ICRS) XWRK = XDATT (IWRK) Do IDCR = ICRS1 - 1, 1, - 1 If (XWRK >= XDATT (IWRKT (IDCR))) Exit IWRKT (IDCR+1) = IWRKT (IDCR) End Do IWRKT (IDCR+1) = IWRK IWRK1 = IWRKT (ICRS1) XWRK1 = XDATT (IWRK1) Else If (ICRS2 < NORD) Then IWRKT (ICRS1) = IWRKT (ICRS) IWRK1 = IWRKT (ICRS1) XWRK1 = XDATT (IWRK1) Endif End If ICRS1 = MIN (NORD, ICRS1 + 1) ICRS2 = MIN (NORD, ICRS2 + 1) End Do Else Do ICRS = ISTRT (IDON1), IDON + 6 If (XDATT (IWRKT (ICRS)) >= XWRK1) Exit IWRK = IWRKT (ICRS) XWRK = XDATT (IWRK) Do IDCR = ICRS1 - 1, 1, - 1 If (XWRK >= XDATT (IWRKT (IDCR))) Exit IWRKT (IDCR+1) = IWRKT (IDCR) End Do IWRKT (IDCR+1) = IWRK IWRK1 = IWRKT (ICRS1) XWRK1 = XDATT (IWRK1) End Do End If End Do ires_med = IWRK1 Return Else ires_med = IMED7 Return End If Else ! If (NLEQ > NMED) ! Not enough high values XWRK1 = -XHUGE NORD = NLEQ - NEQU - NMED + 1 IDON1 = 0 ICRS1 = 1 ICRS2 = 0 Do IDON = 1, NTRI, 7 IDON1 = IDON1 + 1 If (ICRS2 < NORD) Then ! Do ICRS = IDON, IENDT (IDON1) If (XDATT(IWRKT (ICRS)) > XWRK1) Then IWRK = IWRKT (ICRS) XWRK = XDATT (IWRK) IDCR = ICRS1 - 1 Do IDCR = ICRS1 - 1, 1, - 1 If (XWRK <= XDATT(IWRKT (IDCR))) Exit IWRKT (IDCR+1) = IWRKT (IDCR) End Do IWRKT (IDCR+1) = IWRK IWRK1 = IWRKT(ICRS1) XWRK1 = XDATT(IWRK1) Else If (ICRS2 < NORD) Then IWRKT (ICRS1) = IWRKT (ICRS) IWRK1 = IWRKT(ICRS1) XWRK1 = XDATT(IWRK1) End If End If ICRS1 = MIN (NORD, ICRS1 + 1) ICRS2 = MIN (NORD, ICRS2 + 1) End Do Else Do ICRS = IENDT (IDON1), IDON, -1 If (XDATT(IWRKT (ICRS)) <= XWRK1) Exit IWRK = IWRKT (ICRS) XWRK = XDATT (IWRK) IDCR = ICRS1 - 1 Do IDCR = ICRS1 - 1, 1, - 1 If (XWRK <= XDATT(IWRKT (IDCR))) Exit IWRKT (IDCR+1) = IWRKT (IDCR) End Do IWRKT (IDCR+1) = IWRK IWRK1 = IWRKT(ICRS1) XWRK1 = XDATT(IWRK1) End Do Endif End Do ! ires_med = IWRK1 Return End If ! END Subroutine r_med Subroutine I_indmed (XDONT, INDM) ! Returns index of median value of XDONT. ! __________________________________________________________ Integer, Dimension (:), Intent (In) :: XDONT Integer, Intent (Out) :: INDM ! __________________________________________________________ Integer :: IDON ! Allocate (IDONT (SIZE(XDONT))) Do IDON = 1, SIZE(XDONT) IDONT (IDON) = IDON End Do ! Call i_med (XDONT, IDONT, INDM) ! Deallocate (IDONT) End Subroutine I_indmed Recursive Subroutine i_med (XDATT, IDATT, ires_med) ! Finds the index of the median of XDONT using the recursive procedure ! described in Knuth, The Art of Computer Programming, ! vol. 3, 5.3.3 - This procedure is linear in time, and ! does not require to be able to interpolate in the ! set as the one used in INDNTH. It also has better worst ! case behavior than INDNTH, but is about 30% slower in ! average for random uniformly distributed values. ! __________________________________________________________ Integer, Dimension (:), Intent (In) :: XDATT Integer, Dimension (:), Intent (In) :: IDATT Integer, Intent (Out) :: ires_med ! __________________________________________________________ ! Integer, Parameter :: XHUGE = HUGE (XDATT) Integer :: XWRK, XWRK1, XMED7, XMAX, XMIN ! Integer, Dimension (7*(((Size (IDATT)+6)/7+6)/7)) :: ISTRT, IENDT, IMEDT Integer, Dimension (7*((Size(IDATT)+6)/7)) :: IWRKT Integer :: NTRI, NMED, NORD, NEQU, NLEQ, IMED, IDON, IDON1 Integer :: IDEB, ITMP, IDCR, ICRS, ICRS1, ICRS2, IMAX, IMIN Integer :: IWRK, IWRK1, IMED1, IMED7, NDAT ! NDAT = Size (IDATT) NMED = (NDAT+1) / 2 IWRKT = IDATT ! ! If the number of values is small, then use insertion sort ! If (NDAT < 35) Then ! ! Bring minimum to first location to save test in decreasing loop ! IDCR = NDAT If (XDATT (IWRKT (1)) < XDATT (IWRKT (IDCR))) Then IWRK = IWRKT (1) Else IWRK = IWRKT (IDCR) IWRKT (IDCR) = IWRKT (1) Endif XWRK = XDATT (IWRK) Do ITMP = 1, NDAT - 2 IDCR = IDCR - 1 IWRK1 = IWRKT (IDCR) XWRK1 = XDATT (IWRK1) If (XWRK1 < XWRK) Then IWRKT (IDCR) = IWRK XWRK = XWRK1 IWRK = IWRK1 Endif End Do IWRKT (1) = IWRK ! ! Sort the first half, until we have NMED sorted values ! Do ICRS = 3, NMED XWRK = XDATT (IWRKT (ICRS)) IWRK = IWRKT (ICRS) IDCR = ICRS - 1 Do If (XWRK >= XDATT (IWRKT(IDCR))) Exit IWRKT (IDCR+1) = IWRKT (IDCR) IDCR = IDCR - 1 End Do IWRKT (IDCR+1) = IWRK End Do ! ! Insert any value less than the current median in the first half ! XWRK1 = XDATT (IWRKT (NMED)) Do ICRS = NMED+1, NDAT XWRK = XDATT (IWRKT (ICRS)) IWRK = IWRKT (ICRS) If (XWRK < XWRK1) Then IDCR = NMED - 1 Do If (XWRK >= XDATT (IWRKT(IDCR))) Exit IWRKT (IDCR+1) = IWRKT (IDCR) IDCR = IDCR - 1 End Do IWRKT (IDCR+1) = IWRK XWRK1 = XDATT (IWRKT (NMED)) End If End Do ires_med = IWRKT (NMED) Return End If ! ! Make sorted subsets of 7 elements ! This is done by a variant of insertion sort where a first ! pass is used to bring the smallest element to the first position ! decreasing disorder at the same time, so that we may remove ! remove the loop test in the insertion loop. ! IMAX = 1 IMIN = 1 XMAX = XDATT (IWRKT(IMAX)) XMIN = XDATT (IWRKT(IMIN)) DO IDEB = 1, NDAT-6, 7 IDCR = IDEB + 6 If (XDATT (IWRKT(IDEB)) < XDATT (IWRKT(IDCR))) Then IWRK = IWRKT(IDEB) Else IWRK = IWRKT (IDCR) IWRKT (IDCR) = IWRKT(IDEB) Endif XWRK = XDATT (IWRK) Do ITMP = 1, 5 IDCR = IDCR - 1 IWRK1 = IWRKT (IDCR) XWRK1 = XDATT (IWRK1) If (XWRK1 < XWRK) Then IWRKT (IDCR) = IWRK IWRK = IWRK1 XWRK = XWRK1 Endif End Do IWRKT (IDEB) = IWRK If (XWRK < XMIN) Then IMIN = IWRK XMIN = XWRK End If Do ICRS = IDEB+1, IDEB+5 IWRK = IWRKT (ICRS+1) XWRK = XDATT (IWRK) IDON = IWRKT(ICRS) If (XWRK < XDATT(IDON)) Then IWRKT (ICRS+1) = IDON IDCR = ICRS IWRK1 = IWRKT (IDCR-1) XWRK1 = XDATT (IWRK1) Do If (XWRK >= XWRK1) Exit IWRKT (IDCR) = IWRK1 IDCR = IDCR - 1 IWRK1 = IWRKT (IDCR-1) XWRK1 = XDATT (IWRK1) End Do IWRKT (IDCR) = IWRK EndIf End Do If (XWRK > XMAX) Then IMAX = IWRK XMAX = XWRK End If End Do ! ! Add-up alternatively MAX and MIN values to make the number of data ! an exact multiple of 7. ! IDEB = 7 * (NDAT/7) NTRI = NDAT If (IDEB < NDAT) Then ! Do ICRS = IDEB+1, NDAT XWRK1 = XDATT (IWRKT (ICRS)) IF (XWRK1 > XMAX) Then IMAX = IWRKT (ICRS) XMAX = XWRK1 End If IF (XWRK1 < XMIN) Then IMIN = IWRKT (ICRS) XMIN = XWRK1 End If End Do IWRK1 = IMAX Do ICRS = NDAT+1, IDEB+7 IWRKT (ICRS) = IWRK1 If (IWRK1 == IMAX) Then IWRK1 = IMIN Else NMED = NMED + 1 IWRK1 = IMAX End If End Do ! Do ICRS = IDEB+2, IDEB+7 IWRK = IWRKT (ICRS) XWRK = XDATT (IWRK) Do IDCR = ICRS - 1, IDEB+1, - 1 If (XWRK >= XDATT (IWRKT(IDCR))) Exit IWRKT (IDCR+1) = IWRKT (IDCR) End Do IWRKT (IDCR+1) = IWRK End Do ! NTRI = IDEB+7 End If ! ! Make the set of the indices of median values of each sorted subset ! IDON1 = 0 Do IDON = 1, NTRI, 7 IDON1 = IDON1 + 1 IMEDT (IDON1) = IWRKT (IDON + 3) End Do ! ! Find XMED7, the median of the medians ! Call i_med (XDATT, IMEDT(1:IDON1), IMED7) XMED7 = XDATT (IMED7) ! ! Count how many values are not higher than (and how many equal to) XMED7 ! This number is at least 4 * 1/2 * (N/7) : 4 values in each of the ! subsets where the median is lower than the median of medians. For similar ! reasons, we also have at least 2N/7 values not lower than XMED7. At the ! same time, we find in each subset the index of the last value < XMED7, ! and that of the first > XMED7. These indices will be used to restrict the ! search for the median as the Kth element in the subset (> or <) where ! we know it to be. ! IDON1 = 1 NLEQ = 0 NEQU = 0 Do IDON = 1, NTRI, 7 IMED = IDON+3 If (XDATT (IWRKT (IMED)) > XMED7) Then IMED = IMED - 2 If (XDATT (IWRKT (IMED)) > XMED7) Then IMED = IMED - 1 Else If (XDATT (IWRKT (IMED)) < XMED7) Then IMED = IMED + 1 Endif Else If (XDATT (IWRKT (IMED)) < XMED7) Then IMED = IMED + 2 If (XDATT (IWRKT (IMED)) > XMED7) Then IMED = IMED - 1 Else If (XDATT (IWRKT (IMED)) < XMED7) Then IMED = IMED + 1 Endif Endif If (XDATT (IWRKT (IMED)) > XMED7) Then NLEQ = NLEQ + IMED - IDON IENDT (IDON1) = IMED - 1 ISTRT (IDON1) = IMED Else If (XDATT (IWRKT (IMED)) < XMED7) Then NLEQ = NLEQ + IMED - IDON + 1 IENDT (IDON1) = IMED ISTRT (IDON1) = IMED + 1 Else ! If (XDATT (IWRKT (IMED)) == XMED7) NLEQ = NLEQ + IMED - IDON + 1 NEQU = NEQU + 1 IENDT (IDON1) = IMED - 1 Do IMED1 = IMED - 1, IDON, -1 If (XDATT (IWRKT (IMED1)) == XMED7) Then NEQU = NEQU + 1 IENDT (IDON1) = IMED1 - 1 Else Exit End If End Do ISTRT (IDON1) = IMED + 1 Do IMED1 = IMED + 1, IDON + 6 If (XDATT (IWRKT (IMED1)) == XMED7) Then NEQU = NEQU + 1 NLEQ = NLEQ + 1 ISTRT (IDON1) = IMED1 + 1 Else Exit End If End Do Endif IDON1 = IDON1 + 1 End Do ! ! Carry out a partial insertion sort to find the Kth smallest of the ! large values, or the Kth largest of the small values, according to ! what is needed. ! ! If (NLEQ - NEQU + 1 <= NMED) Then If (NLEQ < NMED) Then ! Not enough low values IWRK1 = IMAX XWRK1 = XDATT (IWRK1) NORD = NMED - NLEQ IDON1 = 0 ICRS1 = 1 ICRS2 = 0 IDCR = 0 Do IDON = 1, NTRI, 7 IDON1 = IDON1 + 1 If (ICRS2 < NORD) Then Do ICRS = ISTRT (IDON1), IDON + 6 If (XDATT (IWRKT (ICRS)) < XWRK1) Then IWRK = IWRKT (ICRS) XWRK = XDATT (IWRK) Do IDCR = ICRS1 - 1, 1, - 1 If (XWRK >= XDATT (IWRKT (IDCR))) Exit IWRKT (IDCR+1) = IWRKT (IDCR) End Do IWRKT (IDCR+1) = IWRK IWRK1 = IWRKT (ICRS1) XWRK1 = XDATT (IWRK1) Else If (ICRS2 < NORD) Then IWRKT (ICRS1) = IWRKT (ICRS) IWRK1 = IWRKT (ICRS1) XWRK1 = XDATT (IWRK1) Endif End If ICRS1 = MIN (NORD, ICRS1 + 1) ICRS2 = MIN (NORD, ICRS2 + 1) End Do Else Do ICRS = ISTRT (IDON1), IDON + 6 If (XDATT (IWRKT (ICRS)) >= XWRK1) Exit IWRK = IWRKT (ICRS) XWRK = XDATT (IWRK) Do IDCR = ICRS1 - 1, 1, - 1 If (XWRK >= XDATT (IWRKT (IDCR))) Exit IWRKT (IDCR+1) = IWRKT (IDCR) End Do IWRKT (IDCR+1) = IWRK IWRK1 = IWRKT (ICRS1) XWRK1 = XDATT (IWRK1) End Do End If End Do ires_med = IWRK1 Return Else ires_med = IMED7 Return End If Else ! If (NLEQ > NMED) ! Not enough high values XWRK1 = -XHUGE NORD = NLEQ - NEQU - NMED + 1 IDON1 = 0 ICRS1 = 1 ICRS2 = 0 Do IDON = 1, NTRI, 7 IDON1 = IDON1 + 1 If (ICRS2 < NORD) Then ! Do ICRS = IDON, IENDT (IDON1) If (XDATT(IWRKT (ICRS)) > XWRK1) Then IWRK = IWRKT (ICRS) XWRK = XDATT (IWRK) IDCR = ICRS1 - 1 Do IDCR = ICRS1 - 1, 1, - 1 If (XWRK <= XDATT(IWRKT (IDCR))) Exit IWRKT (IDCR+1) = IWRKT (IDCR) End Do IWRKT (IDCR+1) = IWRK IWRK1 = IWRKT(ICRS1) XWRK1 = XDATT(IWRK1) Else If (ICRS2 < NORD) Then IWRKT (ICRS1) = IWRKT (ICRS) IWRK1 = IWRKT(ICRS1) XWRK1 = XDATT(IWRK1) End If End If ICRS1 = MIN (NORD, ICRS1 + 1) ICRS2 = MIN (NORD, ICRS2 + 1) End Do Else Do ICRS = IENDT (IDON1), IDON, -1 If (XDATT(IWRKT (ICRS)) <= XWRK1) Exit IWRK = IWRKT (ICRS) XWRK = XDATT (IWRK) IDCR = ICRS1 - 1 Do IDCR = ICRS1 - 1, 1, - 1 If (XWRK <= XDATT(IWRKT (IDCR))) Exit IWRKT (IDCR+1) = IWRKT (IDCR) End Do IWRKT (IDCR+1) = IWRK IWRK1 = IWRKT(ICRS1) XWRK1 = XDATT(IWRK1) End Do Endif End Do ! ires_med = IWRK1 Return End If ! END Subroutine i_med end module m_indmed