PROGRAM WCSUM
      IMPLICIT INTEGER(A-Z)
      DOUBLE PRECISION TIMEB,TIMEC,A(300,300),Z,SUMA(300),ZB,RHALF,
     1       A1(300,300),BUP2(300),DELMIN,DELTA,RSUM(20),ZXB,DELK,RMAX,
     1       RFIRST,B1(300,300),BAP,MINK,KSUM
      INTEGER N(20),X(300),XB(300),S(300),T(300),XREOR(300)
      OPEN(1,FILE='AMAT.DAT')
      OPEN(3,FILE='RESULT')
C
C  ##################################################################
C  5/18/2004 - BRANCH AND BOUND ALGORITHM FOR MINIMIZING THE
C              WITHIN CLUSTER SUM OF DISSIMILARITIES
C              -- NEAREST NEIGHBOR SEPARATION
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  -- NEAREST NEIGHBOR SEPARATION REORDERING OF THE OBJECTS
C  ##################################################################
C
      MAXHLP = 250
      READ(1,*) E1                      ! Read number of objects
      WRITE(*,*) 'TYPE 1 FOR HALF MATRIX OR TYPE 2 FOR FULL MATRIX'
      READ(*,*) ITYPE
      IF(ITYPE.EQ.2) THEN
        READ(1,*) ((A1(I,J),J=1,E1),I=1,E1)
      ELSE
        DO J = 2,E1
          READ(1,*) (A1(I,J),I=1,J-1)
        END DO
        DO J = 2,E1
          DO I = 1,J-1
            A1(J,I) = A1(I,J)
          END DO
        END DO
      END IF
      WRITE(*,*) ' PLEASE INPUT NUMBER OF CLUSTERS 2 TO 10'
      READ(*,*) C
      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
C ################ RELABELING ALGORITHM ########################
C
C      DO I = 1,E1            ! UNCOMMENT THESE 4 LINES IF NO
C        S(I) = I             ! REORDERING OF OBJECTS IS DESIRED
C      END DO
C      GO TO 586
      IFIRST = 0
      ILAST = E1+1
      RHALF = FLOAT(E1)/2.
      WRITE(*,*) RHALF
 580  RMAX=999999.
      DO 501 I = 1,E1-1
        IF(T(I).EQ.1) GO TO 501
        DO 502 J = 1+I,E1
          IF(T(J).EQ.1) GO TO 502
          IF(A1(I,J).LT.RMAX) THEN
            RMAX=A1(I,J)
            ISEL=I
            JSEL=J
          END IF
 502    CONTINUE
 501  CONTINUE
      IFIRST=IFIRST+1
      S(IFIRST)=ISEL
      ILAST=ILAST-1
      S(ILAST)=JSEL
      T(ISEL)=1
      T(JSEL)=1
      RFIRST = IFIRST
      WRITE(*,*) IFIRST,S(IFIRST),ILAST,S(ILAST)
      IF(RFIRST.EQ.RHALF) GO TO 505
      IF(RFIRST.GT.RHALF-.50001) THEN
        IFIRST = IFIRST + 1
        DO I = 1,E1
          IF(T(I).EQ.0) S(IFIRST) = I
        END DO
        GO TO 505
      END IF
      GO TO 580
C
 505  DO I = 1,E1
        DO J = 1,E1
          B1(I,J)=A1(S(I),S(J))
        END DO
      END DO
      DO I = 1,E1
        DO J = 1,E1
          A1(I,J)= B1(I,J)
        END DO
      END DO
C
 586  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)
          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
      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) SUMA(M)=SUMA(M)+A(P,J)
      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  ZB = 0.
      DO 310 I = 1,C-Q
        ZB = ZB + SUMA(I)
 310  CONTINUE
      BAP=0
      DO J = P+1,E                ! For each unassigned object
        MINK=999999
        DO K = 1,C
          KSUM=0
          DO I = 1,P
            IF(X(I).EQ.K) KSUM=KSUM+A(I,J)
          END DO
          IF(KSUM.LT.MINK) MINK=KSUM
        END DO
        BAP=BAP+MINK
      END DO
      IF(Z.LE.ZB+BUP2(P)+BAP) GO TO 500
C
C   ###############################
C   STEP 4: UPDATE INCUMBENT
C   ###############################
C
 400  IF(P.NE.E) GO TO 100
      Z = ZB
      DO I = 1,E
        XB(I)=X(I)
      END DO
C
C   ###############################
C   STEP 5: DETERMINE ACTION AFTER FATHOM
C   ###############################
C
C 500  IF(M.EQ.C.OR.(N(M).EQ.1.AND.N(M+1).EQ.0)) GO TO 700
 500  IF(M.EQ.C.OR.N(M).EQ.1) GO TO 700    ! MODIFIED 1/7/2005
C
C   ###############################
C   STEP 6: FATHOM: BRANCH RIGHT ON GROUP
C   ###############################
C
 600  DO J = 1,P-1                         ! X(P) = 0 DELETED 1/7/2005
        IF(X(J).EQ.M) SUMA(M)=SUMA(M)-A(P,J)
      END DO
      N(M)=N(M)-1
C      IF(N(M).EQ.0) Q = Q+1     ! ADDED 4/20/04, DELETED 1/7/2005
      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) SUMA(M)=SUMA(M)+A(P,J)
      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) SUMA(M)=SUMA(M)-A(P,J)
      END DO
      N(M)=N(M)-1
      P = P-1
      IF(N(M).EQ.0) Q = Q+1
      IF(P.GT.0) THEN         ! MODIFIED 1/7/2005
        M=X(P)
        GO TO 500
      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.
        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)
          END DO
        END DO
        II = E1-E
        DELMIN = 99999999.
        DO 710 K = 1,C
          DELK = RSUM(K)
          DO 711 I = 1,E
            K1 = XB(I)
            IF(K1.NE.K) GO TO 711
            DELK = DELK + A1(II,I+II)
 711      CONTINUE
          DELTA = DELK - RSUM(K) ! /DFLOAT(N(K)+1) - RSUM(K)/DFLOAT(N(K))
          IF(DELTA.LT.DELMIN) DELMIN = DELTA
 710    CONTINUE
        ZXB = Z + DELMIN + .0001
      END IF
 4500 CONTINUE
C
      DO I = 1,E1         ! MAP BACK TO THE ORIGINAL OBJECT ORDERING
       I1 = S(I)
       XREOR(I1) = XB(I)
      END DO
      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) Z
      write(*,825) timec-timeb
      WRITE(*,820) (XREOR(I),I=1,E)
      WRITE(3,828) Z, TIMEC-TIMEB
 828  FORMAT(F15.5,F10.2)
 820  FORMAT(30I3)
 824  FORMAT(' MINIMUM WITHIN CLUSTER SUMS OF DISSIMILARITIES ',F15.5)
 825  format(' TOTAL CPU TIME = ',f16.2)
 830  format(' NUMBER OF OBJECTS ',I3,' Z = ',2F12.5)
 889  STOP
      END