*----------------------------------------------------------------------- * rsorti - passive radix sort * by Andy Allinger, 2011-2019, released to the public domain * This program may be used by any person for any purpose * * origin: Herman Hollerith, 1887 * * NOTE: * 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___________________________ * IX(N) Integer Both Array to sort * IND(N) Integer Out Permutation vector * N Integer In Length of the array * B Integer In Number of significant bits * IW(N) Integer Neither Workspace *----------------------------------------------------------------------- SUBROUTINE RSORTI (IX, IND, N, B, IW) IMPLICIT NONE INTEGER IX, IND, N, B, IW DIMENSION IX(N), IND(N), IW(N) INTEGER BS, I, J, K, ILIM, P1OLD, P0OLD, P1, P0 LOGICAL ODD INTEGER BITSZ ! external function * Validate BS = BITSZ() IF (B .EQ. 0) THEN PRINT *, 'rsorti: Bit size must not be zero' RETURN END IF IF (ABS(B) > BS) THEN PRINT *, 'rsorti: Bit size too large' RETURN END IF DO J = 1, N ! Test for negative values. IF (IX(J) < 0 .AND. B > 0 .AND. B .NE. BS) THEN PRINT *, 'Set argument B negative to sort negative numbers' RETURN END IF END DO DO J = 1, N ! initialize IND(J) = J END DO IF (N < 2) RETURN *----------------------------------------------------------------------- * Alternate between referencing by IND and referencing by IW *----------------------------------------------------------------------- ILIM = 0 IF (B < 0) ILIM = MIN( ABS(B)-1, BS-2 ) ! -31 and -32 treated the same IF (B > 0) ILIM = B - 1 P1 = N+1 P0 = N ODD = .FALSE. DO I = 0, ILIM P1OLD = P1 P0OLD = P0 ! save the value from previous bit P1 = N+1 P0 = 0 ! start a fresh count for next bit IF (ODD) THEN DO J = 1, P0OLD, 1 ! copy data from the zeros IF ( BTEST(IX(IW(J)), I) ) THEN P1 = P1 - 1 IND(P1) = IW(J) ELSE P0 = P0 + 1 IND(P0) = IW(J) END IF END DO DO J = N, P1OLD, -1 ! copy data from the ones IF ( BTEST(IX(IW(J)), I) ) THEN P1 = P1 - 1 IND(P1) = IW(J) ELSE P0 = P0 + 1 IND(P0) = IW(J) END IF END DO ELSE DO J = 1, P0OLD, 1 ! copy data from the zeros IF ( BTEST(IX(IND(J)), I) ) THEN P1 = P1 - 1 IW(P1) = IND(J) ELSE P0 = P0 + 1 IW(P0) = IND(J) END IF END DO DO J = N, P1OLD, -1 ! copy data from the ones IF ( BTEST(IX(IND(J)), I) ) THEN P1 = P1 - 1 IW(P1) = IND(J) ELSE P0 = P0 + 1 IW(P0) = IND(J) END IF END DO END IF ! even or odd i ODD = .NOT. ODD END DO ! next i *----------------------------------------------------------------------- * the sign bit *----------------------------------------------------------------------- IF (B < 0) THEN P1OLD = P1 P0OLD = P0 P1 = N+1 P0 = 0 IF (ODD) THEN DO J = 1, P0OLD, 1 IF ( BTEST(IX(IW(J)), BS-1) ) THEN P0 = P0 + 1 IND(P0) = IW(J) ELSE P1 = P1 - 1 IND(P1) = IW(J) END IF END DO DO J = N, P1OLD, -1 IF ( BTEST(IX(IW(J)), BS-1) ) THEN P0 = P0 + 1 IND(P0) = IW(J) ELSE P1 = P1 - 1 IND(P1) = IW(J) END IF END DO ELSE DO J = 1, P0OLD, 1 IF ( BTEST(IX(IND(J)), BS-1) ) THEN P0 = P0 + 1 IW(P0) = IND(J) ELSE P1 = P1 - 1 IW(P1) = IND(J) END IF END DO DO J = N, P1OLD, -1 IF ( BTEST(IX(IND(J)), BS-1) ) THEN P0 = P0 + 1 IW(P0) = IND(J) ELSE P1 = P1 - 1 IW(P1) = IND(J) END IF END DO END IF ODD = .NOT. ODD END IF *----------------------------------------------------------------------- * Put indices into IND *----------------------------------------------------------------------- K = P1 IF (ODD) THEN DO J = 1, P0, 1 IND(J) = IW(J) END DO DO J = N, P1, -1 IND(K) = IW(J) K = K + 1 END DO * Swap values in the ones portion ELSE DO J = N, (N+P1+2)/2, -1 I = IND(J) IND(J) = IND(K) IND(K) = I K = K + 1 END DO END IF RETURN END ! of rsorti *----------------------------------------------------------------------- * BIT_SIZE replacement for old compilers. *----------------------------------------------------------------------- FUNCTION BITSZ () IMPLICIT NONE INTEGER BITSZ INTEGER I, J J = 0 J = NOT(J) BITSZ = 0 DO I = 1, 1024 J = ISHFT(J, -1) IF (.NOT. BTEST(J, 0)) THEN BITSZ = I RETURN END IF END DO RETURN ! error! END ! of bitsz