*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