      PROGRAM EXCH
      IMPLICIT INTEGER(A-Z)
      DOUBLE PRECISION TIMEA, TIMEB, TLIMIT, Z, GBEST, ZTRIAL
      REAL S1,S2,S3
      INTEGER NR(40),A(200,200),NC(40), S(200),
     1        CMEM(200),RMEM(200),GR(200),GC(200),
     1        TRIAL0(40,40), TRIAL1(40,40),
     1        SUMBEST0(40,40), SUMBEST1(40,40),CENTR(40,200),
     1        CENTC(200,40),IDX(200)
C
C  ################################################################
C     DESTROY-AND-RECONSTRUCT HEURISTIC
C        1. OBJECT TRANSFERS ONLY -- NO EXCHANGES
C        2. RANDOMLY CHOOSE K1/K2 CENTROIDS -- ASSIGN EACH OBJECT TO NEAREST
C        3. SAME AS DRH - BUT RESIDUAL CELLS PERMITTED
C  ################################################################
C
      OPEN(2, FILE='SBT.DAT')
      OPEN(4, FILE='TESTROW')
      OPEN(5, FILE='TESTCOL')
      OPEN(3, FILE='RESULT')
      OPEN(6, FILE='PARTIT')
      OPEN(8, FILE='AMAT')
      OPEN(9, FILE='MAXIT')
C
C  INPUT FORMAT:
C  LINE 1: NUMBER OF ROW AND COLUMN OBJECTS
C  LINE 2: NUMBER OF ONES IN ROW, THEN THE COLUMNS THAT ARE ONES
C
      READ(2,*) RO,CO
      DO 611 I = 1,RO
        DO J = 1,CO
          IDX(J) = 0
        END DO
        READ(2,*) IJNK,IPT,(IDX(J),J=1,IPT)
        DO J = 1,IPT
          IF(IDX(J).NE.0) A(I,IDX(J)) = 1
        END DO
 611  CONTINUE
C
C OUTPUT FORMAT:
C THE SAME MATRIX IS WRITTEN OUT TO AMAT
C
      WRITE(8,*) RO,CO
      DO I = 1,RO
        WRITE(8,622) (A(I,J),J=1,CO)
      END DO
 622  FORMAT(100I2)
C
C THE MAXIMUM POSSIBLE NUMBER OF CLUSTERS IS THE MIN{RO, CO}
      NCMAX = RO
      IF(CO.LT.NCMAX) NCMAX = CO
      IF(NCMAX.GT.20) NCMAX = 20
C
      CALL GETTIM (IHR, IMIN, ISEC, I100)
      CALL GETDAT (IYR, IMON, IDAY)
      TIMEA=DFLOAT(86400*IDAY+3600*IHR+60*IMIN+ISEC)+DFLOAT(I100)/100.
C
C  ################################################################
C  RATHER THAN A FIXED TIME LIMIT, I'LL FIX ITERATIONS
C  ################################################################
C
C  GBEST IS GLOBAL MAXIMUM GEI ACROSS ALL NUMBERS OF CLUSTERS
C  MAXIT IS THE MAXIMUM NUMBER OF ITERATIONS FOR EACH K
C  MAXPERT IS THE MAXIMUM PERTURBATION SIZE
C
      GBEST = -999.
      NBEST = 0
      READ(9,*) MAXIT
      MAXPERT = (CO+RO)/2
C
C   MAXPERT CONTROLS MAXIMUM LEVEL OF DESTRUCTION - REASSIGNMENTS
C   IPERT IS THE CURRENT NUMBER OF REASSIGNMENTS FOR DESTRUCTION
C   INTENSE IS THE NUMBER OF ITERATIONS FOR SEARCH INTENSITY AT A
C           GIVEN DESTRUCTION LEVELS WITH NO IMPROVEMENT
C   MAXIT IS AN UPPER LIMIT ON TOTAL SOLUTIONS GENERATED
C
C   THE OUTER LOOP (8000) WILL RUN THE ALGORITHM FOR DIFFERENT K 
C
      NOIMP = 0
      DO 8000 INC = 2,NCMAX-2
        NOIMP = NOIMP+1
        NREPS = 0
        RC = INC
        CC = RC
        DO I = 1,RC
          NR(I) = 0
        END DO
        DO I = 1,CC
          NC(I) = 0
        END DO
        DO I = 1,RO
          RMEM(I) = 0
        END DO
        DO I = 1,CO
          CMEM(I) = 0
        END DO
        ONES = 0
        DO I = 1,RO
          DO J = 1,CO
            ONES = ONES + A(I,J)
          END DO
        END DO
