Module m_valmed Integer, Parameter :: kdp = selected_real_kind(15) public :: valmed private :: kdp private :: R_valmed, I_valmed, D_valmed interface valmed module procedure d_valmed, r_valmed, i_valmed end interface valmed contains Recursive Function D_valmed (XDONT) Result (res_med) ! Finds 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) :: XDONT Real (kind=kdp) :: res_med ! __________________________________________________________ Real (kind=kdp), Parameter :: XHUGE = HUGE (XDONT) Real (kind=kdp), Dimension (SIZE(XDONT)+6) :: XWRKT Real (kind=kdp) :: XWRK, XWRK1, XMED7 ! Integer, Dimension ((SIZE(XDONT)+6)/7) :: ISTRT, IENDT, IMEDT Integer :: NDON, NTRI, NMED, NORD, NEQU, NLEQ, IMED, IDON, IDON1 Integer :: IDEB, IWRK, IDCR, ICRS, ICRS1, ICRS2, IMED1 ! NDON = SIZE (XDONT) NMED = (NDON+1) / 2 ! write(unit=*,fmt=*) NMED, NDON ! ! If the number of values is small, then use insertion sort ! If (NDON < 35) Then ! ! Bring minimum to first location to save test in decreasing loop ! IDCR = NDON If (XDONT (1) < XDONT (NDON)) Then XWRK = XDONT (1) XWRKT (IDCR) = XDONT (IDCR) Else XWRK = XDONT (IDCR) XWRKT (IDCR) = XDONT (1) Endif Do IWRK = 1, NDON - 2 IDCR = IDCR - 1 XWRK1 = XDONT (IDCR) If (XWRK1 < XWRK) Then XWRKT (IDCR) = XWRK XWRK = XWRK1 Else XWRKT (IDCR) = XWRK1 Endif End Do XWRKT (1) = XWRK ! ! Sort the first half, until we have NMED sorted values ! Do ICRS = 3, NMED XWRK = XWRKT (ICRS) IDCR = ICRS - 1 Do If (XWRK >= XWRKT(IDCR)) Exit XWRKT (IDCR+1) = XWRKT (IDCR) IDCR = IDCR - 1 End Do XWRKT (IDCR+1) = XWRK End Do ! ! Insert any value less than the current median in the first half ! Do ICRS = NMED+1, NDON XWRK = XWRKT (ICRS) If (XWRK < XWRKT (NMED)) Then IDCR = NMED - 1 Do If (XWRK >= XWRKT(IDCR)) Exit XWRKT (IDCR+1) = XWRKT (IDCR) IDCR = IDCR - 1 End Do XWRKT (IDCR+1) = XWRK End If End Do res_med = XWRKT (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. ! DO IDEB = 1, NDON-6, 7 IDCR = IDEB + 6 If (XDONT (IDEB) < XDONT (IDCR)) Then XWRK = XDONT (IDEB) XWRKT (IDCR) = XDONT (IDCR) Else XWRK = XDONT (IDCR) XWRKT (IDCR) = XDONT (IDEB) Endif Do IWRK = 1, 5 IDCR = IDCR - 1 XWRK1 = XDONT (IDCR) If (XWRK1 < XWRK) Then XWRKT (IDCR) = XWRK XWRK = XWRK1 Else XWRKT (IDCR) = XWRK1 Endif End Do XWRKT (IDEB) = XWRK Do ICRS = IDEB+2, IDEB+6 XWRK = XWRKT (ICRS) If (XWRK < XWRKT(ICRS-1)) Then XWRKT (ICRS) = XWRKT (ICRS-1) IDCR = ICRS - 1 XWRK1 = XWRKT (IDCR-1) Do If (XWRK >= XWRK1) Exit XWRKT (IDCR) = XWRK1 IDCR = IDCR - 1 XWRK1 = XWRKT (IDCR-1) End Do XWRKT (IDCR) = XWRK EndIf End Do End Do ! ! Add-up alternatively + and - HUGE values to make the number of data ! an exact multiple of 7. ! IDEB = 7 * (NDON/7) NTRI = NDON If (IDEB < NDON) Then ! XWRK1 = XHUGE Do ICRS = IDEB+1, IDEB+7 If (ICRS <= NDON) Then XWRKT (ICRS) = XDONT (ICRS) Else If (XWRK1 /= XHUGE) NMED = NMED + 1 XWRKT (ICRS) = XWRK1 XWRK1 = - XWRK1 Endif End Do ! Do ICRS = IDEB+2, IDEB+7 XWRK = XWRKT (ICRS) Do IDCR = ICRS - 1, IDEB+1, - 1 If (XWRK >= XWRKT(IDCR)) Exit XWRKT (IDCR+1) = XWRKT (IDCR) End Do XWRKT (IDCR+1) = XWRK 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) = IDON + 3 End Do ! ! Find XMED7, the median of the medians ! XMED7 = D_valmed (XWRKT (IMEDT)) ! ! 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 (XWRKT (IMED) > XMED7) Then IMED = IMED - 2 If (XWRKT (IMED) > XMED7) Then IMED = IMED - 1 Else If (XWRKT (IMED) < XMED7) Then IMED = IMED + 1 Endif Else If (XWRKT (IMED) < XMED7) Then IMED = IMED + 2 If (XWRKT (IMED) > XMED7) Then IMED = IMED - 1 Else If (XWRKT (IMED) < XMED7) Then IMED = IMED + 1 Endif Endif If (XWRKT (IMED) > XMED7) Then NLEQ = NLEQ + IMED - IDON IENDT (IDON1) = IMED - 1 ISTRT (IDON1) = IMED Else If (XWRKT (IMED) < XMED7) Then NLEQ = NLEQ + IMED - IDON + 1 IENDT (IDON1) = IMED ISTRT (IDON1) = IMED + 1 Else ! If (XWRKT (IMED) == XMED7) NLEQ = NLEQ + IMED - IDON + 1 NEQU = NEQU + 1 IENDT (IDON1) = IMED - 1 Do IMED1 = IMED - 1, IDON, -1 If (XWRKT (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 (XWRKT (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 XWRK1 = XHUGE 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 (XWRKT(ICRS) < XWRK1) Then XWRK = XWRKT (ICRS) Do IDCR = ICRS1 - 1, 1, - 1 If (XWRK >= XWRKT(IDCR)) Exit XWRKT (IDCR+1) = XWRKT (IDCR) End Do XWRKT (IDCR+1) = XWRK XWRK1 = XWRKT(ICRS1) Else If (ICRS2 < NORD) Then XWRKT (ICRS1) = XWRKT (ICRS) XWRK1 = XWRKT(ICRS1) Endif End If ICRS1 = MIN (NORD, ICRS1 + 1) ICRS2 = MIN (NORD, ICRS2 + 1) End Do Else Do ICRS = ISTRT (IDON1), IDON + 6 If (XWRKT(ICRS) >= XWRK1) Exit XWRK = XWRKT (ICRS) Do IDCR = ICRS1 - 1, 1, - 1 If (XWRK >= XWRKT(IDCR)) Exit XWRKT (IDCR+1) = XWRKT (IDCR) End Do XWRKT (IDCR+1) = XWRK XWRK1 = XWRKT(ICRS1) End Do End If End Do res_med = XWRK1 Return Else res_med = XMED7 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 (XWRKT(ICRS) > XWRK1) Then XWRK = XWRKT (ICRS) IDCR = ICRS1 - 1 Do IDCR = ICRS1 - 1, 1, - 1 If (XWRK <= XWRKT(IDCR)) Exit XWRKT (IDCR+1) = XWRKT (IDCR) End Do XWRKT (IDCR+1) = XWRK XWRK1 = XWRKT(ICRS1) Else If (ICRS2 < NORD) Then XWRKT (ICRS1) = XWRKT (ICRS) XWRK1 = XWRKT (ICRS1) End If End If ICRS1 = MIN (NORD, ICRS1 + 1) ICRS2 = MIN (NORD, ICRS2 + 1) End Do Else Do ICRS = IENDT (IDON1), IDON, -1 If (XWRKT(ICRS) > XWRK1) Then XWRK = XWRKT (ICRS) IDCR = ICRS1 - 1 Do IDCR = ICRS1 - 1, 1, - 1 If (XWRK <= XWRKT(IDCR)) Exit XWRKT (IDCR+1) = XWRKT (IDCR) End Do XWRKT (IDCR+1) = XWRK XWRK1 = XWRKT(ICRS1) Else Exit End If End Do Endif End Do ! res_med = XWRK1 Return End If ! End Function D_valmed Recursive Function R_valmed (XDONT) Result (res_med) ! Finds 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) :: XDONT Real :: res_med ! __________________________________________________________ Real, Parameter :: XHUGE = HUGE (XDONT) Real, Dimension (SIZE(XDONT)+6) :: XWRKT Real :: XWRK, XWRK1, XMED7 ! Integer, Dimension ((SIZE(XDONT)+6)/7) :: ISTRT, IENDT, IMEDT Integer :: NDON, NTRI, NMED, NORD, NEQU, NLEQ, IMED, IDON, IDON1 Integer :: IDEB, IWRK, IDCR, ICRS, ICRS1, ICRS2, IMED1 ! NDON = SIZE (XDONT) NMED = (NDON+1) / 2 ! write(unit=*,fmt=*) NMED, NDON ! ! If the number of values is small, then use insertion sort ! If (NDON < 35) Then ! ! Bring minimum to first location to save test in decreasing loop ! IDCR = NDON If (XDONT (1) < XDONT (NDON)) Then XWRK = XDONT (1) XWRKT (IDCR) = XDONT (IDCR) Else XWRK = XDONT (IDCR) XWRKT (IDCR) = XDONT (1) Endif Do IWRK = 1, NDON - 2 IDCR = IDCR - 1 XWRK1 = XDONT (IDCR) If (XWRK1 < XWRK) Then XWRKT (IDCR) = XWRK XWRK = XWRK1 Else XWRKT (IDCR) = XWRK1 Endif End Do XWRKT (1) = XWRK ! ! Sort the first half, until we have NMED sorted values ! Do ICRS = 3, NMED XWRK = XWRKT (ICRS) IDCR = ICRS - 1 Do If (XWRK >= XWRKT(IDCR)) Exit XWRKT (IDCR+1) = XWRKT (IDCR) IDCR = IDCR - 1 End Do XWRKT (IDCR+1) = XWRK End Do ! ! Insert any value less than the current median in the first half ! Do ICRS = NMED+1, NDON XWRK = XWRKT (ICRS) If (XWRK < XWRKT (NMED)) Then IDCR = NMED - 1 Do If (XWRK >= XWRKT(IDCR)) Exit XWRKT (IDCR+1) = XWRKT (IDCR) IDCR = IDCR - 1 End Do XWRKT (IDCR+1) = XWRK End If End Do res_med = XWRKT (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. ! DO IDEB = 1, NDON-6, 7 IDCR = IDEB + 6 If (XDONT (IDEB) < XDONT (IDCR)) Then XWRK = XDONT (IDEB) XWRKT (IDCR) = XDONT (IDCR) Else XWRK = XDONT (IDCR) XWRKT (IDCR) = XDONT (IDEB) Endif Do IWRK = 1, 5 IDCR = IDCR - 1 XWRK1 = XDONT (IDCR) If (XWRK1 < XWRK) Then XWRKT (IDCR) = XWRK XWRK = XWRK1 Else XWRKT (IDCR) = XWRK1 Endif End Do XWRKT (IDEB) = XWRK Do ICRS = IDEB+2, IDEB+6 XWRK = XWRKT (ICRS) If (XWRK < XWRKT(ICRS-1)) Then XWRKT (ICRS) = XWRKT (ICRS-1) IDCR = ICRS - 1 XWRK1 = XWRKT (IDCR-1) Do If (XWRK >= XWRK1) Exit XWRKT (IDCR) = XWRK1 IDCR = IDCR - 1 XWRK1 = XWRKT (IDCR-1) End Do XWRKT (IDCR) = XWRK EndIf End Do End Do ! ! Add-up alternatively + and - HUGE values to make the number of data ! an exact multiple of 7. ! IDEB = 7 * (NDON/7) NTRI = NDON If (IDEB < NDON) Then ! XWRK1 = XHUGE Do ICRS = IDEB+1, IDEB+7 If (ICRS <= NDON) Then XWRKT (ICRS) = XDONT (ICRS) Else If (XWRK1 /= XHUGE) NMED = NMED + 1 XWRKT (ICRS) = XWRK1 XWRK1 = - XWRK1 Endif End Do ! Do ICRS = IDEB+2, IDEB+7 XWRK = XWRKT (ICRS) Do IDCR = ICRS - 1, IDEB+1, - 1 If (XWRK >= XWRKT(IDCR)) Exit XWRKT (IDCR+1) = XWRKT (IDCR) End Do XWRKT (IDCR+1) = XWRK 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) = IDON + 3 End Do ! ! Find XMED7, the median of the medians ! XMED7 = R_valmed (XWRKT (IMEDT)) ! ! 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 (XWRKT (IMED) > XMED7) Then IMED = IMED - 2 If (XWRKT (IMED) > XMED7) Then IMED = IMED - 1 Else If (XWRKT (IMED) < XMED7) Then IMED = IMED + 1 Endif Else If (XWRKT (IMED) < XMED7) Then IMED = IMED + 2 If (XWRKT (IMED) > XMED7) Then IMED = IMED - 1 Else If (XWRKT (IMED) < XMED7) Then IMED = IMED + 1 Endif Endif If (XWRKT (IMED) > XMED7) Then NLEQ = NLEQ + IMED - IDON IENDT (IDON1) = IMED - 1 ISTRT (IDON1) = IMED Else If (XWRKT (IMED) < XMED7) Then NLEQ = NLEQ + IMED - IDON + 1 IENDT (IDON1) = IMED ISTRT (IDON1) = IMED + 1 Else ! If (XWRKT (IMED) == XMED7) NLEQ = NLEQ + IMED - IDON + 1 NEQU = NEQU + 1 IENDT (IDON1) = IMED - 1 Do IMED1 = IMED - 1, IDON, -1 If (XWRKT (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 (XWRKT (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 XWRK1 = XHUGE 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 (XWRKT(ICRS) < XWRK1) Then XWRK = XWRKT (ICRS) Do IDCR = ICRS1 - 1, 1, - 1 If (XWRK >= XWRKT(IDCR)) Exit XWRKT (IDCR+1) = XWRKT (IDCR) End Do XWRKT (IDCR+1) = XWRK XWRK1 = XWRKT(ICRS1) Else If (ICRS2 < NORD) Then XWRKT (ICRS1) = XWRKT (ICRS) XWRK1 = XWRKT(ICRS1) Endif End If ICRS1 = MIN (NORD, ICRS1 + 1) ICRS2 = MIN (NORD, ICRS2 + 1) End Do Else Do ICRS = ISTRT (IDON1), IDON + 6 If (XWRKT(ICRS) >= XWRK1) Exit XWRK = XWRKT (ICRS) Do IDCR = ICRS1 - 1, 1, - 1 If (XWRK >= XWRKT(IDCR)) Exit XWRKT (IDCR+1) = XWRKT (IDCR) End Do XWRKT (IDCR+1) = XWRK XWRK1 = XWRKT(ICRS1) End Do End If End Do res_med = XWRK1 Return Else res_med = XMED7 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 (XWRKT(ICRS) > XWRK1) Then XWRK = XWRKT (ICRS) IDCR = ICRS1 - 1 Do IDCR = ICRS1 - 1, 1, - 1 If (XWRK <= XWRKT(IDCR)) Exit XWRKT (IDCR+1) = XWRKT (IDCR) End Do XWRKT (IDCR+1) = XWRK XWRK1 = XWRKT(ICRS1) Else If (ICRS2 < NORD) Then XWRKT (ICRS1) = XWRKT (ICRS) XWRK1 = XWRKT (ICRS1) End If End If ICRS1 = MIN (NORD, ICRS1 + 1) ICRS2 = MIN (NORD, ICRS2 + 1) End Do Else Do ICRS = IENDT (IDON1), IDON, -1 If (XWRKT(ICRS) > XWRK1) Then XWRK = XWRKT (ICRS) IDCR = ICRS1 - 1 Do IDCR = ICRS1 - 1, 1, - 1 If (XWRK <= XWRKT(IDCR)) Exit XWRKT (IDCR+1) = XWRKT (IDCR) End Do XWRKT (IDCR+1) = XWRK XWRK1 = XWRKT(ICRS1) Else Exit End If End Do Endif End Do ! res_med = XWRK1 Return End If ! End Function R_valmed Recursive Function I_valmed (XDONT) Result (res_med) ! Finds 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) :: XDONT Integer :: res_med ! __________________________________________________________ Integer, Parameter :: XHUGE = HUGE (XDONT) Integer, Dimension (SIZE(XDONT)+6) :: XWRKT Integer :: XWRK, XWRK1, XMED7 ! Integer, Dimension ((SIZE(XDONT)+6)/7) :: ISTRT, IENDT, IMEDT Integer :: NDON, NTRI, NMED, NORD, NEQU, NLEQ, IMED, IDON, IDON1 Integer :: IDEB, IWRK, IDCR, ICRS, ICRS1, ICRS2, IMED1 ! NDON = SIZE (XDONT) NMED = (NDON+1) / 2 ! write(unit=*,fmt=*) NMED, NDON ! ! If the number of values is small, then use insertion sort ! If (NDON < 35) Then ! ! Bring minimum to first location to save test in decreasing loop ! IDCR = NDON If (XDONT (1) < XDONT (NDON)) Then XWRK = XDONT (1) XWRKT (IDCR) = XDONT (IDCR) Else XWRK = XDONT (IDCR) XWRKT (IDCR) = XDONT (1) Endif Do IWRK = 1, NDON - 2 IDCR = IDCR - 1 XWRK1 = XDONT (IDCR) If (XWRK1 < XWRK) Then XWRKT (IDCR) = XWRK XWRK = XWRK1 Else XWRKT (IDCR) = XWRK1 Endif End Do XWRKT (1) = XWRK ! ! Sort the first half, until we have NMED sorted values ! Do ICRS = 3, NMED XWRK = XWRKT (ICRS) IDCR = ICRS - 1 Do If (XWRK >= XWRKT(IDCR)) Exit XWRKT (IDCR+1) = XWRKT (IDCR) IDCR = IDCR - 1 End Do XWRKT (IDCR+1) = XWRK End Do ! ! Insert any value less than the current median in the first half ! Do ICRS = NMED+1, NDON XWRK = XWRKT (ICRS) If (XWRK < XWRKT (NMED)) Then IDCR = NMED - 1 Do If (XWRK >= XWRKT(IDCR)) Exit XWRKT (IDCR+1) = XWRKT (IDCR) IDCR = IDCR - 1 End Do XWRKT (IDCR+1) = XWRK End If End Do res_med = XWRKT (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. ! DO IDEB = 1, NDON-6, 7 IDCR = IDEB + 6 If (XDONT (IDEB) < XDONT (IDCR)) Then XWRK = XDONT (IDEB) XWRKT (IDCR) = XDONT (IDCR) Else XWRK = XDONT (IDCR) XWRKT (IDCR) = XDONT (IDEB) Endif Do IWRK = 1, 5 IDCR = IDCR - 1 XWRK1 = XDONT (IDCR) If (XWRK1 < XWRK) Then XWRKT (IDCR) = XWRK XWRK = XWRK1 Else XWRKT (IDCR) = XWRK1 Endif End Do XWRKT (IDEB) = XWRK Do ICRS = IDEB+2, IDEB+6 XWRK = XWRKT (ICRS) If (XWRK < XWRKT(ICRS-1)) Then XWRKT (ICRS) = XWRKT (ICRS-1) IDCR = ICRS - 1 XWRK1 = XWRKT (IDCR-1) Do If (XWRK >= XWRK1) Exit XWRKT (IDCR) = XWRK1 IDCR = IDCR - 1 XWRK1 = XWRKT (IDCR-1) End Do XWRKT (IDCR) = XWRK EndIf End Do End Do ! ! Add-up alternatively + and - HUGE values to make the number of data ! an exact multiple of 7. ! IDEB = 7 * (NDON/7) NTRI = NDON If (IDEB < NDON) Then ! XWRK1 = XHUGE Do ICRS = IDEB+1, IDEB+7 If (ICRS <= NDON) Then XWRKT (ICRS) = XDONT (ICRS) Else If (XWRK1 /= XHUGE) NMED = NMED + 1 XWRKT (ICRS) = XWRK1 XWRK1 = - XWRK1 Endif End Do ! Do ICRS = IDEB+2, IDEB+7 XWRK = XWRKT (ICRS) Do IDCR = ICRS - 1, IDEB+1, - 1 If (XWRK >= XWRKT(IDCR)) Exit XWRKT (IDCR+1) = XWRKT (IDCR) End Do XWRKT (IDCR+1) = XWRK 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) = IDON + 3 End Do ! ! Find XMED7, the median of the medians ! XMED7 = I_valmed (XWRKT (IMEDT)) ! ! 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 (XWRKT (IMED) > XMED7) Then IMED = IMED - 2 If (XWRKT (IMED) > XMED7) Then IMED = IMED - 1 Else If (XWRKT (IMED) < XMED7) Then IMED = IMED + 1 Endif Else If (XWRKT (IMED) < XMED7) Then IMED = IMED + 2 If (XWRKT (IMED) > XMED7) Then IMED = IMED - 1 Else If (XWRKT (IMED) < XMED7) Then IMED = IMED + 1 Endif Endif If (XWRKT (IMED) > XMED7) Then NLEQ = NLEQ + IMED - IDON IENDT (IDON1) = IMED - 1 ISTRT (IDON1) = IMED Else If (XWRKT (IMED) < XMED7) Then NLEQ = NLEQ + IMED - IDON + 1 IENDT (IDON1) = IMED ISTRT (IDON1) = IMED + 1 Else ! If (XWRKT (IMED) == XMED7) NLEQ = NLEQ + IMED - IDON + 1 NEQU = NEQU + 1 IENDT (IDON1) = IMED - 1 Do IMED1 = IMED - 1, IDON, -1 If (XWRKT (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 (XWRKT (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 XWRK1 = XHUGE 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 (XWRKT(ICRS) < XWRK1) Then XWRK = XWRKT (ICRS) Do IDCR = ICRS1 - 1, 1, - 1 If (XWRK >= XWRKT(IDCR)) Exit XWRKT (IDCR+1) = XWRKT (IDCR) End Do XWRKT (IDCR+1) = XWRK XWRK1 = XWRKT(ICRS1) Else If (ICRS2 < NORD) Then XWRKT (ICRS1) = XWRKT (ICRS) XWRK1 = XWRKT(ICRS1) Endif End If ICRS1 = MIN (NORD, ICRS1 + 1) ICRS2 = MIN (NORD, ICRS2 + 1) End Do Else Do ICRS = ISTRT (IDON1), IDON + 6 If (XWRKT(ICRS) >= XWRK1) Exit XWRK = XWRKT (ICRS) Do IDCR = ICRS1 - 1, 1, - 1 If (XWRK >= XWRKT(IDCR)) Exit XWRKT (IDCR+1) = XWRKT (IDCR) End Do XWRKT (IDCR+1) = XWRK XWRK1 = XWRKT(ICRS1) End Do End If End Do res_med = XWRK1 Return Else res_med = XMED7 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 (XWRKT(ICRS) > XWRK1) Then XWRK = XWRKT (ICRS) IDCR = ICRS1 - 1 Do IDCR = ICRS1 - 1, 1, - 1 If (XWRK <= XWRKT(IDCR)) Exit XWRKT (IDCR+1) = XWRKT (IDCR) End Do XWRKT (IDCR+1) = XWRK XWRK1 = XWRKT(ICRS1) Else If (ICRS2 < NORD) Then XWRKT (ICRS1) = XWRKT (ICRS) XWRK1 = XWRKT (ICRS1) End If End If ICRS1 = MIN (NORD, ICRS1 + 1) ICRS2 = MIN (NORD, ICRS2 + 1) End Do Else Do ICRS = IENDT (IDON1), IDON, -1 If (XWRKT(ICRS) > XWRK1) Then XWRK = XWRKT (ICRS) IDCR = ICRS1 - 1 Do IDCR = ICRS1 - 1, 1, - 1 If (XWRK <= XWRKT(IDCR)) Exit XWRKT (IDCR+1) = XWRKT (IDCR) End Do XWRKT (IDCR+1) = XWRK XWRK1 = XWRKT(ICRS1) Else Exit End If End Do Endif End Do ! res_med = XWRK1 Return End If ! End Function I_valmed end module m_valmed