PROGRAM t_median
IMPLICIT NONE
INTEGER :: n
REAL, ALLOCATABLE :: x(:)
REAL :: xmed
INTERFACE
SUBROUTINE median(x, n, xmed)
IMPLICIT NONE
INTEGER, INTENT(IN) :: n
REAL, INTENT(IN OUT), DIMENSION(:) :: x
REAL, INTENT(OUT) :: xmed
END SUBROUTINE median
END INTERFACE
DO
WRITE(*, *)'Enter n: '
READ(*, *) n
ALLOCATE( x(n) )
CALL RANDOM_NUMBER(x)
CALL median(x, n, xmed)
WRITE(*, 900) x(1), x(n), xmed
900 FORMAT(' First & last = ', 2f10.4, ' Median = ', f10.4/)
DEALLOCATE( x )
END DO
STOP
END PROGRAM t_median
SUBROUTINE median(x, n, xmed)
! Find the median of X(1), ... , X(N), using as much of the quicksort
! algorithm as is needed to isolate it.
! N.B. On exit, the array X is partially ordered.
! Latest revision - 26 November 1996
IMPLICIT NONE
INTEGER, INTENT(IN) :: n
REAL, INTENT(IN OUT), DIMENSION(:) :: x
REAL, INTENT(OUT) :: xmed
! Local variables
REAL :: temp, xhi, xlo, xmax, xmin
LOGICAL :: odd
INTEGER :: hi, lo, nby2, nby2p1, mid, i, j, k
nby2 = n / 2
nby2p1 = nby2 + 1
odd = .true.
! HI & LO are position limits encompassing the median.
IF (n == 2 * nby2) odd = .false.
lo = 1
hi = n
IF (n < 3) THEN
IF (n < 1) THEN
xmed = 0.0
RETURN
END IF
xmed = x(1)
IF (n == 1) RETURN
xmed = 0.5*(xmed + x(2))
RETURN
END IF
! Find median of 1st, middle & last values.
10 mid = (lo + hi)/2
xmed = x(mid)
xlo = x(lo)
xhi = x(hi)
IF (xhi < xlo) THEN ! Swap xhi & xlo
temp = xhi
xhi = xlo
xlo = temp
END IF
IF (xmed > xhi) THEN
xmed = xhi
ELSE IF (xmed < xlo) THEN
xmed = xlo
END IF
! The basic quicksort algorithm to move all values <= the sort key (XMED)
! to the left-hand end, and all higher values to the other end.
i = lo
j = hi
50 DO
IF (x(i) >= xmed) EXIT
i = i + 1
END DO
DO
IF (x(j) <= xmed) EXIT
j = j - 1
END DO
IF (i < j) THEN
temp = x(i)
x(i) = x(j)
x(j) = temp
i = i + 1
j = j - 1
! Decide which half the median is in.
IF (i <= j) GO TO 50
END IF
IF (.NOT. odd) THEN
IF (j == nby2 .AND. i == nby2p1) GO TO 130
IF (j < nby2) lo = i
IF (i > nby2p1) hi = j
IF (i /= j) GO TO 100
IF (i == nby2) lo = nby2
IF (j == nby2p1) hi = nby2p1
ELSE
IF (j < nby2p1) lo = i
IF (i > nby2p1) hi = j
IF (i /= j) GO TO 100
! Test whether median has been isolated.
IF (i == nby2p1) RETURN
END IF
100 IF (lo < hi - 1) GO TO 10
IF (.NOT. odd) THEN
xmed = 0.5*(x(nby2) + x(nby2p1))
RETURN
END IF
temp = x(lo)
IF (temp > x(hi)) THEN
x(lo) = x(hi)
x(hi) = temp
END IF
xmed = x(nby2p1)
RETURN
! Special case, N even, J = N/2 & I = J + 1, so the median is
! between the two halves of the series. Find max. of the first
! half & min. of the second half, then average.
130 xmax = x(1)
DO k = lo, j
xmax = MAX(xmax, x(k))
END DO
xmin = x(n)
DO k = i, hi
xmin = MIN(xmin, x(k))
END DO
xmed = 0.5*(xmin + xmax)
RETURN
END SUBROUTINE median