C ALGORITHM 479 COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN COMM. ACM, VOL. 17, NO. 06, C P. 321. C TO CLUSTER A POINT SET USING THIS ALGORITHM, TWO THINGS A 10 C NEED TO BE DONE. (1) BUILD THE MINIMAL SPANNING TREE BY A 20 C CALLING GROW, AND (2) DELETE ITS INCONSISTENT BRANCHES BY A 30 C CALLING CLUSTR. ONCE STEP (1) HAS BEEN DONE, STEP (2) CAN A 40 C BE REPEATED OVER AND OVER WITH DIFFERENT PARAMETERS. A 50 C SEE THE BEGINNINGS OF GROW AND CLUSTR FOR EXPLANATIONS OF A 60 C THE PARAMETERS. A 70 C CURRENTLY, THE ARRAYS ARE DIMENSIONED TO HANDLE UP TO 100 A 80 C POINTS. TO CHANGE THIS, SIMPLY CHANGE THE SIZE OF THE A 90 C ARRAYS MST, NIT, AND UI IN GROW AS DIRECTED BELOW THEIR A 100 C DECLARATIONS. ALSO, CHANGE THE LENGTHS OF A 110 C THE ARRAYS EDGE ST, EDGE PT, AVE, SQ, AND NUMNEI AS A 120 C DIRECTED IN THE SUBROUTINE CLUSTR. IN ADDITION, IF THE A 130 C PARAMETER D IN CLUSTR WILL BE LARGER THAN 18, CHANGE THE A 140 C LENGTHS OF THE ARRAYS NEIG ST AND NEIG PT AS DIRECTED. A 150 SUBROUTINE GROW(DATA, DIMEN, NUMPTS, PRINT) A 160 INTEGER DIMEN, NUMPTS DIMENSION DATA(1) LOGICAL PRINT C THIS SUBROUTINE COMPUTES THE MINIMAL SPANNING TREE OF THE C COMPLETE GRAPH ON THE NUM PTS POINTS IN ARRAY DATA . C EACH POINT IS A VECTOR WITH DIMEN COMPONENTS STORED IN C CONTIGUOUS LOCATIONS IN THE ARRAY DATA. SPECIFICALLY, C DATA( (K-1)*DIMEN +I ) IS THE I-TH COMPONENT OF THE K-TH C VECTOR. THE ARRAY DATA MAY CONTAIN NUMBERS IN EITHER C INTEGER OR FLOATING POINT FORMAT AS LONG AS THE FORMAT IS C CONSISTENT WITH THE TYPE SPECIFICATION OF THE PARAMETERS C IN THE FUNCTION DIST. C IF THE PARAMETER PRINT HAS THE VALUE .TRUE., THEN A C A DESCRIPTION OF THE MINIMAL SPANNING TREE IS PRINTED ON C UNIT 6. EACH NODE IS LABELED WITH AN INTEGER INDICATING C ITS RELATIVE POSITION IN THE ARRAY DATA. INTEGER DIM, N, MST(800), LOC(1), NBR(1), NXT(1) REAL WT(1) EQUIVALENCE (MST,LOC,NBR,WT,NXT) COMMON DIM, N, MST INTEGER LASTPT, FREE, PT C MST (ALIAS LOC, NBR, WT, NXT) IS A DESCRIPTION OF THE C MINIMAL SPANNING TREE. IT CONTAINS ONE LIST FOR EACH NODE. C THE POINTERS TO THE HEADS OF THESE LISTS ARE STORED IN THE C FIRST N=NUM PTS LOCATIONS OF MST AND GO BY THE NAME MST. C THE FIRST ELEMENT OF EACH LIST CONSISTS OF FOUR FIELDS C STORED IN CONTIGUOUS WORDS OF MST. EACH FIELD IS CALLED BY C A NAME WHICH IS AN ALIAS OF MST. C FIELD 1: LOCATION IN DATA OF THE NODE (LOC) C FIELD 2: NAME OF NEIGHBORING NODE (NBR) C FIELD 3: WEIGHT OF THIS BRANCH (WT) C FIELD 4: POINTER TO NEXT NEIGHBOR OR END MARK=0 (NXT) C EACH ADDITIONAL ELEMENT OF THE LIST CONSISTS OF THREE C FIELDS. FIELD 1 ABOVE IS OMITTED. C THE LENGTH OF THE ARRAY MST MUST BE AT LEAST 8*N . C THE MINIMAL SPANNING TREE IS COMPUTED USING THE ALGORITHM C OF PRIM AND DIJKSTRA AS IMPLEMENTED BY WHITNEY (CACM 15, C APR 1972). C EACH COLUMN OF NIT IS A PAIR (NIT(1,I),NIT(2,I),I=1,NITP) C DENOTING A NODE NOT (YET) IN THE TREE AND ITS NEAREST C NEIGHBOR IN THE CURRENT TREE. UI(I) IS THE LENGTH OF THE C EDGE (NIT(1,I),NIT(2,I)). THE LENGTH OF THE ARRAY UI AND C THE NUMBER OF COLUMNS OF NIT CANNOT BE LESS THAN N. INTEGER NIT(2,100) REAL UI(100) DIM = DIMEN N = NUMPTS C COMPUTE MINIMAL SPANNING TREE USING ALGORITHM OF WHITNEY C INITIALIZE NODE LABEL ARRAYS AND SET UP LIST FOR NODE N=KP NITP = N - 1 KP = N KPDATA = (KP-1)*DIM + 1 DO 10 I=1,NITP IDATA = (I-1)*DIM + 1 NIT(1,I) = I UI(I) = DIST(DATA(IDATA),DATA(KPDATA),DIM) NIT(2,I) = KP 10 CONTINUE FREE = N + 1 MST(KP) = FREE LOC(FREE) = (KP-1)*DIM + 1 FREE = FREE + 1 NXT(FREE+2) = 0 C UPDATE LABEL OF NODES NOT YET IN TREE 20 KPDATA = (KP-1)*DIM + 1 DO 30 I=1,NITP IDATA = (NIT(1,I)-1)*DIM + 1 D = DIST(DATA(IDATA),DATA(KPDATA),DIM) IF (UI(I).LE.D) GO TO 30 UI(I) = D NIT(2,I) = KP 30 CONTINUE C FIND NODE OUTSIDE TREE NEAREST TO TREE UK = UI(1) DO 40 I=1,NITP IF (UI(I).GT.UK) GO TO 40 UK = UI(I) K = I 40 CONTINUE C ADD NEW EDGE TO MST C ADD NEIGHBOR TO LIST OF NODE NIT(2,K) C CHANGE END OF LIST MARK TO POINT TO NEXT NEIGHBOR PT = LASTPT(NIT(2,K)) NXT(PT) = FREE C ENTER NAME OF NEIGHBOR NBR(FREE) = NIT(1,K) C ENTER WEIGHT OF THIS BRANCH (OFFSET PICKS UP WT FIELD) WT(FREE+1) = UI(K) C PUT IN END OF LIST MARK (OFFSET PICKS UP POINTER FIELD) NXT(FREE+2) = 0 FREE = FREE + 3 C NEW NODE--CREATE ITS NEIGHBOR LIST C SET UP HEAD POINTER NODE = NIT(1,K) MST(NODE) = FREE C ENTER LOCATION OF THIS NODE IN DATA LOC(FREE) = (NODE-1)*DIM + 1 C ENTER NAME OF NEIGHBORING NODE (OFFSET PICKS UP NBR FIELD) NBR(FREE+1) = NIT(2,K) C ENTER WEIGHT OF THIS BRANCH (OFFSET PICKS UP WT FIELD) WT(FREE+2) = UI(K) C ENTER END OF LIST MARK (OFFSET PICKS UP POINTER FIELD) NXT(FREE+3) = 0 FREE = FREE + 4 KP = NIT(1,K) C DELETE NEW TREE NODE FROM ARRAY NIT UI(K) = UI(NITP) NIT(1,K) = NIT(1,NITP) NIT(2,K) = NIT(2,NITP) NITP = NITP - 1 C THE MST IS FINISHED WHEN IT CONTAINS ALL NODES IF (NITP.NE.0) GO TO 20 IF (PRINT) CALL PRTREE RETURN END SUBROUTINE CLUSTR(D, FACTOR, SPREAD, C, CLEN, PRINT, MINPTS) INTEGER D, CLEN, C(CLEN) REAL FACTOR, SPREAD LOGICAL PRINT C THIS SUBROUTINE FINDS THE CLUSTERS OF A POINT SET USING C A MINIMAL SPANNING TREE CLUSTERING METHOD OF ZAHN. THE C MINIMAL SPANNING TREE, COMPUTED BY SUBROUTINE GROW, IS C STORED IN BLANK COMMON. C THE ZAHN ALGORITHM FINDS CLUSTERS BY DELETING INCONSISTENT C EDGES FROM THE MINIMAL SPANNING TREE, AN INCONSISTENT EDGE C BEING ONE WHOSE WEIGHT IS SIGNIFICANTLY LARGER THAN THE C AVERAGE WEIGHT OF NEARBY EDGES. C NEARBY MEANS CONNECTED TO THE EDGE IN QUESTION BY A C PATH CONTAINING D OR FEWER EDGES. C SIGNIFICANTLY LARGER MEANS C WEIGHT .GT. FACTOR * AVERAGE C AND WEIGHT .GT. AVERAGE + SPREAD * STANDARD DEVIATION C WHERE THE AVERAGE AND STANDARD DEVIATION ARE COMPUTED ON C THE WEIGHTS OF NEARBY EDGES. C THE OUTPUT VECTOR C DESCRIBES THE CLUSTERS DETERMINED. C IT IS ARRANGED IN BLOCKS, EACH BLOCK DESCRIBING ONE C CLUSTER. THE FIRST ELEMENT IN EACH BLOCK IS THE NUMBER C OF NODES IN THE CLUSTER. THE REMAINING ELEMENTS ARE THE C LABELS OF THE NODES IN THE CLUSTER, THE LABEL INDICATING C THE RELATIVE POSITION OF THE NODE IN THE ARRAY DATA. THE C FIRST BLOCK STARTS AT C(2). C C(1) IS THE NUMBER OF CLUSTERS FOUND BY THE ALGORITHM. C THE VALUE OF C LEN SHOULD BE THE TRUE SIZE OF C THE ARRAY C. IT IS USED TO PREVENT INVALID SUBSCRIPTS. C IF C LEN IS ZERO, THE ARRAY C WILL NOT BE USED. C IF THE PARAMETER PRINT HAS THE VALUE .TRUE., CLUSTERS C ARE PRINTED OUT ON UNIT 6. INTEGER EDGEST(101), EDGELN, EDGEPT(101) REAL AVE(100), SQ(100), SUPPWT, W INTEGER NUMNEI(100) INTEGER NEIGST(20), NEIGLN, NEIGPT(20) C THE ARRAY EDGE ST (EDGE STACK) IS A STACK OF NODES USED TO C DIRECT THE SEARCH THROUGH THE TREE FOR INCONSISTENT EDGES. C ITS LENGTH (EDGE LN) CAN GROW AS LARGE AS ONE MORE THAN C THE NUMBER OF NODES IN THE TREE. C THE ARRAY EDGE PT (EDGE POINTERS) IS A STACK OF POINTERS C TO THE NEXT UNEXAMINED NEIGHBORING NODE OF THE NODE IN THE C SAME POSITION IN EDGE ST. THUS THE LENGTH OF EDGE PT IS C ALWAYS THE SAME AS THAT OF EDGE ST. C THE ARRAY NEIG ST (NEIGHBOR STACK) IS A STACK OF NODES C USED TO DIRECT THE AVERAGING OF THE WEIGHTS OF NEARBY C EDGES. ITS LENGTH (NEIG LN) CAN GROW AS LARGE AS D+2. C THE ARRAY NEIG PT IS USED IN CONJUNCTION WITH NEIG ST. ITS C LENGTH CAN GROW AS LARGE A D+2. C THE ARRAYS AVE AND SQ ARE USED TO EXPEDITE THE CALCULATION C OF AVERAGE WEIGHTS. SPECIFICALLY, AVE(I) STORES THE SUM OF C THE WEIGHTS OF EDGES EXTENDING FROM THE I-TH NODE AND C SQ(I) STORES THE SUM OF THE SQUARES. SIMILARLY, NUMNEI(I) C STORES THE NUMBER OF NEIGHBORS OF THE I-TH NODE. THUS EACH C OF THESE ARRAYS MUST BE AS LONG AS THE NUMBER OF NODES. INTEGER FINDCN, A, B, DLESS1 INTEGER CLS, INCLS(1), PARENT(1), BAKWRD, BEGCLS EQUIVALENCE (INCLS,EDGEST), (PARENT,EDGEPT) INTEGER CP, OTHEND INTEGER DIM, N, MST(1), LOC(1), NBR(1), NXT(1) REAL WT(1) EQUIVALENCE (MST,LOC,NBR,WT,NXT) COMMON DIM, N, MST IF(MINPTS.LE.N) GOTO 5 C(1) = 0 RETURN 5 IF(PRINT)WRITE(6,99998)D,FACTOR,SPREAD IF(PRINT)WRITE(6,99996)MINPTS 99996 FORMAT(1H ,10X,39HMINIMUM NUMBER OF POINTS PER CLUSTER +IS, I9) DLESS1 = D - 1 C COMPUTATION SECTION C SUM BRANCH WEIGHTS OFF EACH NODE (DEPTH 1) DO 20 NODE=1,N NUMNEI(NODE) = 1 K = MST(NODE) AVE(NODE) = WT(K+2) SQ(NODE) = WT(K+2)**2 K = NXT(K+3) 10 IF (K.EQ.0) GO TO 20 AVE(NODE) = AVE(NODE) + WT(K+1) SQ(NODE) = SQ(NODE) + WT(K+1)**2 NUMNEI(NODE) = NUMNEI(NODE) + 1 K = NXT(K+2) GO TO 10 20 CONTINUE C INITIALIZE EDGE STACK WITH NODE 1 SURROUNDED BY ITS FIRST C TWO NEIGHBORS. SINCE THE TOP TWO ELEMENTS OF THE STACK C INDICATE THE DIRECTION OF TRAVEL ALONG A BRANCH, THE C SEARCH WILL FIRST BE DIRECTED AWAY FROM NODE 1 IN THE C DIRECTION OF ITS FIRST NEIGHBOR. WHEN ALL THE TREE IN THAT C DIRECTION IS SEARCHED, THE SEARCH WILL PROCEDE AWAY FROM C ITS FIRST NEIGHBOR TOWARD NODE 1. C THE EDGE PT STACK IS USED TO KEEP TRACK OF THE NEIGHBORS C OF THE CORRESPONDING NODE IN EDGE ST WHICH HAVE ALREADY C BEEN SEARCHED. EDGE PT(I) POINTS TO THE LOCATION OF C EDGE ST(I+1) IN THE LIST OF NEIGHBORS OF EDGE ST(I) . EDGELN = 3 K = MST(1) EDGEST(2) = LOC(K)/DIM + 1 EDGEST(1) = NBR(K+1) EDGEST(3) = NBR(K+1) EDGEPT(1) = FINDCN(EDGEST(1),EDGEST(2)) EDGEPT(2) = K + 1 EDGEPT(3) = -1 C CLIMB TREE TO NEXT UNTESTED BRANCH 30 CALL CLIMB(EDGEPT, EDGEST, EDGELN, N) IF (EDGELN.LE.2) GO TO 70 C CHECK THE EDGE BETWEEN NODE EDGE ST(EDGE LN -1) AND C NODE EDGE ST(EDGE LN) FOR INCONSISTENCY. A = EDGEST(EDGELN-1) B = EDGEST(EDGELN) C SUM WEIGHTS OF ALL BRANCHES NEARBY BRANCH A--B NEARBY = 0 AV = 0. STDDEV = 0. C INITIALZE NEIG ST TO SUM WEIGHTS HEADING OFF NODE B NEIGLN = 2 NEIGST(1) = A NEIGPT(1) = EDGEPT(EDGELN-1) NEIGST(2) = B NEIGPT(2) = -1 ASSIGN 50 TO OTHEND C GO OUT TO DEPTH D-1 ALONG BRANCHES NOT YET ADDED 40 CALL CLIMB(NEIGPT, NEIGST, NEIGLN, DLESS1) C ADD WEIGHTS OF BRANCHES OFF THE TOP NODE LESS THE WEIGHT C OF THE BRANCH SUPPORTING IT K = NEIGPT(NEIGLN-1) SUPPWT = WT(K+1) K = NEIGST(NEIGLN) AV = AV + AVE(K) - SUPPWT STDDEV = STDDEV + SQ(K) - SUPPWT**2 NEARBY = NEARBY + NUMNEI(K) - 1 C WHEN DEPTH OF STACK RETURNS TO 2, ALL BRANCH WEIGHTS OFF C THIS END HAVE BEEN ADDED IF (NEIGLN.LE.2) GO TO OTHEND, (50,60) NEIGLN = NEIGLN - 1 GO TO 40 C INITIALZE NEIG ST TO SUM WEIGHTS HEADING OFF NODE A 50 NEIGLN = 2 NEIGST(1) = B NEIGPT(1) = FINDCN(B,A) NEIGST(2) = A NEIGPT(2) = -1 ASSIGN 60 TO OTHEND GO TO 40 C TEST BRANCH A--B FOR INCONSISTENCY. 60 AV = AV/FLOAT(NEARBY) STDDEV = SQRT(ABS(STDDEV/FLOAT(NEARBY)-AV**2)) K = EDGEPT(EDGELN-1) W = WT(K+1) EDGELN = EDGELN - 1 IF (W.LE.AV+SPREAD*STDDEV .OR. W.LE.FACTOR*AV) GO TO 30 C BRANCH A--B IS INCONSISTENT. DELETE IT. NBR(K) = -IABS(NBR(K)) K = NEIGPT(1) NBR(K) = -IABS(NBR(K)) GO TO 30 C OUTPUT SECTION C WE COLLECT THE CLUSTERS AS FOLLOWS: 1. START WITH FIRST C NODE. 2. THROW IN ITS NEIGHBORS. 3. THROW IN NEIGHBORS C OF NEIGHBORS UNTIL NO NEW ONES CAN BE FOUND. 4. EACH C TIME A DELETED BRANCH IS ENCOUNTERED, PUT OTHER END IN A C LIST OF UNUSED NODES (AT TOP OF ARRAY IN CLS). 5. WHEN C A FULL CLUSTER IS COLLECTED , OUTPUT IT. 6. START AGAIN C AT STEP 2 WITH A NODE FROM THE LIST OF UNUSED NODES. 70 NUMIN = 0 CLS = 0 CP = 1 K = MST(1) NXTCLS = N INCLS(NXTCLS) = LOC(K)/DIM + 1 PARENT(NXTCLS) = 0 BAKWRD = 0 C START CLUSTER WITH NEXT AVAILABLE UNUSED NODE 80 CLS = CLS + 1 NUMIN = NUMIN + 1 BEGCLS = NUMIN NXTCN = NUMIN NODE = INCLS(NXTCLS) INLIST = PARENT(NXTCLS) INCLS(NUMIN) = NODE NXTCLS = NXTCLS + 1 C LET K POINT TO FIRST NEIGHBOR OF NODE 90 K = MST(NODE) + 1 C ADD NEIGHBOR TO CLUSTER AND RECORD IT ANCESTRY 100 NXTNBR = NBR(K) IF (NXTNBR.LT.0) GO TO 110 IF (NXTNBR.EQ.BAKWRD) GO TO 120 NUMIN = NUMIN + 1 INCLS(NUMIN) = NXTNBR PARENT(NUMIN) = NODE GO TO 120 C THIS NEIGHBOR IS IN A DIFFERENT CLUSTER--ADD TO UNUSED 110 NXTNBR = -NXTNBR IF (NXTNBR.EQ.INLIST) GO TO 120 NXTCLS = NXTCLS - 1 INCLS(NXTCLS) = NXTNBR PARENT(NXTCLS) = NODE C GET NEXT NEIGHBOR 120 K = NXT(K+2) IF (K.NE.0) GO TO 100 C ADD LIST OF NEIGHBORS OF NEXT ELEMENT OF THIS CLUSTER NXTCN = NXTCN + 1 IF (NXTCN.GT.NUMIN) GO TO 130 NODE = INCLS(NXTCN) BAKWRD = PARENT(NXTCN) GO TO 90 C END OF CLUSTER--DO OUTPUT 130 IF((NUMIN-BEGCLS+1).LT.MINPTS)GOTO 150 CALL STORE(NUMIN-BEGCLS+1, C, CP, CLEN) IF (PRINT) WRITE (6,99999) CLS DO 140 I=BEGCLS,NUMIN IF (PRINT) WRITE (6,99997) INCLS(I) CALL STORE(INCLS(I), C, CP, CLEN) 140 CONTINUE GOTO 160 150 CLS = CLS - 1 160 IF (NUMIN.LT.N) GO TO 80 CP = 0 CALL STORE(CLS, C, CP, CLEN) CALL FIXMST RETURN 99999 FORMAT(1H0/8H0CLUSTER, I5, 12H CONSISTS OF) 99998 FORMAT(44H1THE TREE HAS BEEN CLUSTERED SEARCHING TO A , * 8HDEPTH OF, I3/11X, 28HINCONSISTENT EDGES HAVE BEEN, * 27H DETERMINED BY A FACTOR OF , G11.4/11X, 10HAND A SPRE, * 6HAD OF , G11.4, 21H STANDARD DEVIATIONS.) 99997 FORMAT(10X, 4HNODE, I5) END REAL FUNCTION DIST(A, B, N) C 10 INTEGER N REAL A(N), B(N) C THIS FUNCTION COMPUTES THE WEIGHT OF THE BRANCH BETWEEN C NODE A AND NODE B. IT SHOULD BE WRITTEN TO SUIT THE DATA. C THE TYPE DECLARATION OF A AND B SHOULD MATCH THE DATA. C THIS VERSION COMPUTES THE USUAL EUCLIDEAN DISTANCE. DIST = (A(1)-B(1))**2 DO 10 I=2,N DIST = DIST + (A(I)-B(I))**2 10 CONTINUE DIST = SQRT(DIST) RETURN END SUBROUTINE CLIMB(POINTR, STACK, LN, D) D 10 INTEGER POINTR(1), STACK(1), LN, D INTEGER SPACE(2), MST(1), NBR(1), NXT(1) EQUIVALENCE (MST,NBR,NXT) COMMON SPACE, MST C STARTING FROM THE NODE ON TOP OF THE STACK, CLIMB OUT C TO DEPTH D OR TO A TERMINAL NODE, WHICHEVER OCCURS FIRST 10 IF (LN.EQ.D+2) RETURN K = POINTR(LN) IF (K) 20, 30, 40 C SET POINTER TO FIRST NEIGHBOR OF TOP NODE 20 NODE = STACK(LN) POINTR(LN) = MST(NODE) + 1 GO TO 50 C BACK DOWN FROM TERMINAL NODE 30 LN = LN - 1 C CLIMB OUT ON NEXT NEIGHBOR IF POSSIBLE 40 POINTR(LN) = NXT(K+2) IF (POINTR(LN).EQ.0) RETURN C CHECK DIRECTION 50 K = POINTR(LN) NEIGHB = IABS(NBR(K)) IF (NEIGHB.EQ.STACK(LN-1)) GO TO 40 C CLIMB OUT ON NEIGHBORING NODE LN = LN + 1 STACK(LN) = NEIGHB POINTR(LN) = -1 GO TO 10 END INTEGER FUNCTION LASTPT(NODE) E 10 C THE VALUE OF THIS FUNCTION POINTS TO THE END OF THE LIST C OF NEIGHBORS OF NODE. INTEGER SPACE(2), MST(1), NXT(1) EQUIVALENCE (MST,NXT) COMMON SPACE, MST C OFFSET PICKS UP POINTER FIELD LASTPT = MST(NODE) + 3 10 IF (NXT(LASTPT).EQ.0) RETURN LASTPT = NXT(LASTPT) + 2 GO TO 10 END INTEGER FUNCTION FINDCN(A, B) F 10 INTEGER A, B INTEGER SPACE(2), MST(1), NBR(1), NXT(1) EQUIVALENCE (MST,NBR,NXT) COMMON SPACE, MST C THIS FUNCTION LOCATES NODE B IN THE LIST OF NEIGHBORS OF A C OFFSET PICKS UP NEIGHBOR FIELD FINDCN = MST(A) + 1 10 IF (IABS(NBR(FINDCN)).EQ.B) RETURN FINDCN = NXT(FINDCN+2) IF (FINDCN.NE.0) GO TO 10 WRITE (6,99999) B, A 99999 FORMAT(5H0NODE, I3, 26H IS NOT A NEIGHBOR OF NODE, I3) RETURN END SUBROUTINE STORE(VALUE, ARRAY, LOC, N) G 10 INTEGER VALUE, ARRAY(N), LOC, N C THIS SUBROUTINE IS USED TO STORE VALUES INTO THE ARRAY C WHICH IS THE FOURTH PARAMETER OF CLUSTR. IF (N.EQ.0) RETURN LOC = LOC + 1 IF (LOC.GT.N) GO TO 10 ARRAY(LOC) = VALUE RETURN 10 WRITE (6,99999) VALUE 99999 FORMAT(41H THE ARRAY USED TO STORE A DESCRIPTION OF/3H TH, * 30HE CLUSTERS IS NOT LONG ENOUGH /15H ITS NEXT VALUE, * 11H SHOULD BE , I10) RETURN END SUBROUTINE PRTREE H 10 C THE DESCRIPTION OF THE MINIMAL SPANNING TREE PRINTED HERE C LABELS EACH NODE SEQUENTIALLY AS IT OCCURS IN DATA . INTEGER DIM, N, MST(1), LOC(1), NBR(1), NXT(1) REAL WT(1) EQUIVALENCE (MST,LOC,NBR,WT,NXT) COMMON DIM, N, MST DO 20 NODE=1,N WRITE (6,99999) NODE K = MST(NODE) + 1 10 WRITE (6,99998) NBR(K), WT(K+1) K = NXT(K+2) IF (K.NE.0) GO TO 10 20 CONTINUE RETURN 99999 FORMAT(5H0NODE, I3/16H NEIGHBORS ARE) 99998 FORMAT(10X, 4HNODE, I5, 14H AT DISTANCE , G11.4) END SUBROUTINE FIXMST I 10 INTEGER DIM, N, MST(1), NBR(1), NXT(1) EQUIVALENCE (MST,NBR,NXT) COMMON DIM, N, MST DO 20 I=1,N K = MST(I) + 1 10 NBR(K) = IABS(NBR(K)) K = NXT(K+2) IF (K.NE.0) GO TO 10 20 CONTINUE RETURN END