*DECK DPSORT
SUBROUTINE DPSORT (DX, N, IPERM, KFLAG, IER)
C***BEGIN PROLOGUE DPSORT
C***PURPOSE Return the permutation vector generated by sorting a given
C array and, optionally, rearrange the elements of the array.
C The array may be sorted in increasing or decreasing order.
C A slightly modified quicksort algorithm is used.
C***LIBRARY SLATEC
C***CATEGORY N6A1B, N6A2B
C***TYPE DOUBLE PRECISION (SPSORT-S, DPSORT-D, IPSORT-I, HPSORT-H)
C***KEYWORDS NUMBER SORTING, PASSIVE SORTING, SINGLETON QUICKSORT, SORT
C***AUTHOR Jones, R. E., (SNLA)
C Rhoads, G. S., (NBS)
C Wisniewski, J. A., (SNLA)
C***DESCRIPTION
C
C DPSORT returns the permutation vector IPERM generated by sorting
C the array DX and, optionally, rearranges the values in DX. DX may
C be sorted in increasing or decreasing order. A slightly modified
C quicksort algorithm is used.
C
C IPERM is such that DX(IPERM(I)) is the Ith value in the
C rearrangement of DX. IPERM may be applied to another array by
C calling IPPERM, SPPERM, DPPERM or HPPERM.
C
C The main difference between DPSORT and its active sorting equivalent
C DSORT is that the data are referenced indirectly rather than
C directly. Therefore, DPSORT should require approximately twice as
C long to execute as DSORT. However, DPSORT is more general.
C
C Description of Parameters
C DX - input/output -- double precision array of values to be
C sorted. If ABS(KFLAG) = 2, then the values in DX will be
C rearranged on output; otherwise, they are unchanged.
C N - input -- number of values in array DX to be sorted.
C IPERM - output -- permutation array such that IPERM(I) is the
C index of the value in the original order of the
C DX array that is in the Ith location in the sorted
C order.
C KFLAG - input -- control parameter:
C = 2 means return the permutation vector resulting from
C sorting DX in increasing order and sort DX also.
C = 1 means return the permutation vector resulting from
C sorting DX in increasing order and do not sort DX.
C = -1 means return the permutation vector resulting from
C sorting DX in decreasing order and do not sort DX.
C = -2 means return the permutation vector resulting from
C sorting DX in decreasing order and sort DX also.
C IER - output -- error indicator:
C = 0 if no error,
C = 1 if N is zero or negative,
C = 2 if KFLAG is not 2, 1, -1, or -2.
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***ROUTINES CALLED XERMSG
C***REVISION HISTORY (YYMMDD)
C 761101 DATE WRITTEN
C 761118 Modified by John A. Wisniewski to use the Singleton
C quicksort algorithm.
C 870423 Modified by Gregory S. Rhoads for passive sorting with the
C option for the rearrangement of the original data.
C 890619 Double precision version of SPSORT created by D. W. Lozier.
C 890620 Algorithm for rearranging the data vector corrected by R.
C Boisvert.
C 890622 Prologue upgraded to Version 4.0 style by D. Lozier.
C 891128 Error when KFLAG.LT.0 and N=1 corrected by R. Boisvert.
C 920507 Modified by M. McClain to revise prologue text.
C 920818 Declarations section rebuilt and code restructured to use
C IF-THEN-ELSE-ENDIF. (SMR, WRB)
C***END PROLOGUE DPSORT
C .. Scalar Arguments ..
INTEGER IER, KFLAG, N
C .. Array Arguments ..
DOUBLE PRECISION DX(*)
INTEGER IPERM(*)
C .. Local Scalars ..
DOUBLE PRECISION R, TEMP
INTEGER I, IJ, INDX, INDX0, ISTRT, J, K, KK, L, LM, LMT, M, NN
C .. Local Arrays ..
INTEGER IL(21), IU(21)
C .. External Subroutines ..
EXTERNAL XERMSG
C .. Intrinsic Functions ..
INTRINSIC ABS, INT
C***FIRST EXECUTABLE STATEMENT DPSORT
IER = 0
NN = N
IF (NN .LT. 1) THEN
IER = 1
CALL XERMSG ('SLATEC', 'DPSORT',
+ 'The number of values to be sorted, N, is not positive.',
+ IER, 1)
RETURN
ENDIF
C
KK = ABS(KFLAG)
IF (KK.NE.1 .AND. KK.NE.2) THEN
IER = 2
CALL XERMSG ('SLATEC', 'DPSORT',
+ 'The sort control parameter, KFLAG, is not 2, 1, -1, or -2.',
+ IER, 1)
RETURN
ENDIF
C
C Initialize permutation vector
C
DO 10 I=1,NN
IPERM(I) = I
10 CONTINUE
C
C Return if only one value is to be sorted
C
IF (NN .EQ. 1) RETURN
C
C Alter array DX to get decreasing order if needed
C
IF (KFLAG .LE. -1) THEN
DO 20 I=1,NN
DX(I) = -DX(I)
20 CONTINUE
ENDIF
C
C Sort DX only
C
M = 1
I = 1
J = NN
R = .375D0
C
30 IF (I .EQ. J) GO TO 80
IF (R .LE. 0.5898437D0) THEN
R = R+3.90625D-2
ELSE
R = R-0.21875D0
ENDIF
C
40 K = I
C
C Select a central element of the array and save it in location L
C
IJ = I + INT((J-I)*R)
LM = IPERM(IJ)
C
C If first element of array is greater than LM, interchange with LM
C
IF (DX(IPERM(I)) .GT. DX(LM)) THEN
IPERM(IJ) = IPERM(I)
IPERM(I) = LM
LM = IPERM(IJ)
ENDIF
L = J
C
C If last element of array is less than LM, interchange with LM
C
IF (DX(IPERM(J)) .LT. DX(LM)) THEN
IPERM(IJ) = IPERM(J)
IPERM(J) = LM
LM = IPERM(IJ)
C
C If first element of array is greater than LM, interchange
C with LM
C
IF (DX(IPERM(I)) .GT. DX(LM)) THEN
IPERM(IJ) = IPERM(I)
IPERM(I) = LM
LM = IPERM(IJ)
ENDIF
ENDIF
GO TO 60
50 LMT = IPERM(L)
IPERM(L) = IPERM(K)
IPERM(K) = LMT
C
C Find an element in the second half of the array which is smaller
C than LM
C
60 L = L-1
IF (DX(IPERM(L)) .GT. DX(LM)) GO TO 60
C
C Find an element in the first half of the array which is greater
C than LM
C
70 K = K+1
IF (DX(IPERM(K)) .LT. DX(LM)) GO TO 70
C
C Interchange these elements
C
IF (K .LE. L) GO TO 50
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 90
C
C Begin again on another portion of the unsorted array
C
80 M = M-1
IF (M .EQ. 0) GO TO 120
I = IL(M)
J = IU(M)
C
90 IF (J-I .GE. 1) GO TO 40
IF (I .EQ. 1) GO TO 30
I = I-1
C
100 I = I+1
IF (I .EQ. J) GO TO 80
LM = IPERM(I+1)
IF (DX(IPERM(I)) .LE. DX(LM)) GO TO 100
K = I
C
110 IPERM(K+1) = IPERM(K)
K = K-1
IF (DX(LM) .LT. DX(IPERM(K))) GO TO 110
IPERM(K+1) = LM
GO TO 100
C
C Clean up
C
120 IF (KFLAG .LE. -1) THEN
DO 130 I=1,NN
DX(I) = -DX(I)
130 CONTINUE
ENDIF
C
C Rearrange the values of DX if desired
C
IF (KK .EQ. 2) THEN
C
C Use the IPERM vector as a flag.
C If IPERM(I) < 0, then the I-th value is in correct location
C
DO 150 ISTRT=1,NN
IF (IPERM(ISTRT) .GE. 0) THEN
INDX = ISTRT
INDX0 = INDX
TEMP = DX(ISTRT)
140 IF (IPERM(INDX) .GT. 0) THEN
DX(INDX) = DX(IPERM(INDX))
INDX0 = INDX
IPERM(INDX) = -IPERM(INDX)
INDX = ABS(IPERM(INDX))
GO TO 140
ENDIF
DX(INDX0) = TEMP
ENDIF
150 CONTINUE
C
C Revert the signs of the IPERM values
C
DO 160 I=1,NN
IPERM(I) = -IPERM(I)
160 CONTINUE
C
ENDIF
C
RETURN
END