      SUBROUTINE DOT(MM, M, N, A, CLAB, RLAB, TITLE, KL, DMIWRK, IWORK,
     *               IWORK1, DMWORK, WORK, WORK1, CWORK1, CWORK2, IERR,
     *               OUNIT)
C
C<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
C
C   PURPOSE
C   -------
C
C      CREATES A TREE OF CLUSTERS OF CASES FOR CATEGORICAL DATA BY
C      MINIMUM-MUTATION FITS (CHANGES OF VARIABLE VALUES BETWEEN A
C      CLUSTER AND ITS DIRECT ANCESTOR IN THE TREE ARE MINIMIZED.
C
C   DESCRIPTION
C   -----------
C
C   1.  THE MINIMUM MUTATION FIT ALGORITHM SETS UP A BINARY TREE OF
C       CLUSTERS IN POINTER FORM.  EACH CASE IS A SINGLETON CLUSTER AT
C       THE LEAF LEVEL OF THE TREE AND THE OTHER NODES LINK UP THE
C       CASES INTO CLUSTERS.  THE TREE IS BINARY SINCE EACH CLUSTER OF
C       CARDINALITY GREATER THAN 1 IS DIVIDED UP INTO TWO SUBCLUSTERS.
C       THE ANCESTOR OF A CLUSTER IN THE TREE IS THE SMALLEST CLUSTER
C       THAT CONTAINS IT.  VALUES ARE ASSIGNED TO THE ANCESTORS SO THAT
C       A MINIMUM NUMBER OF MUTATIONS (OR CHANGES OF A VARIABLE VALUE
C       BETWEEN A CLUSTER AND ITS ANCESTOR) OCCUR.
C
C   2.  THE DATA SHOULD BE REAL NUMBERS WHICH CORRESPOND TO NO MORE
C       THAN 36 CATAGORIES.  THE INITIAL TREE IS FORMED BY THE JOINING
C       ALGORITHM (SEE SUBROUTINE JOIN).  THE DATA MATRIX IS DESTROYED
C       DURING EXECUTION.
C
C   3.  THE OUTPUT IS WRITTEN ON FORTRAN UNIT OUNIT AND STARTS WITH THE
C       MATRIX OF CATEGORY VALUES.
C
C   4.  NEXT IS THE STANDARD TREE OUTPUT OF THE CLUSTERS FROM JOIN.  IT
C       LISTS THE CASES VERTICALLY AND HAS HORIZONTAL LINES EMANATING
C       FROM EACH CASE.  EACH CLUSTER WILL CORRESPOND TO A VERTICAL
C       LINE BETWEEN TWO HORIZONTAL LINES.  THE CASES BETWEEN AND
C       INCLUDED IN THE HORIZONTAL LINES ARE THE MEMBERS OF THE
C       CLUSTER.  THE DISTANCE FROM THE CASE NAMES TO THE VERTICAL
C       LINES CORRESPOND TO THE DISTANCE AT WHICH THE TWO CLUSTERS WERE
C       JOINED.
C
C   5.  NEXT IS THE DOT DIAGRAM OF THE CLUSTERS.  IT LISTS THE CASES
C       VERTICALLY, AND EACH COLUMN CORRESPONDS TO THE CLUSTER WHOSE
C       NUMBER IS AT THE BOTTOM OF THE DIAGRAM.  IF THE CASE IS A
C       MEMBER OF CLUSTER I, THEN THE ENTRY IN THE I-TH COLUMN OF THAT
C       CASE'S ROW IS THE LETTER O.  OTHERWISE, THE ENTRY IS A DOT.
C
C   6.  THE LAST DIAGRAM IS AN ARRAY WHOSE ROWS ARE THE VARIABLES AND
C       THE COLUMNS ARE THE CLUSTERS NUMBERED BY THE SAME NUMBERS AS IN
C       THE DIAGRAM IN 4.  THE ENTRY FOR CLUSTER I FOR VARIABLE J IS A
C       DOT IF THE VALUE FROM THE MINIMUM MUTATION FIT FOR CLUSTER I IS
C       THE SAME AS THE VALUE FROM THE MINIMUM MUTATION FIT FOR THE
C       DIRECT ANCESTOR OF CLUSTER I.  OTHERWISE, THE ENTRY IS THE
C       MINIMUM MUTATION VALUE FOR CLUSTER I.
C
C   7.  TO RECOVER THE CATEGORY VALUES FOR A CASE, START WITH THE
C       MINIMUM MUTATION VALUES FOR THE LAST CLUSTER AND GO THROUGH THE
C       CLUSTERS IN DECREASING ORDER CHANGING THE VALUES FOR A VARIABLE
C       IF AND ONLY IF THE CASE IS IN THE CLUSTER (IF THE LETTER O IS
C       IN THE COLUMN FOR THE CLUSTER) AND THE MINIMUM MUTATION VALUE
C       FOR THAT VARIABLE IN THAT CLUSTER IS NOT A DOT.  THE VALUE IS
C       THEN CHANGED TO THE MINIMUM MUTATION VALUE.  AFTER ALL CLUSTERS
C       HAVE BEEN PROCESSED, THE VALUES FOR THE VARIABLES WILL BE THE
C       ORIGINAL CATEGORY VALUES FOR THE CASE.
C
C   INPUT PARAMETERS
C   ----------------
C
C   MM    INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         THE FIRST DIMENSION OF THE MATRIX A.  MUST BE AT LEAST M.
C
C   M     INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         THE NUMBER OF CASES.
C
C   N     INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         THE NUMBER OF VARIABLES.
C
C   A     REAL MATRIX WHOSE FIRST DIMENSION MUST BE MM AND WHOSE SECOND
C            DIMENSION MUST BE AT LEAST N (DESTROYED DURING EXECUTION).
C         THE MATRIX OF CATEGORY VALUES.
C
C         A(I,J) IS THE CATEGORY FOR THE J-TH VARIABLE FOR THE I-TH
C                CASE.
C
C   CLAB  VECTOR OF 4-CHARACTER VARIABLES DIMENSIONED AT LEAST M.
C            (UNCHANGED ON OUTPUT).
C         THE LABELS OF THE VARIABLES.
C
C   RLAB  VECTOR OF 4-CHARACTER VARIABLES DIMENSIONED AT LEAST N.
C            (UNCHANGED ON OUTPUT).
C         THE LABELS OF THE CASES.
C
C   TITLE 10-CHARACTER VARIABLE (UNCHANGED ON OUTPUT).
C         TITLE OF THE DATA SET.
C
C   KL    INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         MAXIMUM CATEGORY VALUE.  MUST BE LESS THAN 36.
C
C         THE CATEGORY VALUES MUST BE BETWEEN 1 AND KL.  IF NOT, AN
C         ERROR MESSAGE WILL BE PRINTED AND EXECUTION WILL BE
C         TERMINATED.
C
C   DMIWRK INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         THE FIRST DIMENSION OF THE MATRIX IWORK. MUST BE AT LEAST 2*M.
C
C   IWORK INTEGER MATRIX WHOSE FIRST DIMENSION MUST BE DMIWRK AND SECOND
C            DIMENSION MUST BE AT LEAST M.
C         WORK MATRIX.
C
C   IWORK1 INTEGER VECTOR DIMENSIONED AT LEAST 2*M.
C         WORK VECTOR.
C
C   DMWORK INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         THE FIRST DIMENSION OF THE MATRIX WORK.  MUST BE AT LEAST M.
C
C   WORK  REAL MATRIX WHOSE FIRST DIMENSION MUST BE DMWORK AND WHOSE
C            SECOND DIMENSION MUST BE AT LEAST N.
C         WORK MATRIX.
C
C   WORK1 REAL VECTOR DIMENSIONED AT LEAST 2*M.
C         WORK VECTOR.
C
C   CWORK1 VECTOR OF 4-CHARACTER VARIABLES DIMENSIONED AT LEAST 2*M+1.
C         WORK VECTOR.
C
C   CWORK2 VECTOR OF 1-CHARACTER VARIABLES DIMENSIONED AT LEAST 2*M+10.
C         WORK VECTOR.
C
C   OUNIT INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         UNIT NUMBER FOR OUTPUT.
C
C   OUTPUT PARAMETER
C   ----------------
C
C   IERR  INTEGER SCALAR.
C         ERROR FLAG.
C
C         IERR = 0, NO ERRORS WERE DETECTED DURING EXECUTION
C
C         IERR = 1, EITHER THE FIRST AND LAST CASES OR THE CLUSTER
C                   DIAMETER FOR A CLUSTER IS OUT OF BOUNDS.  THE
C                   CLUSTER AND ITS VALUES ARE PRINTED ON UNIT OUNIT.
C                   EXECUTION WILL CONTINUE WITH QUESTIONABLE RESULTS
C                   FOR THAT CLUSTER.  ERROR FLAG SET IN JOINING
C                   ROUTINE.
C
C         IERR = 2, A CLUSTER BOUNDARY IS ILLEGAL.  THE CLUSTER AND THE
C                   BOUNDARY ARE PRINTED ON UNIT OUNIT.  EXECUTION WILL
C                   CONTINUE WITH QUESTIONABLE RESULTS FOR THAT CLUSTER.
C
C         IERR = 3, TWO CLUSTERS OVERLAP.  THE NUMBERS OF THE TWO
C                   CLUSTERS ARE PRINTED ON UNIT OUNIT.  EXECUTION WILL
C                   CONTINUE WITH QUESTIONABLE RESULTS FOR THAT CLUSTER.
C
C         IERR = 4, A VALUE OF THE VARIABLE IS GREATER THAN THE MAXIMUM
C                   VALUE.  THE VALUE AND THE MAXIMUM VALUE IS PRINTED
C                   ON UNIT OUNIT AND EXECUTION IS TERMINATED.  INCREASE
C                   THE VALUE OF KL OR CHANGE THE DATA VALUE.
C
C   REFERENCES
C   ----------
C
C     HARTIGAN, J. A. (1975).  CLUSTERING ALGORITHMS, JOHN WILEY &
C        SONS, INC., NEW YORK.  PAGES 233-250.
C
C     HARTIGAN, J. A. (1972).  "MINIMUM MUTATION FITS TO A GIVEN TREE"
C        BIOMETRICS.  VOL. 29, PAGES 53-65.
C
C<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
C
      INTEGER DMIWRK, DMWORK, OUNIT
      DIMENSION A(MM,*), IWORK1(*), WORK1(*), IWORK(DMIWRK,*),
     *          WORK(DMWORK,*)
      CHARACTER*4 CLAB(*), RLAB(*), CWORK1(*)
      CHARACTER*1 AA(36), CWORK2(*)
      CHARACTER*10 TITLE
      DATA AA/'.','1','2','3','4','5','6','7','8','9','A','B',
     *        'C','D','E','F','G','H','I','J','K','L','M','N',
     *        'O','P','Q','R','S','T','U','V','W','X','Y','Z'/
C
      IERR = 0
      KT = 2 * M - 1
C
C     ELIMINATE DECIMAL PARTS OF THE CATEGORY VALUES
C
      DO 10 I = 1 , M
         DO 10 J = 1 , N
  10        A(I,J) = REAL(INT(A(I,J)))
      IF (OUNIT .GT. 0) WRITE(OUNIT,1)
  1   FORMAT('1')
      CALL OUT(MM,M,N,A,CLAB,RLAB,TITLE,OUNIT)
C
C     COPY DATA
C
      DO 30 I = 1 , M
         DO 20 J = 1 , N
  20        WORK(I,J) = A(I,J)
  30     CWORK1(I) = RLAB(I)
C
C     GET TREE USING JOIN
C
      CALL JOIN(MM,M,N,A,RLAB,TITLE,DMIWRK,IWORK,WORK1,CWORK1(M+1),IERR,
     *          OUNIT)
C
C     PERMUTE TABLE TO AGREE WITH THE ORDERING PRODUCED BY JOIN
C     RESTORE DATA DESTROYED BY JOIN
C
      DO 60 I = 1 , M
         DO 50 K = 1 , M
            IF (CWORK1(K) .NE. CWORK1(M+I+1)) GOTO 50
            DO 40 J = 1 , N
  40           A(I,J) = WORK(K,J)
            RLAB(I) = CWORK1(K)
            GOTO 60
  50     CONTINUE
  60  CONTINUE
      CALL NBTOMT(M+1,M,DMIWRK,IWORK,KT,IWORK1,IERR,OUNIT)
      IF (OUNIT .GT. 0) WRITE(OUNIT,1)
C
C     CHECK RANGE OF A
C
      DO 70 I=1,M
         DO 70 J=1,N
            IF(A(I,J).LT.1..OR.A(I,J).GT.KL) THEN
            IF (OUNIT .GT. 0) WRITE(OUNIT,2) I,J,KL
    2       FORMAT(I5,' ,',I5,' TH DATA VALUE IS NOT IN THE RANGE (1,',
     *             I2,')')
            IERR = 4
            RETURN
            ENDIF
  70  CONTINUE
C
C     COUNT NUMBER IN EACH CLUSTER
C
      DO 80 I=1,KT
         IWORK(I,1)=0
   80    IF(I.LE.M) IWORK(I,1)=1
      DO 90 I=1,KT-1
         J=IWORK1(I)
         IWORK(J,1)=IWORK(J,1)+IWORK(I,1)
   90 CONTINUE
      IF (OUNIT .GT. 0) WRITE(OUNIT,3)
    3 FORMAT('0TREE')
C
C     PRINT TREE, A LINE AT A TIME.
C
      DO 110 I=1,M
         J=KT
         IWORK(J,1)=IWORK(J,1)-1
         CWORK2(KT)=AA(25)
         DO 100 L=2,KT
            LL=KT-L+1
            IF(IWORK(LL,1).GT.0) IWORK(LL,1)=-IWORK(LL,1)
            CWORK2(LL)=AA(1)
            IF(IWORK1(LL).EQ.J.AND.IWORK(LL,1).NE.0) THEN
               J=LL
               IWORK(LL,1)=-IWORK(LL,1)-1
               CWORK2(J)=AA(25)
            ENDIF
  100    CONTINUE
         IF (OUNIT .GT. 0) WRITE(OUNIT,7) J,RLAB(J),(CWORK2(L),L=1,KT)
  110 CONTINUE
C
C     GO THROUGH VARIABLES, ONE AT A TIME
C
      IKT = KT / 10
      INTKT = INT(KT)
      IF (OUNIT .GT. 0) THEN
         WRITE(OUNIT,4)((I,J=1,10),I=1,IKT)
    4    FORMAT(22X,4I1,1X,17(5I1,1X))
         WRITE(OUNIT,5)((I,I=0,9),J=1,IKT+1)
    5    FORMAT(10X,4I1,1X,19(5I1,1X))
         DO 120 III = 1 , INTKT
  120       CWORK2(III) = ' '
         IKT = IKT * 10 + 10 + (IKT+1) * 2
         DO 130 III = 1 , IKT
  130       CWORK2(III) = '-'
         WRITE(OUNIT,6) (CWORK2(III),III=1,IKT)
    6    FORMAT(10X,120(A1))
      ENDIF
      DO 160 JJ=1,N
C
C     FIND THE MINIMUM MUTATION FIT FOR EACH VARIABLE
C
         DO 140 I=1,M
  140       WORK1(I)=A(I,JJ)
         CALL MMFIT(M,WORK1,KT,IWORK1,KL,DMIWRK,IWORK,IERR,OUNIT)
         IF (IERR .EQ. 4) RETURN
         DO 150 I=1,KT
C
C     REPLACE WITH A ZERO IF THE MINIMUM MUTATION VALUES CORRESPOND FOR
C     A CLUSTER AND ITS DIRECT ANCESTOR
C
            IF(I.NE.KT) THEN
               II=IWORK1(I)
               IF(WORK1(I).EQ.WORK1(II)) WORK1(I)=0.
            ENDIF
            J=INT(WORK1(I)+1)
            IF(J.LT.1) J=1
            IF(J.GT.35) J=36
            CWORK2(I)=AA(J)
  150    CONTINUE
         IF (OUNIT .GT. 0) WRITE(OUNIT,7) JJ,CLAB(JJ),(CWORK2(I),I=1,KT)
    7    FORMAT(I4,A4,3X,20(3A1,1X,2A1))
  160 CONTINUE
      RETURN
      END
      SUBROUTINE OUT(MM, M, N, A, CLAB, RLAB, TITLE, OUNIT)
C
C<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
C
C   PURPOSE
C   -------
C
C      PRINTS A DATA MATRIX WITH THE APPROPRIATE LABELS FOR THE CASES
C      AND VARIABLES WHICH IS OUTPUT FROM MANY OF THE SUBROUTINES FROM
C      THE CLUSTER SUBLIBRARY
C
C   DESCRIPTION
C   -----------
C
C   1.  THE NAME OF THE ARRAY IS PRINTED OUT FOLLOWED BY THE MATRIX OF
C       DATA VALUES WITH APPROPRIATE ROW AND COLUMN LABELS ON FORTRAN
C       UNIT OUNIT.
C
C   INPUT PARAMETERS
C   ----------------
C
C   MM    INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         THE FIRST DIMENSION OF THE MATRIX A.  MUST BE AT LEAST M.
C
C   M     INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         THE NUMBER OF CASES.
C
C   N     INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         THE NUMBER OF VARIABLES.
C
C   A     REAL MATRIX WHOSE FIRST DIMENSION MUST BE MM AND SECOND
C            DIMENSION MUST BE AT LEAST N  (UNCHANGED ON OUTPUT).
C         THE MATRIX OF DATA VALUES.
C
C         A(I,J) IS THE VALUE FOR THE J-TH VARIABLE FOR THE I-TH CASE.
C
C   CLAB  VECTOR OF 4-CHARACTER VARIABLES DIMENSIONED AT LEAST N.
C            (UNCHANGED ON OUTPUT).
C         THE LABELS OF THE VARIABLES.
C
C   RLAB  VECTOR OF 4-CHARACTER VARIABLES DIMENSIONED AT LEAST M.
C            (UNCHANGED ON OUTPUT).
C         THE LABELS OF THE CASES.
C
C   TITLE 10-CHARACTER VARIABLE (UNCHANGED ON OUTPUT).
C         TITLE OF THE DATA SET.
C
C   OUNIT INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         UNIT NUMBER FOR OUTPUT.
C
C   REFERENCE
C   ---------
C
C     HARTIGAN, J. A. (1975).  CLUSTERING ALGORITHMS, JOHN WILEY &
C        SONS, INC., NEW YORK.  PAGES 40-41.
C
C<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
C
      DIMENSION A(MM,*)
      INTEGER OUNIT
      CHARACTER*4 CLAB(*), RLAB(*)
      CHARACTER*10 TITLE
C
      IF (OUNIT .LE. 0) RETURN
      WRITE(OUNIT,1) TITLE
      WRITE(OUNIT,2)(CLAB(I),I=1,N)
      DO 10 I=1,M
   10    WRITE(OUNIT,3) RLAB(I), (A(I,J),J=1,N)
    1 FORMAT(' ARRAY ', A10, /)
    2 FORMAT(7X, 11(4X, A4, 3X))
    3 FORMAT(1X, A4, 11(F10.4,1X))
      RETURN
      END
      SUBROUTINE NBTOMT(M, K, DMNB, NB, KT, MT, IERR, OUNIT)
C
C<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
C
C   PURPOSE
C   -------
C
C      CONVERTS A REGULAR TREE DESCRIPTION OF OBJECTS TO A LINKED LIST
C      DESCRIPTION
C
C   DESCRIPTION
C   -----------
C
C   1.  THE LINKED LIST IS AN ARRAY OF SIZE BETWEEN M AND 2*M.  EACH
C       ENTRY IN THE ARRAY CORRESPONDS TO A NODE IN THE TREE AND THE
C       VALUE FOR NODE I IS THE NODE THAT IS THE DIRECT ANCESTOR OF
C       NODE I.
C
C   INPUT PARAMETERS
C   ----------------
C
C   M     INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         THE NUMBER OF OBJECTS.
C
C   K     INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         THE NUMBER OF CLUSTERS.
C
C   DMNB  INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         THE FIRST DIMENSION OF THE MATRIX NB.  MUST BE AT LEAST 3.
C
C   NB    INTEGER MATRIX WHOSE FIRST DIMENSION MUST BE DMNB AND SECOND
C            DIMENSION MUST BE AT LEAST K (UNCHANGED ON OUTPUT).
C         STORES REGULAR DESCRIPTION OF TREE.
C
C          NB(1,K) = FIRST OBJECT IN CLUSTER K
C          NB(2,K) = LAST OBJECT IN CLUSTER K
C          NB(3,K) = CLUSTER DIAMETER
C
C   KT    INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         THE NUMBER OF NODES IN TREE.
C
C   OUNIT INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         UNIT NUMBER FOR ERROR MESSAGES.
C
C   OUTPUT PARAMETERS
C   -----------------
C
C   MT    INTEGER VECTOR DIMENSIONED AT LEAST KT.
C         THE VECTOR OF POINTERS FOR THE TREE.
C
C         MT(I) IS THE ANCESTOR NODE FOR NODE I.
C
C   IERR  INTEGER SCALAR.
C         ERROR FLAG.
C
C         IERR = 0, NO ERRORS WERE DETECTED DURING EXECUTION
C
C         IERR = 2, A CLUSTER BOUNDARY IS ILLEGAL.  THE CLUSTER AND THE
C                   BOUNDARY ARE PRINTED ON UNIT OUNIT.  EXECUTION WILL
C                   CONTINUE WITH QUESTIONABLE RESULTS FOR THAT CLUSTER.
C
C         IERR = 3, TWO CLUSTERS OVERLAP.  THE NUMBERS OF THE TWO
C                   CLUSTERS ARE PRINTED ON UNIT 6.  EXECUTION WILL
C                   CONTINUE WITH QUESTIONABLE RESULTS FOR THAT CLUSTER.
C
C   REFERENCE
C   ---------
C
C     HARTIGAN, J. A. (1975).  CLUSTERING ALGORITHMS, JOHN WILEY &
C        SONS, INC., NEW YORK.  PAGE 166.
C
C<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
C
      INTEGER DMNB, OUNIT
      DIMENSION NB(DMNB,*), MT(*)
C
C     CHECK NB ARRAY
C
      DO 10 I=1,K
         NB(1,I) = NB(1,I) - 1
         NB(2,I) = NB(2,I) - 1
         IF(NB(1,I).LT.1) THEN
            IF (OUNIT .GT. 0) WRITE(OUNIT,1) I
            IERR = 2
         ENDIF
         IF(NB(2,I).LT.NB(1,I)) THEN
            IF (OUNIT .GT. 0) WRITE(OUNIT,2) I
            IERR = 2
         ENDIF
         IF(NB(2,I).GT.M) THEN
            IF (OUNIT .GT. 0) WRITE(OUNIT,3) I, M
            IERR = 2
         ENDIF
  10  CONTINUE
    1 FORMAT(' CLUSTER',I5,' HAS BOUNDARY LESS THAN 1')
    2 FORMAT(' CLUSTER',I5,' HAS FIRST BOUNDARY EXCEEDING SECOND BOUNDAR
     *Y')
    3 FORMAT(' CLUSTER',I5,' HAS BOUNDARY GREATER THAN',I5)
C
C     CHECK NB ARRAY FOR OVERLAPS
C
      DO 20 I=1,K
         DO 20 J=1,I
            I1=NB(1,I)
            I2=NB(2,I)
            J1=NB(1,J)
            J2=NB(2,J)
            IF((I1-J1)*(I1-J2)*(I2-J1)*(I2-J2).LT.0) THEN
               IF (OUNIT .GT. 0) WRITE(OUNIT,4) I, J
               IERR = 3
            ENDIF
   20 CONTINUE
    4 FORMAT(' CLUSTERS',I5,' AND',I5,' OVERLAP')
C
C     CONSTRUCT MT
C
      DO 30 I=1,KT
   30    MT(I)=0
      DO 70 I=M,KT
C
C     FIND CLUSTER UNASSIGNED, WITH MINIMUM NUMBER OF ELEMENTS
C
         IM=I
         MIN=M
         DO 40 J=1,K
            IF(NB(1,J).GE.0) THEN
               MM=NB(2,J)-NB(1,J)+1
               IF(MM.LE.MIN) THEN
                  IM=J
                  MIN=MM
               ENDIF
            ENDIF
   40    CONTINUE
         IF(MIN.NE.M) THEN
C
C     FIND SMALLEST CLUSTER INCLUDED IN I, NOT YET ASSIGNED
C
            JL=NB(1,IM)
            JU=NB(2,IM)
            NB(1,IM)=-JL
            DO 60 J=JL,JU
               L=J
   50          IF(MT(L).NE.0) THEN
                  L=MT(L)
                  IF(MT(L).EQ.L) GO TO 60
                  GO TO 50
               ENDIF
   60          MT(L)=I
         ENDIF
   70 CONTINUE
      MT(KT)=KT
      DO 80 I=1,K
   80    IF(NB(1,I).LT.0) NB(1,I)=-NB(1,I)
      RETURN
      END
      SUBROUTINE JOIN(MM, M, N, A, RLAB, TITLE, DMIWRK, IWORK, WORK,
     *                TLAB, IERR, OUNIT)
C
C<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
C
C   PURPOSE
C   -------
C
C      USES A GENERAL JOINING ALGORITHM TO FORM AND OUTPUT A TREE OF
C      CLUSTERS OF CASES
C
C   DESCRIPTION
C   -----------
C
C   1.  THE VARIABLES MUST BE SCALED SIMILARLY (CLUSTER SUBROUTINE
C       STAND CAN BE USED TO STANDARDIZE THE VARIABLES.)
C
C   2.  INITIALLY, EACH OBJECT IS ASSUMED TO BE A CLUSTER.  THEN THE
C       TWO CLOSEST CLUSTERS ARE AMALGAMATED TO FORM A NEW CLUSTER.
C       THIS PROCESS OF COMBINING THE CLOSEST TWO CLUSTERS CONTINUES
C       UNTIL ONLY ONE CLUSTER REMAINS.  A TREE IS MAINTAINED TO SHOW
C       ALL CLUSTERS AT EACH STEP.
C
C   3.  THE OUTPUT IS THE STANDARD TREE WRITTEN ON FORTRAN UNIT OUNIT.
C       IT LISTS THE CASES VERTICALLY AND HAS HORIZONTAL LINES
C       EMANATING FROM EACH CASE.  EACH CLUSTER WILL CORRESPOND TO A
C       VERTICAL LINE BETWEEN TWO HORIZONTAL LINES.  THE CASES BETWEEN
C       AND INCLUDED IN THE HORIZONTAL LINES ARE THE MEMBERS OF THE
C       CLUSTER.  THE DISTANCE FROM THE CASE NAMES TO THE VERTICAL
C       LINES CORRESPOND TO THE CLUSTER DIAMETER OR THE DISTANCE AT
C       WHICH THE TWO CLUSTERS WERE JOINED.
C
C   INPUT PARAMETERS
C   ----------------
C
C   MM    INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         THE FIRST DIMENSION OF THE MATRIX A.  MUST BE AT LEAST M.
C
C   M     INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         THE NUMBER OF CASES.
C
C   N     INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         THE NUMBER OF VARIABLES.
C
C   A     REAL MATRIX WHOSE FIRST DIMENSION MUST BE MM AND WHOSE SECOND
C            DIMENSION MUST BE AT LEAST N (DESTROYED DURING EXECUTION).
C         THE MATRIX OF DISTANCES BETWEEN CASES.
C
C         A(I,J) IS THE VALUE FOR THE J-TH VARIABLE FOR THE I-TH CASE.
C
C   RLAB  VECTOR OF 4-CHARACTER VARIABLES DIMENSIONED AT LEAST N.
C            (UNCHANGED ON OUTPUT).
C         THE LABELS OF THE CASES.
C
C   TITLE 10-CHARACTER VARIABLE (UNCHANGED ON OUTPUT).
C         TITLE OF THE DATA SET.
C
C   DMIWRK INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         THE FIRST DIMENSION OF THE MATRIX IWORK.  MUST BE AT LEAST 4.
C
C   IWORK INTEGER MATRIX WHOSE FIRST DIMENSION MUST BE DMIWRK AND SECOND
C            DIMENSION MUST BE AT LEAST M.
C         WORK MATRIX.
C
C   WORK  REAL VECTOR DIMENSIONED AT LEAST M.
C         WORK VECTOR.
C
C   TLAB  VECTOR OF 4-CHARACTER VARIABLES DIMENSIONED AT LEAST M+1.
C         WORK VECTOR.
C
C   OUNIT INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         UNIT NUMBER FOR OUTPUT.
C
C   OUTPUT PARAMETER
C   ----------------
C
C   IERR  INTEGER SCALAR.
C         ERROR FLAG.
C
C         IERR = 0, NO ERRORS WERE DETECTED DURING EXECUTION
C
C         IERR = 1, EITHER THE FIRST AND LAST CASES OR THE CLUSTER
C                   DIAMETER FOR A CLUSTER IS OUT OF BOUNDS.  THE
C                   CLUSTER AND ITS VALUES ARE PRINTED ON UNIT OUNIT.
C                   EXECUTION WILL CONTINUE WITH QUESTIONABLE RESULTS
C                   FOR THAT CLUSTER.
C
C   REFERENCES
C   ----------
C
C     HARTIGAN, J. A. (1975).  CLUSTERING ALGORITHMS, JOHN WILEY &
C        SONS, INC., NEW YORK.  PAGES 216-232.
C
C     HARTIGAN, J. A. (1975) PRINTER GRAPHICS FOR CLUSTERING. JOURNAL OF
C        STATISTICAL COMPUTATION AND SIMULATION. VOLUME 4,PAGES 187-213.
C
C<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
C
      INTEGER DMIWRK, OUNIT
      DIMENSION A(MM,*), IWORK(DMIWRK,*), WORK(*)
      CHARACTER*4 RLAB(*), TLAB(*)
      CHARACTER*10 TITLE
C
      IERR = 0
      DO 10 I=1,M
         TLAB(I) = RLAB(I)
   10    IWORK(3,I)=I
C
C     FIND CLOSEST PAIR
C
      DO 40 K=1,M-1
         DMIN=R1MACH(2)
         IMIN=2
         JMIN=3
         DO 30 I=1,M
            IF(IWORK(3,I).GE.0) THEN
               DO 20 J=I,M
                  IF(IWORK(3,J).GE.0.AND.J.NE.I) THEN
                     Z=DIST(N,A(I,1),A(J,1),MM,2.)
                     IF(Z.LT.DMIN) THEN
                        IMIN=I
                        JMIN=J
                        DMIN=Z
                     ENDIF
                  ENDIF
   20          CONTINUE
            ENDIF
   30    CONTINUE
C
C     FORM A NEW CLUSTER BY AMALGAMATING ROWS OF CLOSEST PAIR
C
         I=IMIN
         IWORK(1,K)=I
         II=IWORK(3,JMIN)
         WORK(K)=DMIN
         IWORK(2,K)=II
         IWORK(3,JMIN)=-IWORK(3,I)
         IWORK(3,I)=II
         CALL AMALG(N,A(IMIN,1),A(JMIN,1),MM)
   40 CONTINUE
C
C     REORDER OBJECTS, ASSIGN LABELS
C
      J=M+1
      KC=M-1
      IF(M.EQ.2) RETURN
      L=IWORK(2,KC)
      DO 50 JJ=1,M
         J=M-JJ+2
         LL=-IWORK(3,L)
         IWORK(3,L)=J
         L=LL
   50 CONTINUE
      DO 60 K=1,KC
         L=IWORK(1,K)
         IWORK(1,K)=IWORK(3,L)
         L=IWORK(2,K)
   60    IWORK(2,K)=IWORK(3,L)
      DO 70 I=1,M
         J=IWORK(3,I)
   70    TLAB(J)=RLAB(I)
C
C     DEFINE LAST TWO CLUSTERS, AND CLUSTER DIAMETERS
C
      DMAX=0.
      DO 80 I=1,KC
   80    IF(WORK(I).GT.DMAX) DMAX=WORK(I)
      IF (OUNIT .GT. 0) THEN
         WRITE(OUNIT,1)
    1    FORMAT('1')
         WRITE(OUNIT,2) DMAX
    2    FORMAT(' MAXIMUM JOINING DISTANCE =',F20.6)
      ENDIF
      DO 90 I=1,KC
   90    IWORK(3,I)=(WORK(I)/DMAX)*100
      DO 100 J=1,3
         DO 100 K=KC,M
  100       IWORK(J,K)=IWORK(J,KC)
C
C     PRODUCE OUTPUT
C
      TLAB(1) = TITLE
      CALL TREE1(M+1,M-1,DMIWRK,IWORK,TLAB,IERR,OUNIT)
      RETURN
      END
      SUBROUTINE MMFIT(M, VAR, K, MT, NVAL, DMIWRK, IWORK, IERR, OUNIT)
C
C<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
C
C   PURPOSE
C   -------
C
C      FINDS THE MINIMUM MUTATION FIT OF A VARIABLE TO A TREE
C
C   INPUT PARAMETERS
C   ----------------
C
C   M     INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         THE NUMBER OF OBJECTS.
C
C   VAR   REAL VECTOR DIMENSIONED AT LEAST K (CHANGED ON OUTPUT).
C         VAR(I) FOR I=1,...,M IS THE CATEGORY VALUE FOR CASE I FOR THE
C            VARIABLE.
C
C   K     INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         THE NUMBER OF NODES IN THE TREE.
C
C   MT    INTEGER VECTOR DIMENSIONED AT LEAST K (CHANGED ON OUTPUT).
C         MT(I) STORES THE NODE THAT IS THE DIRECT ANCESTOR OF NODE I IN
C            THE TREE.  (UPDATED DURING ROUTINE).
C
C   NVAL  INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         MAXIMUM VALUE OF VAR.
C
C   DMIWRK INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         THE FIRST DIMENSION OF THE MATRIX IWORK.  MUST BE AT LEAST K.
C
C   IWORK INTEGER MATRIX WHOSE FIRST DIMENSION MUST BE DMIWRK AND SECOND
C            DIMENSION MUST BE AT LEAST NVAL.
C         WORK MATRIX.
C
C   OUNIT INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         UNIT NUMBER FOR ERROR MESSAGES.
C
C   OUTPUT PARAMETER
C   ----------------
C
C   VAR   REAL VECTOR DIMENSIONED AT LEAST K.
C         VAR(I) FOR I=1,...,K IS THE FITTED VALUE FOR NODE I.
C
C   IERR  INTEGER SCALAR.
C         ERROR FLAG.
C
C         IERR = 0, NO ERRORS WERE DETECTED DURING EXECUTION
C
C         IERR = 4, A VALUE OF THE VARIABLE IS GREATER THAN THE MAXIMUM
C                   VALUE.  INCREASE THE VALUE OF NVAL.
C
C   REFERENCES
C   ----------
C
C     HARTIGAN, J. A. (1975).  CLUSTERING ALGORITHMS, JOHN WILEY &
C        SONS, INC., NEW YORK.  PAGE 247.
C
C     HARTIGAN, J. A. (1972).  "MINIMUM MUTATION FITS TO A GIVEN TREE"
C        BIOMETRICS.  VOL. 29, PAGES 53-65.
C
C<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
C
      INTEGER DMIWRK, OUNIT
      DIMENSION VAR(*), IWORK(DMIWRK,*), MT(*)
C
C     CHECK RANGE OF VAR
C
      DO 10 I=1,M
         IF(VAR(I).LT.1..OR.VAR(I).GT.NVAL) THEN
            IF (OUNIT .GT. 0) WRITE(OUNIT,1) I
            IERR=4
         ENDIF
   10 CONTINUE
    1 FORMAT(I5,' TH VALUE OF VARIABLE OUT OF RANGE')
      IF(IERR.NE.0) RETURN
C
C     OPTIMAL ASSIGNMENT TO NODES IGNORING REST OF TREE
C
      DO 20 I=M+1,K
   20    VAR(I)=0.
      DO 30 I=1,K
         DO 30 J=1,NVAL
   30       IWORK(I,J)=0
      DO 40 I=1,M
         J=INT(VAR(I))
   40    IWORK(I,J)=1
      DO 80 I=1,K
         IMAX=0
         DO 50 J=1,NVAL
   50       IF(IWORK(I,J).GT.IMAX) IMAX=IWORK(I,J)
         DO 60 J=1,NVAL
            IF(IWORK(I,J).LT.IMAX) IWORK(I,J)=0
   60       IF(IWORK(I,J).NE.0) IWORK(I,J)=1
         IF(I.NE.K) THEN
            II=MT(I)
            DO 70 J=1,NVAL
   70          IWORK(II,J)=IWORK(II,J)+IWORK(I,J)
         ENDIF
   80 CONTINUE
C
C     FIND OPTIMAL ASSIGNMENT OVER WHOLE TREE
C
      DO 90 I=1,K
         II=K-I+1
         IT=MT(II)
         DO 90 J=1,NVAL
   90       IF((VAR(II).EQ.0..OR.INT(VAR(IT)).EQ.J).AND.IWORK(II,J).EQ.1
     *         ) VAR(II)=J
      RETURN
      END
      REAL FUNCTION DIST(M, X, Y, N, P)
C
C<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
C
C   PURPOSE
C   -------
C
C      COMPUTES THE EUCLIDEAN DISTANCE BETWEEN EITHER TWO CASES OR TWO
C      VARIABLES
C
C   DESCRIPTION
C   -----------
C
C   1.  THE PARAMETERS X AND Y ARE ASSUMED TO BE EITHER A VECTOR OR A
C       STARTING LOCATION IN A MATRIX.  IF THEY ARE VECTORS OR THE
C       DISTANCE BETWEEN TWO COLUMNS OF A MATRIX ARE REQUIRED, SET N=1.
C       IF THE DISTANCE BETWEEN TWO ROWS OF A MATRIX ARE REQUIRED, SET
C       N EQUAL TO THE LEADING DIMENSION OF THE MATRIX AS DECLARED IN
C       THE MAIN PROGRAM.
C
C   2.  P IS A PARAMETER DECIDING WHICH DISTANCE FUNCTION TO USE.  IF P
C       IS GREATER THAN ZERO, THE MINKOWSKI DISTANCE IS COMPUTED,
C       RAISING THE DIFFERENCES TO THE P-TH POWER.  OTHERWISE, THE
C       DISTANCE RETURNED WILL BE THE PROPORTION OF CORRESPONDING
C       ELEMENTS OF X AND Y WHOSE DIFFERENCE IS GREATER THAN THE
C       ABSOLUTE VALUE OF P.
C
C   3.  MISSING VALUES ARE DENOTED BY 99999.  AND ARE IGNORED BY THE
C       DISTANCE CALCULATIONS.
C
C   INPUT PARAMETERS
C   ----------------
C
C   M     INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         THE NUMBER OF ELEMENTS.
C
C   X     REAL VECTOR DIMENSIONED AT LEAST (M-1)*N+1 (UNCHANGED ON
C            OUTPUT).
C         THE STARTING LOCATION OF THE FIRST VECTOR OF ELEMENTS.
C
C   Y     REAL VECTOR DIMENSIONED AT LEAST (M-1)*N+1 (UNCHANGED ON
C            OUTPUT).
C         THE STARTING LOCATION OF THE SECOND VECTOR OF ELEMENTS.
C
C   N     INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         THE SKIP FACTOR BETWEEN ELEMENTS.
C
C         IF THE DISTANCE BETWEEN COLUMNS OF A MATRIX IS DESIRED
C            OR IF X AND Y ARE INDIVIDUAL VECTORS, SET N = 1.
C
C         IF THE DISTANCE BETWEEN ROWS OF A MATRIX IS DESIRED,
C            SET N = DIM, WHERE DIM IS THE LEADING DIMENSION OF THE
C            ARRAY THE ELEMENTS X AND Y COME FROM.
C
C   P     REAL SCALAR (UNCHANGED ON OUTPUT).
C         P .GT. 0 IMPLIES THE P-TH POWER DISTANCE.
C            I.E. SUM EACH DIFFERENCE RAISED TO THE P-TH POWER, THEN
C            RETURN THE P-TH ROOT OF THE SUM AS THE DISTANCE
C         P .LE. 0 RETURNS THE PROPORTION OF DISTANCES GREATER THAN THE
C            ABSOLUTE VALUE OF P.
C
C   OUTPUT PARAMETER
C   ----------------
C
C   DIST  REAL SCALAR.
C         THE DISTANCE BETWEEN THE TWO VECTORS.
C
C   REFERENCE
C   ---------
C
C     HARTIGAN, J. A. (1975).  CLUSTERING ALGORITHMS, JOHN WILEY &
C        SONS, INC., NEW YORK.  PAGE 58-65, 69.
C
C<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
C
      DIMENSION X(*),Y(*)
      DD=0.
      DP=0.
      DO 10 I=1,M
         II = (I-1) * N + 1
         IF(X(II).NE.99999..AND.Y(II).NE.99999.) THEN
            DP=DP+1
            DIF=ABS(X(II)-Y(II))
            IF(P.GT.0.) DD=DD+DIF**P
            IF(P.LE.0..AND.DIF+P.GT.0.) DD=DD+1
         ENDIF
   10 CONTINUE
      IF(DP.NE.0.) THEN
         IF(P.GT.0.) DD=(DD/DP)**(1./P)
         IF(P.LE.0.) DD=DD/DP
      ENDIF
      DIST=DD
      RETURN
      END
      REAL FUNCTION R1MACH(I)
C***BEGIN PROLOGUE  R1MACH
C***DATE WRITTEN   790101   (YYMMDD)
C***REVISION DATE  860825   (YYMMDD)
C***CATEGORY NO.  R1
C***KEYWORDS  MACHINE CONSTANTS
C***AUTHOR  FOX, P. A., (BELL LABS)
C           HALL, A. D., (BELL LABS)
C           SCHRYER, N. L., (BELL LABS)
C***PURPOSE  Returns single precision machine dependent constants
C***DESCRIPTION
C
C     This is the CMLIB version of R1MACH, the real machine
C     constants subroutine originally developed for the PORT library.
C
C     R1MACH can be used to obtain machine-dependent parameters
C     for the local machine environment.  It is a function
C     subroutine with one (input) argument, and can be called
C     as follows, for example
C
C          A = R1MACH(I)
C
C     where I=1,...,5.  The (output) value of A above is
C     determined by the (input) value of I.  The results for
C     various values of I are discussed below.
C
C  Single-Precision Machine Constants
C  R1MACH(1) = B**(EMIN-1), the smallest positive magnitude.
C  R1MACH(2) = B**EMAX*(1 - B**(-T)), the largest magnitude.
C  R1MACH(3) = B**(-T), the smallest relative spacing.
C  R1MACH(4) = B**(1-T), the largest relative spacing.
C  R1MACH(5) = LOG10(B)
C***REFERENCES  FOX, P.A., HALL, A.D., SCHRYER, N.L, *FRAMEWORK FOR
C                 A PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHE-
C                 MATICAL SOFTWARE, VOL. 4, NO. 2, JUNE 1978,
C                 PP. 177-188.
C***ROUTINES CALLED  XERROR
C***END PROLOGUE  R1MACH
C
      INTEGER SMALL(2)
      INTEGER LARGE(2)
      INTEGER RIGHT(2)
      INTEGER DIVER(2)
      INTEGER LOG10(2)
C
      REAL RMACH(5)
C
      EQUIVALENCE (RMACH(1),SMALL(1))
      EQUIVALENCE (RMACH(2),LARGE(1))
      EQUIVALENCE (RMACH(3),RIGHT(1))
      EQUIVALENCE (RMACH(4),DIVER(1))
      EQUIVALENCE (RMACH(5),LOG10(1))
C
C
C     MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T
C     3B SERIES, MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T
C     PC 7300), SUN SPARCSTATIONS, SILICON GRAPHCS WORKSTATIONS, HP
C     9000 WORKSTATIONS, IBM RS/6000 WORKSTATIONS, AND 8087 BASED 
C     MICROS (E.G. IBM PC AND AT&T 6300).
C
C === MACHINE = ATT.3B
C === MACHINE = ATT.6300
C === MACHINE = ATT.7300
C === MACHINE = HP.9000
C === MACHINE = IBM.PC
C === MACHINE = IBM.RS6000
C === MACHINE = IEEE.MOST-SIG-BYTE-FIRST
C === MACHINE = IEEE.LEAST-SIG-BYTE-FIRST
C === MACHINE = SGI
C === MACHINE = SUN
C === MACHINE = 68000
C === MACHINE = 8087
C      DATA SMALL(1) /     8388608 /
C      DATA LARGE(1) /  2139095039 /
C      DATA RIGHT(1) /   864026624 /
C      DATA DIVER(1) /   872415232 /
C      DATA LOG10(1) /  1050288283 /
C
C     MACHINE CONSTANTS FOR SUN WORKSTATIONS.  f77 WITH -r8 OPTION.
C     MACHINE CONSTANTS FOR IBM RS/6000 WORKSTATIONS WITH -qautodbl=dblpad.
C
C === MACHINE = IBM.RS6000.XLF-WITH-AUTODBL-OPTION
C === MACHINE = SUN.F77-WITH-R8-OPTION
C === MACHINE = SGI.ORIGIN.F77-WITH-R8-D16-OPTIONS
C      DATA RMACH(1) / 2.2250738585072E-307 /
C      DATA RMACH(2) / 1.7976931348623E308  /
C      DATA RMACH(3) / 1.1102230246252E-16  /
C      DATA RMACH(4) / 2.2204460492503E-16  /
C      DATA RMACH(5) / 0.30102999566398     /
C
C     MACHINE CONSTANTS FOR AMDAHL MACHINES.
C
C === MACHINE = AMDAHL
C      DATA SMALL(1) /    1048576 /
C      DATA LARGE(1) / 2147483647 /
C      DATA RIGHT(1) /  990904320 /
C      DATA DIVER(1) / 1007681536 /
C      DATA LOG10(1) / 1091781651 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM.
C
C === MACHINE = BURROUGHS.1700
C      DATA RMACH(1) / Z400800000 /
C      DATA RMACH(2) / Z5FFFFFFFF /
C      DATA RMACH(3) / Z4E9800000 /
C      DATA RMACH(4) / Z4EA800000 /
C      DATA RMACH(5) / Z500E730E8 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 5700/6700/7700 SYSTEMS.
C
C === MACHINE = BURROUGHS.5700
C === MACHINE = BURROUGHS.6700
C === MACHINE = BURROUGHS.7700
C      DATA RMACH(1) / O1771000000000000 /
C      DATA RMACH(2) / O0777777777777777 /
C      DATA RMACH(3) / O1311000000000000 /
C      DATA RMACH(4) / O1301000000000000 /
C      DATA RMACH(5) / O1157163034761675 /
C
C     MACHINE CONSTANTS FOR THE CONVEX C1, C2, C3 SERIES (NATIVE MODE)
C
C === MACHINE = CONVEX
C      DATA RMACH(1) / 2.9387360E-39 /
C      DATA RMACH(2) / 1.7014117E+38 /
C      DATA RMACH(3) / 5.9604645E-08 /
C      DATA RMACH(4) / 1.1920929E-07 /
C      DATA RMACH(5) / 3.0102999E-01 /
C
C     MACHINE CONSTANTS FOR THE CONVEX C1, C2, C3 (NATIVE MODE)
C     WITH -P8 OPTION
C
C === MACHINE = CONVEX.P8
C      DATA RMACH(1) / 5.562684646268007E-309 /
C      DATA RMACH(2) / 8.988465674311577E+307 /
C      DATA RMACH(3) / 1.110223024625157E-016 /
C      DATA RMACH(4) / 2.220446049250313E-016 /
C      DATA RMACH(5) / 3.010299956639812E-001 /
C
C     MACHINE CONSTANTS FOR THE CONVEX C1, C2, C3 (IEEE MODE)
C
C === MACHINE = CONVEX.IEEE
C      DATA RMACH(1) / 1.1754945E-38 /
C      DATA RMACH(2) / 3.4028234E+38 /
C      DATA RMACH(3) / 5.9604645E-08 /
C      DATA RMACH(4) / 1.1920929E-07 /
C      DATA RMACH(5) / 3.0102999E-01 /
C
C     MACHINE CONSTANTS FOR THE CONVEX C1, C2, C3 (IEEE MODE)
C     WITH -P8 OPTION
C
C === MACHINE = CONVEX.IEEE.P8
C      DATA RMACH(1) / 2.225073858507202E-308 /
C      DATA RMACH(2) / 1.797693134862315E+308 /
C      DATA RMACH(3) / 1.110223024625157E-016 /
C      DATA RMACH(4) / 2.220446049250313E-016 /
C      DATA RMACH(5) / 3.010299956639812E-001 /
C
C     MACHINE CONSTANTS FOR THE CYBER 170/180 SERIES USING NOS (FTN5).
C
C === MACHINE = CYBER.170.NOS
C === MACHINE = CYBER.180.NOS
C      DATA RMACH(1) / O"00014000000000000000" /
C      DATA RMACH(2) / O"37767777777777777777" /
C      DATA RMACH(3) / O"16404000000000000000" /
C      DATA RMACH(4) / O"16414000000000000000" /
C      DATA RMACH(5) / O"17164642023241175720" /
C
C     MACHINE CONSTANTS FOR THE CDC 180 SERIES USING NOS/VE
C
C === MACHINE = CYBER.180.NOS/VE
C      DATA RMACH(1) / Z"3001800000000000" /
C      DATA RMACH(2) / Z"4FFEFFFFFFFFFFFE" /
C      DATA RMACH(3) / Z"3FD2800000000000" /
C      DATA RMACH(4) / Z"3FD3800000000000" /
C      DATA RMACH(5) / Z"3FFF9A209A84FBCF" /
C
C     MACHINE CONSTANTS FOR THE CYBER 205
C
C === MACHINE = CYBER.205
C      DATA RMACH(1) / X'9000400000000000' /
C      DATA RMACH(2) / X'6FFF7FFFFFFFFFFF' /
C      DATA RMACH(3) / X'FFA3400000000000' /
C      DATA RMACH(4) / X'FFA4400000000000' /
C      DATA RMACH(5) / X'FFD04D104D427DE8' /
C
C     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES.
C
C === MACHINE = CDC.6000
C === MACHINE = CDC.7000
C      DATA RMACH(1) / 00014000000000000000B /
C      DATA RMACH(2) / 37767777777777777777B /
C      DATA RMACH(3) / 16404000000000000000B /
C      DATA RMACH(4) / 16414000000000000000B /
C      DATA RMACH(5) / 17164642023241175720B /
C
C     MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3.
C
C === MACHINE = CRAY.46-BIT-INTEGER
C === MACHINE = CRAY.64-BIT-INTEGER
C      DATA RMACH(1) / 200034000000000000000B /
C      DATA RMACH(2) / 577767777777777777776B /
C      DATA RMACH(3) / 377224000000000000000B /
C      DATA RMACH(4) / 377234000000000000000B /
C      DATA RMACH(5) / 377774642023241175720B /
C
C     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200
C
C     NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING LINE -
C     STATIC RMACH(5)
C
C === MACHINE = DATA_GENERAL.ECLIPSE.S/200
C      DATA SMALL/20K,0/,LARGE/77777K,177777K/
C      DATA RIGHT/35420K,0/,DIVER/36020K,0/
C      DATA LOG10/40423K,42023K/
C
C     ELXSI 6400
C
C === MACHINE = ELSXI.6400
C      DATA SMALL(1) / '00800000'X /
C      DATA LARGE(1) / '7F7FFFFF'X /
C      DATA RIGHT(1) / '33800000'X /
C      DATA DIVER(1) / '34000000'X /
C      DATA LOG10(1) / '3E9A209B'X /
C
C     MACHINE CONSTANTS FOR THE HARRIS 220
C     MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7
C
C === MACHINE = HARRIS.220
C === MACHINE = HARRIS.SLASH6
C === MACHINE = HARRIS.SLASH7
C      DATA SMALL(1),SMALL(2) / '20000000, '00000201 /
C      DATA LARGE(1),LARGE(2) / '37777777, '00000177 /
C      DATA RIGHT(1),RIGHT(2) / '20000000, '00000352 /
C      DATA DIVER(1),DIVER(2) / '20000000, '00000353 /
C      DATA LOG10(1),LOG10(2) / '23210115, '00000377 /
C
C     MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES.
C     MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES.
C
C === MACHINE = HONEYWELL.600/6000
C === MACHINE = HONEYWELL.DPS.8/70
C      DATA RMACH(1) / O402400000000 /
C      DATA RMACH(2) / O376777777777 /
C      DATA RMACH(3) / O714400000000 /
C      DATA RMACH(4) / O716400000000 /
C      DATA RMACH(5) / O776464202324 /
C
C     MACHINE CONSTANTS FOR THE HP 2100
C     3 WORD DOUBLE PRECISION WITH FTN4
C
C === MACHINE = HP.2100.3_WORD_DP
C      DATA SMALL(1), SMALL(2) / 40000B,       1 /
C      DATA LARGE(1), LARGE(2) / 77777B, 177776B /
C      DATA RIGHT(1), RIGHT(2) / 40000B,    325B /
C      DATA DIVER(1), DIVER(2) / 40000B,    327B /
C      DATA LOG10(1), LOG10(2) / 46420B,  46777B /
C
C     MACHINE CONSTANTS FOR THE HP 2100
C     4 WORD DOUBLE PRECISION WITH FTN4
C
C === MACHINE = HP.2100.4_WORD_DP
C      DATA SMALL(1), SMALL(2) / 40000B,       1 /
C      DATA LARGE91), LARGE(2) / 77777B, 177776B /
C      DATA RIGHT(1), RIGHT(2) / 40000B,    325B /
C      DATA DIVER(1), DIVER(2) / 40000B,    327B /
C      DATA LOG10(1), LOG10(2) / 46420B,  46777B /
C
C     MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
C     THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86 AND
C     THE INTERDATA 3230 AND INTERDATA 7/32.
C
C === MACHINE = IBM.360
C === MACHINE = IBM.370
C === MACHINE = XEROX.SIGMA.5
C === MACHINE = XEROX.SIGMA.7
C === MACHINE = XEROX.SIGMA.9
C === MACHINE = SEL.85
C === MACHINE = SEL.86
C === MACHINE = INTERDATA.3230
C === MACHINE = INTERDATA.7/32
C      DATA RMACH(1) / Z00100000 /
C      DATA RMACH(2) / Z7FFFFFFF /
C      DATA RMACH(3) / Z3B100000 /
C      DATA RMACH(4) / Z3C100000 /
C      DATA RMACH(5) / Z41134413 /
C
C     MACHINE CONSTANTS FOR THE INTERDATA 8/32
C     WITH THE UNIX SYSTEM FORTRAN 77 COMPILER.
C
C     FOR THE INTERDATA FORTRAN VII COMPILER REPLACE
C     THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S.
C
C === MACHINE = INTERDATA.8/32.UNIX
C      DATA RMACH(1) / Z'00100000' /
C      DATA RMACH(2) / Z'7EFFFFFF' /
C      DATA RMACH(3) / Z'3B100000' /
C      DATA RMACH(4) / Z'3C100000' /
C      DATA RMACH(5) / Z'41134413' /
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KA OR KI PROCESSOR).
C
C === MACHINE = PDP-10.KA
C === MACHINE = PDP-10.KI
C      DATA RMACH(1) / "000400000000 /
C      DATA RMACH(2) / "377777777777 /
C      DATA RMACH(3) / "146400000000 /
C      DATA RMACH(4) / "147400000000 /
C      DATA RMACH(5) / "177464202324 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C     32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
C
C === MACHINE = PDP-11.32-BIT
C      DATA SMALL(1) /    8388608 /
C      DATA LARGE(1) / 2147483647 /
C      DATA RIGHT(1) /  880803840 /
C      DATA DIVER(1) /  889192448 /
C      DATA LOG10(1) / 1067065499 /
C
C      DATA RMACH(1) / O00040000000 /
C      DATA RMACH(2) / O17777777777 /
C      DATA RMACH(3) / O06440000000 /
C      DATA RMACH(4) / O06500000000 /
C      DATA RMACH(5) / O07746420233 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C     16-BIT INTEGERS  (EXPRESSED IN INTEGER AND OCTAL).
C
C === MACHINE = PDP-11.16-BIT
C      DATA SMALL(1),SMALL(2) /   128,     0 /
C      DATA LARGE(1),LARGE(2) / 32767,    -1 /
C      DATA RIGHT(1),RIGHT(2) / 13440,     0 /
C      DATA DIVER(1),DIVER(2) / 13568,     0 /
C      DATA LOG10(1),LOG10(2) / 16282,  8347 /
C
C      DATA SMALL(1),SMALL(2) / O000200, O000000 /
C      DATA LARGE(1),LARGE(2) / O077777, O177777 /
C      DATA RIGHT(1),RIGHT(2) / O032200, O000000 /
C      DATA DIVER(1),DIVER(2) / O032400, O000000 /
C      DATA LOG10(1),LOG10(2) / O037632, O020233 /
C
C     MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000.
C
C === MACHINE = SEQUENT.BALANCE.8000
C      DATA SMALL(1) / $00800000 /
C      DATA LARGE(1) / $7F7FFFFF /
C      DATA RIGHT(1) / $33800000 /
C      DATA DIVER(1) / $34000000 /
C      DATA LOG10(1) / $3E9A209B /
C
C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES.
C
C === MACHINE = UNIVAC.1100
C      DATA RMACH(1) / O000400000000 /
C      DATA RMACH(2) / O377777777777 /
C      DATA RMACH(3) / O146400000000 /
C      DATA RMACH(4) / O147400000000 /
C      DATA RMACH(5) / O177464202324 /
C
C     MACHINE CONSTANTS FOR THE VAX 11/780
C    (EXPRESSED IN INTEGER AND HEXADECIMAL)
C  *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS***
C
C === MACHINE = VAX.11/780
C      DATA SMALL(1) /       128 /
C      DATA LARGE(1) /    -32769 /
C      DATA RIGHT(1) /     13440 /
C      DATA DIVER(1) /     13568 /
C      DATA LOG10(1) / 547045274 /
C
C  ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSTEMS***
C
C      DATA SMALL(1) / Z00000080 /
C      DATA LARGE(1) / ZFFFF7FFF /
C      DATA RIGHT(1) / Z00003480 /
C      DATA DIVER(1) / Z00003500 /
C      DATA LOG10(1) / Z209B3F9A /
C
C
C***FIRST EXECUTABLE STATEMENT  R1MACH
C
      R1MACH = RMACH(I)
      RETURN
