      PROGRAM DPCLUSRS
      ALLOCATABLE::Q(:,:),NA(:),XI(:,:),NXI(:,:),CLUS(:,:)
      ALLOCATABLE::NB(:),NU(:),NENDPT(:,:),Y(:)
      CHARACTER*80 FNAME,GNAME,INFORM
C****************INPUT SETUP************
      WRITE(*,*) 'MAIN MENU'
      WRITE(*,*) '1  AVERAGE PROXIMITY WITHIN SUBSET (USING TWICE'
      WRITE(*,*) '   THE NUMBER OF OBJECTS)'
      WRITE(*,*) '2  AVERAGE PROXIMITY WITHIN SUBSET (USING THE'
      WRITE(*,*) '   NUMBER OF OBJECT PAIRS)'
      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'
      READ(*,1) NSUB
1     FORMAT(BN,I6)
552   FORMAT(1H )
      WRITE(*,*) 'NUMBER OF ROWS/COLUMNS FOR THE LOWER TRIANGULAR'
      WRITE(*,*) 'MATRIX-- NO DIAGONAL'
      READ(*,1) N
      IERR=0
      IF(N.GE.3) GO TO 503
      WRITE(*,*) 'MATRIX OUT OF RANGE'
      STOP
503   ALLOCATE(Q(N,N),NA(N),XI(NSUB,N),NXI(NSUB,N),CLUS(N,N),STAT=IERR)
      IF(IERR.NE.0) WRITE(*,*) 'ARRAY ALLOCATION ERROR'
      IF(IERR.EQ.0) WRITE(*,*) 'ARRAY ALLOCATION SUCCESSFUL'
      IF(IERR.NE.0) STOP
      ALLOCATE (NB(N),NU(N),NENDPT(2,N),Y(N),STAT=IERR)
      IF(IERR.NE.0) WRITE(*,*) 'ARRAY ALLOCATION ERROR'
      IF(IERR.EQ.0) WRITE(*,*) 'ARRAY ALLOCATION SUCCESSFUL'
      IF(IERR.NE.0) STOP
      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 CONSTRAINED ORDER IS TO BE THE IDENTITY'
      WRITE(*,*) '1 IF IT IS TO BE INPUT'
      READ(*,1) NCON
      DO 6001 I=1,N
      NA(I)=I
6001  CONTINUE
      IF(NCON.EQ.1) WRITE(*,*) 'ENTER ORDER SEPARATED BY BLANKS'
      IF(NCON.EQ.1) READ(*,*) (NA(I),I=1,N)
      WRITE(2,*) 'CONSTRAINT ORDERING'
      WRITE(2,601)(NA(I),I=1,N)
601   FORMAT(20I4)
      WRITE(2,552)
      WRITE(2,*) 'INPUT MATRIX'
      DO 1002 I=2,N
      ISTOP=I-1
      READ(1,INFORM) (Q(I,J),J=1,ISTOP)
      WRITE(2,552)
      WRITE(2,100) (Q(I,J),J=1,ISTOP)
100   FORMAT(10F8.3)
1002  CONTINUE
      DO 1004 I=1,N
      Q(I,I)=0.0
      DO 1004 J=I+1,N
      Q(I,J)=Q(J,I)
1004  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,N
      XI(I,J) = Z
      NXI(I,J) = NZ
1000  CONTINUE
C**************INITIAL SINGLE SUBSET PROCESSING***********
      DO 200 NLB=1,N
      DO 200 NUB=NLB,N
      K=NUB-NLB+1
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=NLB,NUB
      DO 11 J=NLB,NUB
11    SUM=SUM+Q(NA(I),NA(J))
      GO TO (751,752,755) NOPT
751   SUM=(0.5/FLOAT(K))*SUM
      GO TO 755
752   IF(K.EQ.1) GO TO 755
      SUM=(1.0/FLOAT(K*(K-1)))*SUM
      GO TO 755
753   DO 756 I=NLB,NUB
      DO 756 J=NLB,NUB
      IF(SUM.LT.Q(NA(I),NA(J))) SUM=Q(NA(I),NA(J))
756   CONTINUE
      GO TO 755
900   JK=1
      DO 951 I=1,N
      DO 902 J=NLB,NUB
      IF(NA(J).EQ.I) GO TO 951
902   CONTINUE
      NB(JK)=I
      JK=JK+1
951   CONTINUE
      GO TO(613,602,602,903,903,613,613) NOPT-4
