      PROGRAM XCLUST
C***********************************************************************
C*                                                                     *
C*  FORTRAN CODE WRITTEN FOR INCLUSION IN IBM RESEARCH REPORT RC20525, *
C*  'FORTRAN ROUTINES FOR USE WITH THE METHOD OF L-MOMENTS, VERSION 3' *
C*                                                                     *
C*  J. R. M. HOSKING                                                   *
C*  IBM RESEARCH DIVISION                                              *
C*  T. J. WATSON RESEARCH CENTER                                       *
C*  YORKTOWN HEIGHTS                                                   *
C*  NEW YORK 10598, U.S.A.                                             *
C*                                                                     *
C*  VERSION 3     AUGUST 1996                                          *
C*                                                                     *
C*  VERSION 3.02  MARCH 1997                                           *
C*  * Minor change to FORMAT statement 6080                            *
C*                                                                     *
C***********************************************************************
C
C  EXAMPLE PROGRAM FOR CLUSTER ANALYSIS.  THE PROGRAM READS IN
C  ATTRIBUTES FOR A NUMBER OF SITES, TRANSFORMS THE ATTRIBUTES, FORMS
C  CLUSTERS BY WARD'S METHOD, PRINTS INFORMATION ABOUT THE CLUSTERS,
C  AND REFINES THE CLUSTERS USING THE K-MEANS ALGORITHM.
C    THE ANALYSIS FOLLOWS HOSKING AND WALLIS ("REGIONAL FREQUENCY
C  ANALYSIS: AN APPROACH BASED ON L-MOMENTS", CAMBRIDGE UNIV. PRESS,
C  1997, SEC. 9.2).
C
C  PARAMETERS OF PROGRAM:
C  MAXNS  - MAXIMUM NUMBER OF SITES
C  NATMAX - MAXIMUM NUMBER OF ATTRIBUTES
C  INFILE - STREAM NUMBER TO WHICH INPUT FILE IS ATTACHED
C  KOUT   - STREAM NUMBER TO WHICH OUTPUT FILE IS ATTACHED
C
C  ROUTINES USED: CLUAGG,CLUINF,CLUKM,KMNS,OPTRA,QTRAN
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MAXNS=104,NATMAX=10,INFILE=7,KOUT=6)
C
      PARAMETER (NWORK=MAXNS*(MAXNS-1)/2,NRWORK=NWORK)
      CHARACTER*12 ID(MAXNS)
      INTEGER MERGE(2,MAXNS),IWORK(MAXNS),NREC(MAXNS),
     *  IASSGN(MAXNS),LIST(MAXNS),NUM(MAXNS)
      DOUBLE PRECISION X(MAXNS,NATMAX),WGSS(MAXNS),WORK(NWORK),
     *  WEIGHT(NATMAX),XMOM(5,MAXNS),Z(NATMAX),CENT(NATMAX,MAXNS)
      REAL RWORK(NRWORK)
      DATA ZERO/0D0/,ONE/1D0/,THREE/3D0/
C
      N=104
      READ(INFILE,*)
C
C         FOR EACH DATA POINT ...
C
      DO 10 I=1,N
C
C         ... READ THE DATA ...
C
      READ(INFILE,5000)ID(I),XLAT,XLONG,AREA,ELEV
C
C         ... AND COMPUTE THE TRANSFORMED ATTRIBUTES
C
      X(I,1)=DLOG(AREA)
      X(I,2)=DSQRT(ELEV)
      X(I,3)=XLAT
      X(I,4)=XLONG
C
C         END OF LOOP OVER DATA POINTS
C
   10 CONTINUE
C
C         SET WEIGHTS FOR EACH ATTRIBUTE
C
      NATT=4
      WEIGHT(1)=THREE
      WEIGHT(2)=ONE
      WEIGHT(3)=ONE
      WEIGHT(4)=ONE
C
C         FOR EACH ATTRIBUTE ...
C
      DO 40 J=1,NATT
C
C         ... CALCULATE ITS STANDARD DEVIATION ACROSS THE DATA POINTS ...
C
      SUM1=ZERO
      SUM2=ZERO
      DO 20 I=1,N
      SUM1=SUM1+X(I,J)
      SUM2=SUM2+X(I,J)**2
   20 CONTINUE
      SD=DSQRT((SUM2-SUM1*SUM1/N)/(N-ONE))
C
C         ... DIVIDE THE WEIGHT BY THIS STANDARD DEVIATION ...
C
      WEIGHT(J)=WEIGHT(J)/SD
C
C         ... AND APPLY THE WEIGHT TO EACH DATA POINT
C
      DO 30 I=1,N
   30 X(I,J)=X(I,J)*WEIGHT(J)
