PROGRAM BBBI
      IMPLICIT INTEGER(A-Z)
      DOUBLE PRECISION TIMEB,TIMEC,A(150,150),Z,SUMA(150),ZB,
     1       A1(150,150),BUP2(150),DELMIN,DELTA,RSUM(20),ZXB,DELK,
     1       B(150,150),B1(150,150),WA,WB,SUMB(150),ZBA,ZBB,DELKB,
     1       DELTAB,RSUMB(20),ZBESTA,ZBESTB
      INTEGER N(20),X(150),XB(150)
      OPEN(1,FILE='AMAT.DAT')
      OPEN(3,FILE='BMAT.DAT')
C
C  ##################################################################
C  BICRITERION WITHING CLUSTER SUM OF SQUARES CLUSTERING
C            - BRANCH AND BOUND ALGORITHM FOR MINIMIZING WEIGHTED
C              WITHIN CLUSTER SUM OF SQUARES (SUM OF DISTANCES IN
C              CLUSTER DIVIDED BY N) FOR TWO DISTANCE MATRICES
C  -- SOLVES PROBLEMS FROM SIZE C+1 TO N FROM THE BACK
C  -- INCREMENTAL SOLUTION APPROACH ALLOWS TIGHT BOUNDS TO BE RAPIDLY
C  -- PROVIDED.
C  ##################################################################
C
      MAXHLP = 140
      READ(1,*) E1                        ! NUMBER OF OBJECTS
      READ(3,*) E1
      WRITE(*,*) ' TYPE 1 FOR HALF MATRIX OR 2 FOR FULL MATRIX INPUT'
      READ(*,*) ITYPE
      IF(ITYPE.EQ.1) THEN
      DO I = 2,E1
      READ(1,*) (A1(I,J),J=1,I-1)         ! PROXIMITY MATRIX
      READ(3,*) (B1(I,J),J=1,I-1)         ! PROXIMITY MATRIX
      END DO
      DO I = 2,E1
        DO J = 1,I-1
          A1(J,I) = A1(I,J)
          B1(J,I) = B1(I,J)
        END DO
      END DO
      ELSE
      READ(1,*) ((A1(I,J),J=1,E1),I=1,E1)   ! PROXIMITY MATRIX
      READ(3,*) ((B1(I,J),J=1,E1),I=1,E1)   ! PROXIMITY MATRIX
      END IF
      ZXB=99999999.
      CALL GETTIM (IHR, IMIN, ISEC, I100)
      CALL GETDAT (IYR, IMON, IDAY)  
      TIMEB = DFLOAT(86400*IDAY+3600*IHR+60*IMIN+ISEC)+DFLOAT(I100)/100.  
C
      WRITE(*,*) ' INPUT NUMBER OF CLUSTERS '
      READ(*,*) C
      WRITE(*,*) ' INPUT WEIGHTS '
      READ(*,*) WA,WB
      E=C
      DO 4500 IJKL = C+1,E1
        E=E+1
        IF(E.GT.MAXHLP.AND.E.LT.E1) GO TO 4499
        DO I = E1-E+1,E1
          DO J = E1-E+1,E1
            A(I-E1+E,J-E1+E)=A1(I,J)
            B(I-E1+E,J-E1+E)=B1(I,J)
          END DO
        END DO
C
C   ######################
C   STEP 0: INITIALIZE
C   ######################
C
      P=0
      Q=C
      Z = ZXB
      DO K = 1,C
        N(K)=0
      END DO
      DO I = 1,E
        X(I)=0
        suma(i)=0.0d0
        sumb(i)=0.0d0
      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
      DO J = 1,P-1
        IF(X(J).EQ.M) THEN
          SUMA(M)=SUMA(M)+A(P,J)
          SUMB(M)=SUMB(M)+B(P,J)
        END IF
      END DO
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  ZBA = 0.
      ZBB = 0.
      DO 310 I = 1,C-Q
        ZBA = ZBA + SUMA(I) / DFLOAT(N(I))
        ZBB = ZBB + SUMB(I) / DFLOAT(N(I))
 310  CONTINUE
      IF(Z.LE.WA*ZBA+WB*ZBB+BUP2(P)) GO TO 500
