Trailing-Edge
-
PDP-10 Archives
-
decuslib20-05
-
decus/20-0137/clustr/clustr.for
There are 2 other files named clustr.for in the archive. Click here to see a list.
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