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),MAXOLD,MAXNEW,
     1                 RMAX1,RMAX2,MAXIND,rbef
      REAL S1
      INTEGER N(250),X(250),XB(250),XBB(250),V(250,250),CL(250,250),
     1        IOPEN(250),NIC(250),MEMB(250),S(250),NP(250),CNUM(250),
     1        XT(250)
      OPEN(1,FILE='AMAT.DAT')
      OPEN(6,FILE='RESULTS')
      EPS=1.0D-08
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
      WRITE(*,*) 'PLEASE TYPE 1 FOR COMPLETE LINKAGE STARTING SOLUTION '
      WRITE(*,*) 'OR 2 FOR RANDOM STARTING SOLUTION'
      READ(*,*) ISTRT
      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,C
          NIC(I)=0
        END DO
        DO J = 1,E
          MEMB(J)=0
          CNUM(J) = 1
          CL(J,1) = J 
        END DO
C
      IF(ISTRT.EQ.1) THEN     ! BIASED COMPLETE LINKAGE SOLUTION
C
      DO KKK = E,C+1,-1
        RMAX1=9.9D+12
        RMAX2=9.9D+12
        DO 1 II = 1,E-1
          IF(CNUM(II).EQ.0) GO TO 1
          DO 2 JJ = II+1,E
            IF(CNUM(JJ).EQ.0) GO TO 2
            MAXIND = 0.0D0
            DO 3 K = 1,CNUM(II)
              I1 = CL(II,K)
              DO 4 L = 1,CNUM(JJ)
                J1 = CL(JJ,L)
                IF(A(I1,J1).GT.MAXIND) MAXIND = A(I1,J1)
   4          CONTINUE
   3        CONTINUE
            IF(MAXIND.LT.RMAX1) THEN
              RMAX2 = RMAX1
              K2A = K1A
              K2B = K1B
              RMAX1=MAXIND
              K1A = II
              K1B=JJ
            ELSEIF(MAXIND.LT.RMAX2) THEN
              RMAX2=MAXIND
              K2A = II
              K2B=JJ
            END IF
  2       CONTINUE
  1     CONTINUE
        CALL RANDOM(S1)
        IF(S1.GT..8.OR.RMAX2.GT.99999.0d0) THEN
          DO I = CNUM(K1A)+1,CNUM(K1A)+CNUM(K1B)
            CL(K1A,I) = CL(K1B,I-CNUM(K1A))
          END DO
          CNUM(K1A) = CNUM(K1A)+CNUM(K1B)
          CNUM(K1B) = 0
        ELSE
          DO I = CNUM(K2A)+1,CNUM(K2A)+CNUM(K2B)
            CL(K2A,I) = CL(K2B,I-CNUM(K2A))
          END DO
          CNUM(K2A) = CNUM(K2A)+CNUM(K2B)
          CNUM(K2B) = 0
        END IF
      END DO
      ICT = 0
      DO 420 I = 1,E
        IF(CNUM(I).EQ.0) GO TO 420
        ICT=ICT+1
        NIC(ICT) = CNUM(I)
        DO J = 1,NIC(ICT)
          J1 = CL(I,J)
          MEMB(J1) = ICT
        END DO
 420  CONTINUE
C
      ELSE                          ! RANDOM STARTING SOLUTION
C
      DO I = 1,E
        CALL RANDOM(S1)
        KSEL = S1 * FLOAT(C) + 1.
        IF(KSEL.GT.C) KSEL = 1
        MEMB(I) = KSEL
        NIC(KSEL) = NIC(KSEL)+1
      END DO
 451  ITRIG = 0
      DO 450 K = 1,C
        IF(NIC(K).GT.0) GO TO 450
        ITRIG = 1
        CALL RANDOM(S1)
        ISEL = S1 * FLOAT(E) + 1.
        KREM = MEMB(ISEL)
        NIC(KREM) = NIC(KREM)-1
        NIC(K) = NIC(K) + 1
        MEMB(ISEL) = K
 450  CONTINUE
      IF(ITRIG.EQ.1) GO TO 451
C
      END IF
C
      Z = 0.0D0
      DO I = 1,E-1
        DO J = I + 1,E
          IF(MEMB(I).EQ.MEMB(J).AND.A(I,J).GT.Z+EPS) Z = A(I,J)
        END DO
      END DO
      RBEF = Z
c
 4499 IGLOB=0
 4500 ITRIG=0                                   ! RELOCATION
      DO 70 I = 1,E                             ! I IS CANDIDATE TO MOVE
        KPI = MEMB(I)
        MAXOLD = 0.0D0
        DO 71 J = 1,E                           ! CHECK IF OBJECT J IS IN
          KPJ = MEMB(J)
          IF(KPI.NE.KPJ) GO TO 71
C          IF(A(I,J).GT.Z-EPS) GO TO 170
           IF(A(I,J).GT.MAXOLD) MAXOLD=A(I,J)
 71     CONTINUE
C        GO TO 70
 170    DO 72 KK = 1,C                        ! TRY TO MOVE I FROM K TO KK
          IF(KPI.EQ.KK) GO TO 72
          MAXNEW=0.0D0
          DO 73 L = 1,E
            KPL = MEMB(L)
            IF(KPL.NE.KK) GO TO 73
