      PROGRAM DPCLUSTR
      ALLOCATABLE::Q(:,:),XI(:,:,:),NXI(:,:,:,:),CLUS(:,:,:,:),Y(:)
      ALLOCATABLE::NR(:),NC(:),NA(:),NB(:),NU(:),NENDPT(:,:),SCR(:,:)
      COMMON/PARAM/NOPT,N,MOPT,NSUB,NROW,NCOL
      CHARACTER*80 FNAME,GNAME,INFORM
C*********************INPUT SETUP**************
      WRITE(*,*) 'MAIN MENU'
      WRITE(*,*) '1  AVERAGE PROXIMITY WITHIN SUBSET (USING THE'
      WRITE(*,*) '   NUMBER OF ROW PLUS COLUMN OBJECTS OBJECTS)'
      WRITE(*,*) '2  AVERAGE PROXIMITY WITHIN SUBSET'
      WRITE(*,*) '3  SUM OF PROXIMITIES WITHIN SUBSET'
      WRITE(*,*) '4  MAXIMUM PROXIMITY WITHIN SUBSET'
      WRITE(*,*) '5  AVERAGE WITHIN MINUS AVERAGE BETWEEN'
      WRITE(*,*) '6  NUMBER OF PROXIMITIES WITHIN GREATER THAN BETWEEN'
      WRITE(*,*) '7  NORMALIZED OPTION 6'
      WRITE(*,*) '8  COMMON OBJECT INCONSISTENCIES WITHIN VERSUS'
      WRITE(*,*) '   BETWEEN'
      WRITE(*,*) '9  NORMALIZED OPTION 8'
      WRITE(*,*) '10 NEGATIVE OF THE AVERAGE BETWEEN'
      WRITE(*,*) '11 MINIMUM CONNECTIVITY PROXIMITY'
      WRITE(*,*) '12 SPANNING TREE LENGTH'
      WRITE(*,*) 'CHOICE'
      READ(*,1) NOPT
      WRITE(*,552)
      WRITE(*,*) 'OPTIMIZATION METHOD'
      WRITE(*,*) ' 1 MINIMUM OF THE SUM'
      WRITE(*,*) ' 2 MINIMUM OF THE MAXIMUM'
      WRITE(*,*) 'CHOICE?'
      READ(*,1) MOPT
      WRITE(*,552)
      WRITE(*,*) 'MAXIMUM NUMBER OF SUBSETS TO BE USED'
      WRITE(*,*) 'NO LARGER THAN THE MINIMUM OF THE NUMBER OF'
      WRITE(*,*) 'ROWS AND COLUMNS'
      READ(*,1) NSUB
1     FORMAT(BN,I6)
552   FORMAT(1H )
      WRITE(*,*) 'NUMBER OF ROWS'
      READ(*,1) NROW
      WRITE(*,*) 'NUMBER OF COLUMNS'
      READ(*,1) NCOL
      N=NROW+NCOL
      IF(MIN(NROW,NCOL).LT.NSUB) WRITE(*,*) 'TOO MANY SUBSETS ASKED'
      IF(MIN(NROW,NCOL).LT.NSUB) STOP
      IF(NROW.LT.3) GO TO 650
      IF(NCOL.LT.3) GO TO 650
      GO TO 503
650   WRITE(*,*) 'MATRIX OUT OF RANGE'
      STOP
