PROGRAM BB4
      IMPLICIT INTEGER(A-Z)
      DOUBLE PRECISION TIMEB,TIMEC,A(250,250),B(250,250),EPS,DIAMB,Z,
     1                 DBEST1,DBEST2,DMAX,MINK,KSUM(20),DIAMC(20),XX,
     1                 DTEST,BESDIAM(20)
      REAL S1
      INTEGER N(250),X(250),XB(250,20),XBB(150),
     1        IOPEN(250),NIC(250),MEMB(250,250),S(250),NP(250),
     1        XBT(250,20)
      OPEN(1,FILE='AMAT.DAT')
      OPEN(6,FILE='RESULTS')
      EPS=1.0D-07
C
C ################################################################
C     MINIMIZE THE MAXIMUM DIAMETER
C
C 1. RUN THE COMPLETE LINKAGE ALGORITHM (ONE REP)
C 2. USE A CONVERGENT ONE-OPT PROCEDURE TO IMPROVE DIAMETER
C 3. RUN THE BRANCH-AND-BOUND ALGORITHM
C
C ################################################################
C
      READ(1,*) E                        ! NUMBER OF OBJECTS
      WRITE(*,*) 'TYPE 1 FOR HALF MATRIX OR TYPE 2 FOR FULL MATRIX'
      READ(*,*) ITYPE
      IF(ITYPE.EQ.2) THEN
        READ(1,*) ((A(I,J),J=1,E),I=1,E)
      ELSE
        DO J = 2,E
          READ(1,*) (A(I,J),I=1,J-1)
        END DO
        DO J = 2,E
          DO I = 1,J-1
            A(J,I) = A(I,J)
          END DO
        END DO
      END IF
      WRITE(*,*) 'PLEASE INPUT NUMBER OF CLUSTERS 2 TO 20'
      READ(*,*) C
      CALL GETTIM (IHR, IMIN, ISEC, I100)
      CALL GETDAT (IYR, IMON, IDAY)  
      TIMEB = DFLOAT(86400*IDAY+3600*IHR+60*IMIN+ISEC)+DFLOAT(I100)/100.  
C
C  #################################################################
C    FIND THE UPPER BOUND USING BIASED-SAMPLING COMPLETE LINK
C    CLUSTER ANALYSIS
C  #################################################################
C
      DIAMB=999999.0D0
      DO 4000 IJKL = 1,100
      DO I = 1,150
        NIC(I)=0
        DO J = 1,150
          MEMB(I,J)=0
        END DO
      END DO
      NG=E
      ANG=E
      DO I = 1,NG
        NIC(I)=1
        MEMB(I,1)=I
        IOPEN(I)=1
      END DO
 5000 DBEST1=999999.0D0
      DBEST2=999999.0D0
      DO 1 I = 1,NG-1
        DO 2 J = I+1,NG
          IF(IOPEN(I).EQ.0.OR.IOPEN(J).EQ.0) GO TO 2
          DMAX=0.0D0
          DO I1 = 1,NIC(I)
            DO J1 = 1,NIC(J)
              M1=MEMB(I,I1)
              M2=MEMB(J,J1)
              IF(A(M1,M2).GT.DMAX) DMAX=A(M1,M2)
            END DO
          END DO
          IF(DMAX.LT.DBEST1-EPS) THEN
            ISET2=ISET1
            JSET2=JSET1
            DBEST2=DBEST1
            ISET1=I
            JSET1=J
            DBEST1=DMAX
          ELSEIF(DMAX.GT.DBEST1+EPS.AND.DMAX.LT.DBEST2-EPS) THEN
            ISET2=I
            JSET2=J
            DBEST2=DMAX
          END IF
   2    CONTINUE
   1  CONTINUE
      CALL RANDOM(S1)
      IF(S1.LT.0.5.OR.DBEST2.GT.99999.0D0) THEN
        ISET=ISET1
        JSET=JSET1
      ELSE
        ISET=ISET2
        JSET=JSET2
      END IF
      L=NIC(ISET)
      NIC(ISET)=NIC(ISET)+NIC(JSET)
      IOPEN(JSET)=0
      DO K = 1,NIC(JSET)
        L=L+1
        MEMB(ISET,L)=MEMB(JSET,K)
      END DO
      ANG=ANG-1
      IF(ANG.GT.C) GO TO 5000
