************************************************************ * * Robert Renka * Oak Ridge Natl. Lab. * * This subroutine uses an order N*LOG(N) quick sort to sort an integer * array X into increasing order. The algorithm is as follows. IND is * initialized to the ordered sequence of indices 1,...,N, and all * interchanges are applied to IND. X is divided into two portions by * picking a central element T. The first and last elements are compared * with T, and interchanges are applied as necessary so that the three * values are in ascending order. Interchanges are then applied to that * all elements greater than T are in the upper portion of the array and * all elements less than T are in the lower portion. The upper and lower * indices of one of the portions are saved in local arrays, and the * process is repeated iteratively on the other portion. When a portion * is completely sorted, the process begins again by retrieving the * indices bounding another unsorted portion. * * INPUT PARAMETERS - N - Length of the array X. * * X - Vector of length N to be sorted. * * IND - Vector of length .GE. N. * * N and X are not altered by this routine. * * OUTPUT PARAMETER - IND - Sequence of indices 1,...,N * permuted in the same fashion as X * would be. Thus, the ordering on * X is defined by Y(I) = X(IND(I)). * * INTRINSIC FUNCTIONS CALLED BY QSORTI - IFIX, FLOAT * ************************************************************ * * NOTE -- IU and IL must be dimensioned .GE. LOG(N) where LOG has base 2. * ************************************************************ SUBROUTINE QSORTI (X, IND, N) INTEGER N, X(N), IND(N) INTEGER IU(32), IL(32) INTEGER M, I, J, K, L, IJ, IT, ITT, INDX, T REAL R * LOCAL PARAMETERS - * * IU,IL = temporary storage for the upper and lower * indices of portions of the array * M = index for IU and IL * I,J = lower and upper indices of a portion of X * K,L = indices in the range I,...,J * IJ = randomly chosen index between I and J * IT,ITT = temporary storage for interchanges in IND * INDX = temporary index for X * R = pseudo random number for generating IJ * T = central element of X IF (N .LE. 0) RETURN * Initialize IND, M, I, J, and R DO 1 I = 1,N 1 IND(I) = I M = 1 I = 1 J = N R = .375 * top of loop 2 IF (I .GE. J) GO TO 10 IF (R > .5898437) GO TO 3 R = R + .0390625 GO TO 4 3 R = R - .21875 * initialize K 4 K = I * select a central element of X and save it in T IJ = I + IFIX(R*FLOAT(J-I)) IT = IND(IJ) T = X(IT) * If the first element of the array is greater than T, * interchange it with T INDX = IND(I) IF (X(INDX) .LE. T) GO TO 5 IND(IJ) = INDX IND(I) = IT IT = INDX T = X(IT) * initialize L 5 L = J * if the last element of the array is less than T, * interchange it with T INDX = IND(J) IF (X(INDX) .GE. T) GO TO 7 IND(IJ) = INDX IND(J) = IT IT = INDX T = X(IT) * if the first element of the array is greater than T, * interchange it with T INDX = IND(I) IF (X(INDX) .LE. T) GO TO 7 IND(IJ) = INDX IND(I) = IT IT = INDX T = X(IT) GO TO 7 * interchange elements K and L 6 ITT = IND(L) IND(L) = IND(K) IND(K) = ITT * Find an element in the upper part of the array which is * not larger than T 7 L = L - 1 INDX = IND(L) IF (X(INDX) > T) GO TO 7 * Find an element in the lower part of the array which is * not smaller than T 8 K = K + 1 INDX = IND(K) IF (X(INDX) < T) GO TO 8 * if K .LE. L, interchange elements K and L IF (K .LE. L) GO TO 6 * save the upper and lower subscripts of the portion of the * array yet to be sorted IF (L-I .LE. J-K) GO TO 9 IL(M) = I IU(M) = L I = K M = M + 1 GO TO 11 9 IL(M) = K IU(M) = J J = L M = M + 1 GO TO 11 * begin again on another unsorted portion of the array 10 M = M - 1 IF (M .EQ. 0) RETURN I = IL(M) J = IU(M) 11 IF (J-I .GE. 11) GO TO 4 IF (I .EQ. 1) GO TO 2 I = I - 1 * sort elements I+1,...,J. Note that 1 .LE. I < J and J-I < 11. 12 I = I + 1 IF (I .EQ. J) GO TO 10 INDX = IND(I+1) T = X(INDX) IT = INDX INDX = IND(I) IF (X(INDX) .LE. T) GO TO 12 K = I 13 IND(K+1) = IND(K) K = K - 1 INDX = IND(K) IF (T < X(INDX)) GO TO 13 IND(K+1) = IT GO TO 12 END ! of QSORTI ************************************************************************ SUBROUTINE QSORTR (X, IND, N) INTEGER N, IND(N) REAL X(N) INTEGER IU(32), IL(32) INTEGER M, I, J, K, L, IJ, IT, ITT, INDX REAL R, T * LOCAL PARAMETERS - * * IU,IL = temporary storage for the upper and lower * indices of portions of the array * M = index for IU and IL * I,J = lower and upper indices of a portion of X * K,L = indices in the range I,...,J * IJ = randomly chosen index between I and J * IT,ITT = temporary storage for interchanges in IND * INDX = temporary index for X * R = pseudo random number for generating IJ * T = central element of X ************************************************************************ IF (N .LE. 0) RETURN * initialize IND, M, I, J, and R DO 1 I = 1, N 1 IND(I) = I M = 1 I = 1 J = N R = .375 * top of loop 2 IF (I .GE. J) GO TO 10 IF (R > .5898437) GO TO 3 R = R + .0390625 GO TO 4 3 R = R - .21875 * initialize K 4 K = I * select a central element of X and save it in T IJ = I + IFIX(R*FLOAT(J-I)) IT = IND(IJ) T = X(IT) * if the first element of the array is greater than T, interchange it with T INDX = IND(I) IF (X(INDX) .LE. T) GO TO 5 IND(IJ) = INDX IND(I) = IT IT = INDX T = X(IT) * initialize L 5 L = J * if the last element of the array is less than T, interchange it with T INDX = IND(J) IF (X(INDX) .GE. T) GO TO 7 IND(IJ) = INDX IND(J) = IT IT = INDX T = X(IT) * if the first element of the array is greater than T, interchange it with T INDX = IND(I) IF (X(INDX) .LE. T) GO TO 7 IND(IJ) = INDX IND(I) = IT IT = INDX T = X(IT) GO TO 7 * interchange elements K and L 6 ITT = IND(L) IND(L) = IND(K) IND(K) = ITT * find an element in the upper part of the array which is not larger than T 7 L = L - 1 INDX = IND(L) IF (X(INDX) > T) GO TO 7 * find an element in the lower part of the array which is not smaller than T 8 K = K + 1 INDX = IND(K) IF (X(INDX) < T) GO TO 8 * if K .LE. L, interchange elements K and L IF (K .LE. L) GO TO 6 * save the upper and lower subscripts of the portion of the * array yet to be sorted IF (L-I .LE. J-K) GO TO 9 IL(M) = I IU(M) = L I = K M = M + 1 GO TO 11 9 IL(M) = K IU(M) = J J = L M = M + 1 GO TO 11 * begin again on another unsorted portion of the array 10 M = M - 1 IF (M .EQ. 0) RETURN I = IL(M) J = IU(M) 11 IF (J-I .GE. 11) GO TO 4 IF (I .EQ. 1) GO TO 2 I = I - 1 * sort elements I+1,...,J. Note that 1 .LE. I < J and J-I < 11. 12 I = I + 1 IF (I .EQ. J) GO TO 10 INDX = IND(I+1) T = X(INDX) IT = INDX INDX = IND(I) IF (X(INDX) .LE. T) GO TO 12 K = I 13 IND(K+1) = IND(K) K = K - 1 INDX = IND(K) IF (T < X(INDX)) GO TO 13 IND(K+1) = IT GO TO 12 END ! of QSORTR ************************************************************ * * Robert Renka * Oak Ridge Natl. Lab. * * This routine applies a set of permutations to a vector. * * INPUT PARAMETERS - N - Length of A and IP. * * IP - Vector containing the sequence of * Integers 1,...,N permuted in the * same fashion that A is to be permuted * * A - Vector to be permuted. * * N and IP are not altered by this routine. * * OUTPUT PARAMETER - A - Reordered vector reflecting the * permutations defined by IP. * ************************************************************ SUBROUTINE IORDER (A, IP, N) INTEGER N, A(N), IP(N) INTEGER NN, K, J, IPJ, TEMP * LOCAL PARAMETERS - * * NN = Local copy of N * K = Index for IP and for the first element of A[] in a permutation * J = Index for IP and A, J .GE. K * IPJ = IP[J] * TEMP = Temporary storage for A[K] NN = N IF (NN < 2) RETURN K = 1 * Loop on permutations 1 J = K TEMP = A(K) * Apply permutation to A. IP[J] is marked (made negative) * as being included in the permutation. 2 IPJ = IP(J) IP(J) = -IPJ IF (IPJ .EQ. K) GO TO 3 A(J) = A(IPJ) J = IPJ GO TO 2 3 A(J) = TEMP * Search for an unmarked element of IP 4 K = K + 1 IF (K > NN) GO TO 5 IF (IP(K) > 0) GO TO 1 GO TO 4 * All permutations have been applied. Unmark IP. 5 DO 6 K = 1,NN 6 IP(K) = -IP(K) RETURN END ! of IORDER ************************************************************ SUBROUTINE RORDER (A, IP, N) INTEGER N, IP(N) REAL A(N) INTEGER NN, K, J, IPJ REAL TEMP * LOCAL PARAMETERS - * NN = Local copy of N * K = Index for IP and for the first element of A in a permutation * J = Index for IP and A, J .GE. K * IPJ = IP(J) * TEMP = Temporary storage for A(K) ************************************************************************ NN = N IF (NN < 2) RETURN K = 1 * Loop on permutations 1 J = K TEMP = A(K) * Apply permutation to A. IP(J) is marked (made negative) * as being included in the permutation. 2 IPJ = IP(J) IP(J) = -IPJ IF (IPJ .EQ. K) GO TO 3 A(J) = A(IPJ) J = IPJ GO TO 2 3 A(J) = TEMP * Search for an unmarked element of IP 4 K = K + 1 IF (K > NN) GO TO 5 IF (IP(K) > 0) GO TO 1 GO TO 4 * All permutations have been applied. Unmark IP. 5 DO 6 K = 1,NN 6 IP(K) = -IP(K) RETURN END ************************************************************************ SUBROUTINE DORDER (A, IP, N) INTEGER N, IP(N) DOUBLE PRECISION A(N) INTEGER NN, K, J, IPJ DOUBLE PRECISION TEMP NN = N IF (NN .LT. 2) RETURN K = 1 * Loop on permutations 1 J = K TEMP = A(K) * Apply permutation to A. IP(J) is marked (made negative) * as being included in the permutation. 2 IPJ = IP(J) IP(J) = -IPJ IF (IPJ .EQ. K) GO TO 3 A(J) = A(IPJ) J = IPJ GO TO 2 3 A(J) = TEMP * Search for an unmarked element of IP 4 K = K + 1 IF (K .GT. NN) GO TO 5 IF (IP(K) .GT. 0) GO TO 1 GO TO 4 * All permutations have been applied. Unmark IP. 5 DO 6 K = 1,NN 6 IP(K) = -IP(K) RETURN END *----------------------------------------------------------------------- * ISORT - sort integers *----------------------------------------------------------------------- * Sort an array and optionally make the same interchanges in * an auxiliary array. The array may be sorted in increasing * or decreasing order. A slightly modified QUICKSORT * algorithm is used. From the SLATEC library, public domain. * * by Jones, R. E., (SNLA) * Kahaner, D. K., (NBS) * Wisniewski, J. A., (SNLA) * * Description of Parameters * IX - integer array of values to be sorted * N - number of values in integer array IX to be sorted * * References: R. C. Singleton, Algorithm 347, An efficient algorithm * for sorting with minimal storage, Communications of * the ACM, 12, 3 (1969), pp. 185-187. * * REVISION HISTORY (YYMMDD) * 761118 Date written * 810801 Modified by David K. Kahaner. * 890531 Changed all specific intrinsics to generic. (WRB) * 890831 Modified array declarations. (WRB) * 891009 Removed unreferenced statement labels. (WRB) * 891009 REVISION DATE from Version 3.2 * 891214 Prologue converted to Version 4.0 format. (BAB) * 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) * 901012 Declared all variables; changed X,Y to IX,IY. (M. McClain) * 920501 Reformatted the REFERENCES section. (DWL, WRB) * 920519 Clarified error messages. (DWL) * 920801 Declarations section rebuilt and code restructured to use * IF-THEN-ELSE-ENDIF. (RWC, WRB) * 111005 Removed call to XERMSG, simplified. (AJA) *----------------------------------------------------------------------- SUBROUTINE ISORT (IX, N) * .. Scalar Arguments .. INTEGER N * .. Array Arguments .. INTEGER IX(N) * .. Local Scalars .. REAL R INTEGER I, IJ, J, K, L, M, NN, T, TT * .. Local Arrays .. INTEGER IL(32), IU(32) * Begin NN = N IF (NN < 1) THEN PRINT *, 'The number of values to be sorted is not positive.' RETURN END IF * Sort IX only M = 1 I = 1 J = NN R = 0.375E0 20 IF (I .EQ. J) GO TO 60 IF (R .LE. 0.5898437E0) THEN R = R+3.90625E-2 ELSE R = R-0.21875E0 ENDIF 30 K = I * Select a central element of the array and save it in location T IJ = I + INT(FLOAT(J-I)*R) T = IX(IJ) * If first element of array is greater than T, interchange with T IF (IX(I) > T) THEN IX(IJ) = IX(I) IX(I) = T T = IX(IJ) ENDIF L = J * If last element of array is less than than T, interchange with T IF (IX(J) < T) THEN IX(IJ) = IX(J) IX(J) = T T = IX(IJ) * If first element of array is greater than T, interchange with T IF (IX(I) > T) THEN IX(IJ) = IX(I) IX(I) = T T = IX(IJ) ENDIF ENDIF * Find an element in the second half of the array which is smaller than T 40 L = L-1 IF (IX(L) > T) GO TO 40 * Find an element in the first half of the array which is greater than T 50 K = K+1 IF (IX(K) < T) GO TO 50 * Interchange these elements IF (K .LE. L) THEN TT = IX(L) IX(L) = IX(K) IX(K) = TT GO TO 40 ENDIF * Save upper and lower subscripts of the array yet to be sorted IF (L-I > 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 * Begin again on another portion of the unsorted array 60 M = M-1 IF (M .EQ. 0) GO TO 200 I = IL(M) J = IU(M) 70 IF (J-I .GE. 1) GO TO 30 IF (I .EQ. 1) GO TO 20 I = I-1 80 I = I+1 IF (I .EQ. J) GO TO 60 T = IX(I+1) IF (IX(I) .LE. T) GO TO 80 K = I 90 IX(K+1) = IX(K) K = K-1 IF (T < IX(K)) GO TO 90 IX(K+1) = T GO TO 80 200 END ! of ISORT *----------------------------------------------------------------------- * SSORT - sort single precision numbers *----------------------------------------------------------------------- * Sort an array in increasing order. A slightly modified QUICKSORT * algorithm is used. From the SLATEC library, public domain. * * By R.E. Jones, (SNLA) * and J. A. Wisniewski, (SNLA) * * Description of arguments * X[] array of values to be sorted * N number of values in array X to be sorted * * REFERENCES R. C. Singleton, Algorithm 347, An efficient algorithm * for sorting with minimal storage, Communications of * the ACM, 12, 3 (1969), pp. 185-187. * * REVISION HISTORY (YYMMDD) * 761101 Date written * 761118 Modified to use the Singleton quicksort algorithm. (JAW) * 890531 Changed all specific intrinsics to generic. (WRB) * 890831 Modified array declarations. (WRB) * 891009 Removed unreferenced statement labels. (WRB) * 891024 Changed category. (WRB) * 891024 REVISION DATE from Version 3.2 * 891214 Prologue converted to Version 4.0 format. (BAB) * 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) * 901012 Declared all variables; changed X,Y to SX,SY. (M. McClain) * 920501 Reformatted the REFERENCES section. (DWL, WRB) * 920519 Clarified error messages. (DWL) * 920801 Declarations section rebuilt and code restructured to use * IF-THEN-ELSE-ENDIF. (RWC, WRB) * 110517 removed call to XERMSG, simplified (AJA) *----------------------------------------------------------------------- SUBROUTINE SSORT (X, N) * .. Scalar Arguments .. INTEGER N * .. Array Arguments .. REAL X(N) * .. Local Scalars .. REAL R, T, TT INTEGER I, IJ, J, K, L, M, NN * .. Local Arrays .. INTEGER IL(32), IU(32) *----------------------------------------------------------------------- NN = N IF (NN < 1) THEN PRINT *, 'The number of values to be sorted is not positive.' RETURN ENDIF * Sort X only M = 1 I = 1 J = NN R = 0.375E0 20 IF (I .EQ. J) GO TO 60 IF (R .LE. 0.5898437E0) THEN R = R+3.90625E-2 ELSE R = R-0.21875E0 ENDIF 30 K = I * Select a central element of the array and save it in location T IJ = I + INT(FLOAT(J-I)*R) T = X(IJ) * If first element of array is greater than T, interchange with T IF (X(I) > T) THEN X(IJ) = X(I) X(I) = T T = X(IJ) ENDIF L = J * If last element of array is less than than T, interchange with T IF (X(J) < T) THEN X(IJ) = X(J) X(J) = T T = X(IJ) * If first element of array is greater than T, interchange with T IF (X(I) > T) THEN X(IJ) = X(I) X(I) = T T = X(IJ) ENDIF ENDIF * Find an element in the second half of the array which is smaller than T 40 L = L-1 IF (X(L) > T) GO TO 40 * Find an element in the first half of the array which is greater than T 50 K = K+1 IF (X(K) < T) GO TO 50 * Interchange these elements IF (K .LE. L) THEN TT = X(L) X(L) = X(K) X(K) = TT GO TO 40 ENDIF * Save upper and lower subscripts of the array yet to be sorted IF (L-I > 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 * Begin again on another portion of the unsorted array 60 M = M-1 IF (M .EQ. 0) GO TO 200 I = IL(M) J = IU(M) 70 IF (J-I .GE. 1) GO TO 30 IF (I .EQ. 1) GO TO 20 I = I-1 80 I = I+1 IF (I .EQ. J) GO TO 60 T = X(I+1) IF (X(I) .LE. T) GO TO 80 K = I 90 X(K+1) = X(K) K = K-1 IF (T < X(K)) GO TO 90 X(K+1) = T GO TO 80 200 END ! of SSORT *----------------------------------------------------------------------- * rsorti - radix sort and carry along an integer array * by Andy Allinger, 2011-2016, released to the public domain * This program may be used by any person for any purpose * * origin: Herman Hollerith, 1887 * * NOTE: This routine is a sort-and-carry. To obtain a permutation * vector to apply to other arrays, you must first initialize IND to * 1, 2, 3...N * * NOTE ALSO: * if B is positive, then the integers will be interpreted as UNSIGNED * note that FORTRAN does not support unsigned integers, so make sure * you know what you are doing if you use B=32 * if B is negative, the sign bit will be inspected as a final step * So, for the range -256...255 set B = -8 * *__Name____Type_____In/Out_________Description__________________________ * D[N] Integer Both Array, returned sorted * IND[N] Integer Both Array to carry along, meaning * permuted in the same fashion as D * N Integer In Length of the array * B Integer In Number of significant bits * W[2@N] Integer Neither Workspace *----------------------------------------------------------------------- SUBROUTINE RSORTI (D, IND, N, B, W) ! partition workspace IMPLICIT NONE INTEGER D, IND, N, B, W DIMENSION D(N), IND(N), W(2*N) IF (N < 2) RETURN CALL RSI01 (D, IND, N, B, W(1), W(N+1)) RETURN END ! of RSORTI *----------------------------------------------------------------------- SUBROUTINE RSI01 (D, D2, N, B, W, W2) IMPLICIT NONE INTEGER D, D2, N, B, W, W2 DIMENSION D(N), D2(N), W(N), W2(N) INTEGER I, J, K, ILIM, P1OLD, P0OLD, P1, P0, S LOGICAL INW ! data is in W INTEGER BITSZ ! external function * validate IF (B .EQ. 0) STOP 'Bit size must not be zero' IF (ABS(B) > BITSZ(N)) STOP 'Bit size to large' *----------------------------------------------------------------------- * sort starts on the least significant bit, position 0 *----------------------------------------------------------------------- P1 = N+1 P0 = 0 DO J = 1, N * test for negative values IF (D(J) < 0 .AND. B > 0 .AND. B .NE. BITSZ(N)) & STOP 'Set argument B negative to sort negative numbers' IF ( BTEST(D(J), 0) ) THEN P1 = P1 - 1 W(P1) = D(J) W2(P1) = D2(J) ELSE P0 = P0 + 1 W(P0) = D(J) W2(P0) = D2(J) END IF END DO INW = .TRUE. *----------------------------------------------------------------------- * now alternate between putting data into D and putting data into W *----------------------------------------------------------------------- ILIM = 0 IF (B < 0) ILIM = MIN( ABS(B)-1, BITSZ(I)-2 ) IF (B > 0) ILIM = B - 1 DO I = 1, ILIM P1OLD = P1 P0OLD = P0 ! save the value from previous bit P1 = N+1 P0 = 0 ! start a fresh count for next bit * odd I IF (INW) THEN * copy data from the zeros DO J = 1, P0OLD, 1 IF ( BTEST(W(J), I) ) THEN P1 = P1 - 1 D(P1) = W(J) D2(P1) = W2(J) ELSE P0 = P0 + 1 D(P0) = W(J) D2(P0) = W2(J) END IF END DO * copy data from the ones DO J = N, P1OLD, -1 IF ( BTEST(W(J), I) ) THEN P1 = P1 - 1 D(P1) = W(J) D2(P1) = W2(J) ELSE P0 = P0 + 1 D(P0) = W(J) D2(P0) = W2(J) END IF END DO * even I ELSE * copy data from the zeros DO J = 1, P0OLD, 1 IF ( BTEST(D(J), I) ) THEN P1 = P1 - 1 W(P1) = D(J) W2(P1) = D2(J) ELSE P0 = P0 + 1 W(P0) = D(J) W2(P0) = D2(J) END IF END DO * copy data from the ones DO J = N, P1OLD, -1 IF ( BTEST(D(J), I) ) THEN P1 = P1 - 1 W(P1) = D(J) W2(P1) = D2(J) ELSE P0 = P0 + 1 W(P0) = D(J) W2(P0) = D2(J) END IF END DO END IF ! even or odd i INW = .NOT. INW END DO ! next i *----------------------------------------------------------------------- * the sign bit *----------------------------------------------------------------------- IF (B < 0) THEN P1OLD = P1 P0OLD = P0 P1 = N+1 P0 = 0 * if sign bit is set, send to the zero end of the work array IF (INW) THEN DO J = 1, P0OLD, 1 IF ( BTEST(W(J), BITSZ(I)-1) ) THEN P0 = P0 + 1 D(P0) = W(J) D2(P0) = W2(J) ELSE P1 = P1 - 1 D(P1) = W(J) D2(P1) = W2(J) END IF END DO DO J = N, P1OLD, -1 IF ( BTEST(W(J), BITSZ(I)-1) ) THEN P0 = P0 + 1 D(P0) = W(J) D2(P0) = W2(J) ELSE P1 = P1 - 1 D(P1) = W(J) D2(P1) = W2(J) END IF END DO ELSE DO J = 1, P0OLD, 1 IF ( BTEST(D(J), BITSZ(I)-1) ) THEN P0 = P0 + 1 W(P0) = D(J) W2(P0) = D2(J) ELSE P1 = P1 - 1 W(P1) = D(J) W2(P1) = D2(J) END IF END DO DO J = N, P1OLD, -1 IF ( BTEST(D(J), BITSZ(I)-1) ) THEN P0 = P0 + 1 W(P0) = D(J) W2(P0) = D2(J) ELSE P1 = P1 - 1 W(P1) = D(J) W2(P1) = D2(J) END IF END DO END IF INW = .NOT. INW END IF *----------------------------------------------------------------------- * put the data back into the input array *----------------------------------------------------------------------- K = P1 IF (INW) THEN DO J = 1, P0, 1 D(J) = W(J) D2(J) = W2(J) END DO DO J = N, P1, -1 D(K) = W(J) D2(K) = W2(J) K = K + 1 END DO * Swap values in the ones portion ELSE DO J = N, (N+P1+2)/2, -1 S = D(J) D(J) = D(K) D(K) = S S = D2(J) D2(J) = D2(K) D2(K) = S K = K + 1 END DO END IF RETURN END ! of RSI01 *********************************************************************** * msortr - merge sort on real numbers * by Andy Allinger, 2012-2016, released to the public domain * This program may be used by any person for any purpose * * origin: John von Neumann, 1945 * * NOTE: This routine is a sort-and-carry. To obtain a permutation * vector to apply to other arrays, you must first initialize INDA to * 1, 2, 3...N * *__Name_____Type_________I/O______Description_______________________ * A[N] Real Both The array, returned sorted * INDA[N] Integer Both Integer array to carry along, meaning * permuted in the same fashion as A * N Integer In Length of array * RW[3@N] Real Neither Workspace * IW[3@N] Integer Neither Workspace *********************************************************************** SUBROUTINE MSORTR (A, INDA, N, RW, IW) ! partition workspace IMPLICIT NONE INTEGER INDA, N, IW REAL A, RW DIMENSION A(N), INDA(N), RW(3*N), IW(3*N) IF (N < 2) RETURN CALL MSR01 (A, INDA, N, RW(1), RW(N+1), RW(2*N+1), & IW(1), IW(N+1), IW(2*N+1)) RETURN END ! of MSORTR *********************************************************************** SUBROUTINE MSR01 (A, INDA, N, B, C, D, INDB, INDC, INDD) IMPLICIT NONE INTEGER INDA, N, INDB, INDC, INDD REAL A, B, C, D DIMENSION A(N), INDA(N), B(N), C(N), D(N), & INDB(N), INDC(N), INDD(N) * local variables INTEGER g, h, hlim, iA, iB, iC, iD, ilimA, ilimB, ilimC, & ilimD, j, jA, jB, jC, jD, kA, kB, kC, kD, LA, LB, LC, LD, & m, o, o0, indx, indy INTEGER BITSZ ! external function REAL x, y * validate IF (n < 2) RETURN * first pass through data jC = 0 jD = 0 DO h = 2, n, 2 x = A(h-1) indx = indA(h-1) y = A(h) indy = indA(h) IF (BTEST(h, 1)) THEN IF (x > y) THEN jC = jC + 1 C(jC) = y indC(jC) = indy jC = jC + 1 C(jC) = x indC(jC) = indx ELSE jC = jC + 1 C(jC) = x indC(jC) = indx jC = jC + 1 C(jC) = y indC(jC) = indy END IF ELSE IF (x > y) THEN jD = jD + 1 D(jD) = y indD(jD) = indy jD = jD + 1 D(jD) = x indD(jD) = indx ELSE jD = jD + 1 D(jD) = x indD(jD) = indx jD = jD + 1 D(jD) = y indD(jD) = indy END IF END IF END DO ! next h * odd n IF (BTEST(n, 0)) THEN IF (BTEST(n, 1)) THEN jD = jD + 1 D(jD) = A(n) indD(jD) = indA(n) ELSE jC = jC + 1 C(jC) = A(n) indC(jC) = indA(n) END IF END IF LC = jC LD = jD * number of passes through data M = BITSZ(N) DO WHILE(.NOT. BTEST(n-1, m-1)) m = m - 1 END DO *********************************************************************** * sort *********************************************************************** o = 2 LA = 0 LB = 0 ! avoid compiler warnings DO g = 2, m ! begin pass through data o0 = o o = o * 2 hlim = (n-1) / o +1 IF (BTEST(g, 0)) THEN kA = 0 kB = 0 jC = 0 jD = 0 DO h = 1, hlim ilimA = MIN(o0, LA - kA) ilimB = MIN(o0, LB - kB) iA = 0 iB = 0 *********************************************************************** IF (BTEST(h, 0)) THEN ! alternate between output arrays *********************************************************************** IF (iA .GE. ilimA) PRINT *, '!!empty buffer A iA=', iA iA = iA + 1 kA = kA + 1 x = A(kA) indx = indA(kA) IF (iB .GE. ilimB) THEN jC = jC + 1 C(jC) = x indC(jC) = indx GO TO 20 END IF iB = iB + 1 kB = kB + 1 y = B(kB) indy = indB(kB) *********************************************************************** * main comparison statement *********************************************************************** 10 IF (x > y) THEN ! send y first jC = jC + 1 C(jC) = y indC(jC) = indy IF (iB < ilimB) THEN ! take another from B[] iB = iB + 1 kB = kB + 1 y = B(kB) indy = indB(kB) GO TO 10 ELSE jC = jC + 1 C(jC) = x indC(jC) = indx END IF ELSE ! send x first jC = jC + 1 C(jC) = x indC(jC) = indx IF (iA < ilimA) THEN ! take another from A[] iA = iA + 1 kA = kA + 1 x = A(kA) indx = indA(kA) GO TO 10 ELSE jC = jC + 1 C(jC) = y indC(jC) = indy END IF END IF *********************************************************************** * one side of the input is bound to become empty before the other *********************************************************************** 20 DO WHILE (iA < ilimA) iA = iA + 1 kA = kA + 1 x = A(kA) indx = indA(kA) jC = jC + 1 C(jC) = x indC(jC) = indx END DO DO WHILE (iB < ilimB) iB = iB + 1 kB = kB + 1 y = B(kB) indy = indB(kB) jC = jC + 1 C(jC) = y indC(jC) = indy END DO *********************************************************************** ELSE ! even h *********************************************************************** IF (iA .GE. ilimA) PRINT *, 'empty buffer A iA=', iA iA = iA + 1 kA = kA + 1 x = A(kA) indx = indA(kA) IF (iB .GE. ilimB) THEN jD = jD + 1 D(jD) = x indD(jD) = indx GO TO 40 END IF iB = iB + 1 kB = kB + 1 y = B(kB) indy = indB(kB) *********************************************************************** * main comparison statement *********************************************************************** 30 IF (x > y) THEN jD = jD + 1 D(jD) = y indD(jD) = indy IF (iB < ilimB) THEN iB = iB + 1 kB = kB + 1 y = B(kB) indy = indB(kB) GO TO 30 ELSE jD = jD + 1 D(jD) = x indD(jD) = indx END IF ELSE jD = jD + 1 D(jD) = x indD(jD) = indx IF (iA < ilimA) THEN iA = iA + 1 kA = kA + 1 x = A(kA) indx = indA(kA) GO TO 30 ELSE jD = jD + 1 D(jD) = y indD(jD) = indy END IF END IF *********************************************************************** * one side of the input is bound to become empty before the other *********************************************************************** 40 DO WHILE (iA < ilimA) iA = iA + 1 kA = kA + 1 x = A(kA) indx = indA(kA) jD = jD + 1 D(jD) = x indD(jD) = indx END DO DO WHILE (iB < ilimB) iB = iB + 1 kB = kB + 1 y = B(kB) indy = indB(kB) jD = jD + 1 D(jD) = y indD(jD) = indy END DO END IF ! even or odd h END DO ! next h * remember lengths of arrays LC = jC LD = jD *********************************************************************** ELSE ! even g *********************************************************************** kC = 0 kD = 0 jA = 0 jB = 0 DO h = 1, hlim ilimC = MIN(o0, LC - kC) ilimD = MIN(o0, LD - kD) iC = 0 iD = 0 IF (BTEST(h, 0)) THEN ! alternate between output arrays IF (iC .GE. ilimC) PRINT *, '!!empty buffer C iC=', iC iC = iC + 1 kC = kC + 1 x = C(kC) indx = indC(kC) IF (iD .GE. ilimD) THEN jA = jA + 1 A(jA) = x indA(jA) = indx GO TO 60 END IF iD = iD + 1 kD = kD + 1 y = D(kD) indy = indD(kD) *********************************************************************** * main comparison statement *********************************************************************** 50 IF (x > y) THEN jA = jA + 1 A(jA) = y indA(jA) = indy IF (iD < ilimD) THEN iD = iD + 1 kD = kD + 1 y = D(kD) indy = indD(kD) GO TO 50 ELSE jA = jA + 1 A(jA) = x indA(jA) = indx END IF ELSE jA = jA + 1 A(jA) = x indA(jA) = indx IF (iC < ilimC) THEN iC = iC + 1 kC = kC + 1 x = C(kC) indx = indC(kC) GO TO 50 ELSE jA = jA + 1 A(jA) = y indA(jA) = indy END IF END IF *********************************************************************** * one side of the input is bound to become empty before the other *********************************************************************** 60 DO WHILE (iC < ilimC) iC = iC + 1 kC = kC + 1 x = C(kC) indx = indC(kC) jA = jA + 1 A(jA) = x indA(jA) = indx END DO DO WHILE (iD < ilimD) iD = iD + 1 kD = kD + 1 y = D(kD) indy = indD(kD) jA = jA + 1 A(jA) = y indA(jA) = indy END DO *********************************************************************** ELSE ! even h *********************************************************************** IF (iC .GE. ilimC) PRINT *, '!!empty buffer C iC=', iC iC = iC + 1 kC = kC + 1 x = C(kC) indx = indC(kC) IF (iD .GE. ilimD) THEN jB = jB + 1 B(jB) = x indB(jB) = indx GO TO 80 END IF iD = iD + 1 kD = kD + 1 y = D(kD) indy = indD(kD) *********************************************************************** * main comparison statement *********************************************************************** 70 IF (x > y) THEN jB = jB + 1 B(jB) = y indB(jB) = indy IF (iD < ilimD) THEN iD = iD + 1 kD = kD + 1 y = D(kD) indy = indD(kD) GO TO 70 ELSE jB = jB + 1 B(jB) = x indB(jB) = indx END IF ELSE jB = jB + 1 B(jB) = x indB(jB) = indx IF (iC < ilimC) THEN iC = iC + 1 kC = kC + 1 x = C(kC) indx = indC(kC) GO TO 70 ELSE jB = jB + 1 B(jB) = y indB(jB) = indy END IF END IF *********************************************************************** * one side of the input is bound to become empty before the other *********************************************************************** 80 DO WHILE (iC < ilimC) iC = iC + 1 kC = kC + 1 x = C(kC) indx = indC(kC) jB = jB + 1 B(jB) = x indB(jB) = indx END DO DO WHILE (iD < ilimD) iD = iD + 1 kD = kD + 1 y = D(kD) indy = indD(kD) jB = jB + 1 B(jB) = y indB(jB) = indy END DO END IF ! even or odd h END DO ! next h * remember lengths of arrays LA = jA LB = jB END IF ! even or odd g END DO ! next g * move output back into A IF (BTEST(m, 0)) THEN DO j = 1, n A(j) = C(j) indA(j) = indC(j) END DO END IF RETURN END ! of MSR01 *----------------------------------------------------------------------- * bisection search for the first element greater than or equal to a value * a fragment of SEVAL, by Rondall Jones, Sandia National Laboratory * for unlimited distribution *----------------------------------------------------------------------- * Name Type In/Out Description * X[N] REAL in array (in increasing order) * N INTEGER in length of X * V REAL in the value sought * I INTEGER out index where X .GE. V found *----------------------------------------------------------------------- SUBROUTINE RBSGE(X, N, V, I) IMPLICIT NONE INTEGER N, I REAL X, V DIMENSION X(N) INTEGER IL, IR * check input IF (N < 1) RETURN * V is greater than X(N) or less than X(1) IF (V < X(1)) THEN I = 1 RETURN ELSE IF (V > X(N)) THEN I = N RETURN END IF * bisection search IL = 1 IR = N 10 I = (IL + IR) / 2 IF (I .EQ. IL) THEN IF (X(I) .GE. V) THEN I = IL ELSE I = IR END IF RETURN END IF IF (X(I) < V) THEN IL = I ELSE IR = I END IF GO TO 10 END ! of RBSGE *----------------------------------------------------------------------- * bisection search an array for the last element less than or equal to a value * a fragment of SEVAL, by Rondall Jones, Sandia National Laboratory * for unlimited distribution * * note that type is changed to INTEGER and .GE. is changed to .LE. *----------------------------------------------------------------------- * Name Type In/Out Description * J[N] INTEGER in array (in increasing order) * N INTEGER in length of J * V INTEGER in the value sought * I INTEGER out index where J .LE. V found * returns 0 when J[1] > V *----------------------------------------------------------------------- SUBROUTINE IBSLE(J, N, V, I) IMPLICIT NONE INTEGER J, N, V, I DIMENSION J(N) INTEGER IL, IR * check input IF (N < 1) RETURN * V is greater than J[N] or less than J[1] IF (V < J(1)) THEN I = 0 RETURN ELSE IF (V > J(N)) THEN I = N RETURN END IF * bisection search IL = 1 IR = N 10 I = (IL + IR) / 2 IF (I .EQ. IL) THEN IF (J(IR) .LE. V) THEN I = IR ELSE I = IL END IF RETURN END IF IF (J(I) .LE. V) THEN IL = I ELSE IR = I END IF GO TO 10 END ! of IBSLE *----------------------------------------------------------------------- * Bisection search for the first element greater than a value * *___Name_______Type_____In/Out_____Description__________________________ * X(N) REAL In Array (in increasing order) * N INTEGER In Length of X * V REAL In The value sought * I INTEGER Out Index where X > V found *----------------------------------------------------------------------- SUBROUTINE RBSGT (X, N, V, I) IMPLICIT NONE INTEGER N, I REAL X, V DIMENSION X(N) INTEGER IL, IR * check input IF (N < 1) RETURN * V is greater than X(N) or less than X(1) IF (V < X(1)) THEN I = 1 RETURN ELSE IF (V > X(N)) THEN I = N RETURN END IF * bisection search IL = 1 IR = N 10 I = (IL + IR) / 2 IF (I .EQ. IL) THEN IF (X(I) > V) THEN I = IL ELSE I = IR END IF RETURN END IF IF (X(I) .LE. V) THEN IL = I ELSE IR = I END IF GO TO 10 END ! of RBSGT