503   ALLOCATE(Q(NROW,NCOL),XI(NSUB,NROW,NCOL),NR(NROW),STAT=IERR)
      IF(IERR.NE.0) WRITE(*,*) 'ARRAY ALLOCATION ERROR'
      IF(IERR.NE.0) STOP
      IF(IERR.EQ.0) WRITE(*,*) 'ARRAYS ALLOCATED SUCCESSFULLY'
      ALLOCATE(CLUS(NROW,NROW,NCOL,NCOL),NC(NCOL),NA(NROW),STAT=IERR)
      IF(IERR.NE.0) WRITE(*,*) 'ARRAY ALLOCATION ERROR'
      IF(IERR.NE.0) STOP
      IF(IERR.EQ.0) WRITE(*,*) 'ARRAYS ALLOCATED SUCCESSFULLY'
      ALLOCATE(NXI(NSUB,NROW,NCOL,2),NU(N),NENDPT(2,N),Y(N),STAT=IERR)
      IF(IERR.NE.0) WRITE(*,*) 'ARRAY ALLOCATION ERROR'
      IF(IERR.NE.0) STOP
      IF(IERR.EQ.0) WRITE(*,*) 'ARRAYS ALLOCATED SUCCESSFULLY'
      ALLOCATE(NB(NCOL),SCR(N,N),STAT=IERR)
      IF(IERR.NE.0) WRITE(*,*) 'ARRAY ALLOCATION ERROR'
      IF(IERR.NE.0) STOP
      IF(IERR.EQ.0) WRITE(*,*) 'ARRAYS ALLOCATED SUCCESSFULLY'
      WRITE(*,*) 'FILE NAME FOR THE INPUT MATRIX'
      READ(*,'(A)') FNAME
      WRITE(*,552)
      WRITE(*,*) 'FILE NAME FOR THE OUTPUT MATRIX'
      READ(*,'(A)') GNAME
      WRITE(*,552)
      WRITE(*,*) 'PROVIDE FORMAT FOR THE INPUT MATRIX'
      READ(*,'(A)') INFORM
      OPEN(1,FILE=FNAME)
      OPEN(2,FILE=GNAME)
      WRITE(2,101) NOPT
101   FORMAT('HETEROGENEITY OPTION',I4)
      WRITE(2,552)
      IF(MOPT.EQ.1) WRITE(2,*) 'MINIMIZATION OF THE SUM'
      IF(MOPT.EQ.2) WRITE(2,*) 'MINIMIZATION OF THE MAXIMUM'
      WRITE(2,552)
      WRITE(*,*) '0 IF ROW CONSTRAINED ORDER IS TO BE THE IDENTITY'
      WRITE(*,*) '1 IF IT IS TO BE INPUT'
      READ(*,1) NCON
      DO 6001 I=1,NROW
      NR(I)=I
6001  CONTINUE
      IF(NCON.EQ.1) WRITE(*,*) 'ENTER ROW ORDER SEPARATED BY BLANKS'
      IF(NCON.EQ.1) READ(*,*) (NR(I),I=1,NROW)
      WRITE(2,*) 'CONSTRAINT ORDERING FOR ROWS'
      WRITE(2,601)(NR(I),I=1,NROW)
601   FORMAT(20I4)
      WRITE(2,552)
      WRITE(*,*) '0 IF COLUMN CONSTRANED ORDER IS TO BE THE IDENTITY'
      WRITE(*,*) '1 IF IT IS TO BE INPUT'
      READ(*,1) NCON
      DO 6020 I=1,NCOL
6020  NC(I)=I
      IF(NCON.EQ.1) WRITE(*,*)'ENTER COLUMN ORDER SEPARATED BY BLANKS'
      IF(NCON.EQ.1) READ(*,*) (NC(I),I=1,NCOL)
      WRITE(2,*) 'CONSTRAINT ORDERING FOR COLUMNS'
      WRITE(2,601) (NC(I),I=1,NCOL)
      WRITE(2,552)
      WRITE(2,*) 'INPUT MATRIX'
      DO 1002 I=1,NROW
      READ(1,INFORM) (Q(I,J),J=1,NCOL)
      WRITE(2,552)
      WRITE(2,100) (Q(I,J),J=1,NCOL)
100   FORMAT(10F8.3)
1002  CONTINUE
      WRITE(2,552)
      WRITE(*,*) 'PROGRAM RUNNING'
C*****************INITIALIZE ARRAYS************
      Z=1.0E+20
      NZ=0
      DO 1000 I=1,NSUB
      DO 1000 J=1,NROW
      DO 1000 L=1,NCOL
      XI(I,J,L) = Z
      NXI(I,J,L,1) = NZ
      NXI(I,J,L,2) = NZ
1000  CONTINUE
C*******************SINGLE SUBSET PROCESSING**********
      DO 200 NLBR=1,NROW
      DO 200 NUBR=NLBR,NROW
      DO 200 NLBC=1,NCOL
      DO 200 NUBC=NLBC,NCOL
      KR=NUBR-NLBR+1
      KC=NUBC-NLBC+1
      K=KC+KR
C****************SUBSET HETEROGENEITY CALCULATIONS***********
      SUM=0.0
      GO TO (13,13,13,753,900,900,900,900,900,900,600,600) NOPT
13    DO 11 I=NLBR,NUBR
      DO 11 J=NLBC,NUBC
11    SUM=SUM+Q(NR(I),NC(J))
      GO TO (751,752,755) NOPT
