C     SSA1                                                              SA1   1 
C     A GENERAL NONMETRIC TECHNIQUE FOR FINDING THE SMALLEST EUCLIDEAN  SA1   2 
C     SPACE FOR A CONFIGURATION OF POINTS.  (L.GUTTMAN AND J.C.LINGOES).SA1   3 
C     PROGRAMMED IN FORTRAN II FOR THE UNIV. OF MICHIGAN IBM-7090 BY    SA1   4 
C     J.C.LINGOES (9/21/64).  MAJOR CHANGES IN THE PRESENT VERSION IN-  SA1   5 
C     CLUDE OPTIONS FOR:  1) ANALYZING AN INDEFINITE NUMBER OF VARIABLESSA1   6 
C     GIVEN A FIXED CONFIGURATION, 2) MINIMIZING KRUSKAL'S STRESS, AND  SA1   7 
C     3) INPUT OF A CONFIGURATION OF ONE'S CHOICE.  ALGORITHMS EMPLOYED:SA1   8 
C     G-L "SOFT-SQUEEZE", DOUBLE-PHASE FOLLOWED BY SINGLE-PHASE (RANK-  SA1   9 
C     IMAGES) FOR SEMI-STRONG MONOTONICITY, WHICH IS OPTIONALLY FOLLOWEDSA1  10 
C     BY G-L "SOFT-SQUEEZE", SINGLE-PHASE (KRUSKAL'S MONOTONE REGRESSIONSA1  11 
C     VALUES) FOR WEAK MONOTONICITY.                                    SA1  12 
C                                                                       SA1  13 
C     ==================================================================SA1  14 
C     *                                                                *SA1  15 
C     *              - REVISED GUTTMAN-LINGOES PROGRAMS -              *SA1  16 
C     *                                                                *SA1  17 
C     *              FORTRAN IV (LEVEL G) VERSION - 1/1/69             *SA1  18 
C     *                         EDITED BY TIDY*                        *SA1  19 
C     *                                                                *SA1  20 
C     *                                                                *SA1  21 
C     *    *H. M. MURPHY, TIDY, A COMPUTER CODE FOR RENUMBERING AND    *SA1  22 
C     *         EDITING FORTRAN SOURCE PROGRAMS.  AD-642-099,          *SA1  23 
C     *             CLEARINGHOUSE, U.S. DEPT. OF COMMERCE,             *SA1  24 
C     *                  SPRINGFIELD, VIRGINIA 22151                   *SA1  25 
C     *                                                                *SA1  26 
C     ==================================================================SA1  27 
C                                                                       SA1  28 
C     DECK SET-UP FOR G-L(SSA-I) -                                      SA1  29 
C                                                                       SA1  30 
C        1.   SYSTEM ID CARD/S.                                         SA1  31 
C        2.   BINARY PROGRAM.                                           SA1  32 
C        3.   TITLE CARD (PUNCH A 1 IN COLUMN 1 AND ANY BCD TITLE IN COLSA1  33 
C             -UMNS 2 TO 72, WHICH WILL BE PRINTED OUT FOR EACH PAGE OF SA1  34 
C             OUTPUT).                                                  SA1  35 
C        4.   PARAMETER CARD, 11 4-COLUMN FIELDS CONTAINING THE FOLLOW- SA1  36 
C             ING INFORMATION SERIATUM -                                SA1  37 
C             A)  NR=THE NUMBER OF VARIABLES .LE. 100 AND .GE. 3,       SA1  38 
C             B)  MIND=0 OR BLANK IF THE PROGRAM IS TO DETERMINE THE MINSA1  39 
C                 -IMUM NUMBER OF DIMENSIONS FOR THE PROBLEM, OTHERWISE SA1  40 
C                 ANY NUMBER BETWEEN 1 AND 10 PROVIDED MIND .LE. MAXD,  SA1  41 
C                 IN WHICH CASE ALL SOLUTIONS FROM MIND TO MAXD WILL BE SA1  42 
C                 PRINTED OUT, UNLESS K OR STRESS BECOMES .LE. .0001    SA1  43 
C                 FOR A GIVEN M, THE NUMBER OF DIMENSIONS.  UNLESS YOU  SA1  44 
C                 KNOW M, SET MIND=0, IN GENERAL,                       SA1  45 
C             C)  MAXD=THE LARGEST NUMBER OF DIMENSIONS DESIRED .LE.    SA1  46 
C                 (NV-1,10)MIN BUT .GE. 1,                              SA1  47 
C             D)  ISIM=0 OR BLANK FOR DISTANCE COEFFICIENTS OR DISSIMI- SA1  48 
C                 LARITY DATA AND 1 IF SIMILARITY DATA, E.G., CORRELA-  SA1  49 
C                 TIONS,                                                SA1  50 
C             E)  IFD=1 IF DISTANCE MATRIX IS TO BE PRINTED FOR 2 OR    SA1  51 
C                 MORE DIMENSIONS, OTHERWISE SET TO ZERO OR LEAVE BLANK,SA1  52 
C             F)  IFC=1 IF COORDINATES ARE TO BE PUNCHED FOR 2 OR MORE  SA1  53 
C                 DIMENSIONS, OTHERWISE SET TO ZERO OR LEAVE BLANK.     SA1  54 
C                 CARDS WILL BE AUTOMATICALLY PUNCHED IF MORE THAN 100  SA1  55 
C                 ITERATIONS ARE REQUIRED FOR CONVERGENCE OF 2 OR MORE  SA1  56 
C                 DIMENSIONS.  THESE CARDS CAN BE USED FOR SUBSEQUENT   SA1  57 
C                 INPUT TO CONTINUE THE ITERATIONS,                     SA1  58 
C             G)  IFGLK=1 IF KRUSKAL'S STRESS IS TO BE MINIMIZED, OTHER-SA1  59 
C                 WISE SET TO ZERO OR LEAVE BLANK,                      SA1  60 
C             H)  IFCONF=1 IF A CONFIGURATION IS TO BE INPUT FOR CON-   SA1  61 
C                 TINUED ITERATIONS OR FOR ADDING POINTS TO A FIXED CON-SA1  62 
C                 FIGURATION, OTHERWISE LEAVE BLANK OR SET TO ZERO.     SA1  63 
C                 MAXD MUST BE THE NUMBER OF DIMENSIONS INPUT,          SA1  64 
C             I)  IFFIX=1 IF INPUT CONFIGURATION IS TO REMAIN FIXED AND SA1  65 
C                 ADDITIONAL POINTS ARE TO BE FITTED TO THIS SPACE,     SA1  66 
C                 OTHERWISE LEAVE BLANK OR SET TO ZERO.  MIND=MAXD.     SA1  67 
C                 TO ACCOMPLISH ANALYSIS OF MORE VARIABLES THAN 100:  1)SA1  68 
C                 ANALYZE SOME REPRESENTATIVE SAMPLE OF VARIABLES (S'S),SA1  69 
C                 SPECIFYING IFC=1;  2) USE THE OUTPUT SOLUTION AS IN-  SA1  70 
C                 PUT FOR IFFIX=1.  NR IS THE ORIGINAL NUMBER OF VARIA- SA1  71 
C                 BLES.  ITEMS E) AND G) ARE DISABLED WHEN IFFIX=1.     SA1  72 
C                 A CONDITIONAL APPROACH IS MADE TO PRESERVE THE ORDER  SA1  73 
C                 FOR JUST THE ARRAY OF VALUES INPUT FOR EACH VARIABLE. SA1  74 
C                 THE COEFFICIENT OF FIT APPLIES ONLY TO THESE VALUES   SA1  75 
C                 VIS-A-VIS THE ORIGINAL SET OF POINTS,                 SA1  76 
C             J)  IFSR=1 IF EITHER ITEM 7. OR 8. IS TO BE GENERATED BY ASA1  77 
C                 SUBROUTINE (WHICH THE USER SUBSTITUTES FOR THE DUMMY  SA1  78 
C                 SUBROUTINES PROVIDED), OTHERWISE SET TO ZERO OR LEAVE SA1  79 
C                 BLANK,                                                SA1  80 
C             K)  EPS=0 OR BLANK IF TIED BLOCKS ARE NOT TO BE FORMED,   SA1  81 
C                 OTHERWISE INSERT IN F-NOTATION (WITH DECIMAL POINT    SA1  82 
C                 PUNCHED) THE CATEGORY WIDTH FOR COEFFICIENTS WHICH ARESA1  83 
C                 TO BE TIED.                                           SA1  84 
C        5.   FORMAT CARD (DESCRIBING IN F-NOTATION WHERE THE DATA      SA1  85 
C             APPEARS ON THE CARDS).                                    SA1  86 
C        6.   IF IFCONF=1, PUNCH CONFIGURATION TO BE INPUT (OR USE OUT- SA1  87 
C             PUT CONFIGURATION) ACCORDING TO FORMAT: 10F8.3.  THERE    SA1  88 
C             SHOULD BE NR SETS OF CARDS, EACH SET OF WHICH SHOULD HAVE SA1  89 
C             MAXD COORDINATES, OTHERWISE OMIT THIS SET WHEN IFCONF=0.  SA1  90 
C        7.   DATA (PUNCH LOWER-HALF OF THE SQUARE-SYMMETRIC MATRIX WITHSA1  91 
C             -OUT THE DIAGONAL ELEMENTS, STARTING A NEW ROW ON A NEW   SA1  92 
C             CARD.  IN TOTAL YOU SHOULD HAVE NR-1 SETS OF CARDS WITH   SA1  93 
C             1 ELEMENT IN THE FIRST SET FOR THE SECOND ROW, 2 ELEMENTS SA1  94 
C             IN THE SECOND SET FOR THE THIRD ROW, ... , AND NR-1 COEF- SA1  95 
C             FICIENTS IN THE LAST OR NR-1'ST SET).  ALL ROWS MUST BE   SA1  96 
C             LEFT-ADJUSTED.  OMIT WHEN IFFIX=1.                        SA1  97 
C        8.   FOR EACH VARIABLE TO BE ADDED (WHEN IFFFIX=1) PUNCH NR    SA1  98 
C             COEFFICIENTS ON AS MANY CARDS AS NECESSARY TO SATISFY NR, SA1  99 
C             ACCORDING TO THE SAME FORMAT AS IN ITEM 5.  WHEN IFFIX=0, SA1 100 
C             OMIT THIS SET OF CARDS.                                   SA1 101 
C        9.   REPEAT ITEMS 3-8 FOR ADDITIONAL RUNS.  ONLY 1 RUN OF A    SA1 102 
C             FIXED CONFIGURATION CAN BE MADE AT A TIME AND THIS SHOULD SA1 103 
C             APPEAR AS THE LAST JOB RUN IN A DECK OF JOBS.             SA1 104 
C                                                                       SA1 105 
C     *** USERS OF THIS PROGRAM ARE EXPECTED TO PROPERLY CREDIT SOURCE  SA1 106 
C     FROM REFERENCES LISTED BELOW.  IF IFGLK IS SET TO 1, THEN A REFER-SA1 107 
C     ENCE TO KRUSKAL'S PAIR OF 1964 PSYCHOMETRIKA PAPERS SHOULD ALSO BESA1 108 
C     MADE. ***                                                         SA1 109 
C                                                                       SA1 110 
C     *** REFERENCES - GUTTMAN, L. A GENERAL NONMETRIC TECHNIQUE FOR    SA1 111 
C                                  FINDING THE SMALLEST COORDINATE SPACESA1 112 
C                                  FOR A CONFIGURATION OF POINTS.  PSY- SA1 113 
C                                  CHOMETRIKA, 1968, 33, 469-506.       SA1 114 
C                      LINGOES, J. C.  NEW COMPUTER DEVELOPMENTS IN PAT-SA1 115 
C                                  TERN ANALYSIS AND NONMETRIC TECH-    SA1 116 
C                                  NIQUES.  IN - USES OF COMPUTERS IN   SA1 117 
C                                  PSYCHOLOGICAL RESEARCH.  GAUTHIER-   SA1 118 
C                                  VILLARS, PARIS, 1966, 1-22.          SA1 119 
C                      LINGOES, J. C.  AN IBM-7090 PROGRAM FOR GUTTMAN- SA1 120 
C                                  LINGOES SMALLEST SPACE ANALYSIS - I. SA1 121 
C                                  BEHAV. SCI., 1965,10,183-184.        SA1 122 
C                      LINGOES, J.C., ROSKAM, E.E.C.I., & GUTTMAN, L.   SA1 123 
C                                  AN EMPIRICAL STUDY OF TWO MULTIDIMEN-SA1 124 
C                                  SIONAL SCALING ALGORITHMS.  MULTIV.  SA1 125 
C                                  BEHAV. RES., 1969, 4,                SA1 126 
C                                                                       SA1 127 
C                                                                       SA1 128 
      DIMENSION RHO(100,100), RHO1(100,100), C(100,100), EVAL(10000), X(SA1 129 
     1101,10), X2(101,11), INDI(5050), INDJ(5050), DIST(5050), PHI(11), SA1 130 
     2NOTIES(1000), FMT(18), PROX(5050), CMEAN(10)                      SA1 131 
      EQUIVALENCE (RHO1,C), (RHO,EVAL), (DIST,PROX)                     SA1 132 
C                                                                       SA1 133 
C     *** SUBROUTINES CALLED:  MXOUT, FIT, SORT, EIGEN, AND PLOT ***    SA1 134 
C     *** DUMMY SUBROUTINES:  COR1 & COR2 ***                           SA1 135 
C                                                                       SA1 136 
C     *** TAPE ASSIGNMENTS                                              SA1 137 
      ITAPE=5                                                           SA1 138 
      JTAPE=6                                                           SA1 139 
C     *** SUBROUTINE DIMENSIONING                                       SA1 140 
      MD=100                                                            SA1 141 
      ND=11                                                             SA1 142 
      LD=MD+1                                                           SA1 143 
C     READ IN TITLE, PARAMETERS, AND FORMAT                             SA1 144 
1     WRITE (JTAPE,166)                                                 SA1 145 
      READ (ITAPE,167)                                                  SA1 146 
      READ (ITAPE,168) NR,MIND,MAXD,ISIM,IFD,IFC,IFGLK,IFCONF,IFFIX,IFSRSA1 147 
     1,EPS                                                              SA1 148 
      READ (ITAPE,169) (FMT(J),J=1,18)                                  SA1 149 
C     INITIALIZATION                                                    SA1 150 
      ISW=0                                                             SA1 151 
      NIT=5                                                             SA1 152 
      IFK=0                                                             SA1 153 
      ITER=100                                                          SA1 154 
      NRP1=NR+1                                                         SA1 155 
      NRM1=NR-1                                                         SA1 156 
      FNR=NR                                                            SA1 157 
      IF (IFCONF) 2,7,2                                                 SA1 158 
C     INPUT CONFIGURATION                                               SA1 159 
2     DO 3 I=1,NR                                                       SA1 160 
3     READ (ITAPE,170) (X(I,J),J=1,MAXD)                                SA1 161 
C     SET COORDINATE MEANS TO ZERO                                      SA1 162 
      DO 6 J=1,MAXD                                                     SA1 163 
      CMEAN(J)=0.                                                       SA1 164 
      DO 4 I=1,NR                                                       SA1 165 
4     CMEAN(J)=CMEAN(J)+X(I,J)                                          SA1 166 
      CMEAN(J)=CMEAN(J)/FNR                                             SA1 167 
      DO 5 K=1,NR                                                       SA1 168 
5     X(K,J)=X(K,J)-CMEAN(J)                                            SA1 169 
6     CONTINUE                                                          SA1 170 
7     IF (IFFIX) 8,15,8                                                 SA1 171 
C     INITIALIZE FOR FIXED-CONFIGURATION OPTION                         SA1 172 
8     NEL=NR                                                            SA1 173 
      FNR=NRP1                                                          SA1 174 
      IND=NR                                                            SA1 175 
      MIND=MAXD                                                         SA1 176 
      WRITE (JTAPE,167)                                                 SA1 177 
      WRITE (JTAPE,171) MAXD,(MM,MM=1,MAXD)                             SA1 178 
      WRITE (JTAPE,177)                                                 SA1 179 
      WRITE (JTAPE,173)                                                 SA1 180 
C     INPUT COEFFICIENTS FOR ADDED VARIABLE                             SA1 181 
9     IF (IFSR) 10,11,10                                                SA1 182 
10    CALL COR2 (NR,IFSR,ISW,MD,FMT,PROX,C)                             SA1 183 
      GO TO 12                                                          SA1 184 
11    READ (ITAPE,FMT) (PROX(J),J=1,NR)                                 SA1 185 
12    DO 13 J=1,MAXD                                                    SA1 186 
13    X(101,J)=0.                                                       SA1 187 
      DO 14 J=1,NEL                                                     SA1 188 
14    EVAL(J)=PROX(J)                                                   SA1 189 
      GO TO 21                                                          SA1 190 
C     INPUT SIMILARITIES/DISSIMILARITIES                                SA1 191 
15    IF (IFSR) 16,17,16                                                SA1 192 
16    CALL COR1 (NR,MD,FMT,PROX,C,IFSR,PROX(LD),EVAL)                   SA1 193 
      GO TO 19                                                          SA1 194 
17    DO 18 I=2,NR                                                      SA1 195 
      IP1=I-1                                                           SA1 196 
18    READ (ITAPE,FMT) (C(I,J),J=1,IP1)                                 SA1 197 
19    JJ=0                                                              SA1 198 
      DO 20 I=1,NRM1                                                    SA1 199 
      IP1=I+1                                                           SA1 200 
      C(I,I)=0.                                                         SA1 201 
      DO 20 J=IP1,NR                                                    SA1 202 
      JJ=JJ+1                                                           SA1 203 
      PROX(JJ)=C(J,I)                                                   SA1 204 
20    EVAL(JJ)=PROX(JJ)                                                 SA1 205 
      C(NR,NR)=0.                                                       SA1 206 
      NEL=JJ                                                            SA1 207 
      NELP1=NEL+1                                                       SA1 208 
      FNEL=NEL                                                          SA1 209 
      LFACT=MIND                                                        SA1 210 
C     PRINT OUT INPUT COEFFICIENTS                                      SA1 211 
      CALL MXOUT (RHO1,NR,0,MD)                                         SA1 212 
C     SORT SUBSCRIPTS OF INDI ACCORDING TO VALUES IN DIST AND DIREC-    SA1 213 
C     TION OF ISIM                                                      SA1 214 
21    CALL SORT (ISW,NEL,PROX,ISIM,INDI)                                SA1 215 
      IF (IFFIX) 24,22,24                                               SA1 216 
C     SUBSTITUTE RANKS FOR SIMILARITIES/DISSIMILARITIES                 SA1 217 
22    S=0.                                                              SA1 218 
      DO 23 J=1,NEL                                                     SA1 219 
      MM=INDI(J)                                                        SA1 220 
      S=S+1.                                                            SA1 221 
23    PROX(MM)=S                                                        SA1 222 
C     CHECK FOR TIES IN INPUT VALUES                                    SA1 223 
24    NOTIES(1)=0                                                       SA1 224 
      J=1                                                               SA1 225 
      II=-1                                                             SA1 226 
      K=J                                                               SA1 227 
      NT=0                                                              SA1 228 
25    LL=1                                                              SA1 229 
      MM=INDI(J)                                                        SA1 230 
      D1=PROX(MM)                                                       SA1 231 
26    K=K+1                                                             SA1 232 
      NN=INDI(K)                                                        SA1 233 
      IF (ABS(EVAL(MM)-EVAL(NN))-EPS) 27,27,28                          SA1 234 
27    LL=LL+1                                                           SA1 235 
      D1=D1+PROX(NN)                                                    SA1 236 
      IF (K-NEL) 26,29,26                                               SA1 237 
28    IF (LL-1) 29,32,29                                                SA1 238 
29    II=II+2                                                           SA1 239 
      NT=NT+2                                                           SA1 240 
C     *** IF NOTIES IS RE-DIMENSIONED, ALSO ADJUST FOLLOWING TEST       SA1 241 
      IF (NT.GT.1000) GO TO 163                                         SA1 242 
      NOTIES(II)=LL                                                     SA1 243 
      NOTIES(II+1)=J                                                    SA1 244 
      IF (IFFIX) 32,30,32                                               SA1 245 
30    D1=D1/FLOAT(LL)                                                   SA1 246 
      JJ=J                                                              SA1 247 
      DO 31 KK=1,LL                                                     SA1 248 
      MM=INDI(JJ)                                                       SA1 249 
      PROX(MM)=D1                                                       SA1 250 
31    JJ=JJ+1                                                           SA1 251 
32    J=K                                                               SA1 252 
      IF (J-NEL) 25,33,25                                               SA1 253 
33    IF (IFFIX+IFCONF) 47,34,47                                        SA1 254 
C     FORM C-MATRIX BASED ON RANKS                                      SA1 255 
34    II=0                                                              SA1 256 
      DO 35 I=1,NRM1                                                    SA1 257 
      IP1=I+1                                                           SA1 258 
      DO 35 J=IP1,NR                                                    SA1 259 
      II=II+1                                                           SA1 260 
      RHO1(I,J)=1.-PROX(II)/FNEL                                        SA1 261 
35    RHO1(J,I)=RHO1(I,J)                                               SA1 262 
      DO 38 I=1,NR                                                      SA1 263 
      RHO1(I,I)=FNR                                                     SA1 264 
      DO 37 J=1,NR                                                      SA1 265 
      IF (I-J) 36,37,36                                                 SA1 266 
36    RHO1(I,I)=RHO1(I,I)-RHO1(I,J)                                     SA1 267 
37    CONTINUE                                                          SA1 268 
38    CONTINUE                                                          SA1 269 
C     OBTAIN INITIAL CONFIGURATION                                      SA1 270 
C     CALL HOUSEHOLDER SUBROUTINE                                       SA1 271 
      KK=MAXD+1                                                         SA1 272 
      CALL EIGEN (C,RHO,NR,PHI,KK,MD,X(1,1),X(1,4),X(1,5),X(1,6),X(1,7),SA1 273 
     1X(1,8),1)                                                         SA1 274 
      DO 39 J=2,KK                                                      SA1 275 
      PHI(J-1)=PHI(J)                                                   SA1 276 
      DO 39 I=1,NR                                                      SA1 277 
39    X(I,J-1)=RHO(I,J)                                                 SA1 278 
C     DETERMINE WHETHER TO GO UP OR DOWN                                SA1 279 
      IF (MIND) 47,40,47                                                SA1 280 
40    DIM=FNR/2.                                                        SA1 281 
      KK=0                                                              SA1 282 
      DO 42 J=1,MAXD                                                    SA1 283 
      IF (PHI(J)-DIM) 43,41,41                                          SA1 284 
41    KK=KK+1                                                           SA1 285 
42    CONTINUE                                                          SA1 286 
43    IF (KK) 45,44,45                                                  SA1 287 
44    MAXD=1                                                            SA1 288 
      GO TO 46                                                          SA1 289 
45    MAXD=KK                                                           SA1 290 
46    MIND=MAXD                                                         SA1 291 
      ITER=25                                                           SA1 292 
47    M=MIND-1                                                          SA1 293 
C     INITIALIZATION                                                    SA1 294 
48    M=M+1                                                             SA1 295 
      ASSIGN 61 TO N1                                                   SA1 296 
      IF (M.EQ.1) NIT=1                                                 SA1 297 
49    STRLST=1.                                                         SA1 298 
      AVST=0.                                                           SA1 299 
      IFBU=0                                                            SA1 300 
      NN=0                                                              SA1 301 
      II=5                                                              SA1 302 
50    ITCT=0                                                            SA1 303 
51    NN=NN+1                                                           SA1 304 
C     CALCULATE DISTANCES                                               SA1 305 
52    IF (IFFIX) 53,56,53                                               SA1 306 
53    DO 55 I=1,NR                                                      SA1 307 
      DIST(I)=0.                                                        SA1 308 
      DO 54 J=1,M                                                       SA1 309 
54    DIST(I)=DIST(I)+(X(I,J)-X(101,J))**2                              SA1 310 
55    DIST(I)=SQRT(DIST(I))                                             SA1 311 
      GO TO 59                                                          SA1 312 
56    JJ=0                                                              SA1 313 
      DO 58 I=1,NRM1                                                    SA1 314 
      IP1=I+1                                                           SA1 315 
      DO 58 J=IP1,NR                                                    SA1 316 
      JJ=JJ+1                                                           SA1 317 
      DIST(JJ)=0.                                                       SA1 318 
      DO 57 K=1,M                                                       SA1 319 
57    DIST(JJ)=DIST(JJ)+(X(I,K)-X(J,K))**2                              SA1 320 
58    DIST(JJ)=SQRT(DIST(JJ))                                           SA1 321 
      IF (IFK) 63,59,63                                                 SA1 322 
59    IF (ITCT) 94,60,94                                                SA1 323 
60    GO TO N1, (61,62)                                                 SA1 324 
C     SORT DISTANCES LOW TO HIGH                                        SA1 325 
61    CALL SORT (0,NEL,DIST,0,INDJ)                                     SA1 326 
      ASSIGN 62 TO N1                                                   SA1 327 
      GO TO 63                                                          SA1 328 
62    CALL SORT (1,NEL,DIST,0,INDJ)                                     SA1 329 
C     IF TIES EXIST IN INPUT, OPTIMALLY PERMUTE INDICES OF INDI TO      SA1 330 
C     CONFORM WITH ORDER OF CORRESPONDING DISTANCES.  PRIMARY APPROACH. SA1 331 
63    IF (NOTIES(1)) 64,66,64                                           SA1 332 
64    DO 65 J=1,NT,2                                                    SA1 333 
      KK=NOTIES(J)                                                      SA1 334 
      LL=NOTIES(J+1)                                                    SA1 335 
      CALL SORT (1,KK,DIST,0,INDI(LL))                                  SA1 336 
65    CONTINUE                                                          SA1 337 
66    IF (IFFIX) 69,67,69                                               SA1 338 
67    IF (IFK) 68,69,68                                                 SA1 339 
C     OBTAIN KRUSKAL'S MONOTONE REGRESSION BEST-FIT VALUES              SA1 340 
68    CALL FIT (NEL,DIST,INDI,EVAL(NELP1),EVAL,INDJ)                    SA1 341 
      GO TO 71                                                          SA1 342 
C     CELL-WISE PERMUTE D TO D* (RANK-IMAGES)                           SA1 343 
69    DO 70 I=1,NEL                                                     SA1 344 
70    EVAL(INDI(I))=DIST(INDJ(I))                                       SA1 345 
C     CALCULATE NORMALIZED PHI/STRESS**2                                SA1 346 
71    STRESS=0.                                                         SA1 347 
      D1=0.                                                             SA1 348 
      DO 72 I=1,NEL                                                     SA1 349 
      D1=D1+DIST(I)**2                                                  SA1 350 
72    STRESS=STRESS+DIST(I)*EVAL(I)                                     SA1 351 
      STRESS=1.-STRESS/D1                                               SA1 352 
      IF (STRESS.GT.STRLST) GO TO 78                                    SA1 353 
C     SAVE BEST CONFIGURATION                                           SA1 354 
      IF (IFFIX) 73,75,73                                               SA1 355 
73    DO 74 J=1,MAXD                                                    SA1 356 
74    X2(101,J)=X(101,J)                                                SA1 357 
      GO TO 77                                                          SA1 358 
75    DO 76 I=1,NR                                                      SA1 359 
      DO 76 J=1,M                                                       SA1 360 
76    X2(I,J)=X(I,J)                                                    SA1 361 
77    STRLST=STRESS                                                     SA1 362 
C     TEST FOR TERMINATION                                              SA1 363 
78    IF (STRESS.LT..000001) GO TO 118                                  SA1 364 
      IF (NN.GE.6) GO TO 79                                             SA1 365 
      AVST=AVST+STRESS                                                  SA1 366 
      PHI(NN)=STRESS                                                    SA1 367 
      GO TO 90                                                          SA1 368 
79    IF (II.EQ.5) II=0                                                 SA1 369 
      II=II+1                                                           SA1 370 
      IF (IFK.EQ.1.AND.((5.*STRESS)/AVST).GT..995) GO TO 83             SA1 371 
      IF (STRESS.LE.PHI(II)) GO TO 87                                   SA1 372 
      IFBU=IFBU+1                                                       SA1 373 
C     RESTORE BEST CONFIGURATION                                        SA1 374 
80    IF (IFFIX) 81,83,81                                               SA1 375 
81    DO 82 J=1,MAXD                                                    SA1 376 
82    X(101,J)=X2(101,J)                                                SA1 377 
      GO TO 85                                                          SA1 378 
83    DO 84 I=1,NR                                                      SA1 379 
      DO 84 J=1,M                                                       SA1 380 
84    X(I,J)=X2(I,J)                                                    SA1 381 
85    STRESS=STRLST                                                     SA1 382 
      IF (IFBU.EQ.2) GO TO 118                                          SA1 383 
      IF (IFK) 118,86,118                                               SA1 384 
86    II=II-1                                                           SA1 385 
      GO TO 50                                                          SA1 386 
87    IF (NN-ITER) 89,88,89                                             SA1 387 
88    IFBU=2                                                            SA1 388 
      GO TO 80                                                          SA1 389 
89    AVST=AVST-PHI(II)+STRESS                                          SA1 390 
      PHI(II)=STRESS                                                    SA1 391 
90    IF (IFK) 94,91,94                                                 SA1 392 
91    IF (IFBU-1) 94,92,94                                              SA1 393 
C     SWITCH TO SINGLE-PHASE ALGORITHM (RANK-IMAGES)                    SA1 394 
92    DO 93 J=1,NEL                                                     SA1 395 
93    EVAL(INDJ(J))=(DIST(INDI(J))+EVAL(INDJ(J)))/2.                    SA1 396 
94    IF (IFFIX) 95,104,95                                              SA1 397 
C     DETERMINE CORRECTIONS FOR ADDED VARIABLE                          SA1 398 
95    DO 100 I=1,NR                                                     SA1 399 
      IF (EVAL(I)+DIST(I)) 97,96,97                                     SA1 400 
96    X2(I,11)=0.                                                       SA1 401 
      GO TO 100                                                         SA1 402 
97    IF (DIST(I)) 99,98,99                                             SA1 403 
98    X2(I,11)=1.-EVAL(I)/.0001                                         SA1 404 
      GO TO 100                                                         SA1 405 
99    X2(I,11)=1.-EVAL(I)/DIST(I)                                       SA1 406 
100   CONTINUE                                                          SA1 407 
      X2(101,11)=FNR                                                    SA1 408 
      DO 101 I=1,NR                                                     SA1 409 
101   X2(101,11)=X2(101,11)-X2(I,11)                                    SA1 410 
C     MODIFY COORDINATES OF ADDED VARIABLE                              SA1 411 
      DO 103 J=1,MAXD                                                   SA1 412 
      X2(101,J)=0.                                                      SA1 413 
      DO 102 I=1,NR                                                     SA1 414 
102   X2(101,J)=X2(101,J)+X(I,J)*X2(I,11)                               SA1 415 
      X2(101,J)=(X2(101,J)+X(101,J)*X2(101,11))/FNR                     SA1 416 
103   X(101,J)=X2(101,J)                                                SA1 417 
      GO TO 116                                                         SA1 418 
C     COMPUTE C-MATRIX                                                  SA1 419 
104   JJ=0                                                              SA1 420 
      DO 109 I=1,NRM1                                                   SA1 421 
      IP1=I+1                                                           SA1 422 
      DO 109 J=IP1,NR                                                   SA1 423 
      JJ=JJ+1                                                           SA1 424 
      IF (DIST(JJ)+EVAL(JJ)) 106,105,106                                SA1 425 
105   C(I,J)=0.                                                         SA1 426 
      GO TO 109                                                         SA1 427 
106   IF (DIST(JJ)) 108,107,108                                         SA1 428 
107   C(I,J)=1.-EVAL(JJ)/.0001                                          SA1 429 
      GO TO 109                                                         SA1 430 
108   C(I,J)=1.-EVAL(JJ)/DIST(JJ)                                       SA1 431 
109   C(J,I)=C(I,J)                                                     SA1 432 
      DO 112 I=1,NR                                                     SA1 433 
      C(I,I)=FNR                                                        SA1 434 
      DO 111 J=1,NR                                                     SA1 435 
      IF (I-J) 110,111,110                                              SA1 436 
110   C(I,I)=C(I,I)-C(I,J)                                              SA1 437 
111   CONTINUE                                                          SA1 438 
112   CONTINUE                                                          SA1 439 
C     APPLY CORRECTIONS TO X                                            SA1 440 
      DO 115 K=1,M                                                      SA1 441 
      DO 114 I=1,NR                                                     SA1 442 
      X2(I,11)=0.                                                       SA1 443 
      DO 113 J=1,NR                                                     SA1 444 
113   X2(I,11)=X2(I,11)+X(J,K)*C(I,J)                                   SA1 445 
114   X2(I,11)=X2(I,11)/FNR                                             SA1 446 
      DO 115 L=1,NR                                                     SA1 447 
115   X(L,K)=X2(L,11)                                                   SA1 448 
116   IF (IFK+IFBU) 51,117,51                                           SA1 449 
C     TEST FOR TERMINATION OF PHASE-I ITERATIONS                        SA1 450 
117   ITCT=ITCT+1                                                       SA1 451 
      IF (NIT-ITCT) 52,50,52                                            SA1 452 
C     COMPUTE FINAL DISTANCES                                           SA1 453 
118   IF (IFFIX) 164,119,164                                            SA1 454 
119   WRITE (JTAPE,167)                                                 SA1 455 
      JJ=0                                                              SA1 456 
      D1=0.                                                             SA1 457 
      DO 121 I=1,NRM1                                                   SA1 458 
      IP1=I+1                                                           SA1 459 
      DO 121 J=IP1,NR                                                   SA1 460 
      JJ=JJ+1                                                           SA1 461 
      DIST(JJ)=0.                                                       SA1 462 
      DO 120 K=1,M                                                      SA1 463 
120   DIST(JJ)=DIST(JJ)+(X(I,K)-X(J,K))**2                              SA1 464 
      D1=D1+DIST(JJ)                                                    SA1 465 
121   DIST(JJ)=SQRT(DIST(JJ))                                           SA1 466 
      IF (IFK) 123,122,123                                              SA1 467 
122   CALL FIT (NEL,DIST,INDI,EVAL(NELP1),EVAL,INDJ)                    SA1 468 
      GO TO 125                                                         SA1 469 
123   CALL SORT (0,NEL,DIST,0,INDJ)                                     SA1 470 
      DO 124 J=1,NEL                                                    SA1 471 
124   EVAL(INDI(J))=DIST(INDJ(J))                                       SA1 472 
125   S=0.                                                              SA1 473 
      DO 126 J=1,NEL                                                    SA1 474 
126   S=S+DIST(J)*EVAL(J)                                               SA1 475 
      IF (IFK) 128,127,128                                              SA1 476 
127   S=SQRT(1.-S/D1)                                                   SA1 477 
      STRESS=SQRT(1.-(1.-STRESS)**2)                                    SA1 478 
      WRITE (JTAPE,171) M,(MM,MM=1,M)                                   SA1 479 
      GO TO 129                                                         SA1 480 
128   S=SQRT(1.-(S/D1)**2)                                              SA1 481 
      STRESS=SQRT(STRESS)                                               SA1 482 
      WRITE (JTAPE,172) M,(MM,MM=1,M)                                   SA1 483 
129   WRITE (JTAPE,173)                                                 SA1 484 
      IF (M-1) 130,135,130                                              SA1 485 
C     PERFORM PRINCIPAL AXIS ROTATION                                   SA1 486 
130   DO 132 I=1,M                                                      SA1 487 
      DO 132 J=I,M                                                      SA1 488 
      C(I,J)=0.                                                         SA1 489 
      DO 131 K=1,NR                                                     SA1 490 
131   C(I,J)=C(I,J)+X(K,I)*X(K,J)                                       SA1 491 
132   C(J,I)=C(I,J)                                                     SA1 492 
C     CALL HOUSEHOLDER SUBROUTINE                                       SA1 493 
      CALL EIGEN (C,RHO,M,PHI,M,MD,X2(1,1),X2(1,4),X2(1,5),X2(1,6),X2(1,SA1 494 
     17),X2(1,8),0)                                                     SA1 495 
      DO 134 I=1,NR                                                     SA1 496 
      DO 133 J=1,M                                                      SA1 497 
      X2(I,J)=0.                                                        SA1 498 
      DO 133 K=1,M                                                      SA1 499 
133   X2(I,J)=X2(I,J)+X(I,K)*RHO(K,J)                                   SA1 500 
      DO 134 L=1,M                                                      SA1 501 
134   X(I,L)=X2(I,L)                                                    SA1 502 
C     NORMALIZE COORDINATES TO LIE IN RANGE OF +1 TO -1                 SA1 503 
135   DIM=0.                                                            SA1 504 
      DO 137 K=1,M                                                      SA1 505 
      RHO(1,K)=0.                                                       SA1 506 
      RHO(2,K)=0.                                                       SA1 507 
      DO 136 I=1,NR                                                     SA1 508 
      RHO(1,K)=AMIN1(RHO(1,K),X(I,K))                                   SA1 509 
136   RHO(2,K)=AMAX1(RHO(2,K),X(I,K))                                   SA1 510 
137   DIM=AMAX1(DIM,(RHO(2,K)-RHO(1,K)))                                SA1 511 
      DIM=2./DIM                                                        SA1 512 
      DO 138 K=1,M                                                      SA1 513 
      DO 138 I=1,NR                                                     SA1 514 
138   X2(I,K)=DIM*(X(I,K)-RHO(1,K))-1.                                  SA1 515 
C     COMPUTE CENTRALITY INDEX, I.E., DISTANCE FROM TRUE ORIGIN         SA1 516 
      MM=M+1                                                            SA1 517 
      DO 140 I=1,NR                                                     SA1 518 
      X2(I,MM)=0.                                                       SA1 519 
      DO 139 K=1,M                                                      SA1 520 
139   X2(I,MM)=X2(I,MM)+X(I,K)**2                                       SA1 521 
      X2(I,MM)=DIM*(SQRT(X2(I,MM)))                                     SA1 522 
140   CONTINUE                                                          SA1 523 
C     PRINT OUT FINAL CONFIGURATION                                     SA1 524 
      DO 141 I=1,NR                                                     SA1 525 
141   WRITE (JTAPE,174) I,(X2(I,K),K=1,MM)                              SA1 526 
      IF (IFK) 143,142,143                                              SA1 527 
142   WRITE (JTAPE,175) STRESS,NN,S                                     SA1 528 
      GO TO 144                                                         SA1 529 
143   WRITE (JTAPE,176) STRESS,NN,S                                     SA1 530 
144   IFGLK=-IFGLK                                                      SA1 531 
      IF (IFGLK) 145,147,146                                            SA1 532 
C     SWITCH TO SINGLE-PHASE ALGORITHM (KRUSKAL'S MONOTONE REGRESSION)  SA1 533 
145   IFK=1                                                             SA1 534 
      GO TO 49                                                          SA1 535 
C     SWITCH TO DOUBLE-PHASE ALGORITHM (RANK-IMAGES)                    SA1 536 
146   IFK=0                                                             SA1 537 
147   IF (LFACT.EQ.0.AND.M.GT.1) GO TO 156                              SA1 538 
148   IF (M-1) 149,162,149                                              SA1 539 
C     CALL ON PLOT SUBROUTINE                                           SA1 540 
149   CALL PLOT (X2,RHO(1,1),RHO(1,2),RHO(1,3),RHO(1,4),RHO(1,5),NR,M,LDSA1 541 
     1,ND)                                                              SA1 542 
      IF (NN.LT.ITER) GO TO 150                                         SA1 543 
      WRITE (JTAPE,180)                                                 SA1 544 
      GO TO 151                                                         SA1 545 
150   IF (IFC) 151,153,151                                              SA1 546 
C     PUNCH COORDINATES                                                 SA1 547 
151   DO 152 I=1,NR                                                     SA1 548 
152   PUNCH 170, (X2(I,J),J=1,M)                                        SA1 549 
153   IF (IFD) 154,162,154                                              SA1 550 
C     PRINT DISTANCES                                                   SA1 551 
154   JJ=0                                                              SA1 552 
      DO 155 I=1,NRM1                                                   SA1 553 
      C(I,I)=0.                                                         SA1 554 
      IP1=I+1                                                           SA1 555 
      DO 155 J=IP1,NR                                                   SA1 556 
      JJ=JJ+1                                                           SA1 557 
155   C(J,I)=DIST(JJ)*DIM                                               SA1 558 
      C(NR,NR)=0.                                                       SA1 559 
      CALL MXOUT (C,NR,1,MD)                                            SA1 560 
      GO TO 162                                                         SA1 561 
C     DETERMINE HOW MANY DIMENSIONS TO DROP WHEN MIND=0                 SA1 562 
156   KK=1                                                              SA1 563 
      MM=M                                                              SA1 564 
      DO 157 J=2,MM                                                     SA1 565 
      IF (PHI(J)/PHI(1).LE..3) GO TO 158                                SA1 566 
      KK=KK+1                                                           SA1 567 
157   CONTINUE                                                          SA1 568 
158   M=KK-1                                                            SA1 569 
      MAXD=KK                                                           SA1 570 
      IF (ITER-100) 159,160,159                                         SA1 571 
159   ITER=ITER+25                                                      SA1 572 
160   IF (MM-KK) 48,161,48                                              SA1 573 
161   M=M+1                                                             SA1 574 
      GO TO 148                                                         SA1 575 
C     TEST FOR END                                                      SA1 576 
162   IF (STRESS.LE..0001) GO TO 1                                      SA1 577 
      IF (M-MAXD) 48,1,48                                               SA1 578 
163   WRITE (JTAPE,179)                                                 SA1 579 
      GO TO 1                                                           SA1 580 
C     PRINT COORDINATES OF ADDED VARIABLE                               SA1 581 
164   IND=IND+1                                                         SA1 582 
      DO 165 J=1,MAXD                                                   SA1 583 
165   X(101,J)=X(101,J)+CMEAN(J)                                        SA1 584 
      WRITE (JTAPE,174) IND,(X(101,J),J=1,MAXD)                         SA1 585 
      ISW=1                                                             SA1 586 
      STRESS=SQRT(1.-(1.-STRESS)**2)                                    SA1 587 
      WRITE (JTAPE,178) STRESS,NN                                       SA1 588 
      IF (IFC.GT.0) PUNCH 170, (X(101,J),J=1,MAXD)                      SA1 589 
      GO TO 9                                                           SA1 590 
C     *** FORMAT STATEMENTS ***                                         SA1 591 
C                                                                       SA1 592 
166   FORMAT ('1')                                                      SA1 593 
167   FORMAT (72H                                                       SA1 594 
     1                 )                                                SA1 595 
168   FORMAT (10I4,F4.0)                                                SA1 596 
169   FORMAT (18A4)                                                     SA1 597 
170   FORMAT (10F8.3)                                                   SA1 598 
171   FORMAT ('0GUTTMAN-LINGOES'' SMALLEST SPACE COORDINATES FOR M =',I3SA1 599 
     1,' (SEMI-STRONG MONOTONICITY).'/'0DIMENSION',10I10)               SA1 600 
172   FORMAT ('0KRUSKAL-GUTTMAN-LINGOES'' SMALLEST SPACE COORDINATES FORSA1 601 
     1 M =',I3,' (WEAK MONOTONICITY).'/'0DIMENSION',10I10)              SA1 602 
173   FORMAT (1H ,130(1H-)/9H0VARIABLE)                                 SA1 603 
174   FORMAT (I5,5X2P11F10.3)                                           SA1 604 
175   FORMAT ('0GUTTMAN-LINGOES'' COEFFICIENT OF ALIENATION =',F8.5,' INSA1 605 
     1',I4,' ITERATIONS.'/' KRUSKAL'' STRESS =',F8.5)                   SA1 606 
176   FORMAT ('0KRUSKAL''S STRESS =',F8.5,' IN',I4,' ITERATIONS.'/' GUTTSA1 607 
     1MAN-LINGOES'' COEFFICIENT OF ALIENATION =',F8.5)                  SA1 608 
177   FORMAT ('+',111X,'G-L''S K  ITERATIONS')                          SA1 609 
178   FORMAT ('+',109XF9.5,I8)                                          SA1 610 
179   FORMAT ('1***** YOU HAVE MORE THAN 500 TIED BLOCKS.  EITHER USE EPSA1 611 
     1S PARAMETER TO REDUCE OR CHANGE DIMENSION STATEMENT FOR NOTIES.') SA1 612 
180   FORMAT ('1***** FURTHER ITERATIONS MAY BE NEEDED IN CURRENT DIMENSSA1 613 
     1IONALITY.')                                                       SA1 614 
      END                                                               SA1 615-
C     COR1                                                              CR1   1 
C     SUBROUTINE TO COMPUTE PRODUCT MOMENT CORRELATIONS (IFSR IN SSA    CR1   2 
C     MUST BE SET EQUAL TO 1) OR COVARIANCES (SET IFSR TO -1), WHICH ARECR1   3 
C     NORMALIZED ON LARGEST VARIANCE FOR A SET OF N VARIABLES.  DATA SETCR1   4 
C     -UP: 1) PREPARE A CARD WITH A "1" IN COLUMN 8 IF THERE ARE MISSINGCR1   5 
C     DATA (OTHERWISE LEAVE FIELD BLANK).  IF MISSING DATA, THEN SPECIFYCR1   6 
C     CODE (WITH DECIMAL POINT PUNCHED) IN COLUMNS 9-16.  THE SAME CODE CR1   7 
C     IS USED FOR ALL MISSING OBSERVATIONS; 2) PUNCH N VARIABLES PER OB-CR1   8 
C     SERVATION ACCORDING TO FORMAT FMT ON AS MANY CARDS AS NECESSARY TOCR1   9 
C     SATISFY N, LEAVING COLUMN 1 BLANK;  AND, 3) FOLLOW LAST DATA SET  CR1  10 
C     WITH T TRAILER CARDS (WHERE T=NUMBER OF CARDS/OBSERVATION), HAVINGCR1  11 
C     A "9" IN COLUMN 1.                                                CR1  12 
C                                                                       CR1  13 
      SUBROUTINE COR1 (N,MD,FMT,D,R,IFCOV,FN,SX)                        CR1  14 
      DIMENSION FMT(1), D(1), R(MD,MD), FN(1), SX(1)                    CR1  15 
C                                                                       CR1  16 
C     *** TAPE ASSIGNMENTS:                                             CR1  17 
      KTAPE=3                                                           CR1  18 
      LTAPE=4                                                           CR1  19 
      ITAPE=5                                                           CR1  20 
      JTAPE=6                                                           CR1  21 
      CALL REWIND (KTAPE)                                               CR1  22 
      CALL REWIND (LTAPE)                                               CR1  23 
      READ (ITAPE,36) MISS,CODE                                         CR1  24 
      IF (MISS) 1,12,1                                                  CR1  25 
1     DO 2 J=1,N                                                        CR1  26 
      SX(J)=0.                                                          CR1  27 
2     FN(J)=0.                                                          CR1  28 
      NN=0                                                              CR1  29 
3     READ (ITAPE,FMT) ITYPE,(D(J),J=1,N)                               CR1  30 
      IF (ITYPE-9) 4,7,4                                                CR1  31 
4     NN=NN+1                                                           CR1  32 
      DO 6 J=1,N                                                        CR1  33 
      IF (D(J)-CODE) 5,6,5                                              CR1  34 
5     FN(J)=FN(J)+1.                                                    CR1  35 
      SX(J)=SX(J)+D(J)                                                  CR1  36 
6     CONTINUE                                                          CR1  37 
      WRITE (KTAPE) (D(J),J=1,N)                                        CR1  38 
      GO TO 3                                                           CR1  39 
7     CALL REWIND (KTAPE)                                               CR1  40 
      DO 8 J=1,N                                                        CR1  41 
8     SX(J)=SX(J)/FN(J)                                                 CR1  42 
      WRITE (JTAPE,32) (SX(J),J=1,N)                                    CR1  43 
      WRITE (JTAPE,35) (FN(J),J=1,N)                                    CR1  44 
      ITYPE=0                                                           CR1  45 
      DO 11 I=1,NN                                                      CR1  46 
      READ (KTAPE) (D(J),J=1,N)                                         CR1  47 
      DO 10 J=1,N                                                       CR1  48 
      IF (D(J)-CODE) 10,9,10                                            CR1  49 
9     D(J)=SX(J)                                                        CR1  50 
10    CONTINUE                                                          CR1  51 
11    WRITE (LTAPE) ITYPE,(D(J),J=1,N)                                  CR1  52 
      ITYPE=9                                                           CR1  53 
      WRITE (LTAPE) ITYPE,(D(J),J=1,N)                                  CR1  54 
      CALL REWIND (LTAPE)                                               CR1  55 
12    DO 13 I=1,N                                                       CR1  56 
      SX(I)=0.                                                          CR1  57 
      DO 13 J=1,I                                                       CR1  58 
13    R(I,J)=0.                                                         CR1  59 
      CALL REWIND (KTAPE)                                               CR1  60 
      NN=0                                                              CR1  61 
14    IF (MISS) 15,16,15                                                CR1  62 
15    READ (LTAPE) ITYPE,(D(J),J=1,N)                                   CR1  63 
      GO TO 17                                                          CR1  64 
16    READ (ITAPE,FMT) ITYPE,(D(J),J=1,N)                               CR1  65 
17    IF (ITYPE-9) 18,21,18                                             CR1  66 
18    NN=NN+1                                                           CR1  67 
      DO 19 I=1,N                                                       CR1  68 
      DO 19 J=1,I                                                       CR1  69 
19    R(I,J)=R(I,J)+D(I)*D(J)                                           CR1  70 
      DO 20 J=1,N                                                       CR1  71 
20    SX(J)=SX(J)+D(J)                                                  CR1  72 
      GO TO 14                                                          CR1  73 
21    WRITE (JTAPE,33) N,NN                                             CR1  74 
      FNS=NN                                                            CR1  75 
      DO 22 I=1,N                                                       CR1  76 
      D(I)=SX(I)/FNS                                                    CR1  77 
      FN(I)=SQRT(ABS(R(I,I)/FNS-D(I)**2))                               CR1  78 
22    WRITE (JTAPE,34) I,SX(I),D(I),R(I,I),FN(I)                        CR1  79 
      DO 23 I=1,N                                                       CR1  80 
      DO 23 J=1,I                                                       CR1  81 
23    R(I,J)=(R(I,J)-(SX(I)*SX(J))/FNS)/FNS                             CR1  82 
      NM1=N-1                                                           CR1  83 
      IF (IFCOV) 24,28,28                                               CR1  84 
24    D(1)=0.                                                           CR1  85 
      DO 26 I=1,N                                                       CR1  86 
      IF (D(1)-R(I,I)) 25,26,26                                         CR1  87 
25    D(1)=R(I,I)                                                       CR1  88 
26    CONTINUE                                                          CR1  89 
      DO 27 I=1,NM1                                                     CR1  90 
      IP1=I+1                                                           CR1  91 
      DO 27 J=IP1,N                                                     CR1  92 
27    R(I,J)=R(J,I)/D(1)                                                CR1  93 
      GO TO 31                                                          CR1  94 
28    DO 29 I=1,N                                                       CR1  95 
      DO 29 J=1,I                                                       CR1  96 
29    R(I,J)=R(I,J)/(FN(I)*FN(J))                                       CR1  97 
      DO 30 I=1,NM1                                                     CR1  98 
      IP1=I+1                                                           CR1  99 
      DO 30 J=IP1,N                                                     CR1 100 
30    R(I,J)=R(J,I)                                                     CR1 101 
31    RETURN                                                            CR1 102 
C     *** FORMAT STATEMENTS ***                                         CR1 103 
C                                                                       CR1 104 
32    FORMAT (31H0MISSING DATA MEANS BY VARIABLE/(1P7E18.7))            CR1 105 
33    FORMAT (1H /'0',16X19HNO. OF VARIABLES = I5/'0',16X22HNO. OF OBSERCR1 106 
     1VATIONS = I5/13H0    VAR. NO.,14X3HSUM,21X4HMEAN,18X11HSUM SQUARESCR1 107 
     2,15X9HSTD. DEV.)                                                  CR1 108 
34    FORMAT (1H0,I9,1P4E25.7)                                          CR1 109 
35    FORMAT (38H0MISSING DATA SAMPLE SIZES BY VARIABLE/(7F18.0))       CR1 110 
36    FORMAT (I8,F8.0)                                                  CR1 111 
      END                                                               CR1 112-
C     COR1                                                              CR1   1 
C     SUBROUTINE TO GENERATE MATRIX OF NORMALIZED AGREEMENT SCORES.     CR1   2 
C     PUNCH DATA WITH FMT, I.E., (I1,...), LEAVING COL. 1 BLANK ON DATA CR1   3 
C     CARD/S (HAVING NR VALUES).  A TRAILER SET FOLLOWS LAST DATA SET   CR1   4 
C     WITH A "9" IN COL. 1.  PROGRAMMED BY LINGOES (1/1/69).            CR1   5 
C                                                                       CR1   6 
      SUBROUTINE COR1 (NR,MD,FMT,REC,C,IFSR,A,B)                        CR1   7 
      DIMENSION FMT(1), REC(1), C(MD,MD), A(1), B(1)                    CR1   8 
C                                                                       CR1   9 
C     *** TAPE ASSIGNMENT:                                              CR1  10 
      ITAPE=5                                                           CR1  11 
C     INITIALIZE                                                        CR1  12 
      NRM1=NR-1                                                         CR1  13 
      JJ=0                                                              CR1  14 
      DO 1 I=2,NR                                                       CR1  15 
      IM1=I-1                                                           CR1  16 
      DO 1 J=1,IM1                                                      CR1  17 
1     C(I,J)=0.                                                         CR1  18 
C     INPUT AN OBSERVATION                                              CR1  19 
2     READ (ITAPE,FMT) I,(REC(J),J=1,NR)                                CR1  20 
      IF (I.EQ.9) GO TO 4                                               CR1  21 
C     ACCUMULATE AGREEMENT SCORES                                       CR1  22 
      JJ=JJ+1                                                           CR1  23 
      DO 3 I=1,NRM1                                                     CR1  24 
      IM1=I+1                                                           CR1  25 
      DO 3 J=IM1,NR                                                     CR1  26 
      IF (REC(I).EQ.REC(J)) C(J,I)=C(J,I)+1.                            CR1  27 
3     CONTINUE                                                          CR1  28 
      GO TO 2                                                           CR1  29 
4     FN=JJ                                                             CR1  30 
C     NORMALIZE AGREEMENT SCORES                                        CR1  31 
      DO 5 I=1,NRM1                                                     CR1  32 
      IM1=I+1                                                           CR1  33 
      DO 5 J=IM1,NR                                                     CR1  34 
5     C(J,I)=C(J,I)/FN                                                  CR1  35 
      RETURN                                                            CR1  36 
      END                                                               CR1  37-
C     COR2                                                              CR2   1 
C     SUBROUTINE TO GENERATE PRODUCT MOMENT CORRELATIONS BETWEEN AN     CR2   2 
C     ADDED VARIABLE/SUBJECT AND NR FIXED VARIABLES/SUBJECTS.  IFSR IN  CR2   3 
C     SSA-I MUST BE SET EQUAL TO THE NUMBER OF OBSERVATIONS, I.E., NS   CR2   4 
C     .LE. 100 AND NR MUST BE .LE. 97.  DATA SET-UP: 1) NS SETS OF CARDSCR2   5 
C     HAVING NR ELEMENTS PER SET REPRESENTING SCORE MATRIX FOR FIXED    CR2   6 
C     VARIABLES; 2) FORMAT CARD DESCRIBING OBSERVATIONS OF ADDED VARI-  CR2   7 
C     ABLE; AND, 3) SETS OF ADDED VARIABLES HAVING NS ELEMENTS REPRE-   CR2   8 
C     SENTING OBSERVATIONS.  PROGRAMMED BY LINGOES (1/1/69).            CR2   9 
C                                                                       CR2  10 
      SUBROUTINE COR2 (NR,NS,ISW,MD,FMT,PROX,C)                         CR2  11 
      DIMENSION PROX(1), FMT(1), C(MD,MD)                               CR2  12 
C                                                                       CR2  13 
C     *** TAPE ASSIGNMENTS:                                             CR2  14 
      ITAPE=5                                                           CR2  15 
      KTAPE=3                                                           CR2  16 
      CALL REWIND (KTAPE)                                               CR2  17 
      IF (ISW) 6,1,6                                                    CR2  18 
1     DO 2 J=1,NS                                                       CR2  19 
2     READ (ITAPE,FMT) (C(I,J),I=1,NR)                                  CR2  20 
      NRP1=NR+1                                                         CR2  21 
      NRP2=NRP1+1                                                       CR2  22 
      NRP3=NRP2+1                                                       CR2  23 
      FNS=NS                                                            CR2  24 
      DO 4 I=1,NR                                                       CR2  25 
      C(NRP1,I)=0.                                                      CR2  26 
      C(NRP2,I)=0.                                                      CR2  27 
      DO 3 J=1,NS                                                       CR2  28 
      C(NRP1,I)=C(NRP1,I)+C(I,J)                                        CR2  29 
3     C(NRP2,I)=C(NRP2,I)+C(I,J)**2                                     CR2  30 
4     C(NRP2,I)=SQRT(ABS(FNS*C(NRP2,I)-C(NRP1,I)**2))                   CR2  31 
      DO 5 I=1,NR                                                       CR2  32 
5     WRITE (KTAPE) (C(I,J),J=1,NS)                                     CR2  33 
      WRITE (KTAPE) (C(NRP1,I),I=1,NR)                                  CR2  34 
      WRITE (KTAPE) (C(NRP2,I),I=1,NR)                                  CR2  35 
      READ (ITAPE,11) (FMT(J),J=1,18)                                   CR2  36 
      GO TO 8                                                           CR2  37 
6     DO 7 I=1,NR                                                       CR2  38 
7     READ (KTAPE) (C(I,J),J=1,NS)                                      CR2  39 
      READ (KTAPE) (C(NRP1,I),I=1,NR)                                   CR2  40 
      READ (KTAPE) (C(NRP2,I),I=1,NR)                                   CR2  41 
8     READ (ITAPE,FMT) (C(NRP3,J),J=1,NS)                               CR2  42 
      DO 10 I=1,NR                                                      CR2  43 
      SY=0.                                                             CR2  44 
      SYSQ=0.                                                           CR2  45 
      PROX(I)=0.                                                        CR2  46 
      DO 9 J=1,NS                                                       CR2  47 
      SY=SY+C(NRP3,J)                                                   CR2  48 
      SYSQ=SYSQ+C(NRP3,J)**2                                            CR2  49 
9     PROX(I)=PROX(I)+C(I,J)*C(NRP3,J)                                  CR2  50 
10    PROX(I)=(FNS*PROX(I)-C(NRP1,I)*SY)/(C(NRP2,I)*SQRT(ABS(FNS*SYSQ-SYCR2  51 
     1**2)))                                                            CR2  52 
      RETURN                                                            CR2  53 
C     *** FORMAT STATEMENT ***                                          CR2  54 
C                                                                       CR2  55 
11    FORMAT (18A4)                                                     CR2  56 
      END                                                               CR2  57-
C     FIT                                                               FIT   1 
C     COMPUTE MONOTONE REGRESSION VALUES ACCORDING TO KRUSKAL'S ALGO-   FIT   2 
C     RITHM.  ADAPTED FROM A VERSION OF M-D-SCAL WRITTEN BY ROSKAM, BOR-FIT   3 
C     GERS, AND HORSTEN OF THE UNIVERSITY OF NIJMEGEN - 11/1/68.  J.C.L.FIT   4 
C                                                                       FIT   5 
      SUBROUTINE FIT (NEL,DIST,INDI,DHAT,SUMD,IBLK)                     FIT   6 
      DIMENSION DIST(1), INDI(1), DHAT(1), SUMD(1), IBLK(1)             FIT   7 
C                                                                       FIT   8 
C     PERMUTE D TO ORDER OF P                                           FIT   9 
      DO 1 J=1,NEL                                                      FIT  10 
1     DHAT(J)=DIST(INDI(J))                                             FIT  11 
      JJ=0                                                              FIT  12 
      NB=0                                                              FIT  13 
2     II=JJ+1                                                           FIT  14 
      JJ=II                                                             FIT  15 
      IF (II.GT.NEL) GO TO 13                                           FIT  16 
      NB=NB+1                                                           FIT  17 
      SUMD(NB)=DHAT(II)                                                 FIT  18 
      IBLK(NB)=1                                                        FIT  19 
      IUP=0                                                             FIT  20 
      IDN=0                                                             FIT  21 
3     IF (JJ-NEL) 5,4,5                                                 FIT  22 
4     IUP=1                                                             FIT  23 
      GO TO 11                                                          FIT  24 
5     IF (DHAT(JJ+1).GT.DHAT(JJ)) GO TO 4                               FIT  25 
      JJ=JJ+1                                                           FIT  26 
      IBLK(NB)=IBLK(NB)+1                                               FIT  27 
      SUMD(NB)=SUMD(NB)+DHAT(JJ)                                        FIT  28 
      AVER=SUMD(NB)/FLOAT(IBLK(NB))                                     FIT  29 
      DO 6 J=II,JJ                                                      FIT  30 
6     DHAT(J)=AVER                                                      FIT  31 
      IDN=0                                                             FIT  32 
7     IF (II-1) 9,8,9                                                   FIT  33 
8     IDN=1                                                             FIT  34 
      GO TO 11                                                          FIT  35 
9     IF (DHAT(II).GT.DHAT(II-1)) GO TO 8                               FIT  36 
      II=II-IBLK(NB-1)                                                  FIT  37 
      SUMD(NB-1)=SUMD(NB-1)+SUMD(NB)                                    FIT  38 
      IBLK(NB-1)=IBLK(NB-1)+IBLK(NB)                                    FIT  39 
      NB=NB-1                                                           FIT  40 
      AVER=SUMD(NB)/FLOAT(IBLK(NB))                                     FIT  41 
      DO 10 J=II,JJ                                                     FIT  42 
10    DHAT(J)=AVER                                                      FIT  43 
      IUP=0                                                             FIT  44 
      GO TO 3                                                           FIT  45 
11    IF (IUP*IDN) 2,12,2                                               FIT  46 
12    IF (IUP) 7,3,7                                                    FIT  47 
C     PERMUTE D-HAT TO ORDER OF P                                       FIT  48 
13    DO 14 J=1,NEL                                                     FIT  49 
14    SUMD(INDI(J))=DHAT(J)                                             FIT  50 
      RETURN                                                            FIT  51 
      END                                                               FIT  52-
C     SORT                                                              SRT   1 
C     THIS SUBROUTINE IS A N-LOG-N LIST SORT USING SHELL'S ALGORITHM.   SRT   2 
C     NUMBERS ARE ORDERED LOW TO HIGH (ISIM=0) OR HIGH TO LOW (ISIM=1)  SRT   3 
C     WITHOUT DISTURBING THE NUMBERS TO BE ORDERED.  PROGRAMMED BY      SRT   4 
C     LINGOES - 11/1/68, U. OF M. COMPUTING CENTER.                     SRT   5 
C                                                                       SRT   6 
      SUBROUTINE SORT (ISW,NEL,DIST,ISIM,INDEX)                         SRT   7 
      DIMENSION DIST(1), INDEX(1)                                       SRT   8 
C                                                                       SRT   9 
      IF (ISW) 3,1,3                                                    SRT  10 
1     DO 2 J=1,NEL                                                      SRT  11 
2     INDEX(J)=J                                                        SRT  12 
3     M=NEL                                                             SRT  13 
4     M=M/2                                                             SRT  14 
      IF (M) 5,12,5                                                     SRT  15 
5     KK=NEL-M                                                          SRT  16 
      J=1                                                               SRT  17 
6     I=J                                                               SRT  18 
7     IPM=I+M                                                           SRT  19 
      IF (ISIM) 9,8,9                                                   SRT  20 
8     IF (DIST(INDEX(I)).GT.DIST(INDEX(IPM))) GO TO 11                  SRT  21 
      GO TO 10                                                          SRT  22 
9     IF (DIST(INDEX(I)).LT.DIST(INDEX(IPM))) GO TO 11                  SRT  23 
10    J=J+1                                                             SRT  24 
      IF (J.GT.KK) GO TO 4                                              SRT  25 
      GO TO 6                                                           SRT  26 
11    LL=INDEX(I)                                                       SRT  27 
      INDEX(I)=INDEX(IPM)                                               SRT  28 
      INDEX(IPM)=LL                                                     SRT  29 
      I=I-M                                                             SRT  30 
      IF (I.LT.1) GO TO 10                                              SRT  31 
      GO TO 7                                                           SRT  32 
12    RETURN                                                            SRT  33 
      END                                                               SRT  34-
C     EIGEN                                                             EGN   1 
C     EIGENVALUES AND NORMALIZED EIGENVECTORS OF A REAL SYMMETRIC MATRIXEGN   2 
C     PROGRAMMED BY GARBOW, ARGONNE, 1965 AND MODIFIED BY LINGOES, U OF EGN   3 
C     M, 1966 USING HOUSEHOLDER'S TRIDIAGONALIZATION PROCEDURE AND      EGN   4 
C     INVERSE ITERATIONS TO OBTAIN EIGENVECTORS.  COMPLETE MULTIPLICITY EGN   5 
C     OF EIGENSYSTEM IS DETERMINED.  IF NZ=1 VECTORS ARE NORMALIZED.    EGN   6 
C     EIGENVALUES ARE RETURNED IN VALU AND NORMALIZED EIGENVECTORS ARE  EGN   7 
C     STORED IN B.  NSUB IS ORDER OF MATRICES A AND B AND MSUB IS THE   EGN   8 
C     NUMBER OF ROOTS AND VECTORS DESIRED.                              EGN   9 
C                                                                       EGN  10 
      SUBROUTINE EIGEN (A,B,NSUB,VALU,MSUB,MD,T,DIAG,SUPERD,U,INDEX,V,NZEGN  11 
     1)                                                                 EGN  12 
C                                                                       EGN  13 
      DIMENSION A(MD,MD), B(MD,MD), VALU(1), T(MD,3), DIAG(1), SUPERD(1)EGN  14 
     1, U(1), INDEX(1), V(1)                                            EGN  15 
      EQUIVALENCE (I1,T1), (I2,T2), (TEMP,T0), (SUM,MATCH), (I,P), (DIV,EGN  16 
     1SCALAR,TAU), (ANORM2,NORM), (VTEMP,VNORM2,VNORM)                  EGN  17 
C                                                                       EGN  18 
C     INITIALIZATION                                                    EGN  19 
      N=NSUB                                                            EGN  20 
      M=MSUB                                                            EGN  21 
      NP1=N+1                                                           EGN  22 
      NM1=N-1                                                           EGN  23 
      E1=1.E-8                                                          EGN  24 
C     GENERATE IDENTITY MATRIX                                          EGN  25 
      DO 3 I=1,N                                                        EGN  26 
      DO 3 J=1,N                                                        EGN  27 
      IF (I-J) 2,1,2                                                    EGN  28 
1     B(I,J)=1.                                                         EGN  29 
      GO TO 3                                                           EGN  30 
2     B(I,J)=0.                                                         EGN  31 
3     CONTINUE                                                          EGN  32 
C     HOUSEHOLDER SIMILARITY TRANSFORMATION TO CO-DIAGONAL FORM         EGN  33 
C     REDUCE COLUMN OF MATRIX                                           EGN  34 
      DO 14 I=1,NM1                                                     EGN  35 
      IF (I-NM1) 4,13,4                                                 EGN  36 
4     I1=I+1                                                            EGN  37 
      I2=I1+1                                                           EGN  38 
      SUM=0.                                                            EGN  39 
      DO 5 J=I2,N                                                       EGN  40 
5     SUM=SUM+A(J,I)**2                                                 EGN  41 
      IF (SUM) 6,13,6                                                   EGN  42 
6     J=I1                                                              EGN  43 
      TEMP=A(J,I)                                                       EGN  44 
      SUM=SQRT(SUM+TEMP**2)                                             EGN  45 
      A(J,I)=-SIGN(SUM,TEMP)                                            EGN  46 
      U(J)=SQRT(1.+ABS(TEMP)/SUM)                                       EGN  47 
      DIV=SIGN(U(J)*SUM,TEMP)                                           EGN  48 
      DO 7 J=I2,N                                                       EGN  49 
7     U(J)=A(J,I)/DIV                                                   EGN  50 
      SCALAR=0.                                                         EGN  51 
      DO 9 J=I1,N                                                       EGN  52 
      V(J)=0.                                                           EGN  53 
      DO 8 K=I1,N                                                       EGN  54 
8     V(J)=V(J)+A(K,J)*U(K)                                             EGN  55 
      SCALAR=SCALAR+V(J)*U(J)                                           EGN  56 
9     CONTINUE                                                          EGN  57 
      SCALAR=SCALAR/2.                                                  EGN  58 
      DO 10 J=I1,N                                                      EGN  59 
      V(J)=V(J)-SCALAR*U(J)                                             EGN  60 
      DO 10 K=I1,J                                                      EGN  61 
      A(K,J)=A(K,J)-(U(K)*V(J)+U(J)*V(K))                               EGN  62 
      A(J,K)=A(K,J)                                                     EGN  63 
10    CONTINUE                                                          EGN  64 
C     SAVE ROTATION FOR LATER APPLICATION TO CO-DIAGONAL VECTORS        EGN  65 
      DO 12 K=2,N                                                       EGN  66 
      TEMP=0.                                                           EGN  67 
      DO 11 J=I1,N                                                      EGN  68 
11    TEMP=TEMP+U(J)*B(J,K)                                             EGN  69 
      DO 12 J=I1,N                                                      EGN  70 
      B(J,K)=B(J,K)-U(J)*TEMP                                           EGN  71 
12    CONTINUE                                                          EGN  72 
C     MOVE CO-DIAGONAL FORM ELEMENTS FOR ITERATIVE PROCEDURE            EGN  73 
13    J=I                                                               EGN  74 
      DIAG(I)=A(J,I)                                                    EGN  75 
      SUPERD(I)=A(J+1,I)                                                EGN  76 
14    CONTINUE                                                          EGN  77 
      DIAG(N)=A(N,N)                                                    EGN  78 
C     DETERMINE EIGENVALUES FROM STURM CHAIN OF CO-DIAGONAL MINORS      EGN  79 
C     CALCULATE NORM OF MATRIX AND INITIALIZE EIGENVALUE BOUNDS         EGN  80 
      ANORM2=DIAG(1)**2                                                 EGN  81 
      DO 15 L=2,N                                                       EGN  82 
      V(L-1)=SUPERD(L-1)**2                                             EGN  83 
      ANORM2=DIAG(L)**2+V(L-1)+V(L-1)+ANORM2                            EGN  84 
15    CONTINUE                                                          EGN  85 
      ANORM=SQRT(ANORM2)                                                EGN  86 
      DO 16 L=1,M                                                       EGN  87 
      VALU(L)=ANORM                                                     EGN  88 
      U(L)=-ANORM                                                       EGN  89 
16    CONTINUE                                                          EGN  90 
      EPS1=ANORM*E1                                                     EGN  91 
      IF (EPS1) 17,75,17                                                EGN  92 
C     CHOOSE NEW TRIAL VALUE WHILE TESTING BOUNDS FOR CONVERGENCE       EGN  93 
17    DO 35 L=1,M                                                       EGN  94 
      ITER=0                                                            EGN  95 
      VTEMP=EPS1                                                        EGN  96 
18    TAU=(VALU(L)+U(L))/2.                                             EGN  97 
      IF (ITER-10) 20,19,20                                             EGN  98 
19    VTEMP=VTEMP*10.                                                   EGN  99 
      ITER=0                                                            EGN 100 
20    IF (2.*(TAU-U(L))-VTEMP) 35,35,21                                 EGN 101 
C     DETERMINE SIGNS OF PRINCIPAL MINORS                               EGN 102 
21    MATCH=0                                                           EGN 103 
      ITER=ITER+1                                                       EGN 104 
      T2=0.                                                             EGN 105 
      T1=1.                                                             EGN 106 
      DO 30 L1=1,N                                                      EGN 107 
      P=DIAG(L1)-TAU                                                    EGN 108 
      IF (T2) 23,22,23                                                  EGN 109 
22    T1=SIGN(1.,T1)                                                    EGN 110 
23    IF (T1) 25,24,25                                                  EGN 111 
24    T0=-SIGN(1.,T2)                                                   EGN 112 
      T2=0.                                                             EGN 113 
      IF (V(L1-1)) 26,22,26                                             EGN 114 
25    T0=P-V(L1-1)*T2/T1                                                EGN 115 
      T2=1.                                                             EGN 116 
C     COUNT AGREEMENTS IN SIGN (ZERO CONSIDERED POSITIVE)               EGN 117 
26    IF (T0) 29,27,28                                                  EGN 118 
27    T2=T1                                                             EGN 119 
      IF (T2) 29,28,28                                                  EGN 120 
28    MATCH=MATCH+1                                                     EGN 121 
29    T1=T0                                                             EGN 122 
30    CONTINUE                                                          EGN 123 
C     ESTABLISH TIGHTER BOUNDS ON EIGENVALUES                           EGN 124 
      DO 34 L1=L,M                                                      EGN 125 
      IF (L1-MATCH) 33,33,31                                            EGN 126 
31    IF (VALU(L1)-TAU) 18,18,32                                        EGN 127 
32    VALU(L1)=TAU                                                      EGN 128 
      GO TO 34                                                          EGN 129 
33    U(L1)=TAU                                                         EGN 130 
34    CONTINUE                                                          EGN 131 
      GO TO 18                                                          EGN 132 
35    CONTINUE                                                          EGN 133 
C     EIGENVECTORS OF CO-DIAGONAL SYMMETRIC MATRIX -- INVERSE ITERATION EGN 134 
C     CHECK FOR REPEATED VALUE                                          EGN 135 
      DO 68 I=1,M                                                       EGN 136 
      IF (I-2) 37,36,36                                                 EGN 137 
36    IF (VALU(I-1)-VALU(I)-1.E-3) 38,37,37                             EGN 138 
37    I1=-1                                                             EGN 139 
38    I1=I1+1                                                           EGN 140 
C     TRIANGULARIZE CO-DIAGONAL FORM AFTER EIGENVALUE SUBTRACTION       EGN 141 
      DO 43 L=1,N                                                       EGN 142 
      V(L)=EPS1                                                         EGN 143 
      T(L,2)=DIAG(L)-VALU(I)                                            EGN 144 
      IF (L-N) 40,39,40                                                 EGN 145 
39    T(L,3)=0.                                                         EGN 146 
      GO TO 43                                                          EGN 147 
40    T(L,3)=SUPERD(L)                                                  EGN 148 
      IF (T(L,3)) 42,41,42                                              EGN 149 
41    T(L,3)=EPS1                                                       EGN 150 
42    T(L+1,1)=T(L,3)                                                   EGN 151 
43    CONTINUE                                                          EGN 152 
      DO 50 J=1,N                                                       EGN 153 
      T(J,1)=T(J,2)                                                     EGN 154 
      T(J,2)=T(J,3)                                                     EGN 155 
      T(J,3)=0.                                                         EGN 156 
      VTEMP=ABS(T(J,1))                                                 EGN 157 
      IF (J-N) 46,44,46                                                 EGN 158 
44    IF (VTEMP) 50,45,50                                               EGN 159 
45    T(J,1)=EPS1                                                       EGN 160 
      GO TO 50                                                          EGN 161 
46    INDEX(J)=0                                                        EGN 162 
      IF (ABS(T(J+1,1))-VTEMP) 49,49,47                                 EGN 163 
47    INDEX(J)=1                                                        EGN 164 
      DO 48 K=1,3                                                       EGN 165 
      VTEMP=T(J,K)                                                      EGN 166 
      T(J,K)=T(J+1,K)                                                   EGN 167 
      T(J+1,K)=VTEMP                                                    EGN 168 
48    CONTINUE                                                          EGN 169 
49    VTEMP=T(J+1,1)/T(J,1)                                             EGN 170 
      U(J)=VTEMP                                                        EGN 171 
      T(J+1,2)=T(J+1,2)-VTEMP*T(J,2)                                    EGN 172 
      T(J+1,3)=T(J+1,3)-VTEMP*T(J,3)                                    EGN 173 
50    CONTINUE                                                          EGN 174 
      ITER=1                                                            EGN 175 
      IF (I1) 58,51,58                                                  EGN 176 
C     BACK SUBSTITUTE TO OBTAIN EIGENVECTOR                             EGN 177 
51    DO 52 L1=1,N                                                      EGN 178 
      L=NP1-L1                                                          EGN 179 
      V(L)=(V(L)-T(L,2)*V(L+1)-T(L,3)*V(L+2))/T(L,1)                    EGN 180 
52    CONTINUE                                                          EGN 181 
      GO TO (53,58), ITER                                               EGN 182 
C     PERFORM SECOND ITERATION                                          EGN 183 
53    ITER=2                                                            EGN 184 
54    DO 57 L=2,N                                                       EGN 185 
      IF (INDEX(L-1)) 55,56,55                                          EGN 186 
55    VTEMP=V(L-1)                                                      EGN 187 
      V(L-1)=V(L)                                                       EGN 188 
      V(L)=VTEMP                                                        EGN 189 
56    V(L)=V(L)-U(L-1)*V(L-1)                                           EGN 190 
57    CONTINUE                                                          EGN 191 
      GO TO 51                                                          EGN 192 
C     ORTHOGONALIZE VECTOR TO OTHERS ASSOCIATED WITH REPEATED ROOT      EGN 193 
58    IF (I1) 59,62,59                                                  EGN 194 
59    DO 61 L1=1,I1                                                     EGN 195 
      K=I-L1                                                            EGN 196 
      VTEMP=0.                                                          EGN 197 
      DO 60 J=1,N                                                       EGN 198 
60    VTEMP=VTEMP+A(J,K)*V(J)                                           EGN 199 
      DO 61 J=1,N                                                       EGN 200 
61    V(J)=V(J)-A(J,K)*VTEMP                                            EGN 201 
62    GO TO (54,63), ITER                                               EGN 202 
C     NORMALIZE VECTOR TO UNIT LENGTH                                   EGN 203 
63    VNORM2=0.                                                         EGN 204 
      SUM=0.                                                            EGN 205 
      DO 65 L=1,N                                                       EGN 206 
      IF (SUM-ABS(V(L))) 64,65,65                                       EGN 207 
64    SUM=ABS(V(L))                                                     EGN 208 
65    CONTINUE                                                          EGN 209 
      DO 66 L=1,N                                                       EGN 210 
      V(L)=V(L)/SUM                                                     EGN 211 
66    VNORM2=VNORM2+V(L)**2                                             EGN 212 
      VNORM=SQRT(VNORM2)                                                EGN 213 
      DO 67 J=1,N                                                       EGN 214 
67    A(J,I)=V(J)/VNORM                                                 EGN 215 
68    CONTINUE                                                          EGN 216 
C     ROTATION OF CO-DIAGONAL VECTORS INTO MATRIX EIGENVECTORS          EGN 217 
      DO 70 I=1,M                                                       EGN 218 
      DO 69 K=2,N                                                       EGN 219 
      U(K)=0.                                                           EGN 220 
      DO 69 J=2,N                                                       EGN 221 
69    U(K)=U(K)+B(J,K)*A(J,I)                                           EGN 222 
      DO 70 J=2,N                                                       EGN 223 
70    A(J,I)=U(J)                                                       EGN 224 
      IF (NZ.EQ.0) GO TO 73                                             EGN 225 
C     NORMALIZE LENGTH OF VECTORS TO EIGENVALUES AND STORE IN B(I,J)    EGN 226 
      DO 72 J=1,M                                                       EGN 227 
      IF (VALU(J)) 75,75,71                                             EGN 228 
71    VTEMP=SQRT(VALU(J))                                               EGN 229 
      DO 72 I=1,N                                                       EGN 230 
72    B(I,J)=A(I,J)*VTEMP                                               EGN 231 
      GO TO 75                                                          EGN 232 
73    DO 74 J=1,M                                                       EGN 233 
      DO 74 I=1,N                                                       EGN 234 
74    B(I,J)=A(I,J)                                                     EGN 235 
75    RETURN                                                            EGN 236 
      END                                                               EGN 237-
1G-L(SSA-I) OF 13 EXPERIMENTAL (UNREFLECTED) MMPI SCALES - DAMARIN.             
  13   2   2   1   1                                                            
(12F3.2)                                                                        
-54                                                                             
 68-43                                                                          
-36 55-49                                                                       
-01 24-14 42                                                                    
 45-18 64-30-10                                                                 
-15 27-24 36 25-18                                                              
 50-26 72-31-12 69-19                                                           
 19 08 45-01-07 48 06 47                                                        
 59-24 55-08 21 37 11 44 20                                                     
-25 37-21 44 13 00 27-06 29-11                                                  
 69-38 84-36-07 61-16 67 45 64-15                                               
 39 00 55 04 10 50 12 55 51 49 21 57                                            