C
      MCT=0
      NCT=0
      DMAX=0.0D0
      DO 3 I = 1,NG
        IF(IOPEN(I).EQ.0) GO TO 3
        MCT=MCT+1
        NCT=NCT+NIC(I)
        DO J = 1,NIC(I)-1
          DO K = J+1,NIC(I)
            M1=MEMB(I,J)
            M2=MEMB(I,K)
            IF(A(M1,M2).GT.DMAX) DMAX=A(M1,M2)
          END DO
        END DO
 3    CONTINUE
      IF(DMAX.LT.DIAMB) THEN
        DIAMB=DMAX
        DO I = 1,NG
          DO J = 1,C
            XB(I,J)=0
          END DO
        END DO
        NCT=0
        DO 47 I = 1,NG
          IF(IOPEN(I).EQ.0) GO TO 47
          NCT=NCT+1
          DO K = 1,NIC(I)
            KK=MEMB(I,K)
            XB(KK,NCT)=1
          END DO
  47    CONTINUE
      END IF
 4000 CONTINUE
c
 4500 ITRIG=0
      DO 70 I = 1,E                             ! I IS CANDIDATE TO MOVE
        DO K = 1,C                              ! I IS CURRENTLY IN CLUS. KPT
          IF(XB(I,K).EQ.1) KPT=K
        END DO
        DO 71 J = 1,E                           ! CHECK IF OBJECT J IS IN
          DO K = 1,C                            ! SAME CLUSTER AS I
            IF(XB(J,K).EQ.1) JPT=K              ! IF THEY ARE, AND THEY HIT
          END DO                                ! MAX DIAMETER, THEN TRY TO
          IF(KPT.NE.JPT) GO TO 71               ! MOVE I
          IF(A(I,J).LT.DIAMB-EPS) GO TO 71
          DO 72 KK = 1,C                        ! TRY TO MOVE I FROM K TO KK
            IF(K.EQ.KK) GO TO 72
            DO 73 L = 1,E
              IF(XB(L,KK).NE.1) GO TO 73
              IF(A(I,L).GT.DIAMB-EPS) GO TO 72
 73         CONTINUE
            XB(I,KPT)=0
            XB(I,KK)=1
            DIAMB=0.0D0
            DO II = 1,E-1
              DO JJ = II+1,E
                DO LL = 1,C
                  IF(XB(II,LL).EQ.1.AND.XB(JJ,LL).EQ.1) THEN
                    IF(A(II,JJ).GT.DIAMB) DIAMB=A(II,JJ)
                  END IF
                END DO
              END DO
            END DO
            GO TO 4500
 72       CONTINUE
 71     CONTINUE
 70   CONTINUE
C      write(*,*) diamb
C
C  ###########################################################
C    NOW THAT THE HEURISTIC IS COMPLETE, WE RELABEL THE ROWS
C    AND COLUMNS OF THE DISSIMILARITY MATRIX SUCH THAT THOSE
C    VERTICES WITH THE MEXIMUM NUMBER OF VALUES EXCEEDING THE
C    BOUND ARE PLACED FIRST IN THE REORDERING
C  ###########################################################
C
      DO I = 1,E
        S(I)=0
        DO J = 1,E
          IF(A(I,J).GT.DIAMB+EPS) S(I)=S(I)+1
        END DO
      END DO
      DO I = 1,E
        NP(I)=I
      END DO
      DO I = 1,E-1
        DO J = I+1,E
          IF(S(J).GT.S(I)) THEN
            IDUM=S(J)
            IDUM2=NP(J)
            S(J)=S(I)
            NP(J)=NP(I)
            S(I)=IDUM
            NP(I)=IDUM2
          END IF
        END DO
      END DO
      DO I = 1,E
        DO J = 1,E
          B(I,J) = A(NP(I),NP(J))
        END DO
      END DO
      DO I = 1,E
        DO J = 1,E
          A(I,J) = B(I,J)
        END DO
      END DO
C
      DO I = 1,E
        do j = 1,15
          xbt(i,j)=xb(np(i),j)
        end do
      END DO
      DO I = 1,E
        DO J = 1,15
          XB(I,J)=XBT(I,J)
          IF(XB(I,J).EQ.1) XBB(I)=J
        END DO
      END DO
C
C   ######################
C   STEP 0: INITIALIZE
C   ######################
C
      P=0
      Q=C
      Z=999999.
      DO K = 1,C
        N(K)=0
      END DO
      DO I = 1,E
        X(I)=0
      END DO