C            IF(A(I,L).GT.Z-EPS) GO TO 72
            IF(A(I,L).GT.MAXNEW) MAXNEW=A(I,L)
 73       CONTINUE
          IF(MAXOLD.GT.MAXNEW+EPS) then
          MEMB(I)=KK
          KPI=KK
          MAXOLD = MAXNEW
          ITRIG=1
          IGLOB=1
          Z=0.0D0
          DO II = 1,E-1
            DO JJ = II+1,E
              IF(MEMB(II).EQ.MEMB(JJ).AND.A(II,JJ).GT.Z+EPS) Z=A(II,JJ)
            END DO
          END DO
          END IF
 72     CONTINUE
 70   CONTINUE
      IF(ITRIG.EQ.1)  GO TO 4500
C      go to 178
C
 4600 ITRIG=0                                   ! INTERCHANGE
      DO 75 I = 1,E-1                           ! OBJECT 1
        KPI = MEMB(I)
        DO 76 J = I+1,E                         ! OBJECT 2
          KPJ = MEMB(J)
          IF(KPI.EQ.KPJ) GO TO 76
          MAXOLD = 0.0D0
          MAXNEW = 0.0D0
          DO 77 L = 1,E
            IF(L.EQ.I.OR.L.EQ.J) GO TO 77
            KPL = MEMB(L)
            IF(KPL.NE.KPI.AND.KPL.NE.KPJ) GO TO 77
            IF(KPL.EQ.KPI) THEN
              IF(A(I,L).GT.MAXOLD+EPS) MAXOLD = A(I,L)
              IF(A(J,L).GT.MAXNEW+EPS) MAXNEW = A(J,L)
            ELSE
              IF(A(J,L).GT.MAXOLD+EPS) MAXOLD = A(J,L)
              IF(A(I,L).GT.MAXNEW+EPS) MAXNEW = A(I,L)
            END IF
 77       CONTINUE
          IF(MAXOLD.GT.MAXNEW+EPS) THEN
            MEMB(I)=KPJ
            MEMB(J)=KPI
            KPI = KPJ
            ITRIG=1
            IGLOB=1
            Z=0.0D0
            DO II = 1,E-1
              DO JJ = II+1,E
             IF(MEMB(II).EQ.MEMB(JJ).AND.A(II,JJ).GT.Z+EPS) Z=A(II,JJ)
              END DO
            END DO
          END IF
 76     CONTINUE
 75   CONTINUE
      IF(ITRIG.EQ.1)  GO TO 4600
      IF(IGLOB.EQ.1) GO TO 4499
C
 178  write(*,79) ijkl,RBEF,z
 79   format(i5,2f18.5)
      IF(Z.LT.DIAMB) THEN
        DIAMB = Z
        DO I = 1,E
          XB(I) = MEMB(I)
        END DO
      END IF
 4000 CONTINUE
      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 k = 1,c
        nic(k) = 0
      end do
      do i = 1,e
        K = XB(I)
        nic(k) = nic(k) + 1
        v(k,nic(k)) = i
      end do
      write(*,183) (nic(k),k=1,c)
 183  format(15i4)
      nsel = 0
      ksel = 0
  515 nsel = nsel+1
  516 ksel = ksel+1
      if(ksel.gt.c) ksel = 1
      if(nic(ksel).eq.0) go to 516
      if(nsel.le.c) then
        call random(s1)
        kv = s1 * nic(ksel)+1.
        np(nsel) = v(ksel,kv)
        do i = kv,nic(ksel)-1
          v(ksel,i) = v(ksel,i+1)
        end do
        nic(ksel) = nic(ksel)-1
      else
        imax = -1
        do k = 1,nic(ksel)
          icrit = 0
          do i = 1,nsel-1
            if(a(np(i),v(ksel,k)).gt.diamb-eps) icrit = icrit + 1
          end do
          if(icrit.gt.imax) then
            imax = icrit
            kv = k
          end if
        end do
        np(nsel) = v(ksel,kv)
        do i = kv,nic(ksel)-1
          v(ksel,i) = v(ksel,i+1)
        end do
        nic(ksel) = nic(ksel)-1
      end if
      if(nsel.lt.e) go to 515
C
      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
 182  DO I = 1,E
        xt(i)=xb(np(i))
      END DO
      DO I = 1,E
        XB(I)=XT(I)
        XBB(I)=XB(I)
      END DO
C
C   ######################
C   STEP 0: INITIALIZE
C   ######################
C
      P=0
      Q=C
      Z=DIAMB
      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 I = 1,P-1             ! THIS IS BETTER THAN THE COMMENTED
        K1 = X(I)              ! BLOCK BECAUSE OF THE POTENTIAL FOR
        DO J = I+1,P           ! Z TO CHANGE
          K2 = X(J)
          IF(K1.EQ.K2.and.A(I,J).GT.Z-EPS) GO TO 500
        END DO
      END DO
c        IF(K1.EQ.M)THEN
c          IF(A(I,P).GT.Z-EPS) GO TO 500
c        END IF
c      END DO
      if(q.gt.0) go to 400
      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
          IF(KSUM(K).LT.Z-EPS) GO TO 310
        END DO
        GO TO 500
 310  CONTINUE
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 I = 1,P-1
        KI = X(I)
        DO J = I+1,P
          KJ=X(J)
          IF(KI.EQ.KJ.AND.A(I,J).GT.Z-EPS) Z=A(I,J)
        END DO
      END DO
      write(*,*) z
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(*,825) diamb
      write(6,825) diamb
      write(*,826) z
      write(6,826) z
      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 OPTIMAL MINIMUM DIAMETER ',F12.5)
 827  format(' THE TOTAL COMPUTATION TIME   ',F12.2)
 889  STOP
      END