751   SUM=(1.0/FLOAT(K))*SUM
      GO TO 755
752   SUM=(1.0/FLOAT(KR*KC))*SUM
      GO TO 755
753   DO 756 I=NLBR,NUBR
      DO 756 J=NLBC,NUBC
      IF(SUM.LT.Q(NR(I),NC(J))) SUM=Q(NR(I),NC(J))
756   CONTINUE
      GO TO 755
900   JK=1
      DO 951 I=1,NROW
      DO 902 J=NLBR,NUBR
      IF(NR(J).EQ.I) GO TO 951
902   CONTINUE
      NA(JK)=I
      JK=JK+1
951   CONTINUE
      JK=1
      DO 201 I=1,NCOL
      DO 202 J=NLBC,NUBC
      IF(NC(J).EQ.I) GO TO 201
202   CONTINUE
      NB(JK)=I
      JK=JK+1
201   CONTINUE
      GO TO(613,602,602,903,903,613,613) NOPT-4
613   SUMW=0.0
      SUMB=0.0
      DO 904 I=NLBR,NUBR
      DO 904 J=NLBC,NUBC
904   SUMW=SUMW+Q(NR(I),NC(J))
590   IF(KC.EQ.NCOL) GO TO 210
      DO 905 I=NLBR,NUBR
      DO 905 J=1,NCOL-KC
905   SUMB=SUMB+Q(NR(I),NB(J))
210   IF(KR.EQ.NROW) GO TO 591
      DO 203 I=NLBC,NUBC
      DO 203 J=1,NROW-KR
203   SUMB=SUMB+Q(NA(J),NC(I))
591   IF(K.EQ.N) SUM=(1.0/FLOAT(KR*KC))*SUMW
      IF(K.EQ.N) GO TO 755
      SUM=(1.0/FLOAT(KR*KC))*SUMW-(1.0/FLOAT(NROW*NCOL-KR*KC))*SUMB
      IF(NOPT.EQ.10) SUM=-(1.0/FLOAT(NROW*NCOL-KR*KC))*SUMB
      GO TO 755
903   IF(KC.EQ.NCOL) GO TO 220
      DO 907 I=NLBR,NUBR
      DO 907 J=NLBC,NUBC
      DO 908 L=1,NCOL-KC
908   IF(Q(NR(I),NC(J)).GT.Q(NR(I),NB(L))) SUM=SUM+1.0
907   CONTINUE
220   IF(KR.EQ.NROW) GO TO 204
      DO 205 I=NLBC,NUBC
      DO 205 J=NLBR,NUBR
      DO 206 L=1,NROW-KR
206   IF(Q(NR(J),NC(I)).GT.Q(NA(L),NC(I))) SUM=SUM+1.0
205   CONTINUE
204   IF(K.EQ.N) GO TO 755
      NN=KC*(NCOL-KC)*KR+KR*(NROW-KR)*KC
      IF(NOPT.EQ.9) SUM=(1.0/FLOAT(NN))*SUM
      GO TO 755
602   IF(KC.EQ.NCOL) GO TO 211
      DO 603 I=NLBR,NUBR
      DO 603 J=NLBC,NUBC
      DO 612 I1=NLBR,NUBR
      DO 612 J1=1,NCOL-KC
      IF(Q(NR(I),NC(J)).GT.Q(NR(I1),NB(J1))) SUM=SUM+1.0
612   CONTINUE
603   CONTINUE
211   IF(KR.EQ.NROW) GO TO 212
      DO 207 I=NLBR,NUBR
      DO 207 J=NLBC,NUBC
      DO 208 I1=NLBC,NUBC
      DO 208 J1=1,NROW-KR
      IF(Q(NR(I),NC(J)).GT.Q(NA(J1),NC(I1))) SUM=SUM+1.0
208   CONTINUE
207   CONTINUE
212   IF(K.EQ.N) GO TO 755
      IF(NOPT.EQ.7) SUM=(1.0/FLOAT(KR*KC*((NROW-KR)+(NCOL-KC))))*SUM
      GO TO 755
600   IF(K.EQ.1) GO TO 755
      DO 250 I=1,K
      SCR(I,I)=Z
      DO 250 J=I+1,K
      SCR(I,J)=Z
      SCR(J,I)=Z
      IF((I.LE.KR).AND.(J.LE.KR)) GO TO 250
      IF((I.GT.KR).AND.(J.GT.KR)) GO TO 250
      SCR(I,J)=Q(NR(NLBR+I-1),NC(NLBC+J-KR-1))
      SCR(J,I)=SCR(I,J)