C
C  ###############################################################
C              GREEDY INITIALIZATION PROCEDURE - PART 1:
C
C    RANDOMLY SELECT K1 OF THE ROW OBJECTS TO PUT ONE OBJECT IN
C    EACH OF THE FIRST K1 ROW CLUSTERS.  THEN, DO THE SAME FOR
C    TO OBTAIN K2 COLUMN OBJECTS TO PUT ONE OBJECT IN EACH OF
C    THE FIRST K2 COLUMN CLUSTERS.  THEN COMPUTE THE OBJECTIVE
C    VALUE.
C  ################################################################
C
        DO K = 1,RC
 98       CALL RANDOM(S1)
          ISEL = S1 * FLOAT(RO) + 1.
          IF(ISEL.GT.RO) ISEL = RO
          IF(RMEM(ISEL).GT.0) GO TO 98
          RMEM(ISEL) = K
          NR(K) = NR(K) + 1
          DO J = 1,CO
            CENTR(K,J) = A(ISEL,J)
          END DO
        END DO
        DO L = 1,CC
 97       CALL RANDOM(S1)
          JSEL = S1 * FLOAT(CO) + 1.
          IF(JSEL.GT.CO) JSEL = CO
          IF(CMEM(JSEL).GT.0) GO TO 97
          CMEM(JSEL) = L
          NC(L) = NC(L) + 1
          DO I = 1,RO
            CENTC(I,L) = A(I,JSEL)
          END DO
        END DO
        DO 300 I = 1,RO
          IF(RMEM(I).NE.0) GO TO 300
          DMIN = 9999999
          DO K = 1,RC
            DCOM = 0
            DO J = 1,CO
              DCOM = DCOM + (A(I,J)-CENTR(K,J))**2
            END DO
            IF(DCOM.LT.DMIN) THEN
              DMIN = DCOM
              KSEL = K
            END IF
          END DO
          RMEM(I) = KSEL
          NR(KSEL) = NR(KSEL) + 1
 300    CONTINUE
        DO 301 J = 1,CO
          IF(CMEM(J).NE.0) GO TO 301
          DMIN = 9999999
          DO K = 1,CC
            DCOM = 0
            DO I = 1,RO
              DCOM = DCOM + (A(I,J)-CENTC(I,K))**2
            END DO
            IF(DCOM.LT.DMIN) THEN
              DMIN = DCOM
              KSEL = K
            END IF
          END DO
          CMEM(J) = KSEL
          NC(KSEL) = NC(KSEL) + 1
 301    CONTINUE
C
 1700   VOIDS = 0
        EXCEPTS = 0
        DO I = 1,RO
          K = RMEM(I)
          DO J = 1,CO
            L = CMEM(J)
            IF(K.EQ.L) THEN
              IF(A(I,J).EQ.0) VOIDS = VOIDS + 1
            ELSE
              IF(A(I,J).EQ.1) EXCEPTS = EXCEPTS + 1
            END IF
          END DO
        END DO
        Z = DFLOAT(ONES-EXCEPTS)/DFLOAT(ONES+VOIDS)
C
C THE INITIAL SOLUTION FOR RC = CC = INC CLUSTERS IS GENERATED
C NOW, THE SUBROUTINE BELOW WILL USE DESTROY RECONSTRUCT TO
C      IMPROVE THE SOLUTION
C
 1100   ITRIG = 0
C
C       ##################################################
C              EVALUATE ALL ROW OBJECT RELOCATIONS
C       ##################################################
C
        DO 600 I = 1,RO
          K = RMEM(I)
          IF(NR(K)+NC(K).EQ.1) GO TO 600
          DO 601 KK = 1,RC
            IF(KK.EQ.K) GO TO 601
            VTRIAL = VOIDS
            ETRIAL = EXCEPTS
            DO 602 J = 1,CO
              L = CMEM(J)
              IF(K.EQ.L) THEN
                IF(A(I,J).EQ.0) VTRIAL = VTRIAL - 1
              ELSE
                IF(A(I,J).EQ.1) ETRIAL = ETRIAL - 1
              END IF
              IF(KK.EQ.L) THEN
                IF(A(I,J).EQ.0) VTRIAL = VTRIAL + 1
              ELSE
                IF(A(I,J).EQ.1) ETRIAL = ETRIAL + 1
              END IF
  602       CONTINUE
            ZTRIAL = DFLOAT(ONES-ETRIAL)/DFLOAT(ONES+VTRIAL)
            IF(ZTRIAL.GT.Z+.00001) THEN
              Z = ZTRIAL
              NR(K) = NR(K) - 1
              NR(KK) = NR(KK) + 1
              VOIDS = VTRIAL
              EXCEPTS = ETRIAL
              RMEM(I) = KK
              ITRIG = 1
              K = KK
C              go to 600
            END IF
 601      CONTINUE
 600    CONTINUE