C
C         END OF LOOP OVER ATTRIBUTES
C
   40 CONTINUE
C
C         WARD'S ALGORITHM
C
      NX=MAXNS
      NW=NWORK
      CALL CLUAGG(3,X,NX,N,NATT,MERGE,WGSS,IWORK,WORK,NW)
      WRITE(KOUT,6000)
      DO 50 I=1,N-1
      WRITE(KOUT,6010)I,N-I,MERGE(1,I),MERGE(2,I),WGSS(I)
   50 CONTINUE
C
C         PRINT INFORMATION ABOUT THE 7-CLUSTER GROUPING
C
      NCLUST=7
      CALL CLUINF(NCLUST,N,MERGE,IASSGN,LIST,NUM)
      WRITE(KOUT,6020)
      WRITE(KOUT,6030)(IASSGN(I),I=1,N)
      WRITE(KOUT,6040)
      IORIG=0
      DO 60 ICL=1,NCLUST
      NN=NUM(ICL)
      WRITE(KOUT,6050)ICL,NN
      WRITE(KOUT,6060)(IABS(LIST(I)),I=IORIG+1,IORIG+NN)
      IORIG=IORIG+NN
   60 CONTINUE
C
C         ADJUST CLUSTERS BY K-MEANS ALGORITHM
C
      MAXIT=10
      CALL CLUKM(X,NX,N,NATT,NCLUST,IASSGN,LIST,NUM,SS,MAXIT,IWORK,
     *  RWORK,NRWORK)
C
C         PRINT INFORMATION ABOUT ADJUSTED CLUSTERS
C
      WRITE(KOUT,6070)SS
      WRITE(KOUT,6020)
      WRITE(KOUT,6030)(IASSGN(I),I=1,N)
      WRITE(KOUT,6040)
      IORIG=0
      DO 70 ICL=1,NCLUST
      NN=NUM(ICL)
      WRITE(KOUT,6050)ICL,NN
      WRITE(KOUT,6060)(IABS(LIST(I)),I=IORIG+1,IORIG+NN)
      IORIG=IORIG+NN
   70 CONTINUE
C
C         FIND CLUSTER CENTERS, IN SPACE OF TRANSFORMED ATTRIBUTES
C
      WRITE(KOUT,6080)
      DO 80 ICL=1,NCLUST
      DO 80 IATT=1,NATT
   80 CENT(IATT,ICL)=ZERO
      ICL=1
      DO 100 I=1,N
      L=LIST(I)
      ISITE=IABS(L)
      DO 90 IATT=1,NATT
      CENT(IATT,ICL)=CENT(IATT,ICL)+X(ISITE,IATT)
   90 CONTINUE
      IF(L.LT.0)ICL=ICL+1
  100 CONTINUE
      DO 110 ICL=1,NCLUST
      NN=NUM(ICL)
      DO 110 IATT=1,NATT
  110 CENT(IATT,ICL)=CENT(IATT,ICL)/NN
C
C         TRANSFORM BACK TO ORIGINAL ATTRIBUTES
C
      DO 120 ICL=1,NCLUST
      Z(1)=DEXP(CENT(1,ICL)/WEIGHT(1))
      Z(2)=(CENT(2,ICL)/WEIGHT(2))**2
      Z(3)=CENT(3,ICL)/WEIGHT(3)
      Z(4)=CENT(4,ICL)/WEIGHT(4)
      WRITE(KOUT,6090)ICL,(Z(J),J=1,NATT)
  120 CONTINUE
C
      STOP
C
 5000 FORMAT(A8,4F8.0)
 6000 FORMAT(' MERGING SEQUENCE FROM WARD''S ALGORITHM'//
     * 8X,'NUMBER OF    MERGED       SUM OF'/
     * ' STAGE  CLUSTERS    CLUSTERS      SQUARES'/)
 6010 FORMAT(1X,I3,I9,I10,I5,F12.2)
 6020 FORMAT(/' ASSIGNMENT OF SITES TO CLUSTERS')
 6030 FORMAT(1X,10I4)
 6040 FORMAT(//' CLUSTER MEMBERSHIP')
 6050 FORMAT(/' CLUSTER',I4,'  HAS',I4,' MEMBERS:')
 6060 FORMAT(1X,10I4)
 6070 FORMAT(///' ADJUSTED CLUSTERS FROM K-MEANS ALGORITHM'/
     *  ' (SUM OF SQUARES =',F12.2,')')
 6080 FORMAT(/' CLUSTER CENTERS'/
     *  '          AREA      ELEV       LAT      LONG')
 6090 FORMAT(1X,I3,6F10.2)
C
      END