250   CONTINUE
      L=0
      DO 2100 I=2,K
      NU(I)=1
2100  Y(I)=SCR(1,I)
3000  DMIN=1.0E+20
      DO 4100 I=2,K
      IF(Y(I).EQ.(-Z)) GO TO 4100
      IF(Y(I).GE.DMIN) GO TO 4100
      DMIN=Y(I)
      IMIN=I
4100  CONTINUE
      L=L+1
      NENDPT(1,L)=IMIN
      NENDPT(2,L)=NU(IMIN)
      IF(L.EQ.(K-1)) GO TO 604
      Y(IMIN)=-Z
      DO 1110 I=2,K
      IF(Y(I).EQ.(-Z)) GO TO 1110
      D1=SCR(I,IMIN)
      IF(Y(I).LE.D1) GO TO 1110
      NU(I)=IMIN
      Y(I)=D1
1110  CONTINUE
      GO TO 3000
604   SUMMAX=0.0
      DO 605 I=1,K-1
      I1=NENDPT(1,I)
      I2=NENDPT(2,I)
      SUM=SUM+SCR(I1,I2)
      IF(SCR(I1,I2).GT.SUMMAX) SUMMAX=SCR(I1,I2)
605   CONTINUE
      IF(NOPT.EQ.11) SUM=SUMMAX
      GO TO 755
755   IF((NLBR.EQ.1).AND.(NLBC.EQ.1)) XI(1,NUBR,NUBC)=SUM
      IF((NLBR.EQ.1).AND.(NLBC.EQ.1)) NXI(1,NUBR,NUBC,1) = 1
      IF((NLBR.EQ.1).AND.(NLBC.EQ.1)) NXI(1,NUBR,NUBC,2) = 1
      CLUS(NLBR,NUBR,NLBC,NUBC)=SUM
200   CONTINUE
C****************PROCESSING THE RECURSION*************
606   DO 608 K=2,NSUB
      WRITE(*,*) 'PROCESSING A STAGE'
      DO 300 I=K-1,NROW-1
      DO 300 I1=I+1,NROW
      DO 300 J=K-1,NCOL-1
      DO 300 J1=J+1,NCOL
      IF(MOPT.EQ.1) TEMP=XI(K-1,I,J)+CLUS(I+1,I1,J+1,J1)
      IF(MOPT.EQ.2) TEMP=MAX(XI(K-1,I,J),CLUS(I+1,I1,J+1,J1))
      IF(TEMP.GE.XI(K,I1,J1)) GO TO 300
      XI(K,I1,J1) = TEMP
      NXI(K,I1,J1,1) = I+1
      NXI(K,I1,J1,2) = J+1
300   CONTINUE
608   CONTINUE
C*********************PRODUCE OUTPUT*************
      DO 320 K=1,NSUB
      WRITE(2,152)
      WRITE(2,152)
152   FORMAT(1H )
      WRITE(2,151) K
151   FORMAT('NUMBER OF SUBSETS', I4)
      NENDR=NROW
      NENDC=NCOL
      NK=K
319   WRITE(2,150) XI(NK,NENDR,NENDC)
150   FORMAT('CUMULATIVE INDEX',F13.4)
      WRITE(2,*) 'SUBSET ADDED: ROW OBJECTS'
      WRITE(2,149) (NR(J),J=NXI(NK,NENDR,NENDC,1),NENDR)
149   FORMAT(25I3)
      WRITE(2,*) '              COLUMN OBJECTS'
      WRITE(2,149) (NC(J),J=NXI(NK,NENDR,NENDC,2),NENDC)
      T= CLUS(NXI(NK,NENDR,NENDC,1),NENDR,NXI(NK,NENDR,NENDC,2),NENDC)
      WRITE(2,148) T
148   FORMAT('WITH INDEX VALUE', F13.4)
      WRITE(2,152)
      NT1=NXI(NK,NENDR,NENDC,1)
      NT2=NXI(NK,NENDR,NENDC,2)
      IF((NT1.EQ.1).AND.(NT2.EQ.1)) GO TO 320
      NENDR=NT1-1
      NENDC=NT2-1
      NK=NK-1
      GO TO 319
320   CONTINUE
      STOP
      END