C
C
C       ##################################################
C              EVALUATE ALL COLUMN OBJECT RELOCATIONS
C       ##################################################
C
        DO 700 J = 1,CO
          L = CMEM(J)
          IF(NC(L)+NR(L).EQ.1) GO TO 700
          DO 701 LL = 1,CC
            IF(LL.EQ.L) GO TO 701
            VTRIAL = VOIDS
            ETRIAL = EXCEPTS
            DO 702 I = 1,RO
              K = RMEM(I)
              IF(K.EQ.L) THEN
                IF(A(I,J).EQ.0) VTRIAL = VTRIAL - 1
              ELSE
                IF(A(I,J).EQ.1) ETRIAL = ETRIAL - 1
              END IF
              IF(K.EQ.LL) THEN
                IF(A(I,J).EQ.0) VTRIAL = VTRIAL + 1
              ELSE
                IF(A(I,J).EQ.1) ETRIAL = ETRIAL + 1
              END IF
  702       CONTINUE
            ZTRIAL = DFLOAT(ONES-ETRIAL)/DFLOAT(ONES+VTRIAL)
            IF(ZTRIAL.GT.Z+.00001) THEN
              Z = ZTRIAL
              NC(L) = NC(L) - 1
              NC(LL) = NC(LL) + 1
              VOIDS = VTRIAL
              EXCEPTS = ETRIAL
              CMEM(J) = LL
              ITRIG = 1
              L = LL
C              go to 700
            END IF
 701      CONTINUE
 700    CONTINUE
C
        IF(ITRIG.EQ.1) GO TO 1100
C
        IF(Z.GT.GBEST+.00001) THEN
          NOIMP = 0
          vbest = voids
          ebest = excepts
          GBEST = Z
          NBEST = 1
          KBEST = RC
          DO I = 1,RO
            GR(I) = RMEM(I)
          END DO
          DO I = 1,CO
            GC(I) = CMEM(I)
          END DO
        ELSEIF(Z.GE.GBEST-.00001) THEN
          NBEST = NBEST + 1
        END IF
C
        NREPS = NREPS + 1
        IF(NREPS.EQ.MAXIT) THEN
          IF(NOIMP.EQ.2) GO TO 8001
          GO TO 8000
        ELSE
          CALL RANDOM(S1)
          NPERT = 2+FLOAT(MAXPERT)*S1
          DO I = 1,NPERT
            CALL RANDOM(S1)
            IF(S1.LT.0.5) THEN
 53           CALL RANDOM(S2)
              IROW = 1+S2*FLOAT(RO)
              IF(IROW.GT.RO) IROW = 1
              ICLUS = RMEM(IROW)
              IF(NR(ICLUS)+NC(ICLUS).EQ.1) GO TO 53
              CALL RANDOM(S3)
              ICLUS2 = 1+S2*FLOAT(RC)
              IF(ICLUS2.GT.RC) ICLUS2 = 1
              IF(ICLUS2.EQ.ICLUS) ICLUS2 = ICLUS2+1
              IF(ICLUS2.GT.RC) ICLUS2 = 1
              RMEM(IROW) = ICLUS2
              NR(ICLUS2) = NR(ICLUS2)+1
              NR(ICLUS) = NR(ICLUS)-1
            ELSE
 54           CALL RANDOM(S2)
              ICOL = 1+S2*FLOAT(CO)
              IF(ICOL.GT.CO) ICOL = 1
              ICLUS = CMEM(ICOL)
              IF(NC(ICLUS)+NR(ICLUS).EQ.1) GO TO 54
              CALL RANDOM(S3)
              ICLUS2 = 1+S2*FLOAT(CC)
              IF(ICLUS2.GT.CC) ICLUS2 = 1
              IF(ICLUS2.EQ.ICLUS) ICLUS2 = ICLUS2+1
              IF(ICLUS2.GT.CC) ICLUS2 = 1
              CMEM(ICOL) = ICLUS2
              NC(ICLUS2) = NC(ICLUS2)+1
              NC(ICLUS) = NC(ICLUS)-1
            END IF
          END DO
          GO TO 1700
        END IF
 8000 CONTINUE
 8001 DO I = 1,RO
        WRITE(4,*) GR(I)
      END DO
      DO J = 1,CO
        WRITE(5,*) GC(J)
      END DO
      WRITE(*,272) (GR(I),I=1,RO)
      WRITE(*,272) (GC(I),I=1,CO)
      WRITE(6,272) (GR(I),I=1,RO)
      WRITE(6,272) (GC(I),I=1,CO)
 272  FORMAT(40I3)
      CALL GETTIM (IHR, IMIN, ISEC, I100)
      CALL GETDAT (IYR, IMON, IDAY)
      TIMEB=DFLOAT(86400*IDAY+3600*IHR+60*IMIN+ISEC)+DFLOAT(I100)/100.
      WRITE(*,61) GBEST, TIMEB-TIMEA,NREPS,NBEST, VBEST, EBEST, KBEST
      WRITE(3,61) GBEST, TIMEB-TIMEA,NREPS,NBEST, VBEST, EBEST, KBEST
  61  FORMAT(' Z = ',F8.5,' TIME = ',F8.2,' RSTS = ',I8,' BESTS ',I8,
     1      ' VOIDS ',I4,' EXCEPTS ',I4,' CELLS ',I4)
 273  FORMAT(F8.5)
      END