C
C   ###############################
C   STEP 1: INCREMENT SEARCH DEPTH
C   ###############################
C
 100  P=P+1
      M=1
      N(M)=N(M)+1
      IF(N(M).EQ.1) Q=Q-1
      X(P)=M
C
C   ###############################
C   STEP 2: FEASIBILITY
C   ###############################
C
 200  IF(E-P.LT.Q) GO TO 500
C
C   ###############################
C   STEP 3: SUBOPTIMALITY
C   ###############################
C
 300  DO K = 1,C
        DIAMC(K) = 0.0D0
      END DO
      DO I = 1,P-1
        K1 = X(I)
        do j = i+1,p
          k2 = x(j)
          if(k1.eq.k2.and.a(i,j).gt.DIAMC(K1)-eps) DIAMC(K1) = A(I,J)
        end do
      end do
      go to 315
      if(q.gt.0) go to 315
      DO 310 J = P+1,E            ! For each unassigned object
        MINK=999999.0D0
        DO K = 1,C
          KSUM(K)=0.0D0
        END DO
        DO I = 1,P
          KI = X(I)
          IF(A(I,J).GT.KSUM(KI)) KSUM(KI)=A(I,J)
        END DO
        DO K = 1,C
          XX = KSUM(K) - DIAMC(K)
          IF(XX.LT.MINK) THEN
            MINK = XX
            KSEL = K
          END IF
        END DO
        IF(MINK.GT.EPS) DIAMC(KSEL) = DIAMC(KSEL) + MINK
 310  CONTINUE
 315  DTEST = 0.0D0
      DO K = 1,C
        DTEST = DTEST + DIAMC(K)
      END DO
      IF(DTEST.GT.Z-EPS) GO TO 500
C
C   ###############################
C   STEP 4: UPDATE INCUMBENT
C   ###############################
C
 400  IF(P.NE.E) GO TO 100
      DO I = 1,E
        XBB(I)=X(I)
      END DO
      Z = 0.0D0
      DO K = 1,C
        BESDIAM(K) = 0.0D0
      END DO
      DO I = 1,P-1
        KI = X(I)
        DO J = I+1,P
          KJ=X(J)
          IF(KI.EQ.KJ.AND.A(I,J).GT.BESDIAM(KI)-EPS) BESDIAM(KI)=A(I,J)
        END DO
      END DO
      DO K = 1,C
        Z = Z + BESDIAM(K)
      END DO

C
C   ###############################
C   STEP 5: DETERMINE ACTION AFTER FATHOM
C   ###############################
C
 500  IF(M.EQ.C.OR.(N(M).EQ.1.AND.N(M+1).EQ.0)) GO TO 700
C
C   ###############################
C   STEP 6: FATHOM: BRANCH RIGHT ON GROUP
C   ###############################
C
 600  X(P)=0
      N(M)=N(M)-1
      M=M+1
      N(M)=N(M)+1
      IF(N(M).EQ.1) Q=Q-1
      X(P)=M
      GO TO 200
C
C   ###############################
C   STEP 7: FATHOM: DEPTH RETRACTION
C   ###############################
C
 700  X(P)=0
      N(M)=N(M)-1
      P = P-1
      IF(N(M).EQ.0) Q = Q+1
      IF(P.GT.0) THEN
        DO K = 1,C
          IF(X(P).EQ.K) THEN
            M=K
            GO TO 500
          END IF
        END DO
      END IF
C
      CALL GETTIM (IHR, IMIN, ISEC, I100)  
      CALL GETDAT (IYR, IMON, IDAY)  
      TIMEC = DFLOAT(86400*IDAY+3600*IHR+60*IMIN+ISEC)+DFLOAT(I100)/100.  
c
      write(*,826) z
      write(6,826) z
      DO K = 1,C
        WRITE(*,828) K,BESDIAM(K)
        WRITE(6,828) K,BESDIAM(K)
      END DO
      write(*,827) timec-timeb
      write(6,827) timec-timeb
      DO I = 1,E
        S(NP(I))=XBB(I)
      END DO
      WRITE(6,660) (S(I),I=1,E)
 660  FORMAT(24I3)
 820  FORMAT(10I4)
 825  format(' HEURISTIC SOLUTION  DIAMETER ',F12.5)
 826  format(' THE MINIMUM SUM OF DIAMETERS ',F12.5)
 827  format(' THE TOTAL COMPUTATION TIME   ',F12.2)
 828  FORMAT(I5,F15.4)
 889  STOP
      END