SORT Subroutine

public subroutine SORT(X, n, Y)

NAME

sort(3f) - [M_datapac:SORT] sort a vector of sample
observations, also return the positions in the original vector

SYNOPSIS

 SUBROUTINE SORT(X,N,Y)

  real,intent(in)    :: x(:)
  integer,intent(in) :: n
  real,intent(out)   :: y(:)

DESCRIPTION

This subroutine sorts (in ascending order) the N elements of the
REAL vector X using the binary sort algorithm and puts
the resulting N sorted values into the REAL vector Y.

OPTIONS

INPUT

 X  The REAL vector of observations to be sorted.
    The input vector X remains unaltered.

 N  The integer number of observations in the vector X.

OUTPUT

 Y  The REAL vector into which the sorted data values from X will be
    placed in ascending order.

EXAMPLES

Sample program:

program demo_sort
use M_datapac, only : sort
implicit none
integer,parameter            :: isz=20
real                         :: aa(isz)
real                         :: bb(isz)
integer                      :: i
   write(*,*)'initializing array with ',isz,' random numbers'
   call random_seed()
   CALL RANDOM_NUMBER(aa)
   aa=aa*450000.0
   bb=real([(i,i=1,isz)])

   call sort(aa,isz,bb) ! sort data

   write(*,*)'checking if real values are sorted(3f)'
   do i=1,isz-1
      if(bb(i).gt.bb(i+1))then
         write(*,*)'Error in sorting reals small to large ',i,bb(i),bb(i+1)
      endif
   enddo
  write(*,'(2(g0,1x))')'ORIGINAL','SORTED',(aa(i),bb(i),i=1,isz)

end program demo_sort

Results:

AUTHOR

The original DATAPAC library was written by James Filliben of the Statistical
Engineering Division, National Institute of Standards and Technology.

MAINTAINER

John Urban, 2022.05.31

REFERENCES

1. CACM MARCH 1969, page 186 (BINARY SORT ALGORITHM BY RICHARD C. SINGLETON).
2. CACM JANUARY 1970, page 54.
3. CACM OCTOBER 1970, page 624.
1. JACM JANUARY 1961, page 41.

LICENSE

CC0-1.0

Arguments

Type IntentOptional Attributes Name
real(kind=wp), dimension(:) :: X
integer :: n
real(kind=wp), dimension(:) :: Y

Source Code

SUBROUTINE SORT(X,N,Y)
REAL(kind=wp) :: amed, hold, tt, X, Y
integer i, il, ip1, iu, j, jmi, jmk, k, l, lmi,   m, mid, n, nm1
DIMENSION X(:), Y(:)
DIMENSION iu(36), il(36)

!     RESTRICTIONS--THE DIMENSIONS OF THE VECTORS IL AND IU
!                   (DEFINED AND USED INTERNALLY WITHIN
!                   THIS SUBROUTINE) DICTATE THE MAXIMUM
!                   ALLOWABLE VALUE OF N FOR THIS SUBROUTINE.
!                   IF IL AND IU EACH HAVE DIMENSION K,
!                   THEN N MAY NOT EXCEED 2**(K+1) - 1.
!                   FOR THIS SUBROUTINE AS WRITTEN, THE DIMENSIONS
!                   OF IL AND IU HAVE BEEN SET TO 36,
!                   THUS THE MAXIMUM ALLOWABLE VALUE OF N IS
!                   APPROXIMATELY 137 BILLION.
!                   SINCE THIS EXCEEDS THE MAXIMUM ALLOWABLE
!                   VALUE FOR AN INTEGER VARIABLE IN MANY COMPUTERS,
!                   AND SINCE A SORT OF 137 BILLION ELEMENTS
!                   IS PRESENTLY IMPRACTICAL AND UNLIKELY,
!                   THEN THERE IS NO PRACTICAL RESTRICTION
!                   ON THE MAXIMUM VALUE OF N FOR THIS SUBROUTINE.
!                   (IN LIGHT OF THE ABOVE, NO CHECK OF THE
!                   UPPER LIMIT OF N HAS BEEN INCORPORATED
!                   INTO THIS SUBROUTINE.)