C
C   ###############################
C   STEP 4: UPDATE INCUMBENT
C   ###############################
C
 400  IF(P.NE.E) GO TO 100
      Z = WA*ZBA+WB*ZBB
      DO I = 1,E
        XB(I)=X(I)
      END DO
      ZBESTA = ZBA
      ZBESTB = ZBB
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
      DO J = 1,P-1
        IF(X(J).EQ.M) THEN
          SUMA(M)=SUMA(M)-A(P,J)
          SUMB(M)=SUMB(M)-B(P,J)
        END IF
      END DO
      N(M)=N(M)-1
      IF(N(M).EQ.0) Q = Q+1     ! ADDED 4/20/04
      M=M+1
      N(M)=N(M)+1
      IF(N(M).EQ.1) Q=Q-1
      X(P)=M
      DO J = 1,P-1
        IF(X(J).EQ.M) THEN
          SUMA(M)=SUMA(M)+A(P,J)
          SUMB(M)=SUMB(M)+B(P,J)
        END IF
      END DO
      GO TO 200
C
C   ###############################
C   STEP 7: FATHOM: DEPTH RETRACTION
C   ###############################
C
 700  X(P)=0
      DO J = 1,P-1
        IF(X(J).EQ.M) THEN
          SUMA(M)=SUMA(M)-A(P,J)
          SUMB(M)=SUMB(M)-B(P,J)
        END IF
      END DO
      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
      WRITE(*, 830) E,Z,zxb
 4499 DO I = E,2,-1
        BUP2(I)=BUP2(I-1)
      END DO
      BUP2(1)=Z
      IF(E.NE.E1) THEN
        DO K = 1,C
          N(K) = 0
          RSUM(K)=0.
          RSUMB(K) = 0.
        END DO
        DO I = 1,E
          K = XB(I)
          N(K) = N(K)+1
        END DO
        DO I = 1,E-1
          K1 = XB(I)
          DO J = I+1,E
            K2 = XB(J)
            IF(K1.EQ.K2) RSUM(K1) = RSUM(K1) + A(I,J)
            IF(K1.EQ.K2) RSUMB(K1) = RSUMB(K1) + B(I,J)
          END DO
        END DO
        II = E1-E
        DELMIN = 99999999.
        DO 710 K = 1,C
          DELK = RSUM(K)
          DELKB = RSUMB(K)
          DO 711 I = 1,E
            K1 = XB(I)
            IF(K1.NE.K) GO TO 711
            DELK = DELK + A1(II,I+II)
            DELKB = DELKB + B1(II,I+II)
 711      CONTINUE
          DELTA = DELK/DFLOAT(N(K)+1) - RSUM(K)/DFLOAT(N(K))
          DELTAB = DELKB/DFLOAT(N(K)+1) - RSUMB(K)/DFLOAT(N(K))
          DELTA = WA*DELTA + WB*DELTAB
          IF(DELTA.LT.DELMIN) DELMIN = DELTA
 710    CONTINUE
        ZXB = Z + DELMIN + .0001
      END IF
 4500 CONTINUE
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(*,824) zbesta,zbestb
      write(*,825) timec-timeb
      WRITE(*,820) (XB(I),I=1,E)
      WRITE(*,*) ZBESTA,ZBESTB
 820  FORMAT(30I3)
 824  FORMAT(' WCSS MATRIX A = ',F16.5,' WCSS MATRIX B = ',F16.5)
 825  format(' TOTAL CPU TIME ',f16.2)
 830  format(' NUMBER OF OBJECTS ',I3,' Z = ',2F12.5)
 889  STOP
      END