dsort.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	dsort (dx, iy, n, kflag)

Function/Subroutine Documentation

dsort()

subroutine dsort	(	double precision, dimension(*), intent(inout)	dx,
		integer, dimension(*), intent(inout)	iy,
		integer, intent(in)	n,
		integer, intent(in)	kflag
	)

c
c    slight change: XERMSG was removed; error messages are
c                   led to the screen; 
c
C***BEGIN PROLOGUE  DSORT
C***PURPOSE  Sort an array and optionally make the same interchanges in
C            an auxiliary array.  The array may be sorted in increasing
C            or decreasing order.  A slightly modified QUICKSORT
C            algorithm is used.
C***LIBRARY   SLATEC
C***CATEGORY  N6A2B
C***TYPE      DOUBLE PRECISION (SSORT-S, DSORT-D, ISORT-I)
C***KEYWORDS  SINGLETON QUICKSORT, SORT, SORTING
C***AUTHOR  Jones, R. E., (SNLA)
C           Wisniewski, J. A., (SNLA)
C***ROUTINES CALLED  XERMSG
C***DESCRIPTION
C
C   DSORT sorts array DX and optionally makes the same interchanges in
C   array IY.  The array DX may be sorted in increasing order or
C   decreasing order.  A slightly modified quicksort algorithm is used.
C
C   Description of Parameters
C      DX - array of values to be sorted   (usually abscissas)
C      IY - array to be (optionally) carried along
C      N  - number of values in array DX to be sorted
C      KFLAG - control parameter
C            =  2  means sort DX in increasing order and carry IY along.
C            =  1  means sort DX in increasing order (ignoring IY)
C            = -1  means sort DX in decreasing order (ignoring IY)
C            = -2  means sort DX in decreasing order and carry IY along.
C
C***REFERENCES  R. C. Singleton, Algorithm 347, An efficient algorithm
C                 for sorting with minimal storage, Communications of
C                 the ACM, 12, 3 (1969), pp. 185-187.
C***REVISION HISTORY  (YYMMDD)
C   761101  DATE WRITTEN
C   761118  Modified to use the Singleton quicksort algorithm.  (JAW)
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890831  Modified array declarations.  (WRB)
C   891009  Removed unreferenced statement labels.  (WRB)
C   891024  Changed category.  (WRB)
C   891024  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   901012  Declared all variables; changed X,Y to DX,IY; changed
C           code to parallel SSORT. (M. McClain)
C   920501  Reformatted the REFERENCES section.  (DWL, WRB)
C   920519  Clarified error messages.  (DWL)
C   920801  Declarations section rebuilt and code restructured to use
C           IF-THEN-ELSE-ENDIF.  (RWC, WRB)
!   100411  changed the dimension of IL and IU from 21 to 31.
!   150514  inserted intent statements
!
!     field IL and IU have the dimension 31. This is log2 of the largest
!     array size to be sorted. If arrays larger than 2**31 in length have
!     to be sorted, this dimension has to be modified accordingly
!
C***END PROLOGUE  DSORT
      implicit none
C     .. Scalar Arguments ..
      integer kflag, n,iy(*),ty,tty
C     .. Array Arguments ..
      double precision dx(*)
C     .. Local Scalars ..
      double precision r, t, tt
      integer i, ij, j, k, kk, l, m, nn
C     .. Local Arrays ..
      integer il(31), iu(31)
C     .. External Subroutines ..
c      EXTERNAL XERMSG
C     .. Intrinsic Functions ..
      intrinsic abs, int
!
      intent(in) n,kflag 
!
      intent(inout) dx,iy
C***FIRST EXECUTABLE STATEMENT  DSORT
      nn = n
      if (nn .lt. 1) then
         write(*,*) '*error in dsort: the number of values to be'
         write(*,*) '       sorted is not positive'
         call exit(201)
      endif
C
      kk = abs(kflag)
      if (kk.ne.1 .and. kk.ne.2) then
         write(*,*) '*error in dsort: the sort control parameter is'
         write(*,*) '       not 2, 1, -1 or -2'
         call exit(201)
      endif
C
C     Alter array DX to get decreasing order if needed
C
      if (kflag .le. -1) then
         do 10 i=1,nn
            dx(i) = -dx(i)
   10    continue
      endif
C
      if (kk .eq. 2) go to 100
C
C     Sort DX only
C
      m = 1
      i = 1
      j = nn
      r = 0.375d0
C
   20 if (i .eq. j) go to 60
      if (r .le. 0.5898437d0) then
         r = r+3.90625d-2
      else
         r = r-0.21875d0
      endif
C
   30 k = i
C
C     Select a central element of the array and save it in location T
C
      ij = i + int((j-i)*r)
      t = dx(ij)
