      MODULE ARRAYS
      SAVE
      ALLOCATABLE::Q(:,:),XI(:,:),Y(:)
      ALLOCATABLE::NA(:),NB(:),NC(:),NXI(:,:),II(:),JJ(:)
      END MODULE ARRAYS
      PROGRAM DPCLUSUN
      USE ARRAYS
      COMMON/PARAM/NOPT,N,MOPT,NSUB,K
      CHARACTER*80 FNAME,GNAME,INFORM
C*************INPUT SETUP*******************
      WRITE(*,*) 'MAIN MENU'
      WRITE(*,*) '1 UNWEIGHTED GRADIENT WITHIN ROWS AND COLUMNS'
      WRITE(*,*) '2 WEIGHTED GRADIENT WITHIN ROWS AND COLUMNS'
      WRITE(*,*) '3 DEFAYS CRITERION'
      WRITE(*,*) '4 PROXIMITY SUM'
      WRITE(*,*) '5 SKEW-SYMMETRIC COORDINATES'
      WRITE(*,*) 'CHOICE'
      READ(*,1) NOPT
      WRITE(*,552)
      IERR=0
      WRITE(*,*) 'MAXIMUM NUMBER OF SUBSETS TO BE USED'
      READ(*,1) NSUB
1     FORMAT(BN,I6)
552   FORMAT(1H )
      WRITE(*,*) 'NUMBER OF ROWS/COLUMNS'
      READ(*,1) N
      IF((N.GE.3).AND.(N.LE.30)) GO TO 503
      WRITE(*,*) 'PROBLEM SIZE IS OUT OF RANGE'
      STOP
503   CALL CONFIGURE
      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(*,*) ' IS THE INPUT MATRIX SYMMETRIC AND READ AS A LOWER'
      WRITE(*,*) ' TRIANGULAR MATRIX WITHOUT DIAGONAL ENTRIES OR AS'
      WRITE(*,*) ' A COMPLETE SQUARE MATRIX '
      WRITE(*,*) '    1 LOWER TRIANGULAR '
      WRITE(*,*) '    2 COMPLETE SQUARE '
      WRITE(*,*) '    3 COMPLETE SQUARE AND SKEW-SYMMETRIC'
      WRITE(*,*) '    4 COMPLETE SQUARE AND SKEW-SYMMETRIC'
      WRITE(*,*) '      BUT USE SIGN INFORMATION ONLY'
      WRITE(*,*) ' CHOICE? '
      READ(*,1) NSYM
      WRITE(*,552)
504   WRITE(*,*) ' PROVIDE FORMAT FOR THE INPUT MATRIX '
      READ(*,'(A)') INFORM
      WRITE(*,552)
      WRITE(*,551)
      WRITE(*,*) INFORM
551   FORMAT(' INPUT FORMAT IS ')
      OPEN(1,FILE=FNAME)
      OPEN(2,FILE=GNAME,STATUS='NEW')
      WRITE(2,7799) NOPT
7799  FORMAT(' OPTIMIZATION MENU OPTION ',I4)
      WRITE(2,552)
      WRITE(2,*) 'INPUT MATRIX '
      IF(NSYM.GE.2) GO TO 1001
      DO 1002 I=2,N
      ISTOP=I-1
      READ(1,INFORM) (Q(I,J),J=1,ISTOP)
      WRITE(2,552)
      WRITE(2,6000) (Q(I,J),J=1,ISTOP)
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)
      GO TO 7779
1001  DO 1006 I=1,N
      READ(1,INFORM) (Q(I,J),J=1,N)
      WRITE(2,552)
      WRITE(2,6000) (Q(I,J),J=1,N)
      Q(I,I)=0.0
1006  CONTINUE
6000  FORMAT(2X,6F13.3)
      IF(NSYM.NE.4) GO TO 7779
      DO 9020 I=1,N
      DO 9020 J=1,N
      IF(Q(I,J).GT.0.0) Q(I,J)=1.0
      IF(Q(I,J).LT.0.0) Q(I,J)=-1.0
