C	WESTERN MICHIGAN UNIVERSITY
C	CLUSTR.FOR (FILENAME ON LIBRARY DECTAPE)
C	CLUSTR, 1.12.1 (CALLING NAME, SUBLST NO.)
C	SINGLE-LINK CLUSTER ANALYSIS PACKAGE
C	THIS PROGRAM IS A COMBINATION OF ROUTINES FROM SEVERAL SOURCES
C	RUSS BARR PROGRAMMEDD THE MAIN PROG. AND ALSO DID SOME
C	 SUBSTANTIAL ADAPTIVE PROGRAMMING ON OTHER PARTS OF THE PACKAGE.
C	LIBRARY DECTAPE PROGS. USED:  USAGE.MAC
C	APLB 10 PROGS. USED:  IO, GETFOR
C	FORWMU PROGS. USED:  ALLCOR, DEVCHG, EXISTS, PRINTS, GES,
C	 DEVCHR, GETPPN, RUNUUO, JOBNUM
C	INTERNAL SUBR. USED:  MNL, READD, SLINK, SHADE, SORT,
C	 SLINK2, DENDRG, DENDHD, TRANSP, SLINK1
C	FUNCTIONS USED:  SLINK3
C	ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
C
	DOUBLE PRECISION IFLNMI,IFLNMO,DEVNAM,FILNM,DEVN

C---------------COMMON /INIO/ AND COMMON /IOB/ ARE SHARED BY
C--------------- BNKLIB.FOR, SUBR. MNL, READD, SHADE
C---------------COMMON /SPEC/ SHARED BY SUBR. MNL, READD, SHADE
	COMMON/INIO/IFTR,IFTW,DEVN(30),FILNM(30),IPP(30),DEST(30)
	COMMON/IOB/LEFBK,IRTBK,IART,MAXPAG,IPAGE,IPAGCT,IDLG,IRSP,
     1		   IDUMMY,JDUMMY
	COMMON/SPEC/IDVO,IDVI,NDEVO,NDEVI

	DIMENSION INPKOD(7)
C---------------SEE ST. 4202 AND 4202+1 FOR CLIM(16).
	DIMENSION CLIM(16),IFMT(80),TITLE(16)
	DATA INPKOD /7,1,4,5,6,2,3/
C
	DIMENSION S(1)
C
C
C*******NDEVI AND NDEVO ARE THE LOGICAL UNIT NUMBERS FOR
C*******THE INPUT AND OUTPUT RESPECTIVELY.
C*******IDLG AND IRSP ARE THE LOGICAL UNIT NUMBERS OF
C*******THE INPUT AND OUTPUT DIALOGUE CHANNELS REPECTIVELY.
C---------------IDLG, IRSP ARE PASSED TO IO(IN APLB10) THRU
C--------------- COMMON /IOB/
C
	IDLG=-1
	IRSP=-4
	NDEVI=4
	NDEVO=6
	IPAGCT = -1
	WRITE(IDLG,96)