613   SUMW=0.0
      SUMB=0.0
      IF(K.EQ.1) GO TO 590
      DO 904 I=NLB,NUB-1
      DO 904 J=I+1,NUB
904   SUMW=SUMW+Q(NA(I),NA(J))
590   IF(K.EQ.N) GO TO 591
      DO 905 I=NLB,NUB
      DO 905 J=1,N-K
905   SUMB=SUMB+Q(NA(I),NB(J))
591   IF(K.EQ.1) SUM=-(1.0/FLOAT(K*(N-K)))*SUMB
      IF(K.EQ.1) GO TO 755
      IF(K.EQ.N) SUM=(2.0/FLOAT(K*(K-1)))*SUMW
      IF(K.EQ.N) GO TO 755
      SUM=(2.0/FLOAT(K*(K-1)))*SUMW-(1.0/FLOAT(K*(N-K)))*SUMB
      IF(NOPT.EQ.10) SUM=-(1.0/FLOAT(K*(N-K)))*SUMB
      GO TO 755
903   IF((K.EQ.1).OR.(K.EQ.N)) GO TO 755
      DO 907 I=NLB,NUB
      DO 907 J=NLB,NUB
      IF(I.EQ.J) GO TO 907
      DO 908 L=1,N-K
908   IF(Q(NA(I),NA(J)).GT.Q(NA(I),NB(L))) SUM=SUM+1.0
907   CONTINUE
      SUM=2.0*SUM
      IF(NOPT.EQ.9) SUM=(1.0/FLOAT(2*K*(K-1)*(N-K)))*SUM
      GO TO 755
602   IF((K.EQ.1).OR.(K.EQ.N)) GO TO 755
      DO 603 I=NLB,NUB
      DO 603 J=NLB,NUB
      IF(I.EQ.J) GO TO 603
      DO 612 I1=NLB,NUB
      DO 612 J1=1,N-K
      IF(Q(NA(I),NA(J)).GT.Q(NA(I1),NB(J1))) SUM=SUM+1.0
612   CONTINUE
603   CONTINUE
      SUM=2.0*SUM
      IF(NOPT.EQ.7) SUM=(1.0/FLOAT(2*K*(K-1)*K*(N-K)))*SUM
      GO TO 755
600   IF(K.EQ.1) GO TO 755
      IF(K.EQ.2) SUM=Q(NA(NLB),NA(NUB))
      IF(K.EQ.2) GO TO 755
      L=0
      DO 2100 I=NLB+1,NUB
      NU(I)=NLB
2100  Y(I)=Q(NA(NLB),NA(I))
3000  DMIN=1.0E+20
      DO 4100 I=NLB+1,NUB
      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=NLB+1,NUB
      IF(Y(I).EQ.(-Z)) GO TO 1110
      D1=Q(NA(I),NA(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+Q(NA(I1),NA(I2))
      IF(Q(NA(I1),NA(I2)).GT.SUMMAX) SUMMAX=Q(NA(I1),NA(I2))
605   CONTINUE
      IF(NOPT.EQ.11) SUM=SUMMAX
      GO TO 755
755   IF(NLB.EQ.1) XI(1,NUB)=SUM
      IF(NLB.EQ.1) NXI(1,NUB) = 1
      CLUS(NLB,NUB)=SUM
200   CONTINUE
C**********************PROCESSING THE RECURSION*************
606   DO 608 K=2,NSUB
      WRITE(*,*) 'PROCESSING A STAGE'
      DO 300 I=K-1,N-1
      DO 300 J1=I+1,N
      IF(MOPT.EQ.1) TEMP=XI(K-1,I)+CLUS(I+1,J1)
      IF(MOPT.EQ.2) TEMP=MAX(XI(K-1,I),CLUS(I+1,J1))
      IF(TEMP.GE.XI(K,J1)) GO TO 300
      XI(K,J1) = TEMP
      NXI(K,J1) = I+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)
      NEND=N
      NK=K
319   WRITE(2,150) XI(NK,NEND)
150   FORMAT('CUMULATIVE INDEX',F13.4)
      WRITE(2,*) 'SUBSET ADDED'
      WRITE(2,149) (NA(J),J=NXI(NK,NEND),NEND)
149   FORMAT(25I3)
      WRITE(2,148) CLUS(NXI(NK,NEND),NEND)
148   FORMAT('WITH INDEX VALUE', F13.4)
      WRITE(2,152)
      IF(NXI(NK,NEND).EQ.1) GO TO 320
      NEND=NXI(NK,NEND)-1
      NK=NK-1
      GO TO 319
320   CONTINUE
      STOP
      END