!     COMMENT--THE SMALLEST ELEMENT OF THE VECTOR X
!              WILL BE PLACED IN THE FIRST POSITION
!              OF THE VECTOR Y,
!              THE SECOND SMALLEST ELEMENT IN THE VECTOR X
!              WILL BE PLACED IN THE SECOND POSITION
!              OF THE VECTOR Y, ETC.
!     COMMENT--IF THE ANALYST DESIRES A SORT 'IN PLACE',
!              THIS IS DONE BY HAVING THE SAME
!              OUTPUT VECTOR AS INPUT VECTOR IN THE CALLING SEQUENCE.
!              THUS, FOR EXAMPLE, THE CALLING SEQUENCE
!              CALL SORT(X,N,X)
!              IS ALLOWABLE AND WILL RESULT IN
!              THE DESIRED 'IN-PLACE' SORT.
!     COMMENT--THE SORTING ALGORTHM USED HEREIN
!              IS THE BINARY SORT.
!              THIS ALGORTHIM IS EXTREMELY FAST AS THE
!              FOLLOWING TIME TRIALS INDICATE.
!              THESE TIME TRIALS WERE CARRIED OUT ON THE
!              UNIVAC 1108 EXEC 8 SYSTEM AT NBS
!              IN AUGUST OF 1974.
!              BY WAY OF COMPARISON, THE TIME TRIAL VALUES
!              FOR THE EASY-TO-PROGRAM BUT EXTREMELY
!              INEFFICIENT BUBBLE SORT ALGORITHM HAVE
!              ALSO BEEN INCLUDED--
!              NUMBER OF RANDOM        BINARY SORT       BUBBLE SORT
!               NUMBERS SORTED
!                N = 10                 .002 SEC          .002 SEC
!                N = 100                .011 SEC          .045 SEC
!                N = 1000               .141 SEC         4.332 SEC
!                N = 3000               .476 SEC        37.683 SEC
!                N = 10000             1.887 SEC      NOT COMPUTED
!     ORIGINAL VERSION--JUNE      1972.
!     UPDATED         --NOVEMBER  1975.
!
!---------------------------------------------------------------------
!
!     CHECK THE INPUT ARGUMENTS FOR ERRORS
!
      IF ( N<1 ) THEN
         WRITE (G_IO,99001)
99001    FORMAT (' ',                                                   &
     &'***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE SORT   SUBROU&
     &TINE IS NON-POSITIVE *****')
         WRITE (G_IO,99002) N
99002    FORMAT (' ','***** THE VALUE OF THE ARGUMENT IS ',I0,' *****')
         RETURN
      ELSE
         IF ( N==1 ) THEN
            WRITE (G_IO,99003)
99003       FORMAT (' ','***** NON-FATAL DIAGNOSTIC--THE SECOND INPUT ARGUMENT TO THE SORT SUBROUTINE HAS THE VALUE 1 *****')
            Y(1) = X(1)
            RETURN
         ELSE
            hold = X(1)
            DO i = 2 , N
               IF ( X(i)/=hold ) GOTO 50
            ENDDO
            WRITE (G_IO,99004) hold
99004       FORMAT (' ',&
            & '***** NON-FATAL DIAGNOSTIC--THE FIRST  INPUT ARGUMENT (A VECTOR) TO THE SORT   SUBROUTINE HAS ALL ELEMENTS =',&
            & E15.8,&
            & ' *****')
            DO i = 1 , N
               Y(i) = X(i)
            ENDDO
            RETURN
         ENDIF
!
!-----START POINT-----------------------------------------------------
!
!     COPY THE VECTOR X INTO THE VECTOR Y
 50      DO i = 1 , N
            Y(i) = X(i)
         ENDDO
!
!     CHECK TO SEE IF THE INPUT VECTOR IS ALREADY SORTED
!
         nm1 = N - 1
         DO i = 1 , nm1
            ip1 = i + 1
            IF ( Y(i)>Y(ip1) ) GOTO 100
         ENDDO
         RETURN
      ENDIF
 100  m = 1
      i = 1
      j = N
 200  IF ( i>=j ) GOTO 400
 300  k = i
      mid = (i+j)/2
      amed = Y(mid)
      IF ( Y(i)>amed ) THEN
         Y(mid) = Y(i)
         Y(i) = amed
         amed = Y(mid)
      ENDIF
      l = j
      IF ( Y(j)<amed ) THEN
         Y(mid) = Y(j)
         Y(j) = amed
         amed = Y(mid)
         IF ( Y(i)>amed ) THEN
            Y(mid) = Y(i)
            Y(i) = amed
            amed = Y(mid)
         ENDIF
      ENDIF
      DO
         l = l - 1
         IF ( Y(l)<=amed ) THEN
            tt = Y(l)
            DO
               k = k + 1
               IF ( Y(k)>=amed ) THEN
                  IF ( k<=l ) THEN
                     Y(l) = Y(k)
                     Y(k) = tt
                     EXIT
                  ELSE
                     lmi = l - i
                     jmk = j - k
                     IF ( lmi<=jmk ) THEN
                        il(m) = k
                        iu(m) = j
                        j = l
                        m = m + 1
                     ELSE
                        il(m) = i
                        iu(m) = l
                        i = k
                        m = m + 1
                     ENDIF
                     GOTO 500
                  ENDIF
               ENDIF
            ENDDO
         ENDIF
      ENDDO
 400  m = m - 1
      IF ( m==0 ) RETURN
      i = il(m)
      j = iu(m)
 500  jmi = j - i
      IF ( jmi>=11 ) GOTO 300
      IF ( i==1 ) GOTO 200
      i = i - 1
      DO
         i = i + 1
         IF ( i==j ) GOTO 400
         amed = Y(i+1)
         IF ( Y(i)>amed ) THEN
            k = i
            DO
               Y(k+1) = Y(k)
               k = k - 1
               IF ( amed>=Y(k) ) THEN
                  Y(k+1) = amed
                  EXIT
               ENDIF
            ENDDO
         ENDIF
      ENDDO
END SUBROUTINE SORT