************************************************************
*
* 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
*++++++++++++++++++++++++ End of file sort.f +++++++++++++++++++++++++++