PROGRAM TRACE_BB
      IMPLICIT INTEGER(A-Z)
      DOUBLE PRECISION MATRIX(20,20), STAT,ACTUAL, PVALUE, TIMEA, TIMEB,
     1                 NFACTORIAL, PREstat, LowerBound, UpperBound,
     1                 largest, smallest, F,INDEX
      INTEGER PERMUTATION(20), POSITION
      LOGICAL FATHOM, RED, FOUND
C
C  BRANCH-AND-BOUND EXACT TEST FOR THE TRACE STATISTIC
C
      OPEN(1, FILE='MATCH.DAT')
      OPEN(3, FILE='MATCH.OUT')
C     Determine the number of objects (N)
      READ(1,*) N
      DO I = 1,N
        READ(1,*) (MATRIX(I,J),J=1,N)
      END DO
C     Timing is everything!
      CALL GETTIM (IHR, IMIN, ISEC, I100)
      CALL GETDAT (IYR, IMON, IDAY)
      TIMEA=DFLOAT(86400*IDAY+3600*IHR+60*IMIN+ISEC)+DFLOAT(I100)/100.

C     Evaluate the actual statistic and store in a variable (Observed)
	ACTUAL = 0
	DO 10 POSITION = 1, N
		ACTUAL = ACTUAL + MATRIX(POSITION, POSITION)
 10	CONTINUE

C	Initialize an "empty" n × 1 array as zeros (permutation)
	DO 20 POSITION = 1, N
		PERMUTATION(POSITION) = 0		
 20	CONTINUE
	

C	Initialize a counter as zero and a pointer to the location in the permutation,
      INDEX = 0
	POSITION = 0

C	Start building and pruning sequences/permutations
 1      DO WHILE ((POSITION.NE.1) .OR. (PERMUTATION(POSITION).LE.N))   
C         Forward Branch
          POSITION = POSITION + 1
          FATHOM = .FALSE.
          DO WHILE (FATHOM.EQV..FALSE.)
C           Right Branch by loading the next available object, into
C           permutation(Position).  If/when permutation(Position) > n,
C           all objects have been considered.
C           NOTE: RED is short for redundancy.
            RED = .TRUE.
            DO WHILE ((RED.EQV..TRUE.).AND.(PERMUTATION(POSITION).LE.N))
              PERMUTATION(POSITION) = PERMUTATION(POSITION) + 1
              RED = .FALSE.
              DO 30 I = 1,(POSITION - 1)
                IF (PERMUTATION(POSITION).EQ.PERMUTATION(I)) THEN
                  RED = .TRUE.
                END IF
 30           CONTINUE
            END DO 
C           Termination
            IF((POSITION.EQ.1).AND.(PERMUTATION(POSITION).GT.N)) GO TO 1
C           Retraction        
            IF((POSITION.GT.1).AND.(PERMUTATION(POSITION).GT.N)) THEN
              PERMUTATION(POSITION) = 0
              POSITION = POSITION - 1
            ELSEIF (POSITION.EQ.N) THEN
C           Complete sequence is ready for evaluation.
C           Evaluate the statistic for the current permutation as a variable, STAT.
              STAT = 0
              DO 40 I = 1,N
                STAT = STAT + MATRIX(I, PERMUTATION(I))
 40           CONTINUE
              IF (STAT.GE.ACTUAL) INDEX = INDEX + 1
            ELSE
              PREstat = 0
              LowerBound = 0
              UpperBound = 0
              DO 60 I = 1, N
                FOUND = .FALSE.
                DO 80 J = 1,POSITION
                  IF(PERMUTATION(J).EQ.I) THEN
                    FOUND = .TRUE.
                    PREstat = PREstat + MATRIX(J, PERMUTATION(J))
                  ENDIF
 80             CONTINUE
                IF(FOUND.EQV..FALSE.) THEN
                  largest = MATRIX(POSITION + 1, I)
                  smallest = MATRIX(POSITION + 1, I)
                  DO K = Position + 1, N
                    IF(largest.LT.MATRIX(K, I)) largest = MATRIX(K, I)
                    IF(smallest.GT.MATRIX(K, I)) smallest = MATRIX(K, I)
                  END DO
                  LowerBound = LowerBound + smallest
                  UpperBound = UpperBound + largest
                END IF
 60           CONTINUE                    
              FATHOM = .TRUE.
              IF((PREstat + LowerBound).GE.ACTUAL) THEN
                FATHOM = .FALSE.
                F = 1
                DO 70 I = 2, (N - POSITION)
                  F = F*I
 70             CONTINUE
                INDEX = INDEX + F
              ELSEIF((PREstat + UpperBound).LT.ACTUAL) THEN
                FATHOM = .FALSE.
              END IF
            END IF
          END DO
        END DO

C	Return p = (Index/factorial(n)) as the p-value.
        NFACTORIAL = 1
        DO 50 I = 2, N
                NFACTORIAL = NFACTORIAL * I
 50	CONTINUE
        PVALUE= INDEX/NFACTORIAL

C	Write it all out
C	Syntax:
C		      write(*, label) list-of-variables
C		label format format-code
      WRITE(3,51) NFACTORIAL, INDEX
      WRITE(3,52) ACTUAL
      WRITE(3,53) PVALUE
      CALL GETTIM (IHR, IMIN, ISEC, I100)
      CALL GETDAT (IYR, IMON, IDAY)
      TIMEB=DFLOAT(86400*IDAY+3600*IHR+60*IMIN+ISEC)+DFLOAT(I100)/100.
      WRITE(3,54) TIMEB-TIMEA
 51     FORMAT(' PERMUTATIONS ',F13.0,' ABOVE OR EQUAL ',F13.0)
 52	FORMAT(' MATCH STAT ',F14.5)
 53	FORMAT(' ONE TAIL P VALUE ',F14.11)
 54	FORMAT(' CPU TIME = ',F14.2)
      END