C
C     If first element of array is greater than T, interchange with T
C
      if (dx(i) .gt. t) then
         dx(ij) = dx(i)
         dx(i) = t
         t = dx(ij)
      endif
      l = j
C
C     If last element of array is less than than T, interchange with T
C
      if (dx(j) .lt. t) then
         dx(ij) = dx(j)
         dx(j) = t
         t = dx(ij)
C
C        If first element of array is greater than T, interchange with T
C
         if (dx(i) .gt. t) then
            dx(ij) = dx(i)
            dx(i) = t
            t = dx(ij)
         endif
      endif
C
C     Find an element in the second half of the array which is smaller
C     than T
C
   40 l = l-1
      if (dx(l) .gt. t) go to 40
C
C     Find an element in the first half of the array which is greater
C     than T
C
   50 k = k+1
      if (dx(k) .lt. t) go to 50
C
C     Interchange these elements
C
      if (k .le. l) then
         tt = dx(l)
         dx(l) = dx(k)
         dx(k) = tt
         go to 40
      endif
C
C     Save upper and lower subscripts of the array yet to be sorted
C
      if (l-i .gt. j-k) then
         il(m) = i
         iu(m) = l
         i = k
         m = m+1
      else
         il(m) = k
         iu(m) = j
         j = l
         m = m+1
      endif
      go to 70
C
C     Begin again on another portion of the unsorted array
C
   60 m = m-1
      if (m .eq. 0) go to 190
      i = il(m)
      j = iu(m)
C
   70 if (j-i .ge. 1) go to 30
      if (i .eq. 1) go to 20
      i = i-1
C
   80 i = i+1
      if (i .eq. j) go to 60
      t = dx(i+1)
      if (dx(i) .le. t) go to 80
      k = i
C
   90 dx(k+1) = dx(k)
      k = k-1
      if (t .lt. dx(k)) go to 90
      dx(k+1) = t
      go to 80
C
C     Sort DX and carry IY along
C
  100 m = 1
      i = 1
      j = nn
      r = 0.375d0
C
  110 if (i .eq. j) go to 150
      if (r .le. 0.5898437d0) then
         r = r+3.90625d-2
      else
         r = r-0.21875d0
      endif
C
  120 k = i
C
C     Select a central element of the array and save it in location T
C
      ij = i + int((j-i)*r)
      t = dx(ij)
      ty = iy(ij)
C
C     If first element of array is greater than T, interchange with T
C
      if (dx(i) .gt. t) then
         dx(ij) = dx(i)
         dx(i) = t
         t = dx(ij)
         iy(ij) = iy(i)
         iy(i) = ty
         ty = iy(ij)
      endif
      l = j
C
C     If last element of array is less than T, interchange with T
C
      if (dx(j) .lt. t) then
         dx(ij) = dx(j)
         dx(j) = t
         t = dx(ij)
         iy(ij) = iy(j)
         iy(j) = ty
         ty = iy(ij)
C
C        If first element of array is greater than T, interchange with T
C
         if (dx(i) .gt. t) then
            dx(ij) = dx(i)
            dx(i) = t
            t = dx(ij)
            iy(ij) = iy(i)
            iy(i) = ty
            ty = iy(ij)
         endif
      endif
C
C     Find an element in the second half of the array which is smaller
C     than T
C
  130 l = l-1
      if (dx(l) .gt. t) go to 130
C
C     Find an element in the first half of the array which is greater
C     than T
C
  140 k = k+1
      if (dx(k) .lt. t) go to 140
C
C     Interchange these elements
C
      if (k .le. l) then
         tt = dx(l)
         dx(l) = dx(k)
         dx(k) = tt
         tty = iy(l)
         iy(l) = iy(k)
         iy(k) = tty
         go to 130
      endif
C
C     Save upper and lower subscripts of the array yet to be sorted
C
      if (l-i .gt. j-k) then
         il(m) = i
         iu(m) = l
         i = k
         m = m+1
      else
         il(m) = k
         iu(m) = j
         j = l
         m = m+1
      endif
      go to 160
C
C     Begin again on another portion of the unsorted array
C
  150 m = m-1
      if (m .eq. 0) go to 190
      i = il(m)
      j = iu(m)
C
  160 if (j-i .ge. 1) go to 120
      if (i .eq. 1) go to 110
      i = i-1
C
  170 i = i+1
      if (i .eq. j) go to 150
      t = dx(i+1)
      ty = iy(i+1)
      if (dx(i) .le. t) go to 170
      k = i
C
  180 dx(k+1) = dx(k)
      iy(k+1) = iy(k)
      k = k-1
      if (t .lt. dx(k)) go to 180
      dx(k+1) = t
      iy(k+1) = ty
      go to 170
C
C     Clean up
C
  190 if (kflag .le. -1) then
         do 200 i=1,nn
            dx(i) = -dx(i)
  200    continue
      endif
      return