96	FORMAT(/,' WMU SINGLE-LINK CLUSTER ANALYSIS',//)
C	CALL USAGEB('CLUSTR')
C
C---------------1 MEANS OUTPUT PROMPTING IS REQUESTED, NDEVO
C--------------- IS DEVICE NO., DEVNAM RETURNED BY SUBR. IO AND
C--------------- IS DEVICE NAME (MUST BE DOUBL PRECISION), IDVO
	CALL IO(1,NDEVO,DEVNAM,IDVO,IFLNMO,IPJ,IPG,IBNK)
98	CALL IO(0,NDEVI,DEVNAM,IDVI,IFLNMI,IPJ,IPG,IBNK)
C
	WRITE(IDLG,116)
116	FORMAT(' ENTER OUTPUT IDENTIFICATION IF DESIRED, ELSE TYPE',
     1	' <RETURN>',/)
	READ(IRSP,20)TITLE
20	FORMAT(16A5)
100	WRITE(IDLG,102)
102	FORMAT(' ENTER TYPE OF INPUT, TYPE <RETURN> FOR HELP',/)
	READ(IRSP,104)INPRD
104	FORMAT(10I)
	IF(INPRD.GT.0.AND.INPRD.LE.7)GO TO 108
	WRITE(IDLG,106)
C*******NOTE: CONTINUATION COLUMN OF FORMAT 106 INDICATES INTERNAL
C*******REPRESENTATION OF INPUT(INPRD).
106	FORMAT(' TYPE:',/,
     7	'  1 FOR RAW DATA, CORRELATION COMPUTED',/,
     1	'  2 FOR RAW DATA, EUCLIDEAN DISTANCES COMPUTED',/,
     4	'  3 FOR RAW DATA, MANHATTEN METRIC COMPUTED',/,
     5	'  4 FOR RAW DATA, MINKOWSKY METRIC COMPUTED',/,
     6	'  5 FOR RAW DATA, CHI SQUARE COMPUTED',/,
     2	'  6 FOR LOWER TRIANGULAR MATRIX OF DISTANCES',/,
     3	'  7 FOR LOWER TRIANGULAR MATRIX OF CORRELATION',/)

	GO TO 100
108	IF(INPRD.EQ.2)POWER=2
	IF(INPRD.EQ.3)POWER=1
	WRITE(IDLG,110)
110	FORMAT(' ENTER CLUSTER LIMITS CONTROL CODE, TYPE <RETURN>',
     1	' FOR HELP',/)
	READ(IRSP,104)KODE
	IF(KODE.GT.0.AND.KODE.LE.4)GO TO 114
	WRITE(IDLG,112)
112	FORMAT(' TYPE:',/,
     1	'  1 TO DIVIDE THE RANKS OF DISTANCES INTO 15 EQUAL INTERVALS',/,
     1	'  2 TO DIVIDE THE ACTUAL DISTANCES INTO 15 EQUAL INTERVALS',/,
     1	'  3 TO SPECIFY THE INTERVAL MAX AND MIN',/,
     1	'  4 TO ENTER THE ENDPOINTS OF THE INTERVALS',/)
	GO TO 108
114	CONTINUE
700	WRITE(IDLG,702)
702	FORMAT(' ENTER METHOD(S) OF OUTPUT, TYPE <RETURN> FOR HELP',/)
	READ(IRSP,104)ITT
	IF(ITT.NE.0)GO TO 703
	WRITE(IDLG,705)
705	FORMAT(' TYPE:',/,'  1 FOR SHADE DIAGRAM',/,
     1	'  2 FOR DENDROGRAM',/,'  3 FOR BOTH',/)
	GO TO 700
703	IF(ITT.LT.1.OR.ITT.GT.3)GO TO 700
	CORFLG=0
	IF(INPRD.EQ.6.OR.INPRD.EQ.7)CORFLG=-1
704	INPUT=INPKOD(INPRD)
C
C*******INPUT=CONTROL CODE FOR FORM OF DATA FOR INPUT.
C*******SEE SUBROUTINE READD.
C*******KODE - SEE COMMENT AFTER STATEMENT 330 OF SUBROUTINE READD.
C
	CALL GETFOR(IDLG,IRSP,IFMT,ISTD,80,2)
	IF(ISTD.EQ.1)IFMT(1)='(10F)'
	GO TO (600,620,620,600,600,600,600),INPUT
600	WRITE(IDLG,602)
602	FORMAT(' ENTER NUMBER OF ROWS AND COLUMNS(I.E. 13,5) ',$)
	READ(IRSP,604)N,M
604	FORMAT(2I)
	IF(M.GT.0.AND.N.GT.0)GO TO 608
	WRITE(IDLG,606)
606	FORMAT(' ?BOTH M AND N MUST BE GREATER THAN ZERO',/)
	GO TO 600
608	IF(INPUT.NE.5)GO TO 620
	WRITE(IDLG,1008)
1008	FORMAT(' ENTER POWER FOR MINKOWSKY METRIC',/)
	READ(IRSP,1006)POWER
1006	FORMAT(2F)
620	GO TO (650,610,610,650,650,650,650),INPUT
610	WRITE(IDLG,2002)
2002	FORMAT(' HOW MANY VARIABLES? ',$)
	READ(IRSP,1004)M
1004	FORMAT(2I)
650	GO TO (665,665,4001,660),KODE
4001	WRITE(IDLG,4002)
4002	FORMAT(' ENTER CLUSTER LIMITS(MIN,MAX) ',$)
	READ(IRSP,1006)DMIN,DMAX
	GO TO 665
660	WRITE(IDLG,4202)
4202	FORMAT(' ENTER END POINTS OF CLUSTERS(8 PER LINE)',/)
	READ(IRSP,4204)(CLIM(I),I=1,16)
4204	FORMAT(8F)
665	IANS=1
	IF(CORFLG.EQ.-1)GO TO 670
2004	WRITE(IDLG,2006)
2006	FORMAT(' TYPE:',/,'  1 FOR CLUSTERING BY VARIABLE',/,
     1	'  2 FOR CLUSTERING BY OBSERVATION',/,'  3 FOR BOTH',/)
	READ(IRSP,2008)IANS
2008	FORMAT(I)
	IF(IANS.GE.1.AND.IANS.LE.3)GO TO 670
	WRITE(IDLG,2010)
2010	FORMAT(' ?IMPROPER RESPONSE',/)
	GO TO 2004
C
C*******DYNAMIC ALLOCATION SECTION
C
670	IF(N.EQ.0)N=1
C*******4 INDICATES 1ST PASS
C*******IANS IS AN OCTAL NUMBER REPRESENTING 3 BINARY BITS
C*******THE BITS CONTROL THE TASKS OPERATIONS IN DECREASING PRECEDENCE
C*******BIT 4	(LEFTMOST)	READ IN DATA.
C*******BIT 1	PROCESS DATA WITHOUT TRANSPOSING.
C*******BIT 2	TRANSPOSE DATA BEFORE PROCESSING.
C*******(IF BOTH 1 AND 2 SET DO BOTH OPERATIONS).
	IANS=IANS.OR.4
C*******SAVE PROPER VALUES OF M
690	CONTINUE
	GO TO 697
	GO TO (695,696,697),IANS-4
695	IF(INPRD.EQ.1)GO TO 6961
6951	K=N
	GO TO 691
696	IF(INPRD.EQ.1)GO TO 6951
6961	K=M
	GO TO 691
697	K=MAX0(M,N)
691	L=K*(K-1)/2
	MAX=K*K+L+K+K+K+K+MAX0(L,4*K)+2*K+M*N
C*******D START(K*K)
	CALL ALLCOR(MAX,IERR,I1,S)
	IF(IERR.EQ.-1.AND.INPRD.NE.1)GO TO 672
	IF(IERR.EQ.-1.AND.INPRD.EQ.1)GO TO 680
C*******STRING START (L)
	I2=I1+K*K
C*******NCODE START (K)
	I3=I2+L
C*******HA START (K)
	I4=I3+K
C*******HB START (K)
	I5=I4+K
C*******KSCR START (K)
	I6=I5+K
C*******XSCR START (MAX0(L,4*K))
	I7=I6+K
C*******ALABEL START (2*K)
	I8=I7+MAX0(L,4*K)
C*******X START (M*N)
	I9=I8+2*K
C*******I10=MAX
	I10=I9+M*N
C*******RESTORE PROPER VALUE OF M
671	CALL MNL(S(I1),S(I2),S(I3),S(I4),S(I5),S(I6),S(I7),S(I8),S(I9),
     1	CLIM,IFMT,TITLE,N,INPUT,KODE,M,POWER,CORFLG,DMIN,DMAX,ITT,IANS)
	GO TO 98
672	WRITE(IDLG,673)
673	FORMAT(' ?NOT ENOUGH CORE AVAILABLE.',/)
	GO TO 98
680	IF(IANS.NE.5)GO TO 672
	WRITE(IDLG,682)
682	FORMAT(' ?NOT ENOUGH ROOM AVAILABLE TO KEEP DATA IN CORE-',
     1	'CONTINUING',/)
	CORFLG=-1
	L=M*(M-1)/2
	MAX=M*M+L+M+M+M+M+4+MAX0(4,M)+2*M+M
	CALL ALLCOR(MAX,IERR,I1,S)
	IF(IERR.EQ.-1)GO TO 672
	I2=I1+M*M
	I3=I2+L
	I4=I3+M
	I5=I4+M
	I6=I5+M
	I7=I6+M
	I8=I7+4*MAX0(4,M)
	I9=I8+2*M
	I10=I9+M
	GO TO 671
	END
C******************************************************************************
C---------------NCODE, ALABEL MODIFIED BY SUBR. MNL.  SPACE
C--------------- PROVIDED BY DYN. ALLOC. FOR D, STRING, NCODE,
C--------------- HA, HB, KSCR, XSCR, ALABEL, X.CLIM, IFMT,
C--------------- TITLE, N, INPUT, KODE, M, POWER, CORFLG, 
C--------------- DMIN, DMAX, ITT, IANS ARE INPUT.  CORFLG
C--------------- APPARENTLY NOT USED.
C---------------NDEVO, IDVO, INPUT THRU COMMON /SPEC/
	SUBROUTINE MNL(D,STRING,NCODE,HA,HB,KSCR,XSCR,ALABEL,X,CLIM,
     1	IFMT,TITLE,N,INPUT,KODE,M,POWER,CORFLG,DMIN,DMAX,ITT,IANS)
C
	DOUBLE PRECISION ALABEL(1),FILNM,DEVN
	COMMON/INIO/IFTR,IFTW,DEVN(30),FILNM(30),IPP(30),DEST(30)
	COMMON/IOB/LEFBK,IRTBK,IALT,MAXPAG,IPAG,IPAGCT,IDLG,IRSP,
     1		   IDUMMY,JDUMMY
	COMMON /SPEC/IDVO,IDVI,NDEVO,NDEVI
C
	DIMENSION D(1),STRING(1),X(1),NCODE(1),CLIM(16),IFMT(80),
     1	TITLE(16)
	DIMENSION HA(1),HB(1),KSCR(1),XSCR(1)
C
C*******DUPLICATED D CAUSES:
C*******EQUIVALENCE (DD(1),IJ(1))
C*******IN SUBR READ
C
2012	CALL READD(INPUT,KODE,POWER,IANS,N,M,STRING,XSCR,HB,X,D,
     1	CLIM,IFMT,TITLE,DMIN,DMAX,D)
C
C*******M GIVES THE NUMBER OF VARIABLES, AND
C*******IT MUST BE GREATER THAN ZERO.
C*******THE CONTENTS OF ARRAY NCODE CONSTITUTE AN
C*******ORDERING OF THESE OBSERVATIONS USUALLY ONE WHICH
C*******CORRESPONDS TO THE ORDER OF OBSERVATIONS IN A DENDROGRAM
C*******OBTAINED BY SOME TECHNIQUE FOR CLUSTER ANALYSIS.
C*******NTIMES EQUALS THE NUMBER OF TIMES THE DISTANCE
C*******MATRIX IS TO BE PLOTTED (USING DIFFERENT PERMUTATIONS
C*******OF THE OBSERVATIONS.)
C
	IF(IANS.EQ.2)GO TO 711
	IF(ITT.EQ.2)GO TO 706
	WRITE(IDLG,602)
602	FORMAT(' DO YOU WISH TO HAVE A SHADE DIAGRAM',/,
     1	' BEFORE CLUSTERING?(YES OR NO) ',$)
	READ(IRSP,710)BANS
706	WRITE(IDLG,708)
708	FORMAT(' DO YOU WISH TO HAVE LABELS IN THE OUTPUT?',
     1	'(YES OR NO) ',$)
	READ(IRSP,710)LANS
710	FORMAT(A3)
711	IF(LANS.NE.'YES')GO TO 716
	WRITE(IDLG,712)M
712	FORMAT(' ENTER THE',I4,' LABELS, ONE PER LINE(10 CHARACTERS',
     1	' MAX.',/)
	READ(IRSP,714)(ALABEL(I),I=1,M)
714	FORMAT(A10)
716	CONTINUE
	IF(BANS.NE.'YES')GO TO 618
	IF(ITT.EQ.2)GO TO 618
	IF(M.LE.100)GO TO 760
	WRITE(NDEVO,758)
758	FORMAT(/,' SHADE DIAGRAM TOO LARGE TOO PRINT',/)
	GO TO 618
760	WRITE(NDEVO,40)TITLE
	WRITE(NDEVO,762)
762	FORMAT(/,' SHADE DIAGRAM BEFORE CLUSTERING',/)
	DO 766 I=1,M
766	NCODE(I)=I
C*******DUPLICATED D CAUSES
C******  EQUIVALENCE (D(1),DD(1))
C*******IN SUBROUTINE SHADE.
	CALL SHADE(ALABEL,D,D,XSCR,NCODE,INPUT,KODE,M,CLIM)
C
618	CONTINUE
	CALL SLINK(KSCR,XSCR,HA,HB,NCODE,STRING,M,DELHAT)
C
	WRITE(NDEVO,40)TITLE
	WRITE(NDEVO,7531)DELHAT
7531	FORMAT(/,' DELTA-ONE-HAT =  ',E,/)
	WRITE(NDEVO,764)
764	FORMAT(/,' SHADE DIAGRAM AFTER CLUSTERING',/)
	IF(INPUT.EQ.2.OR.INPUT.EQ.3)GO TO 753
	IF((IANS.AND.1).EQ.1)WRITE(NDEVO,751)
751	FORMAT(/,' CLUSTERING BY VARIABLE',/)
	IF((IANS.AND.1).EQ.0)WRITE(NDEVO,752)
752	FORMAT(/,' CLUSTERING BY OBSERVATION',/)
753	IF(ITT.NE.1.AND.ITT.NE.3)GO TO 726
	TYPE 9123,M
9123	FORMAT(' AT 2ND SHADE M=',I10,/)
	IF(M.GT.100)GO TO 730
C*******2 ARGUMENTS D GIVES AN EQUIVALENCE(D(M,M),DD(M,M)) IN SHADE
	CALL SHADE(ALABEL,D,D,XSCR,NCODE,INPUT,KODE,M,CLIM)
726	IF(ITT.NE.2.AND.ITT.NE.3)GO TO 750
	IF(M.GT.50.AND.IDVO.EQ.'TTY')GO TO 740
	CALL DENDHD(NCODE,XSCR,M,NDEVO,ALABEL,LANS)
	CALL DENDRG(HA,M,KSCR,XSCR,NDEVO,IDVO,2)
	GO TO 750
C*******THIS ROUTE IS AN EXIT FROM MNL.
730	WRITE(NDEVO,719)
719	FORMAT(/,' SHADE PLOT TOO LARGE TO PRINT',//)
	GO TO 726
740	WRITE(NDEVO,725)
725	FORMAT(/,' DENDROGRAM TOO LARGE TO PRINT',//)
717	DO 718 I=1,(M/6)*6,6
	WRITE(NDEVO,720)(NCODE((I-1)+J),J=1,6)
720	FORMAT(' NUMBER',5(3X,I3,5X),3X,I3)
	WRITE(NDEVO,721)(ALABEL(NCODE((I-1)+J)),J=1,6)
721	FORMAT(/,' LABEL',2X,5(A10,1X),A10)
	WRITE(NDEVO,722)(HA((I-1)+J),J=1,6)
722	FORMAT(/,' HEIGHT',1X,5(F10.3,1X),F10.3)
718	WRITE(NDEVO,724)
724	FORMAT(/)
	IF((M/6)*6.EQ.M)RETURN
	WRITE(NDEVO,720)(NCODE(I),I=(M/6)*6+1,M)
	WRITE(NDEVO,721)(ALABEL(NCODE(I)),I=(M/6)*6+1,M)
	WRITE(NDEVO,722)(HA(I),I=(M/6)*6+1,M)
	WRITE(NDEVO,724)
750	IF(IANS.NE.7)RETURN
	IANS=2
	GO TO 2012
20	FORMAT(80A1)
40	FORMAT(1H1,16A5,//)
50	FORMAT(2I4)
60	FORMAT(1H1)
	END
C******************************************************************************
C---------------IDVO, NDEVO, IDVI, NDEVI ARE INPUT THRU COMMON /SPEC/.
C---------------A, B, STRING, DD, DMAX, DMIN, CLIM, IJ ARE
C--------------- RETURNED. THE REMAINING ARGS. ARE INPUT.  TITLE APPARENTLY
C--------------- NOT USED.  IDLG INPUT THRU COMMON /IOB/
	SUBROUTINE READD(INPUT,KODE,POWER,IANS,N,M,STRING,A,B,X,IJ,
     1	CLIM,IFMT,TITLE,DMIN,DMAX,DD)
C
	DOUBLE PRECISION FILNM,DEVN
	COMMON/INIO/IFTR,IFTW,DEVN(30),FILNM(30),IPP(30),DEST(30)
	COMMON/IOB/LEFBK,IRTBK,IALT,MAXPAG,IPAGE,IPAGCT,IDLG,IRSP,
     1		   IDUMMY,JDUMMY
	COMMON/SPEC/IDVO,IDVI,NDEVO,NDEVI
C
	DIMENSION STRING(1),X(1),NCODE(1),CLIM(16),IFMT(80),
     1	TITLE(16),IJ(1),A(1),B(1),DD(M,M)
C
C*******ADAPTATION FROM "SUBROUTINE SHADE AND RELATED ROUTINE"
C*******BY ROBERT F. LING, UNIVERSITY OF CHICAGO.
C
C*******ROUTINES FOR 'RAW DATA, CORRELATION COMPUTED' AND
C*******'RAW DATA, CHI SQUARE COMPUTED' ADDED AT WMU
C
C*******SUBROUTINE READD SERVES TO READ DATA WHICH MAY APPEAR
C*******IN ANY ONE OF SEVERAL TYPES. AN ELEMENT D(I,J)
C*******REPRESENTS THE RANK OF SOME MEASURE OF THE DEGREE
C*******OF RELATIONSHIP BETWEEN THE ITH AND JTH OBSERVATIONS OF THE
C*******SET,WHERE THE OBSERVATIONS ARE NUMBERED FROM I TO N.
C
C*******UPON ENTRY TO THE SUBROUTINE, ARGUMENT INPUT
C*******CONTAINS A CONTROL CODE WHICH SPECIFIES THE
C*******TYPE OF INPUT DATA WHICH IS TO BE EXPECTED.
C*******IF THE INPUT DATA IS TO BE VARIABLES FOR THE
C*******OBSERVATIONS,VARIABLE M MUST BE READ TO SPECIFY THE
C*******NUMBER OF VARIABLES ASSOCIATED WITH EACH OBSERVATION. AND
C*******THE DISTANCE BETWEEN ANY PAIR OF OBSERVATIONS MUST
C*******BE COMPUTED USING SOME PARTICULAR MEASURE.
C
C*******THE ALLOWABLE CONTROL CODES CORRESPOND TO THOSE
C*******LISTED IN FORMAT 106 IN MAIN PROGRAM.
C
2	FORMAT(80A1)
4	FORMAT(/,' N= ',I6,5X,'M= ',I6)
6	FORMAT(/,1X,'RAW DATA INPUT'//)
7	FORMAT(1X,5E16.7)
8	FORMAT(/,1X,'POWER= ',F6.0)
9	FORMAT(/,' M= ',I6)
10	FORMAT(/,1X,'INPUT DATA - DISTANCES'//)
11	FORMAT(1X)
14	FORMAT(/,' CLASS',5X,'LOWER LIMIT',5X,'UPPER LIMIT',/)
15	FORMAT(1H0,I3,4X,2(3X,G13.6))
16	FORMAT(/,1X,'INPUT DATA - CORRELATIONS'//)
C
C---------------DEFINED FUNCTION
	IJK(I,J)=N*(J-1)+I
C
C*******VARIABLES ARE TO BE READ FOR THE OBSERVATIONS AND
C*******THE DISTANCES BETWEEN PAIRS OF OBSERVATIONS ARE TO
C*******BE COMPUTED. MATRIX X WILL BE USED TO HOLD THE
C*******VARIABLES.
C
C*******N IS THE NUMBER OF OBSERVATIONS,AND M IS THE NUMBER
C*******OF VARIABLES PER OBSERVATION.
C
C*******IFMT IS A VALID FORTRAN IV CODE. THIS CODE SPECIFIES
C*******THE FORMAT FOR THE VARIABLES ASSOCIATED WITH
C*******EACH OBSERVATION OF THE SET.
C
	IF(IANS.NE.2.OR.INPUT.EQ.7.OR.INPUT.EQ.2.OR.INPUT.EQ.3)
     1	GO TO 1111
	IIM=M
	M=N
	N=IIM
1111	IF(IDVO.NE.'TTY')WRITE(NDEVO,2004)IFMT
2004	FORMAT('1FORMAT=',/,5(1X,16A5,/))
	IF(IDVO.NE.'TTY'.AND.(INPUT.EQ.2.OR.INPUT.EQ.3))WRITE(NDEVO,9)M
	IF(IDVO.NE.'TTY'.AND.(INPUT.NE.2.AND.INPUT.NE.3))
     1	WRITE(NDEVO,4)N,M
	IF(IDVO.NE.'TTY')WRITE(NDEVO,6)
	IF((IANS.AND.4).NE.4)GO TO 1070
	IF(IDVI.EQ.'TTY')WRITE(IDLG,106)
	IF(IDVI.NE.'TTY')WRITE(IDLG,107)
106	FORMAT(' ENTER DATA',/)
107	FORMAT(' DATA BEING READ',/)
1070	GO TO(102,200,200,102,102,102,250),INPUT
102	IF((IANS.AND.4).NE.4)GO TO 110
	IF(IDVO.NE.'TTY'.AND.INPUT.EQ.5)WRITE(NDEVO,8)POWER
C
C*******READ THE VARIABLES FOR EACH OBSERVATION.
C
	DO 105 I=1,N
C---------------SEE ST. 16+2 FOR IJK.
	READ(NDEVI,IFMT) (X(IJK(I,K)),K=1,M)
	IF(IDVO.NE.'TTY')WRITE(NDEVO,7) (X(IJK(I,K)),K=1,M)
105	IF(IDVO.NE.'TTY')WRITE(NDEVO,11)
C*******NOTE: TRANSPOSE IS NORMAL MODE, SO '2'
C*******INHIBITS TRANSPOSE.
110	IF(IANS.EQ.6)GO TO 1110
C*******NOTE N AND M ARE TRANSPOSED IN THE CALL SINCE THE 
C*******SUBR NEEDS STORAGE BY COLUMN NOT ROW
	CALL TRANSP(X,N,M,STRING,IOK)
	IF(IOK.NE.0)GO TO 80
C*******NOW TO REALLY TRANSPOSE M AND N
	IIM=M
	M=N
	N=IIM
C
C*******INITIALIZE FOR CHI SQ.
C
1110	IF(INPUT.NE.6)GO TO 108
	L=0
	DO 109 I=2,N
	B(I)=0
	DO 1090 K=1,M
1090	B(I)=B(I)+X(IJK(I,K))
	DO 109 J=1,I-1
	L=L+1
	A(L)=0
	DO 109 K=1,M
109	A(L)=A(L)+X(IJK(I,K))+X(IJK(J,K))
	B(1)=0
	DO 1091 K=1,M
1091	B(1)=B(1)+X(IJK(1,K))
C
C*******CALCULATE THE DISTANCE BETWEEN PAIRS OF OBSERVATIONS.
C
108	L=0
	DO 120 I=2,N
	IM=I-1
	DO 120 J=1,IM
	L=L+1
	STRING(L)=0.
	IF(INPUT.EQ.6)GO TO 117
C
C*******COMPUTE MINKOWSKY METRIC.
C*******COMPUTE EUCLIDEAN DISTANCES(POWER=2)
C*******COMPUTE MANHATTEN METRIC(POWER=1)
C
115	DO 116 K=1,M
116	STRING(L)=STRING(L)+(ABS(X(IJK(I,K))-X(IJK(J,K))))**POWER
	GO TO 120
C
C*******COMPUTE CHI SQUARE DISTANCES
C
117	DO 118 K=1,M
	EI=(X(IJK(I,K))+X(IJK(J,K)))*B(I)/A(L)
	EJ=(X(IJK(I,K))+X(IJK(J,K)))*B(J)/A(L)
	STRING(L)=STRING(L)+(X(IJK(I,K))-EI)**2./EI+
     1	(X(IJK(J,K))-EJ)**2./EJ
118	CONTINUE
120	CONTINUE
	IIM=M
	M=N
	N=IIM
	GO TO 280
C
C*******A LOWER TRIANGULAR MATRIX OF VALUES TO BE READ.
C
C*******N IS THE NUMBER OF OBSERVATIONS IN THE SET.
C

C*******IFMT IS THE FORMAT CODE WHICH SPECIFIES THE FORMAT
C*******FOR THE VALUES ASSOCIATED WITH EACH ROW OF THE
C*******MATRIX D.
C
200	IF(INPUT.EQ.2)GO TO 210
	IF(IDVO.NE.'TTY')WRITE(NDEVO,16)
	GO TO 215
210	IF(IDVO.NE.'TTY')WRITE(NDEVO,10)
215	M1=M-1
	I2=0
	DO 220 I=1,M1
	I1=I2+1
	I2=I2+I
C
C*******READ THE VALUES FOR THE LOWER TRIANGULAR MATRIX,
C*******ONE ROW AT A TIME.
C
	READ(NDEVI,IFMT)(STRING(L),L=I1,I2)
	IF(IDVO.NE.'TTY')WRITE(NDEVO,7)(STRING(L),L=I1,I2)
220	IF(IDVO.NE.'TTY')WRITE(NDEVO,11)
	L=I2
	GO TO (280,280,235),INPUT
235	CONTINUE
	DO 240 I=1,L
240	STRING(I)=1.-ABS(STRING(I))
	GO TO 280
C
C*******ELEMENTS OF ARRAY SIRING NOW CONTAIN DISTANCES OR
C*******MEASURES OF THE DEGREE OF SIMILARITY BETWEEN PAIRS
C*******OF OBSERVATIONS. ELEMENTS STRING(1),...,STRING(L)
C*******REPRESENT THE PAIR OF OBSERVATIONS (2,1),(3,1),(3,2),...,
C*******(M,1),(M,2),...,(M,M-1).
C
C*******INPUT RAW DATA COLUMNS, CREATE CORRELATION TRIANGLE
C*******AND PASS IT TO BE MADE INTO DISTANCE TRIANGLE.
C
250	IF(IANS.NE.6)GO TO 2506
	DO 260 K=1,N
	READ(NDEVI,IFMT,END=265)(X(IJK(K,J)),J=1,M)
	IF(IDVO.NE.'TTY')WRITE(NDEVO,7)(X(IJK(K,J)),J=1,M)
260	CONTINUE
265	CONTINUE
	CALL TRANSP(X,N,M,STRING,IOK)
	IF(IOK.NE.0)GO TO 80
	MM=M
	M=N
	N=MM
	L=0
	DO 262 I=2,M
	DO 261 J=1,I-1
	L=L+1
261	STRING(L)=0
	A(I)=0
262	B(I)=0
	A(1)=0
	B(1)=0
	DO 268 K=1,N
	DO 266 J=1,M
	A(J)=A(J)+X(IJK(K,J))
266	B(J)=B(J)+X(IJK(K,J))*X(IJK(K,J))
	L=0
	DO 264 I=2,M
	DO 264 J=1,I-1
	L=L+1
264	STRING(L)=STRING(L)+X(IJK(K,I))*X(IJK(K,J))
268	CONTINUE
	GO TO 256
2506	L=0
	MM=M
	IF(IANS.EQ.2)MM=N
	DO 248 I=2,MM
	DO 249 J=1,I-1
	L=L+1
249	STRING(L)=0
	A(I)=0
248	B(I)=0
	A(1)=0
	B(1)=0
	IF(IANS.NE.2)GO TO 251
	CALL TRANSP(X,N,M,STRING,IOK)
	IF(IOK.NE.0)GO TO 80
	IIM=M
	M=N
	N=IIM
251	DO 255 K=1,N
	KA=K
	IF(CORFLG.NE.0)KA=1
	IF(IANS.EQ.2)GO TO 2511
	READ(NDEVI,IFMT,END=256)(X(IJK(KA,J)),J=1,M)
	IF(IDVO.NE.'TTY')WRITE(NDEVO,7)(X(IJK(KA,J)),J=1,M)
2511	DO 252 J=1,M
	A(J)=A(J)+X(IJK(KA,J))
252	B(J)=B(J)+X(IJK(KA,J))*X(IJK(KA,J))
	L=0
	DO 254 I=2,M
	DO 254 J=1,I-1
	L=L+1
254	STRING(L)=STRING(L)+X(IJK(KA,I))*X(IJK(KA,J))
255	CONTINUE
256	L=0
	DO 258 I=2,M
	DO 258 J=1,I-1
	L=L+1
258	STRING(L)=(N*STRING(L)-A(I)*A(J))/
     1	SQRT((N*B(I)-A(I)*A(I))*(N*B(J)-A(J)*A(J)))
	LLL=0
	GO TO 235
C
280	GO TO (300,350,400,420),KODE
C
C*******INITIALIZE ARRAY IJ TO CONTAIN INTEGERS FROM 1 TO L.
C
300	DO 310 I=1,L
310	IJ(I)=I
C
C*******SORT THE ELEMENTS OF STRING AND REARRANGE THE
C*******ELEMENTS OF ARRAY IJ IN THE SAME ORDER AS THAT OF STRING.
C
	CALL SORT(STRING,1,L,IJ,0)
C
C*******REPLACE THE VALUE OF EACH ELEMENT IN STRING BY THE
C*******RANK ASSOCIATED WITH THAT VALUE.
C
	DO 320 I=1,L
	IDUM=IJ(I)
	STRING(IDUM)=I
320	CONTINUE
	L=0
C
C*******PLACE RANK OF EACH ELEMENT BELOW THE DIAGONAL INTO D:
C*******AND SINCE THE MATRIX IS SYMMETRIC,STORE THE SAME
C*******RANK INTO THE CORRESPONDING ELEMENT ABOVE THE DIAGONAL.
C*******ALSO, SET THE DIAGONAL ELEMENTS EQUAL TO RANK ZERO.
C
	DO 330 I=2,M
	IM=I-1
	DO 330 J=1,IM
	L=L+1
	DD(I,J)=STRING(L)
	DD(J,I)=DD(I,J)
330	DD(I,I)=0
C
C*******FIND THE MAXIMUM AND MINIMUM DISTANCE BETWEEN
C*******PAIRS OF OBSERVATIONS. THE MAXIMUM DISTANCE WILL
C*******BE REPRESENTED BY DMAX AND THE MINIMUM DISTANCE BY DMIN.
C
C
C*******THE VARIABLE KODE DETERMINES HOW THE INTERVALS
C*******ARE TO BE MADE. IF KODE EQUALS 1,THEN THE RANKS
C*******OF DISTANCES ARE DIVIDED INTO 15 EQUAL INTERVALS,THE
C*******ENDPOINTS OF WHICH  ARE STORED IN THE VECTOR CLIM.
C*******IF KODE EQUALS 2,THE ACTUAL DISTANCES ARE USED.
C*******THE DISTANCES ARE DIVIDED INTO 15 EQUALLY SPACED
C*******INTERVALS BETWEEN MAX D(I,J) AND MIN D(I,J) IN THE D MATRIX.
C*******IF KODE EQUALS 3,THE THE DISTANCES ARE DIVIDED INTO 15
C*******EQUALLY SPACED INTERVALS BETWEEN VALUES DMAX AND DMIN
C*******WHICH ARE SPECIFIED BY THE USER. IF A DISTANCE
C*******D(I,J) < DMIN,IT IS INCLUDED IN THE FIRST INTERVAL.
C*******IF D(I,J) > DMAX, IT IS INCLUDED INTHE LAST INTERVAL.
C*******IF KODE EQUALS 4, THE ENDPOINTS OF THE PLOTTING
C*******INTERVALS ARE TO BE READ INTO THE VECTOR CLIM FROM
C*******DATA CARDS.
C
	DD(1,1)=0
	GO TO 600
350	DMAX=-1.E20
	DMIN=1.E20
	DO 370 I=1,L
	IF(STRING(I).GE.DMIN) GO TO 360
	DMIN=STRING(I)
	GO TO 370
360	IF(STRING(I).LE.DMAX) GO TO 370
	DMAX=STRING(I)
370	CONTINUE
	IF(INPUT.NE.3)GO TO 440
	DSAVE=DMIN
	DMIN=1.-DMAX
	DMAX=1.-DSAVE
	GO TO 440
400	CONTINUE
	GO TO 440
420	GO TO 500
C
C*******THE ENDPOINTS OF INTERVALS ARE CALCULATED AND
C*******STORED IN THE VECTOR CLIM.
C
440	CLIM(1)=DMIN
	CLIM(16)=DMAX
	DELTA=(DMAX-DMIN)/15.
	DO 450 I=2,15
450	CLIM(I)=CLIM(I-1)+DELTA
500	L=0
C
C*******PLACE THE VALUE OF EACH ELEMENT BELOW THE DIAGONAL INTO D:
C*******AND SINCE THE MATRIX IS SYMMETRIC, STORE THE SAME VALUE
C*******INTO THE CORRESPONDING ELEMENT ABOVE THE DIAGONAL.
C*******ALSO,SET THE DIAGONAL ELEMENTS OF D EQUAL TO ZERO.
C
	DO 520 I=2,M
	IM=I-1
	DO 520 J=1,IM
	L=L+1
	DD(I,J)=STRING(L)
	DD(J,I)=DD(I,J)
520	DD(I,I)=0.
	DD(1,1)=0.
550	IF(IDVO.NE.'TTY')WRITE(NDEVO,14)
	DO 560 I=1,15
560	IF(IDVO.NE.'TTY')WRITE(NDEVO,15)I,CLIM(I),CLIM(I+1)
	IF(IDVO.NE.'TTY'.AND.KODE.EQ.3)WRITE(NDEVO,4003)DMAX,DMIN
4003	FORMAT(/,' MAX='F,'   MIN=',F,/)
600	CONTINUE
	RETURN
80	WRITE(IDLG,81)IOK
81	FORMAT(' ?TRANSPOSE ERROR - I='I10,/)
	CALL EXIT
	END
C******************************************************************************
C---------------X IS RETURNED.  THE REMAINING ARGS. ARE INPUT.
C--------------- D IS EQUIVALENCED TO DD.  SEE COMMENTS IN 
C--------------- SUBR. MNL ST. 766+1.  IDVO, NDEVO, INPUT
C--------------- THRU COMMON /SPEC/
	SUBROUTINE SHADE(ALABEL,D,DD,X,NCODE,INPUT,KODE,N,CLIM)
C
C*******SEE COMMENT IN READD FOR SOURCE
C
	INTEGER D
	DOUBLE PRECISION ALABEL(1),FILNM,DEVN
	COMMON/INIO/IFTR,IFTW,DEVN(30),FILNM(30),IPP(30),DEST(30)
	COMMON/IOB/LEFBK,IRTBK,IALT,MAXPAG,IPAGE,IPAGCT,IDLG,IRSP,
     1		   IDUMMY,JDUMMY
	COMMON/SPEC/IDVO,IDVI,NDEVO,NDEVI
C
	DIMENSION C(15,4),IC(15),DD(N,N)
	DIMENSION D(N,N),X(N,4),NCODE(1),CLIM(16)
C
C*******THE SHADE SUBROUTINE IS A PRINT ROUTINE WHICH
C*******CAN BE USED TO DISPLAY THE RESULTS OF A CLUSTER
C*******ANALYSIS.
C
C*******AN ARRAY ELEMENT D(I,J) CONTAINS A REAL NON-NEGATIVE
C*******VALUE WHICH IS THE RANK OF SOME MEASURE OF THE
C*******DEGREE OF RELATIONSHIP BETWEEN THE ITH AND JTH
C*******OBSERVATIONS OF A SET CONTAINING A MAXIMUM OF M OBSERVATIONS.
C*******N IS RESTRICTED TO BE GREATER THAN ZERO BUT LESS
C*******THAN 101. AN ELEMENT OF NCODE CONTAINS AN OBSERVATION
C*******NUMBER,AND THE CONTENTS OF ARRAY ELEMENTS NCODE(1),
C*******NCODE(2),...,NCODE(N) CONSTITUTE AN ORDERING OF THESE
C*******OBSERVATIONS(USUALLY ONE WHICH CORRESPONDS TO THE ORDER
C*******OF OBSERVATIONS IN A DENDROGRAM OBTAINED BY SOME
C*******TECHNIQUE OF CLUSTER ANALYSIS.)
C*******ARRAY C HOLDS THE CHARACTERS WHICH ARE OVERPRINTED
C*******TO FORM 15 SYMBOLS OF VARIOUS SHADES. THESE
C*******SYMBOLS HAVE RANKS FROM 1 TO 15, AND THE SYMBOL WITH
C*******RANK J IS FORMED BY OVERSTRIKING THE CHARACTERS
C*******CONTAINED IN C(J,1),C(J,2),C(J,3),C(J,4).
C
	DATA C(1,1) /1H./,C(2,1) /1H-/ ,C(3,1) /1H-/ ,
     1     C(1,2) /1H / , C(2,2) /1H /, C(3,2) /1H./ ,
     2     C(1,3) /1H / , C(2,3) /1H / , C(3,3) /1H / ,
     3     C(1,4) /1H / , C(2,4) /1H / , C(3,4) /1H /
	DATA C(4,1) /1H=/ , C(5,1) /1H+/ , C(6,1) /1H=/ ,
     1     C(4,2) /1H / , C(5,2) /1H / , C(6,2) /1HI/ ,
     2     C(4,3) /1H / , C(5,3) /1H / , C(6,3) /1H / ,
     3     C(4,4) /1H / , C(5,4) /1H / , C(6,4) /1H /
	DATA C(7,1) /1H*/ , C(8,1) /1HX/ , C(9,1) /1HX/ ,
     1     C(7,2) /1H / , C(8,2) /1H / , C(9,2) /1HI/ ,
     2     C(7,3) /1H / , C(8,3) /1H / , C(9,3) /1H / ,
     3     C(7,4)/1H / , C(8,4) /1H / , C(9,4) /1H /
	DATA C(10,1) /1HX/ , C(11,1) /1HH/ , C(12,1) /1HH/ ,
     1     C(10,2)  /1H-/ , C(11,2) /1HE/ , C(12,2) /1HE/ ,
     2     C(10,3) /1HI/ , C(11,3) /1H./ , C(12,3) /1H=/ ,
     3     C(10,4) /1H / , C(11,4) /1H / , C(12,4) /1H /
	DATA C(13,1) /1HH/ , C(14,1) /1HH/ , C(15,1) /1HH/ ,
     1     C(13,2) /1HE/ , C(14,2) /1HE/ , C(15,2) /1HE/ ,
     2     C(13,3) /1HT/ , C(14,3) /1HX/ , C(15,3)/1HX/ ,
     3     C(13,4) /1H / , C(14,4) /1H / , C(15,4) /1HT/
	DELTA=FLOAT(N*(N-1)/2)/14.999
	DO 140 I=1,N
	IDUM=NCODE(I)
	J=1
10	IF(J.GE.I) GO TO 70
	JDUM=NCODE(J)
	GO TO (20,30,30,30),KODE
20	JRANK=15-INT(DD(IDUM,JDUM)/DELTA)
	GO TO 50
C
C*******THE DIAGONAL ELEMENTS ARE ALWAYS TAKEN TO HAVE
C*******SYMBOL RANK OF 15(DARKEST SHADE) AND NEED NOT
C*******BE DEFINED IN ARRAY D.
C
30	DO 40 K=2,16
	IF(INPUT.NE.3) GO TO 35
	IF(1.-DD(IDUM,JDUM) .GT.CLIM(K)) GO TO 40
	JRANK=K-1
	GO TO 50
35	IF(DD(IDUM,JDUM) .GT.CLIM(K))GO TO 40
	JRANK=17-K
	GO TO 50
40	CONTINUE
	JRANK=1
	IF(INPUT.NE.3) GO TO 50
	JRANK=15
50	DO 60 K=1,4
60	X(J,K)=C(JRANK,K)
	J=J+1
	GO TO 10
70	DO 80 K=1,4
80	X(J,K)=C(15,K)
90	WRITE(NDEVO,160)ALABEL(NCODE(I)),NCODE(I)
	IF(N.LE.50.AND.IDVO.NE.'TTY') GO TO 110
	DO 100 K=1,4
100	WRITE(NDEVO,170) (X(J,K),J=1,I)
	GO TO 140
110	DO 120 K=1,4
120	WRITE (NDEVO,170) (X(J,K),X(J,K),J=1,I)
	WRITE(NDEVO,180)
	DO 130 K=1,4
130	WRITE(NDEVO,170) (X(J,K),X(J,K),J=1,I)
140	CONTINUE
	WRITE(NDEVO,190)
	DO 150 K=1,4
150	WRITE(NDEVO,200) (C(J,K),J=1,15)
	DO 155 I=1,15
	IC(I)=16-I
	IF(INPUT.NE.3)GO TO 155
	IC(I)=I
155	CONTINUE
	WRITE(NDEVO,210) (IC(I),I=1,15)
	RETURN
160	FORMAT(1X,A10,I4)
170	FORMAT(1H+,16X,100A1)
180	FORMAT(1H )
190	FORMAT(8H0SYMBOLS)
200	FORMAT(1H+,10X,15(A1,3X))
210	FORMAT(8H  CLASS ,2X,15(I2,2X))
	END
C******************************************************************************
C---------------ALL ARGS. ARE INPUT. A, MA ARE MODIFIED.
	SUBROUTINE SORT(A,II,JJ,MA,KSW)
C
C*******SEE COMMENT IN READD FOR SOURCE
C
	DIMENSION A(1),MA(1),IU(16),IL(16)
	INTEGER A
C
C*******THIS ROUTINE IS A VARIATION OF ALGORITHM 347 BY
C*******R.C. SINGLETON (CF. CACM 12,3, PP.185-187.)
C*******THIS ROUTINE SORTS ARRAY A INTO INCREASING ORDER, FROM
C*******A(II) TO A(JJ), AND ALSO REARRANGES THE CORRESPONDING
C*******ELEMENTS OF ARRAY MA IF KSW=0. ORDERING IS BY INTEGER
C*******SUBTRACTION, AND THEREFORE FLOATING NUMBERS MUST BE IN
C*******NORMALIZED FORM. ARRAYS IU(K) AND IL(K) PERMIT
C*******SORTING UP TO 2**(K+1)-1 ELEMENTS.
C
	M=1
	I=II
	J=JJ
5	IF(I.GE.J) GO TO 70
10	K=I
	IJ=(J+I)/2
	T=A(IJ)
	IF(KSW.EQ.0) MT=MA(IJ)
	IF(A(I).LE.T) GO TO 20
	A(IJ)=A(I)
	A(I)=T
	T=A(IJ)
	IF(KSW.NE.0) GO TO 20
	MA(IJ)=MA(I)
	MA(I)=MT
	MT=MA(IJ)
20	L=J
	IF(A(J).GE.T)GO TO 40
	A(IJ)=A(J)
	A(J)=T
	T=A(IJ)
	IF(KSW.NE.0)GO TO 25
	MA(IJ)=MA(J)
	MA(J)=MT
	MT=MA(IJ)
25	CONTINUE
	IF(A(I).LE.T)GO TO 40
	A(IJ)=A(I)
	A(I)=T
	T=A(IJ)
	IF(KSW.NE.0)GO TO 28
	MA(IJ)=MA(I)
	MA(I)=MT
	MT=MA(IJ)
28	CONTINUE
	GO TO 40
30	A(L)=A(K)
	A(K)=TT
	IF(KSW.NE.0) GO TO 40
	MA(L)=MA(K)
	MA(K)=MTT
40	L=L-1
	IF(A(L).GT.T) GO TO 40
	TT=A(L)
	IF(KSW.EQ.0)MTT=MA(L)
50	K=K+1
	IF(A(K).LT.T) GO TO 50
	IF(K.LE.L) GO TO 30
	IF(L-I.LE.J-K) GO TO 60
	IL(M)=I
	IU(M)=L
	I=K
	M=M+1
	GO TO 80
60	IL(M)=K
	IU(M)=J
	J=L
	M=M+1
	GO TO 80
70	M=M-1
	IF(M.EQ.0)RETURN
	I=IL(M)
	J=IU(M)
80	IF(J-I.GE.11)GO TO 10
	IF(I.EQ.II)GO TO 5
	I=I-1
90	I=I+1
	IF(I.EQ.J)GO TO 70
	T=A(I+1)
	IF(KSW.EQ.0)MT=MA(I+1)
	IF(A(I).LE.T)GO TO 90
	K=I
100	A(K+1)=A(K)
	IF(KSW.EQ.0)MA(K+1)=MA(K)
	K=K-1
	IF(T.LT.A(K))GO TO 100
	A(K+1)=T
	IF(KSW.EQ.0)MA(K+1)=MT
	GO TO 90
	END
C******************************************************************************
C---------------NA, NB, HA, HB, NCODE, DELHAT ARE RETURNED.
C--------------- STRING, M INPUT
	SUBROUTINE SLINK(NA,NB,HA,HB,NCODE,STRING,M,DELHAT)
C*******ADAPTED BY WMU FROM "SLINK AN OPTIMALLY EFFICIENT ALGORITHM
C*******FOR THE SINGLE LINK CLUSTR METHOD" BY R. SIBSON OF
C*******KING'S COLLEGE RESEARCH CENTRE, CAMBRIDGE
C
	DIMENSION NA(1),NB(1),HA(1),HB(1),NCODE(1),STRING(1)
	TOP=99999999.0
	SIZE=0.0
	NA(1)=1
	HA(1)=TOP	
	K=0
	DO 1 I=2,M
	I1=I-1
	NA(I)=I
	HA(I)=TOP
	DO 13 J=1,I1
	K=K+1
13	HB(J)=STRING(K)
	DO 2 J=1,I1
	IF(HB(J))3,2,2
3	HB(J)=TOP-1.0
2	SIZE=SIZE+HB(J)
1	CALL SLINK1(NA,HA,HB,I1)
	CALL SLINK2(NA,NB,HA,HB,M)
	DELHAT=(SIZE-SLINK3(NA,HA,M))/SIZE
	DO 4 I=1,M
	J=NB(I)
4	NCODE(J)=I
	RETURN
	END
C******************************************************************************
C---------------ALL ARGS. ARE INPUT.  NA, HA, HB, ARE MODIFIED.
	SUBROUTINE SLINK1(NA,HA,HB,I1)
C
C*******SEE ROUTINE SLINK FOR SOURCE
C
	DIMENSION NA(1),HA(1),HB(1)
	DO 1 J=1,I1
	NEXT=NA(J)
	IF(HA(J)-HB(J))2,3,3
2	H=HB(J)
	IF(HB(NEXT)-H)1,1,4
3	H=HA(J)
	NA(J)=I1+1
	HA(J)=HB(J)
	IF(HB(NEXT)-H)1,1,4
4	HB(NEXT)=H
1	CONTINUE
	DO 5 J=1,I1
	NEXT=NA(J)
	IF(HA(J)-HA(NEXT))5,6,6
6	NA(J)=I1+1
5	CONTINUE
	RETURN
	END
C******************************************************************************
C---------------NA, HA, M ARE INPUT.  NB, HB, ARE RETURNED.
	SUBROUTINE SLINK2(NA,NB,HA,HB,M)
C
C*******SEE ROUTINE SLINK FOR SOURCE
C
	DIMENSION NA(1),NB(1),HA(1),HB(1)
	NB(M)=M
	DO 1 N=2,M
	H=HA(M+1-N)
	NEXT=M
2	NOW=NEXT
	NEXT=NB(NOW)
	IF(H-HA(NEXT))3,2,2
3	NB(NOW)=M+1-N
1	NB(M+1-N)=NEXT
	NEXT=NB(M)
4	NOW=NEXT
	N=NA(NOW)
	IF(NB(NOW))5,999,6
6	NEXT=NB(NOW)
	IF(NB(N))7,999,8
8	NB(NOW)=NB(N)
	NB(N)=-NOW
	IF(NEXT-M)4,11,999
7	NOSE=-NB(N)
	NB(NOW)=NB(NOSE)
	NB(N)=-NOW
	NA(NOW)=NOSE
	IF(NEXT-M)4,11,999
5	NOSE=-NB(NOW)
	NEXT=NB(NOSE)
	IF(NB(N))9,999,10
10	NB(NOSE)=NB(N)
	NB(N)=-NOSE
	IF(NEXT-M)4,11,999
9	NA(NOW)=-NB(N)
	NB(N)=-NOSE
	N=NA(NOW)
	NB(NOSE)=NB(N)
	IF(NEXT-M)4,11,999
11	NEXT=-NB(M)
	DO 12 N=1,M
	HB(N)=HA(N)
	NB(N)=NEXT
12	NEXT=NA(NEXT)
	DO 13 N=1,M
	NOW=NB(N)
	NA(N)=NOW
13	HA(N)=HB(NOW)
	DO 14 N=1,M
	NOW=NA(N)
14	NB(NOW)=N
	RETURN
999	STOP 0
	END
C***********************************************************************
C---------------NA, HA, NOBJ ARE INPUT.  NA APPARENTLY NOT USED.
	FUNCTION SLINK3(NA,HA,NOBJ)
	DIMENSION NA(1),HA(1)
	NOBJ1=NOBJ-1
	SLINK3=0
	DO 1 I=1,NOBJ1
	DO 2 J=1,I
	N1=I+1-J
	IF(HA(N1)-HA(I))2,2,3
2	CONTINUE
	N1=0
3	I1=I+1
	DO 4 J=I1,NOBJ
	IF(HA(J)-HA(I))4,1,1
4	CONTINUE
1	SLINK3=SLINK3+(I-N1)*(J-I)*HA(I)
	RETURN
	END
C******************************************************************************
C---------------K, X ARE RETURNED.  THE OTHER ARGS. ARE INPUT.
	SUBROUTINE DENDRG(H,N,K,X,IO,IDVO,ITYPE)
C
C*******PROGRAM TO GENERATE EITHER A CLUSTER TREE OR A DENDROGRAM FROM A
C*******CLUSTERED LIST OF MAGNITUDES OF LENGTH N. H(N) IS TAKEN AS THE
C*******INFINITY HEIGHT.
C
C*******WRITTEN BY RUSS BARR III - WMU - KALAMAZOO.
C*******DATE: 15 OCT 73
C*******DATE:  8 JUL 76 - MOD TO PRINT ANY SIZE DENDROGRAM
C
C*******ARGUMENTS:
C
C*******H	- HEIGHT VECTOR(LENGTH N)
C*******N	- LENGTH OF ARRAYS H AND K
C*******K	- SCRATCH ARRAY FOR POSITIONING(LENGTH N)
C*******X	- SCRATCH ARRAY FOR PRINTING(LENGTH 2*N-1)
C*******IO	- OUTPUT CHANNEL
C*******IDVO	- OUTPUT DEVICE
C*******ITYPE	= 1 FOR BASIC TREE
C*******= 2 FOR DENROGRAM
C
	DIMENSION H(1),K(1),X(1)
C*******INITIALIZE
	FORM=' '
	IF(IDVO.NE.'TTY')FORM='*'
	BLANK=' '
	AYE='I'
	PLUS='+'
	DASH='-'
	HL=-1
	NA=2*N-1
	BIG=H(N)
	IF=(N-1)/60
	DO 100 I=1,N-1
	DO 98 J=I+1,N
98	IF(H(I).EQ.H(J))H(J)=H(J)+1.E-6
C*******FIND BIGGEST IN H(1) THRU H(N-1)
	IF(H(I).GT.BIG)BIG=H(I)
C*******SET UP K ARRAY
100	K(I)=I*2-1
	K(N)=2*N-1
C*******STORE TRUE VALUE H(N)
	TEMP=H(N)
C*******MAKE H(N) BIGGER THAN THE BIGGEST
	H(N)=BIG+1
	DO 118 I=1,N-1
C*******CLEAR THE PRINT VECTOR
	DO 102 J=1,NA
102	X(J)=BLANK
C*******EXTEND BRANCHES IF NOT YET JOINED
	DO 104 J=1,N
104	IF(H(J).GT.HL)X(K(J))=AYE
C*******WRITE INTERMEDIATE LINE
	NAA=NA
	IF(IF.NE.0)GO TO 107
	DO 105 NAA=NA,1,-1
105	IF(X(NAA).NE.BLANK)GO TO 107
107	IF(NAA.LE.120)GO TO 2107
	DO 3107 J1=1,IF
	IU=J1+19
3107	WRITE(IU,3106)(X(J),J=J1*120+1,MIN0((J1+1)*120,NAA))
3106	FORMAT(121A1)
2107	WRITE(IO,106)FORM,BLANK,(X(J),J=1,MIN0(NAA,120))
106	FORMAT(A1,10X,121A1)
	HN=H(N)
	L=1
C*******FIND NEXT HIGHEST JUNCTION
	DO 108 J=1,N-1
	IF(H(J).LE.HL.OR.H(J).GT.HN)GO TO 108
	HN=H(J)
	L=J
108	CONTINUE
C*******FIND MATCHING BRANCH
	DO 110 J=L+1,N
110	IF(H(J).GT.H(L))GO TO 112
C*******JOIN THEM
112	DO 114 M=K(L),K(J)
114	X(M)=DASH
C*******NEW BRANCH FROM MIDDLE IF DENDOGRAM
	IF(ITYPE.EQ.2)K(J)=(K(J)+K(L))/2
C*******BRANCH SPRINGS FROM HERE EITHER CASE
	X(K(J))=PLUS
	LASTK=K(J)
C*******WRITE JUNCTION LINE
	NAA=NA
	IF(IF.NE.0)GO TO 117
	DO 115 NAA=NA,1,-1
115	IF(X(NAA).NE.BLANK)GO TO 117
117	IF(IF.EQ.0)GO TO 2117
	DO 3117 J1=1,IF
	IU=J1+19
3117	WRITE(IU,3106)(X(J),J=J1*120+1,MIN0((J1+1)*120,NAA))
2117	WRITE(IO,116)FORM,H(L),BLANK,(X(J),J=1,MIN0(NAA,120))
116	FORMAT(A1,F10.3,121A1)
C*******SAVE HEIGHT OF JUNCTION
	HL=H(L)
118	CONTINUE
	IF(IF.EQ.0)GO TO 123
	DO 3121 J1=1,IF
	IU=J1+19
	REWIND(IU)
	WRITE(IO,3199)
3199	FORMAT('1',/)
3120	READ(IU,3106,END=3121)(X(J),J=J1*120+1,MIN0((J1+1)*120,NAA))
	WRITE(IO,3106)FORM,(X(J),J=J1*120+1,MIN0((J1+1)*120,NAA))
3123	FORMAT(A1,1X,120A1)
	GO TO 3120
3121	CLOSE(UNIT=IU,DISPOSE='DELETE')
3122	CONTINUE
123	WRITE(IO,124)
124	FORMAT(//)
C*******RESTORE H(N)
	H(N)=TEMP
	RETURN
	END
C******************************************************************************
C---------------X IS RETURNED, THE OTHER ARGS. ARE INPUT.
	SUBROUTINE DENDHD(H,X,N,IO,ALABEL,LANS)
C
C*******PROGRAM TO GENERATE HEADER FOR DENDROGRAM GENERATED BY
C*******PROGRAM - DENDRG
C
C*******WRITTEN BY RUSS BARR III - WMU - KALAMAZOO
C*******DATE: 15 OCT 73
C
C*******ARGUMENTS:
C
C*******H	- ARRAY OF 3 DIGIT NUMBERS FROM CLUSTR SUBR.
C*******X	- SCRATCH ARRAY
C*******N	- LENGTH OF ARRAYS H AND X
C*******IO	- FORTRAN OUTPUT DEVICE NUMBER
C*******ALABEL	- CONTAINS LABELS FOR DENDRG
C*******LANS	= 'YES' IF LABELS REQUESTED
C*******# 'YES' IF NOT REQUESTED
C
	INTEGER H,X,DUMB
	DIMENSION H(1),X(1),ND(10)
	DOUBLE PRECISION ALABEL(1),BLABEL(2,2),BLANK
	DATA ((BLABEL(I,J),J=1,2),I=1,2)/
     1	'CASE      ','NUMBER    ','CASE      ','LABEL     '/
	DATA BLANK/'          '/
C
	IF=(N-1)/60
	IF(IF.EQ.0)GO TO 102
	DO 101 J1=1,IF
	IU=J1+19
101	OPEN(UNIT=IU,DEVICE='DSK')
102	WRITE(IO,10)
10	FORMAT('1DENDROGRAM',/)
	DO 5 J=1,3
	DO 1 I=1,N
	ENCODE(3,2,DUMB)H(I)
2	FORMAT(I3)
	DECODE(J,21,DUMB)(ND(K),K=1,J)
1	X(I)=ND(J)
	IF(IF.EQ.0)GO TO 108
	DO 115 J1=1,IF
	IU=J1+19
115	WRITE(IU,114)(X(J2),J2=J1*60+1,MIN0((J1+1)*60,N))
114	FORMAT(60(A1,1X))
108	IF(J.LE.2)GO TO 6
	WRITE(IO,4)BLANK,(X(I),I=1,MIN0(N,60))
	GO TO 5
6	WRITE(IO,4)BLABEL(1,J),(X(I),I=1,MIN0(N,60))
5	CONTINUE
4	FORMAT(1X,A10,1X,60(A1,1X))
	WRITE(IO,8)
8	FORMAT(1X)
	IF(IF.EQ.0)GO TO 119
	DO 118 J2=1,IF
	IU=J2+19
118	WRITE(IU,8)
119	IF(LANS.NE.'YES')RETURN
	DO 22 J=1,10
	DO 20 I=1,N
	DECODE(J,21,ALABEL(H(I)))(ND(K),K=1,J)
20	X(I)=ND(J)
21	FORMAT(10A1)
	IF(IF.EQ.0)GO TO 123
	DO 122 J1=1,IF
	IU=J1+19
122	WRITE(IU,114)(X(J2),J2=J1*60+1,MIN0((J1+1)*60,N))
123	IF(J.LE.2)GO TO 23
	WRITE(IO,4)BLANK,(X(I),I=1,MIN0(N,60))
	GO TO 22
23	WRITE(IO,4)BLABEL(2,J),(X(I),I=1,MIN0(N,60))
22	CONTINUE
	WRITE(IO,8)
	RETURN
	END
C******************************************************************************
C---------------A, M, N, MOVE INPUT;  IOK RETURNED.  A, MOVE MODIFIED.
C
C	REPRINTING PRIVILEGES WERE GRANTED BY
C	PERMISSION OF THE ASSOCIATION FOR COM-
C	PUTING MACHINERY, BUT NOT FOR PROFIT.
C
	SUBROUTINE TRANSP(A,M,N,MOVE,IOK)
C*******IN-SITU TRANSPOSITION OF A RECT. MATRIX[F1]
C*******ALGORITHM 380 - CACM MAY 70
C
C*******ADAPTED BY RUSS BARR AT WMU OCT 1973
C*******SIMPLIFIED TO PASS FEWER ARGS(.NOT.(MN.AND.IWRK))
C
C*******ARGUMENTS:
C
C*******A	- 1 DIMENSIONAL ARRAY, LENGTH MN=M*N(MATRIX
C*******	TO BE TRANSPOSED).
C*******M	- NUMBER OF COLUMNS(LENGTH OF ROWS)
C*******N	- NUMBER OF ROWS(LENGTH OF COLUMNS)
C*******MOVE	- 1 DIMENSIONAL ARRAY SCRATCH AREA, LENGTH IWRK=(M+N)/2
C*******IOK	- ERROR CODE	=  0 NORMAL RETURN
C*******			= >0 = FINAL VALUE WHEN SEARCH COMPLETED
C*******			(SOME LOOPS HAVE NOT BEEN MOVED)
	DIMENSION A(1),MOVE(1)
C*******CHECK ARGS AND INIT.
	IF(M.LT.2.OR.N.LT.2)GO TO 60
	MN=M*N
	IWRK=(M+N)/2
	IF(M.EQ.N)GO TO 70
	NCOUNT=2
	M2=M-2
	DO 10 I=1,IWRK
10	MOVE(I)=0
	IF(M2.LT.1)GO TO 12
C*******COUNT NUMBER,NCOUNT,OF SINGLE POINTS.
	DO 11 IA=1,M2
	IB=IA*(N-1)/(M-1)
	IF(IA*(N-1).NE.IB*(M-1))GO TO 11
	NCOUNT=NCOUNT+1
	I=IA*N+IB
	IF(I.GT.IWRK)GO TO 11
	MOVE(I)=1
11	CONTINUE
C*******SET INITIAL VALUES FOR SEARCH
12	K=MN-1
	KMI=K-1
	MAX=MN
	I=1
C*******AT LEAST ONE LOOP MUST BE REARRANGED
	GO TO 30
C*******SEARCH FOR LOOPS TO REARRANGE
20	MAX=K-I
	I=I+1
	KMI=K-I
	IF(I.GT.MAX)GO TO 90
	IF(I.GT.IWRK)GO TO 21
	IF(MOVE(I).LT.1)GO TO 30
	GO TO 20
C
C*******ORIGINAL CODE FOLLOWS
C
C********21	IF(I.EQ.M*I-K*(I/N))GO TO 20
C********	I1=I
C********22	I2=M*I1-K*(I1/N)
C********	IF(I2.LE.I.OR.I2.GE.MAX)GO TO 23
C********	I1=I2
C********	GO TO 22
C
C*******PATCH FROM JAN 72 CACM PAGE 49 FOLLOWS
C
21	I2=M*I-K*(I/N)
	IF(I2.LE.I.OR.I2.GE.MAX)GO TO 20
22	I2=M*I2-K*(I2/N)
	IF(I2.GT.I.AND.I2.LT.MAX)GO TO 22
C*******END OF PATCH
123	IF(I2.NE.I)GO TO 20
C*******REARRANGE ELEMENTS OF A LOOP
30	I1=I
31	B=A(I1+1)
32	I2=M*I1-K*(I1/N)
	IF(I1.LE.IWRK)MOVE(I1)=2
33	NCOUNT=NCOUNT+1
	IF(I2.EQ.I.OR.I2.GE.KMI)GO TO 35
34	A(I1+1)=A(I2+1)
	I1=I2
	GO TO 32
35	IF(MAX.EQ.KMI.OR.I2.EQ.I)GO TO 41
	MAX=KMI
	GO TO 34
C*******TEST FOR SYMMETRIC PAIR OF LOOPS
41	A(I1+1)=B
	IF(NCOUNT.GE.MN)GO TO 60
	IF(I2.EQ.MAX.OR.MAX.EQ.KMI)GO TO 20
	MAX=KMI
	I1=MAX
	GO TO 31
C*******NORMAL RETURN
60	IOK=0
	RETURN
C*******IF MATRIX IS SQUARE EXCHANGE A(I,J) AND A(J,I)
70	N1=N-1
	DO 71 I=1,N1
	J1=I+1
	DO 71 J=J1,N
	I1=I+(J-1)*N
	I2=J+(I-1)*M
	B=A(I1)
	A(I1)=A(I2)
	A(I2)=B
71	CONTINUE
	GO TO 60
C*******ERROR RETURNS
90	IOK=I
91	RETURN
	END