9020  CONTINUE
7779  WRITE(2,552)
      WRITE(*,*) 'PROGRAM RUNNING'
      CALL DYMAIN
      STOP
      END
      SUBROUTINE CONFIGURE
      USE ARRAYS
      COMMON/PARAM/NOPT,N,MOPT,NSUB,K
      ALLOCATE (Q(N,N),XI(NSUB,2**N),Y(N), NXI(NSUB,2**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 (NA(N),NB(N),NC(N),II(N),JJ(N),STAT=IERR)
      IF(IERR.NE.0) WRITE(*,*) 'ARRAY ALLOCATION ERROR'
      IF(IERR.NE.0) STOP
      IF(IERR.EQ.0) WRITE(*,*) 'ARRAY ALLOCATION SUCCESSFUL'
      RETURN
      END SUBROUTINE CONFIGURE
      SUBROUTINE DYMAIN
      USE ARRAYS
      COMMON/PARAM/NOPT,N,MOPT,NSUB,K
      LOGICAL MTC
C****************INITIALIZE ARRAYS*************
      Z=0.0
      NZ=0
      NSUM=(2**N)-1
      DO 1000 I=1,NSUB
      DO 1000 J=1,NSUM
      XI(I,J) = Z
      NXI(I,J) = NZ
1000  CONTINUE
C************INITIAL SINGLE SUBSET PROCESSING*********
      K=0
4     K=K+1
      IF(K.GT.N) GO TO 5
      WRITE(*,3012) K
3012  FORMAT('STAGE ONE PROCESSING SUBSETS OF SIZE',I4)
      MTC=.FALSE.  
6     GO TO 30
30    IF(MTC) GO TO 40
      M2 = 0
      NH=K
      GO TO 50
40    IF(M2.LT.(N-NH)) NH=0
      NH=NH+1
      M2=NA(K+1-NH)
50    DO 51 J=1,NH
51    NA(K+J-NH) = M2+J
      MTC=NA(1).NE.(N-K+1)
      INDEX=0
      DO 10 I=1,K
10    INDEX=INDEX+2**(NA(I)-1)
C*************SUBSET HETEROGENEITY CALCULATIONS***********
      GO TO (6050,6050,6001,6050,6001) NOPT
6050  XI(1,INDEX)=0.0
      NXI(1,INDEX)=INDEX
      IF(MTC) GO TO 6
      GO TO 4
6001  JK=1
      DO 951 I=1,N
      DO 902 J=1,K
      IF(NA(J).EQ.I) GO TO 951
902   CONTINUE
      NB(JK)=I
      JK=JK+1
951   CONTINUE
      IF((N-K).EQ.0) GO TO 6050
      SUM=0.0
      DO 6002 I=1,K
      DO 6002 J=1,N-K
6002  SUM=SUM+Q(NA(I),NB(J))
      SUM=(SUM*SUM)/FLOAT(K)
      XI(1,INDEX)=SUM
      NXI(1,INDEX) = INDEX
      IF(MTC) GO TO 6
      GO TO 4
C*************PROCESSING THE RECURSION***********
5     DO 610 K=2,NSUB
      WRITE(*,*) 'PROCESSING A STAGE',K
      DO 60 II1=0,1
      DO 60 JJ1=0,MOD(II1+1,2)
      II(1)=II1
      JJ(1)=JJ1
      DO 61 II2=0,1
      DO 61 JJ2=0,MOD(II2+1,2)
      II(2)=II2
      JJ(2)=JJ2
      DO 62 II3=0,1
      DO 62 JJ3=0,MOD(II3+1,2)
      II(3)=II3
      JJ(3)=JJ3
      IF(N.EQ.3) CALL EVAL
      IF(N.EQ.3) GO TO 62
      DO 63 II4=0,1
      DO 63 JJ4=0,MOD(II4+1,2)
      II(4)=II4
      JJ(4)=JJ4
      IF(N.EQ.4) CALL EVAL
      IF(N.EQ.4) GO TO 63
      DO 64 II5=0,1
      DO 64 JJ5=0,MOD(II5+1,2)
      II(5)=II5
      JJ(5)=JJ5
      IF(N.EQ.5) CALL EVAL
      IF(N.EQ.5) GO TO 64
      DO 65 II6=0,1
      DO 65 JJ6=0,MOD(II6+1,2)
      II(6)=II6
      JJ(6)=JJ6
      IF(N.EQ.6) CALL EVAL
      IF(N.EQ.6) GO TO 65
      DO 66 II7=0,1
      DO 66 JJ7=0,MOD(II7+1,2)
      II(7)=II7
      JJ(7)=JJ7
      IF(N.EQ.7) CALL EVAL
      IF(N.EQ.7) GO TO 66
      DO 67 II8=0,1
      DO 67 JJ8=0,MOD(II8+1,2)
      II(8)=II8
      JJ(8)=JJ8
      IF(N.EQ.8) CALL EVAL
      IF(N.EQ.8) GO TO 67
      DO 68 II9=0,1
      DO 68 JJ9=0,MOD(II9+1,2)
      II(9)=II9
      JJ(9)=JJ9
      IF(N.EQ.9) CALL EVAL
      IF(N.EQ.9) GO TO 68
      DO 69 II10=0,1
      DO 69 JJ10=0,MOD(II10+1,2)
      II(10)=II10
      JJ(10)=JJ10
      IF(N.EQ.10) CALL EVAL
      IF(N.EQ.10) GO TO 69
      DO 70 II11=0,1
      DO 70 JJ11=0,MOD(II11+1,2)
      II(11)=II11
      JJ(11)=JJ11
      IF(N.EQ.11) CALL EVAL
      IF(N.EQ.11) GO TO 70
      DO 71 II12=0,1
      DO 71 JJ12=0,MOD(II12+1,2)
      II(12)=II12
      JJ(12)=JJ12
      IF(N.EQ.12) CALL EVAL
      IF(N.EQ.12) GO TO 71
      DO 72 II13=0,1
      DO 72 JJ13=0,MOD(II13+1,2)
      II(13)=II13
      JJ(13)=JJ13
      IF(N.EQ.13) CALL EVAL
      IF(N.EQ.13) GO TO 72
      DO 73 II14=0,1
      DO 73 JJ14=0,MOD(II14+1,2)
      II(14)=II14
      JJ(14)=JJ14
      IF(N.EQ.14) CALL EVAL
      IF(N.EQ.14) GO TO 73
      DO 74 II15=0,1
      DO 74 JJ15=0,MOD(II15+1,2)
      II(15)=II15
      JJ(15)=JJ15
      IF(N.EQ.15) CALL EVAL
      IF(N.EQ.15) GO TO 74
      DO 75 II16=0,1
      DO 75 JJ16=0,MOD(II16+1,2)
      II(16)=II16
      JJ(16)=JJ16
      IF(N.EQ.16) CALL EVAL
      IF(N.EQ.16) GO TO 75
      DO 76 II17=0,1
      DO 76 JJ17=0,MOD(II17+1,2)
      II(17)=II17
      JJ(17)=JJ17
      IF(N.EQ.17) CALL EVAL
      IF(N.EQ.17) GO TO 76
      DO 77 II18=0,1
      DO 77 JJ18=0,MOD(II18+1,2)
      II(18)=II18
      JJ(18)=JJ18
      IF(N.EQ.18) CALL EVAL
      IF(N.EQ.18) GO TO 77
      DO 78 II19=0,1
      DO 78 JJ19=0,MOD(II19+1,2)
      II(19)=II19
      JJ(19)=JJ19
      IF(N.EQ.19) CALL EVAL
      IF(N.EQ.19) GO TO 78
      DO 79 II20=0,1
      DO 79 JJ20=0,MOD(II20+1,2)
      II(20)=II20
      JJ(20)=JJ20
      IF(N.EQ.20) CALL EVAL
      IF(N.EQ.20) GO TO 79
      DO 80 II21=0,1
      DO 80 JJ21=0,MOD(II21+1,21)
      II(21)=II21
      JJ(21)=JJ21
      IF(N.EQ.21) CALL EVAL
      IF(N.EQ.21) GO TO 80
      DO 81 II22=0,1
      DO 81 JJ22=0,MOD(II22+1,2)
      II(22)=II22
      JJ(22)=JJ22
      IF(N.EQ.22) CALL EVAL
      IF(N.EQ.22) GO TO 81
      DO 82 II23=0,1
      DO 82 JJ23=0,MOD(II23+1,2)
      II(23)=II23
      JJ(23)=JJ23
      IF(N.EQ.23) CALL EVAL
      IF(N.EQ.23) GO TO 82
      DO 83 II24=0,1
      DO 83 JJ24=0,MOD(II24+1,2)
      II(24)=II24
      JJ(24)=JJ24
      IF(N.EQ.24) CALL EVAL
      IF(N.EQ.24) GO TO 83
      DO 84 II25=0,1
      DO 84 JJ25=0,MOD(II25+1,2)
      II(25)=II25
      JJ(25)=JJ25
      IF(N.EQ.25) CALL EVAL
      IF(N.EQ.25) GO TO 84
      DO 85 II26=0,1
      DO 85 JJ26=0,MOD(II26+1,2)
      II(26)=II26
      JJ(26)=JJ26
      IF(N.EQ.26) CALL EVAL
      IF(N.EQ.26) GO TO 85
      DO 86 II27=0,1
      DO 86 JJ27=0,MOD(II27+1,2)
      II(27)=II27
      JJ(27)=JJ27
      IF(N.EQ.27) CALL EVAL
      IF(N.EQ.27) GO TO 86
      DO 87 II28=0,1
      DO 87 JJ28=0,MOD(II28+1,2)
      II(28)=II28
      JJ(28)=JJ28
      IF(N.EQ.28) CALL EVAL
      IF(N.EQ.28) GO TO 87
      DO 88 II29=0,1
      DO 88 JJ29=0,MOD(II29+1,2)
      II(29)=II29
      JJ(29)=JJ(29)
      IF(N.EQ.29) CALL EVAL
      IF(N.EQ.29) GO TO 88
      DO 89 II30=0,1
      DO 89 JJ30=0,MOD(II30+1,2)
      II(30)=II30
      JJ(30)=JJ30
      IF(N.EQ.30) CALL EVAL
      IF(N.EQ.30) GO TO 89
      STOP
89    CONTINUE
88    CONTINUE
87    CONTINUE
86    CONTINUE
85    CONTINUE
84    CONTINUE
83    CONTINUE
82    CONTINUE
81    CONTINUE
80    CONTINUE
79    CONTINUE
78    CONTINUE
77    CONTINUE
76    CONTINUE
75    CONTINUE
74    CONTINUE
73    CONTINUE
72    CONTINUE
71    CONTINUE
70    CONTINUE
69    CONTINUE
68    CONTINUE
67    CONTINUE
66    CONTINUE
65    CONTINUE
64    CONTINUE
63    CONTINUE
62    CONTINUE
61    CONTINUE
60    CONTINUE
610   CONTINUE
C**************PRODUCE OUTPUT***********
100   DO 150 K=1,NSUB
      WRITE(2,152)
      WRITE(2,152)
      WRITE(2,151) K
151   FORMAT( 'NUMBER OF SUBSETS', I4)
152   FORMAT(1H )
      INDEXA=NSUM
      DO 153 KK=1,K
      NSTAGE=K-KK+1
      INDEXB=NXI(NSTAGE,INDEXA)
      DO 780 I=1,N
780   NA(I)=0
      NSPOT=1
      DO 154 I=1,N
      NTEST=2**(N-I)
      IF(NTEST.GT.INDEXB) GO TO 154
      NPRES=N-I+1
      NA(NSPOT)=NPRES
      NSPOT=NSPOT+1
155   FORMAT(I4)
      INDEXB=INDEXB-2**(N-I)
154   CONTINUE
      WRITE(2,152)
      WRITE(2,781) (NA(I),I=1,NSPOT-1)
      IF((INDEXA-NXI(NSTAGE,INDEXA)).EQ.0) CHANGE=XI(NSTAGE,INDEXA)
      IF((INDEXA-NXI(NSTAGE,INDEXA)).GT.0)
     &CHANGE=XI(NSTAGE,INDEXA)-XI(NSTAGE-1,INDEXA-NXI(NSTAGE,INDEXA))
781   FORMAT('SUBSET ADDED',25I3)
      WRITE(2,782) CHANGE,XI(NSTAGE,INDEXA)
782   FORMAT('INDEX FOR SUBSET',F13.4,'   CUMULATIVE INDEX',F13.4)
      INDEXA=INDEXA-NXI(NSTAGE,INDEXA)
      DO 6025 I=1,NSPOT-1
      DO 6026 J=1,N
      IF(NA(I).EQ.J) NB(J)=NSTAGE
6026  CONTINUE
6025  CONTINUE
      NC(NSTAGE)=NSPOT-1
153   CONTINUE
      WRITE(2,152)
      WRITE(2,*) 'SUMMARY SUBSET MEMBERSHIP FOR THE OBJECTS'
      WRITE(2,6027) (NB(I),I=1,N)
6027  FORMAT(20I3)
      WRITE(2,152)
      IF(((NOPT.EQ.1).OR.(NOPT.EQ.2)).OR.(NOPT.EQ.4)) GO TO 150
      DO 6028 I=1,K
6028  Y(I)=0.0
      DO 6029 I=1,K
      DO 6030 I1=1,N
      DO 6030 J1=1,N
      IF(NB(I1).NE.I) GO TO 6030
      IF(NB(J1).LT.I) Y(I)=Y(I)+Q(I1,J1)
      IF((NB(J1).GT.I).AND.(NOPT.EQ.5)) Y(I)=Y(I)+Q(I1,J1)
      IF((NB(J1).GT.I).AND.(NOPT.EQ.3)) Y(I)=Y(I)-Q(I1,J1)
6030  CONTINUE
6029  CONTINUE
      DO 6031 I=1,K
6031  Y(I)=Y(I)/FLOAT(N*NC(I))
      WRITE(2,*) 'COORDINATES FOR THE SUBSETS'
      WRITE(2,6032) (Y(I),I=1,K)
6032  FORMAT(10F8.3)
      WRITE(2,152)
      SUM=0.0
      DO 6033 I=1,N
      DO 6033 J=1,N
      IF(I.EQ.J) GO TO 6033
      IF(NOPT.EQ.5) SUM=SUM+(Q(I,J)-Y(NB(I))+Y(NB(J)))**2
      IF(NOPT.EQ.3) SUM=SUM+(Q(I,J)-ABS(Y(NB(I))-Y(NB(J))))**2
6033  CONTINUE
      WRITE(2,6034) SUM
6034  FORMAT('RESIDUAL SUM OF SQUARES',F15.3)
      WRITE(2,152)
150   CONTINUE
      RETURN
      END
      SUBROUTINE EVAL
      USE ARRAYS
      COMMON/PARAM/NOPT,N,MOPT,NSUB,K
      INDEXA=0
      INDEXB=0
      IA=0
      IB=0
      IC=0
      DO 10401 I=1,N
      IF(II(I).EQ.1) IA=IA+1
      IF(II(I).EQ.1) NA(IA)=I
      IF(JJ(I).EQ.1) IB=IB+1
      IF(JJ(I).EQ.1) NB(IB)=I
      IF((II(I).EQ.0).AND.(JJ(I).EQ.0)) IC=IC+1
      IF((II(I).EQ.0).AND.(JJ(I).EQ.0)) NC(IC)=I
10401 CONTINUE
      IF((IA.LT.(K-1)).OR.(IB.EQ.0)) RETURN
      SUM=0.0
      GO TO (6011,6012,6013,6014,6013) NOPT
6014  DO 6015 I=1,IA
      DO 6015 J=1,IB
      SUM=SUM+Q(NA(I),NB(J))
6015  CONTINUE
      GO TO 6016
6013  SUMA=0.0
      SUMB=0.0
      DO 6017 I=1,IB
      DO 6017 J=1,IA
      SUMA=SUMA+Q(NB(I),NA(J))
6017  CONTINUE
      IF(IC.EQ.0) GO TO 6019
      DO 6018 I=1,IB
      DO 6018 J=1,IC
      SUMB=SUMB+Q(NB(I),NC(J))
6018  CONTINUE
6019  IF(NOPT.EQ.5) SUM=((SUMA+SUMB)**2)/FLOAT(IB)
      IF(NOPT.EQ.3) SUM=((SUMA-SUMB)**2)/FLOAT(IB)
      GO TO 6016
6012  SUMA=0.0
      SUMB=0.0
      IF(IC.EQ.0) GO TO 6021
      DO 6020 I=1,IA
      DO 6020 J=1,IB
      DO 6020 KKK=1,IC
      SUMA=SUMA+Q(NA(I),NC(KKK))-Q(NA(I),NB(J))
      SUMB=SUMB+Q(NA(I),NC(KKK))-Q(NB(J),NC(KKK))
6020  CONTINUE
6021  SUM=SUMA+SUMB
      GO TO 6016
6011  SUMA=0.0
      SUMB=0.0
      IF(IC.EQ.0) GO TO 6023
      DO 6022 I=1,IA
      DO 6022 J=1,IB
      DO 6022 KKK=1,IC
      IF((Q(NA(I),NC(KKK))-Q(NA(I),NB(J))).GT.(0.0)) SUMA=SUMA+1.0
      IF((Q(NA(I),NC(KKK))-Q(NA(I),NB(J))).LT.(0.0)) SUMA=SUMA-1.0
      IF((Q(NA(I),NC(KKK))-Q(NB(J),NC(KKK))).GT.(0.0)) SUMB=SUMB+1.0
      IF((Q(NA(I),NC(KKK))-Q(NB(J),NC(KKK))).LT.(0.0)) SUMB=SUMB-1.0
6022  CONTINUE
6023  SUM=SUMA+SUMB
      GO TO 6016
6016  DO 10301 I=1,N
      INDEXA=INDEXA+(2**(I-1))*II(I)
      INDEXB=INDEXB+(2**(I-1))*JJ(I)
10301 CONTINUE
      INDEXAB=INDEXA+INDEXB
      TEMP=XI(K-1,INDEXA)+SUM
      IF(TEMP.LT.XI(K,INDEXAB)) RETURN
      XI(K,INDEXAB)=TEMP
      NXI(K,INDEXAB)=INDEXB
      RETURN
      END