C
      END
      SUBROUTINE TREE1(M, K, DMNB, NB, RL, IERR, OUNIT)
C
C<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
C
C   PURPOSE
C   -------
C
C      PRINTS A TREE OF OBJECTS
C
C   DESCRIPTION
C   -----------
C
C   1.  THERE EXISTS AN ORDERING OF THE OBJECTS SUCH THAT EVERY CLUSTER
C       CONSISTS OF A SET OF OBJECTS CONTIGUOUS IN THE ORDER.  THE
C       ALGORITHM IS GIVEN ON PAGE 156 OF THE FIRST REFERENCE.  THE
C       OBJECTS ARE STORED IN THE VECTOR RL IN SUCH A SORTED ORDER.
C
C   2.  THIS ORDERING OF THE OBJECTS ALLOWS THE CLUSTERS TO BE NAMED BY
C       GIVING THE LOCATION OF THE FIRST AND LAST OBJECTS IN THE ARRAY
C       FOR EACH CLUSTER.  THE FIRST TWO COLUMNS OF THE NB ARRAY STORE
C       THESE FIRST AND LAST OBJECTS AND THE THIRD COLUMN STORE THE
C       DIAMETER OF THE CLUSTER, OR THE DISTANCE BETWEEN THE FIRST AND
C       LAST OBJECTS.
C
C   3.  FOR EACH OBJECT, SEARCH EACH CLUSTER PLACING A MARK AT THE
C       DIAMETER DISTANCE OF EACH CLUSTER THE OBJECT BELONGS TO.  THE
C       MARK IS A COMMA IF THE OBJECT IS ON THE BORDER OF A CLUSTER OR
C       A QUOTE MARK IF THE OBJECT IS IN THE INTERIOR OF A CLUSTER.
C       THE CLUSTERS CAN BE RECOVERED FROM THE DIAGRAM BY LINKING THE
C       QUOTE MARKS TO THE COMMAS ON BOTH SIDES OF THEM.
C
C   INPUT PARAMETERS
C   ----------------
C
C   M     INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         THE NUMBER OF OBJECTS.
C
C   K     INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         THE NUMBER OF CLUSTERS.
C
C   DMNB  INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         LEADING DIMENSION OF MATRIX NB.  MUST BE AT LEAST 3.
C
C   NB    REAL MATRIX WHOSE FIRST DIMENSION MUST BE DMNB AND SECOND
C            DIMENSION MUST BE AT LEAST K (UNCHANGED ON OUTPUT).
C         THE MATRIX DEFINING THE BOUNDARIES OF THE CLUSTERS.
C
C         NB(1,I) IS 1 + THE FIRST OBJECT IN CLUSTER I
C         NB(2,I) IS 1 + THE LAST OBJECT IN CLUSTER I
C         NB(3,I) IS THE CLUSTER DIAMETER
C
C   RL    VECTOR OF 4-CHARACTER VARIABLES DIMENSIONED AT LEAST M+1
C            (UNCHANGED ON OUTPUT).
C         RL(1) IS THE LABEL OF THE TREE.
C         RL(I+1) IS THE LABEL OF THE I-TH OBJECT.
C
C   OUNIT INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         UNIT NUMBER FOR OUTPUT.
C
C   OUTPUT PARAMETER
C   ----------------
C
C   IERR  INTEGER SCALAR.
C         ERROR FLAG.
C
C         IERR = 0, NO ERRORS WERE DETECTED DURING EXECUTION
C
C         IERR = 1, EITHER THE FIRST AND LAST CASES OR THE CLUSTER
C                   DIAMETER FOR A CLUSTER IS OUT OF BOUNDS.  THE
C                   CLUSTER ITS VALUES ARE PRINTED ON UNIT OUNIT.
C                   EXECUTION WILL CONTINUE WITH QUESTIONABLE RESULTS
C                   FOR THAT CLUSTER.
C
C   REFERENCES
C   ----------
C
C     HARTIGAN, J. A. (1975).  CLUSTERING ALGORITHMS, JOHN WILEY &
C        SONS, INC., NEW YORK.  PAGE 165.
C
C     HARTIGAN, J. A. (1975) PRINTER GRAPHICS FOR CLUSTERING. JOURNAL OF
C        STATISTICAL COMPUTATION AND SIMULATION. VOLUME 4,PAGES 187-213.
C
C<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
C
      INTEGER DMNB, OUNIT
      DIMENSION NB(DMNB,*)
      CHARACTER*4 RL(*)
      CHARACTER*1 A,B,C,D,E,AA(101)
      DATA A,B,C,D,E/',','-','''','/',' '/
C
C     CHECK NB ARRAY FOR ILLEGAL VALUES
C
      IF (OUNIT .LE. 0) RETURN
      DO 10 KK=1,K
         I=NB(1,KK)
         J=NB(2,KK)
         L=NB(3,KK)
         IF(I.LT.1.OR.I.GT.J.OR.J.GT.M.OR.L.LT.0.OR.L.GT.100) THEN
            WRITE(OUNIT,1) KK,I,J,L
            IERR = 1
         ENDIF
   10 CONTINUE
    1 FORMAT(' ERROR IN CLUSTER',I4,' FIRST OBJ = ',I3,' LAST OBJ = '
     *,I3,' DIAM = ',I3)
C
C     PRINT TREE, CASE BY CASE
C
      DO 50 I=1,M
         DO 20 L=1,101
   20      AA(L)=E
         DO 40 KK=1,K
C
C     DETERMINE IF CLUSTER KK CONTAINS CASE I
C
            IF(NB(1,KK).LE.I.AND.NB(2,KK).GE.I) THEN
               J=NB(3,KK)+1
               AA(J)=C
               IF(NB(1,KK).EQ.I.OR.NB(2,KK).EQ.I) THEN
                  DO 30 L=1,J
   30                IF(AA(L).EQ.E) AA(L)=B
                  IF(NB(2,KK).EQ.I) AA(J)=C
                  IF(NB(1,KK).EQ.I) AA(J)=A
               ENDIF
            ENDIF
   40    CONTINUE
         WRITE(OUNIT,2) I,RL(I),(AA(J),J=1,101)
    2    FORMAT(I5,4X,A4,1X,101A1)
   50 CONTINUE
      RETURN
      END
      SUBROUTINE AMALG(M, X, Y, N)
C
C<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
C
C   PURPOSE
C   -------
C
C      AMALGAMATES TWO OBJECTS
C
C   DESCRIPTION
C   -----------
C
C   1.  THE PARAMETERS X AND Y ARE ASSUMED TO BE EITHER A VECTOR OR A
C       STARTING LOCATION IN A MATRIX.  IF THEY ARE VECTORS OR THE
C       AMALGAMATION OF TWO COLUMNS OF A MATRIX ARE REQUIRED, SET N=1.
C       IF THE AMALGAMATION OF TWO ROWS OF A MATRIX ARE REQUIRED, SET N
C       EQUAL TO THE LEADING DIMENSION OF THE MATRIX AS DECLARED IN THE
C       MAIN PROGRAM.
C
C   2.  BOTH OBJECTS WILL BE SET TO THE ARITHMETIC MEAN OF THE
C       CORRESPONDING ELEMENTS OF BOTH OBJECTS.
C
C   INPUT PARAMETERS
C   ----------------
C
C   M     INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         THE NUMBER OF ELEMENTS.
C
C   X     REAL VECTOR DIMENSIONED AT LEAST M*N (CHANGED ON OUTPUT).
C         THE STARTING LOCATION OF THE FIRST VECTOR OF ELEMENTS.
C
C   Y     REAL VECTOR DIMENSIONED AT LEAST M*N (CHANGED ON OUTPUT).
C         THE STARTING LOCATION OF THE SECOND VECTOR OF ELEMENTS.
C
C   N     INTEGER SCALAR (UNCHANGED ON OUTPUT).
C         THE SKIP FACTOR BETWEEN DESIRED ELEMENTS.
C
C         IF THE DISTANCE BETWEEN COLUMNS OF A MATRIX IS DESIRED,
C            SET N = 1.
C         IF THE DISTANCE BETWEEN ROWS OF A MATRIX IS DESIRED,
C            SET N = DIM, WHERE DIM IS THE LEADING DIMENSION OF THE
C            ARRAY THE ELEMENTS X AND Y COME FROM.
C
C   OUTPUT PARAMETERS
C   -----------------
C
C      X AND Y HOLD THE AMALGAMATED VECTORS.
C
C   REFERENCE
C   ---------
C
C     HARTIGAN, J. A. (1975).  CLUSTERING ALGORITHMS, JOHN WILEY &
C        SONS, INC., NEW YORK.  PAGE 232.
C
C<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
C
      DIMENSION X(*), Y(*)
C
      DO 10 I=1,M
         II = (I-1) * N + 1
         X(II)=(X(II)+Y(II))/2.
   10    Y(II)=X(II)
      RETURN
      END
