!DEC$ ATTRIBUTES DLLEXPORT :: DSORTILU SUBROUTINE DSORTilu (DX, IX, N, KFLAG, IERR) 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***MODIFIED by K.Lipnikov (http://humpty.inm.ras.ru/~kos) C***DESCRIPTION C C DSORT sorts array DX and optionally makes the same interchanges in C array IX. 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 IX - 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 IX along. C = 1 means sort DX in increasing order (ignoring IX) C = -1 means sort DX in decreasing order (ignoring IX) C = -2 means sort DX in decreasing order and carry IX along. C C IERR - the error's indicator C = 0 - Okay C != 0 - errors 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***ROUTINES CALLED XERMSG 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,IX; 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) C***END PROLOGUE DSORT C .. Scalar Arguments .. INTEGER KFLAG, N, IERR C .. Array Arguments .. INTEGER DX(*) DOUBLE PRECISION IX(*) C .. Local Scalars .. DOUBLE PRECISION R INTEGER T, TT DOUBLE PRECISION ITX,ITTX INTEGER I, IJ, J, K, KK, L, M, NN C .. Local Arrays .. INTEGER IL(21), IU(21) C .. Intrinsic Functions .. INTRINSIC ABS, INT C***FIRST EXECUTABLE STATEMENT DSORT IERR = 0 NN = N IF (NN .LT. 1) THEN IERR = 1201 RETURN ENDIF C KK = ABS(KFLAG) IF (KK.NE.1 .AND. KK.NE.2) THEN IERR = 1202 RETURN ENDIF C 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 IX 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) ITX = IX(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) IX(IJ) = IX(I) IX(I) = ITX ITX = IX(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) IX(IJ) = IX(J) IX(J) = ITX ITX = IX(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) IX(IJ) = IX(I) IX(I) = ITX ITX = IX(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 ITTX = IX(L) IX(L) = IX(K) IX(K) = ITTX 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) ITX = IX(I+1) IF (DX(I) .LE. T) GO TO 170 K = I C 180 DX(K+1) = DX(K) IX(K+1) = IX(K) K = K-1 IF (T .LT. DX(K)) GO TO 180 DX(K+1) = T IX(K+1) = ITX 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 END