C kyst.f
C The authors of this software are Joseph B Kruskal, Forrest Young,
C and Judith B Seery.

C Copyright (c) 1993 by AT&T.
C Permission to use, copy, modify, and distribute this software for any
C purpose without fee is hereby granted, provided that this entire
C notice is included in all copies of any software which is or includes
C a copy or modification of this software and in all copies of the
C supporting documentation for such software.
C THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMP-
C LIED WARRANTY.  IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
C REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
C OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.

C This software comes from the FIRST MDS Package of AT&T Bell Laboratories

C The manual of how to use kyst2a is available in several different formats:
C see mds/kyst2a_manual.  There are few differences between kyst and kyst2a.

C For explanation of the method this software implements, see
C "Multidimensional Scaling" (1978) by Joseph B. Kruskal and Myron Wish,
C Sage Publications: Beverly Hills, Calif.,

C "Multidimensional Scaling and Other Methods for Discovering Structure"
C (1977) by Joseph B. Kruskal, pp 296-339 in
C "Statistical Methods for Digital Computers" by Kurt Enslein,
C Anthony Ralston, and Herbert S. Wilf, Wiley: New York,

C "Multidimensional Scaling by Optimizing Goodness of Fit to a Nonmetric
C Hypothesis" (1964) by Joseph B. Kruskal in Psychometrika 29, pp 1-27,

C "Nonmetric Multidimensional Scaling: A Numerical Method" (1964) by
C Joseph B. Kruskal in Psychometrika 29(2) pp 115-129.
C----+----@----+----@----+----@----+----@----+----@----+----@----+----@

$      FORTREX STAB
CMAIN                                                                   MAIN   1
C *********************     KYST             JANUARY,1973      *********MAIN   2
C                                                                       MAIN   3
C  KRUSKAL-YOUNG-SHEPARD-TORGERSON MULTIDIMENSIONAL SCALING PROGRAM     MAIN   4
C                                                                       MAIN   5
C                                                                       MAIN   6
C        WRITTEN BY                                                     MAIN   7
C                                                                       MAIN   8
C               DR. J. B. KRUSKAL                                       MAIN   9
C               BELL TELEPHONE LABORATORIES                             MAIN  10
C               MURRAY HILL, N. J.                                      MAIN  11
C                                                                       MAIN  12
C               DR. F. W. YOUNG                                         MAIN  13
C               PSYCHOMETRIC LABORATORY                                 MAIN  14
C               UNIVERSITY OF NORTH CAROLINA                            MAIN  15
C               CHAPEL HILL, N. C.                                      MAIN  16
C        ASSISTED BY                                                    MAIN  17
C               JUDITH B. SEERY                                         MAIN  18
C               BELL TELEPHONE LABORATORIES                             MAIN  19
C               MURRAY HILL, N. J.                                      MAIN  20
C                                                                       MAIN  21
C                                                                       MAIN  22
C                                                                       MAIN  23
C **********************************************************************MAIN  24
C                                                                       MAIN  25
C GENERAL REMARKS.                                                      MAIN  26
C                                                                       MAIN  27
C                                                                       MAIN  28
C KYST INCLUDES THE FOLLOWING SUBROUTINES                               MAIN  29
C                                                                       MAIN  30
C     BLDA             EQSOLV          PLOT           SERCH             MAIN  31
C     BSEC1            FITM            REGR           SGEV              MAIN  32
C     CCACT            FITMS           RM1POW         SORT              MAIN  33
C     CONFIG           FITP            RPOWER         ST03              MAIN  34
C     DATAPR           INICON          RROOT          WTRAN             MAIN  35
C     DFLT             NEWSTP          RUNIFV         XITEM             MAIN  36
C     DTRAN            NRMLZ           SBK            NVIT1             MAIN  37
C                                                                       MAIN  38
C ALL ROUTINES ARE WRITTEN IN FORTRAN IV, AND PUNCHED IN THE IBMEL 029  MAIN  39
C     CHARACTER SET                                                     MAIN  40
C                                                                       MAIN  41
C NO USE IS MADE OF SPECIAL OR NON-STANDARD SOFTWARE.                   MAIN  42
C                                                                       MAIN  43
C ALL INPUT AND OUTPUT IS ONTO FORTRAN LOGICAL UNITS WITH NUMBERS       MAIN  44
C     CONTROLLED BY THESE VARIABLES  LREAD,LPRINT,LPUNCH,LSCRAT.        MAIN  45
C     UNIT NUMBERS 5,6,7,8  HAVE BEEN USED RESPECTIVELY. TO CHANGE THESEMAIN  46
C     ASSIGNMENTS, CHANGE THE VALUES FOR THE FOUR VARIABLES SET IN THE  MAIN  47
C     BLOCK DATA SUBPROGRAM.                                            MAIN  48
C                                                                       MAIN  49
C THE SCRATCH UNIT IS USED IN A MINOR WAY BY CCACT ONLY.  MANY          MAIN  50
C     INSTALLATIONS WILL HAVE ALTERNATE METHODS OF DOING THE SAME THING.MAIN  51
C                                                                       MAIN  52
C KYST USES THREE MACHINE DEPENDENT CONSTANTS  PRECSN (RELATIVE MACHINE MAIN  53
C     PRECISION), XMAG (LEAST POSITIVE MACHINE NUMBER), XMAXN (MINIMUM  MAIN  54
C     OF LARGEST POSITIVE MACHINE NUMBER AND MINUS LARGEST NEGATIVE     MAIN  55
C     MACHINE NUMBER).  VALUES 1.5E-8, 1.0E-38, 1.0E38 HAVE BEEN USED   MAIN  56
C     RESPECTIVELY.  TO CHANGE THESE ASSIGNMENTS, CHANGE THE VALUES     MAIN  57
C     FOR THESE VARIABLES SET IN THE BLOCK DATA SUBPROGRAM.             MAIN  58
C                                                                       MAIN  59
C **********************************************************************MAIN  60
C                                                                       MAIN  61
C                                                                       MAIN  62
C MAIN PROGRAM     FOR KYST                  JANUARY,1973               MAIN  63
C WRITTEN BY J.KRUSKAL                                                  MAIN  64
C MODIFIED FOR KYST BY J.KRUSKAL AND J.SEERY JANUARY,1973               MAIN  65
       DIMENSION   DATA(1800), IJ(1800), DIST(1800), DHAT(1800)         MAIN  66
      DIMENSION                WW(1800),LDIST(1800)                     MAIN  67
       DIMENSION    X(100,6), GR(100,6),  GL(100,6)                     MAIN  68
      DIMENSION CTITLE(80),FMAT(80),TITLE(80),LBLOCK(1800)              MAIN  69
       INTEGER    GRNO(101), NOGRPS, LSPLIT,         SPLITH             MAIN  70
       INTEGER GRSDIS(100), SDSWIT,SDSWT1                               MAIN  71
       REAL    GRSTRS(100),PCOEFF(5)                                    MAIN  72
      REAL ITEM(101),PTID(100),PMAT(1800,2),RVEC(1800),XNUM(6)          MAIN  73
      DIMENSION  STPL(10),DIMSV(10)                                     MAIN  74
      INTEGER DUMMY(25),RTYPE                                           MAIN  75
C                                                                       MAIN  76
C                                                                       MAIN  77
      COMMON /ACCUR/ PRECSN,XMAG,XMAXN                                  MAIN  78
      COMMON /KYST1/ DATA,WW,IJ,X                                       MAIN  79
      COMMON /KYST2/ GR,GL,RVEC,PCOEFF,DUMMY,PMAT                       MAIN  80
       COMMON  /MDSCL1/  LREAD, LPRINT, LPUNCH, LSCRAT                  MAIN  81
       COMMON  /MDSCL2/                                                 MAIN  82
     1         LDIMX, LDIMN, LDIMD, CUTOFF, STRMIN                      MAIN  83
     2 ,       SFGRMN, COSAVW, ACSAVW, IFIRST, MATSWP                   MAIN  84
     3 ,       SDSWIT,  LCSWIT, LFITSW, R, NOIT                         MAIN  85
     4 ,       SRATST, LSCH, LPUNSW, LSPL, LRANDC                       MAIN  86
     5 ,       LDAPRT, LDIPRT, LREG,   LHIPRT, NUDASW                   MAIN  87
     6 ,        LPOLSP, LCONSW, LNFIXZ, DCON1, DCON2                    MAIN  88
     7 ,       DCON3, DCON4, DCON5, WCON1, WCON2                        MAIN  89
     8 ,       WCON3, WCON4, WCON5, NOITIN, LPLCON                      MAIN  90
     9 ,       LPLSCT, LCOORS                                           MAIN  91
      COMMON /METRIC/ RTYPE,RM1,RECR,RR                                 MAIN  92
      COMMON /PLTCHR/ PTID,ITEM                                         MAIN  93
C                                                                       MAIN  94
       EQUIVALENCE (DIST,LDIST) , (LBLOCK,RVEC)                         MAIN  95
      EQUIVALENCE (PMAT(1,1),DIST(1)), (PMAT(1,2),DHAT(1))              MAIN  96
      EQUIVALENCE (XNUM(1),ONE),(XNUM(2),TWO),(XNUM(3),THREE),          MAIN  97
     .            (XNUM(4),FOUR),(XNUM(5),FIVE),(XNUM(6),SIX)           MAIN  98
       EXTERNAL                    WTRAN       , DTRAN                  MAIN  99
C                                                                       MAIN 100
      DATA  AAA,BEE,CEE,DEE,EEE,FFF /1HA,1HB,1HC,1HD,1HE,1HF/           MAIN 101
      DATA  ONE,TWO,THREE,FOUR,FIVE,SIX /1H1,1H2,1H3,1H4,1H5,1H6/       MAIN 102
      DATA KPACK,IJPACK,NPREVZ /10000,100,0/ , AIE/1H-/                 MAIN 103
C                                                                       MAIN 104
C      INITIALIZE PARAMETERS                                            MAIN 105
C                                                                       MAIN 106
 100   CONTINUE                                                         MAIN 107
C                                                                       MAIN 108
C               ALPHABETICAL ORDER                                      MAIN 109
       ACSAVW=0.66                                                      MAIN 110
       COSAVW=0.66                                                      MAIN 111
       CUTOFF=-1.23E+20                                                 MAIN 112
      DCON1=0.0                                                         MAIN 113
      DCON2=1.0                                                         MAIN 114
      DCON3=1.0                                                         MAIN 115
      DCON4=0.0                                                         MAIN 116
      DCON5=1.0                                                         MAIN 117
       GRNO(1) = 1                                                      MAIN 118
       IFIRST=1                                                         MAIN 119
      ISTC=0                                                            MAIN 120
      LCONSW=4                                                          MAIN 121
      LCOORS=2                                                          MAIN 122
       LCSWIT=1                                                         MAIN 123
       LDAPRT = 1                                                       MAIN 124
       LDIMD=1                                                          MAIN 125
       LDIMN=2                                                          MAIN 126
       LDIMX=2                                                          MAIN 127
       LDIPRT = 1                                                       MAIN 128
       LFITSW=1                                                         MAIN 129
       LHIPRT = 2                                                       MAIN 130
       LNFIXZ = 0                                                       MAIN 131
       LPOLSP = 100                                                     MAIN 132
      LPLCON=2                                                          MAIN 133
      LPLSCT=1                                                          MAIN 134
       LPUNSW = 2                                                       MAIN 135
       LRANDC=-101                                                      MAIN 136
      LREG=0                                                            MAIN 137
       LSCH=1                                                           MAIN 138
       LSPL=302                                                         MAIN 139
       MATSWP=101                                                       MAIN 140
       NOIT=50                                                          MAIN 141
      NOITIN=1                                                          MAIN 142
       NUDASW=1                                                         MAIN 143
       R=2.0                                                            MAIN 144
      SDSIN=1.0                                                         MAIN 145
       SDSWIT = 10                                                      MAIN 146
       SFGRMN=0.0                                                       MAIN 147
       SRATST=0.999                                                     MAIN 148
       STRESS=1.0                                                       MAIN 149
       STRMIN=0.01                                                      MAIN 150
      WCON1=0.0                                                         MAIN 151
      WCON2=1.0                                                         MAIN 152
      WCON3=1.0                                                         MAIN 153
      WCON4=0.0                                                         MAIN 154
      WCON5=1.0                                                         MAIN 155
       WRITE(LPRINT, 17)                                                MAIN 156
 17   FORMAT(1H1)                                                       MAIN 157
C                                                                       MAIN 158
C      CONTROL CARD READ                                                MAIN 159
C                                                                       MAIN 160
 1000  LCSWIT=1                                                         MAIN 161
       CALL CCACT                                                       MAIN 162
       GO TO (1000, 1100, 1165, 1182, 1200, 2000, 9000, 1190 ), LCSWIT  MAIN 163
C                                                                       MAIN 164
C      DATA READ                                                        MAIN 165
C                                                                       MAIN 166
 1100  CONTINUE                                                         MAIN 167
       NREPL1=1                                                         MAIN 168
       NREPL3 = 1                                                       MAIN 169
       LSPLIT=LSPL/100                                                  MAIN 170
       SPLITH = MOD(LSPL,100)                                           MAIN 171
       MATSW = MOD(MATSWP,100)                                          MAIN 172
       LBLKDS = MATSWP/100                                              MAIN 173
      LPOSW = MOD(LPOLSP,100)                                           MAIN 174
      LFITRM = LPOLSP/100                                               MAIN 175
      IF(LREG.EQ.0) GO TO 1104                                          MAIN 176
      SDSWIT=LREG                                                       MAIN 177
      IF(LREG.LT.10) SDSWIT = SDSWIT+LPOSW                              MAIN 178
1104  CONTINUE                                                          MAIN 179
      SDSWIT = SDSWIT+100*LFITRM                                        MAIN 180
C                                                                       MAIN 181
       IF(NUDASW.NE.1) MLASTD=MM                                        MAIN 182
       IF(NUDASW.NE.1) GO TO 1106                                       MAIN 183
       NUDASW=2                                                         MAIN 184
       MLASTD=0                                                         MAIN 185
      MMIN=0                                                            MAIN 186
       NOGRPS=0                                                         MAIN 187
1106   M=MLASTD                                                         MAIN 188
       MA = M+1                                                         MAIN 189
C                                                                       MAIN 190
      READ(LREAD,10)   TITLE                                            MAIN 191
      WRITE(LPRINT,11) TITLE                                            MAIN 192
       IF(MATSW.GE.4) GO TO 1102                                        MAIN 193
 1101  READ     (LREAD,12)     NPART,NREPL2,NREPL4                      MAIN 194
       WRITE    (LPRINT,13)    NPART,NREPL2,NREPL4                      MAIN 195
       NROWS = NPART                                                    MAIN 196
       NCOLS = NPART                                                    MAIN 197
       GO TO 1103                                                       MAIN 198
 1102  READ     (LREAD,12)     NROWS,NCOLS,NREPL2,NREPL4                MAIN 199
       WRITE    (LPRINT,13)    NROWS,NCOLS,NREPL2,NREPL4                MAIN 200
       NPART = NROWS+NCOLS                                              MAIN 201
C                                                                       MAIN 202
 1103  N=NPART                                                          MAIN 203
       IF(NREPL4.NE.0) NREPL3 = NREPL4                                  MAIN 204
       IF(NREPL2.NE.0) NREPL1 = NREPL2                                  MAIN 205
       IF(LBLKDS.EQ.2) N = NPART*NREPL3                                 MAIN 206
       IF(NREPL2.EQ.0  .OR.  NREPL4.EQ.0  )  WRITE  (LPRINT,16)         MAIN 207
 1105 READ(LREAD,10)     FMAT                                           MAIN 208
      WRITE(LPRINT,11)   FMAT                                           MAIN 209
C                                                                       MAIN 210
       DO 1162 NREPL=1,NREPL3                                           MAIN 211
       IJBLKD=0                                                         MAIN 212
       IF(LBLKDS.EQ.2) IJBLKD =(NREPL-1)*NPART                          MAIN 213
       MLASTG=M                                                         MAIN 214
       IA = 1                                                           MAIN 215
       IB = NROWS                                                       MAIN 216
       IF(MATSW.EQ.2) IA = IFIRST                                       MAIN 217
       IF(MATSW.EQ.3) IB = NROWS-(IFIRST-1)                             MAIN 218
C                                                                       MAIN 219
       DO 1160 I = IA,IB                                                MAIN 220
       MLASTR=M                                                         MAIN 221
       LTEMP = NCOLS                                                    MAIN 222
       IF(MATSW.EQ.2) LTEMP = I-IFIRST+1                                MAIN 223
       IF(MATSW.EQ.3) LTEMP = NCOLS-I-IFIRST+2                          MAIN 224
       MB = MA + LTEMP * NREPL1 - 1                                     MAIN 225
       ITRUE = I+IJBLKD                                                 MAIN 226
       IF(MATSW.EQ.4) ITRUE = ITRUE+NCOLS                               MAIN 227
C                                                                       MAIN 228
 1130  READ           (LREAD, FMAT) (DATA(MP),MP=MA,MB)                 MAIN 229
C                                                                       MAIN 230
       DO 1150 MP=MA,MB                                                 MAIN 231
       IF( DATA(MP)-CUTOFF )  1150, 1150, 1140                          MAIN 232
 1140  CONTINUE                                                         MAIN 233
       M = M + 1                                                        MAIN 234
       DATA(M)=DATA(MP)                                                 MAIN 235
       WW(M) = 1.0                                                      MAIN 236
       J=((MP-MA)/NREPL1)+1                                             MAIN 237
       IF(MATSW.EQ.3) J = J + (I-1) + (IFIRST-1)                        MAIN 238
       J = J+IJBLKD                                                     MAIN 239
       IF(MATSW.EQ.5) J = J + NROWS                                     MAIN 240
      IJ(M)=IJPACK*(ITRUE-1)+J-1                                        MAIN 241
       LDIST(M) = MP                                                    MAIN 242
 1150  CONTINUE                                                         MAIN 243
C                                                                       MAIN 244
       IF(LSPLIT.EQ.1 .AND. M.GT.MLASTR) GO TO 1155                     MAIN 245
       GO TO 1160                                                       MAIN 246
 1155  NOGRPS=NOGRPS+1                                                  MAIN 247
       GRNO(NOGRPS+1) = M+1                                             MAIN 248
       GRSDIS(NOGRPS)=SDSWIT                                            MAIN 249
 1160  MA = M + 1                                                       MAIN 250
C                                                                       MAIN 251
       IF(LSPLIT.EQ.2 .AND. M.GT.MLASTG) GO TO 1161                     MAIN 252
       GO TO 1162                                                       MAIN 253
 1161  NOGRPS=NOGRPS+1                                                  MAIN 254
       GRNO(NOGRPS+1) = M+1                                             MAIN 255
       GRSDIS(NOGRPS)=SDSWIT                                            MAIN 256
 1162  CONTINUE                                                         MAIN 257
C                                                                       MAIN 258
       IF(LSPLIT.EQ.3 .AND. M.GT.MLASTD) GO TO 1163                     MAIN 259
       IF(LSPLIT.EQ.4 .AND. SPLITH.EQ.2 .AND. M.GT.MLASTD) GO TO 1163   MAIN 260
       GO TO 1164                                                       MAIN 261
 1163  NOGRPS=NOGRPS+1                                                  MAIN 262
       GRSDIS(NOGRPS)=SDSWIT                                            MAIN 263
C                                                                       MAIN 264
 1164  CONTINUE                                                         MAIN 265
       GRNO(NOGRPS+1) = M+1                                             MAIN 266
       IF(LSPLIT.LE.3) SPLITH=2                                         MAIN 267
       IF(LSPLIT.EQ.4) SPLITH=1                                         MAIN 268
      IF((NOGRPS.EQ.1).OR.(MMIN.EQ.0)) MMIN=M                           MAIN 269
C                                                                       MAIN 270
       LSPL = 100 * LSPLIT + SPLITH                                     MAIN 271
C                                                                       MAIN 272
       MM = M                                                           MAIN 273
       GO TO 1000                                                       MAIN 274
C                                                                       MAIN 275
C      WEIGHTS READ                                                     MAIN 276
C                                                                       MAIN 277
 1165  CONTINUE                                                         MAIN 278
       M=MLASTD                                                         MAIN 279
       MA=M+1                                                           MAIN 280
       DO 1181 NREPL=1,NREPL3                                           MAIN 281
       IA = 1                                                           MAIN 282
       IB = NROWS                                                       MAIN 283
       IF(MATSW.EQ.2) IA = IFIRST                                       MAIN 284
       IF(MATSW.EQ.3) IB = NROWS-(IFIRST-1)                             MAIN 285
       DO 1180 I = IA,IB                                                MAIN 286
       LTEMP = NCOLS                                                    MAIN 287
       IF(MATSW.EQ.2) LTEMP = I-IFIRST+1                                MAIN 288
       IF(MATSW.EQ.3) LTEMP = NCOLS-I-IFIRST+2                          MAIN 289
       MB = MA + (LTEMP-1)*NREPL1                                       MAIN 290
 1172  READ           (LREAD,FMAT)(WW(MP),MP=MA,MB)                     MAIN 291
       DO 1177 MP=MA,MB                                                 MAIN 292
       IF(LDIST(M+1).GT.MP) GO TO 1177                                  MAIN 293
       M = M + 1                                                        MAIN 294
       WW(M)=WW(MP)                                                     MAIN 295
 1177  CONTINUE                                                         MAIN 296
 1180  MA = M + 1                                                       MAIN 297
 1181  CONTINUE                                                         MAIN 298
       GO TO 1000                                                       MAIN 299
C                                                                       MAIN 300
C      WEIGHT FORMATION BY WFUNCTION                                    MAIN 301
C                                                                       MAIN 302
 1182  CONTINUE                                                         MAIN 303
      WRITE(LPRINT,18) WCON1,WCON2,WCON3,WCON4,WCON5                    MAIN 304
 18   FORMAT(1H0,10X,69H** WEIGHTS FORMED ACCORDING TO RULE  WW(I,J)=S+TMAIN 305
     .*((A+B*DATA(I,J))**P)/,15X, 9H WHERE A=,F15.7,5X,2HB=,F15.7,5X,   MAIN 306
     .2HP=,F15.7,5X,2HS=,F15.7,5X,2HT=,F15.7)                           MAIN 307
       MA=MLASTD+1                                                      MAIN 308
       DO 1185 M=MA,MM                                                  MAIN 309
       TEMP1=DATA(M)                                                    MAIN 310
       WW(M) = WTRAN(TEMP1)                                             MAIN 311
 1185  CONTINUE                                                         MAIN 312
       GO TO 1000                                                       MAIN 313
C                                                                       MAIN 314
C     DATA TRANSFORMATION BY DFUNCTION                                  MAIN 315
C                                                                       MAIN 316
 1190 CONTINUE                                                          MAIN 317
      WRITE(LPRINT,19) DCON1,DCON2,DCON3,DCON4,DCON5                    MAIN 318
 19   FORMAT(1H0,10X,73H** DATA TRANSFORMED ACCORDING TO RULE  DATA(I,J)MAIN 319
     .=S+T*((A=B*DATA(I,J))**P),/,15X, 9H WHERE A=,F15.7,5X,2HB=,       MAIN 320
     .F15.7,5X,2HP=,F15.7,5X,2HS=,F15.7,5X,2HT=,F15.7)                  MAIN 321
      MA=MLASTD+1                                                       MAIN 322
      DO 1195 M=MA,MM                                                   MAIN 323
      DATA(M)=DTRAN(DATA(M))                                            MAIN 324
 1195 CONTINUE                                                          MAIN 325
      GO TO 1000                                                        MAIN 326
C                                                                       MAIN 327
C      CONFIGURATION READ                                               MAIN 328
C                                                                       MAIN 329
 1200  LCONSW=2                                                         MAIN 330
      READ(LREAD,10)  CTITLE                                            MAIN 331
      WRITE(LPRINT,11)CTITLE                                            MAIN 332
       READ           (LREAD, 12) NCON, LDIMCO                          MAIN 333
       WRITE            (LPRINT, 13) NCON, LDIMCO                       MAIN 334
      READ(LREAD,10)     FMAT                                           MAIN 335
      WRITE(LPRINT,11)   FMAT                                           MAIN 336
       DO 1210 I=1,NCON                                                 MAIN 337
       READ           (LREAD, FMAT) (X(I,L),L=1,LDIMCO)                 MAIN 338
      WRITE(LPRINT,4008) I,(X(I,L),L=1,LDIMCO)                          MAIN 339
 4008 FORMAT(1X,I2,10F7.3)                                              MAIN 340
 1210  CONTINUE                                                         MAIN 341
       GO TO 1000                                                       MAIN 342
C                                                                       MAIN 343
C      SOME INPUT FORMATS                                               MAIN 344
C                                                                       MAIN 345
 10   FORMAT (80A1)                                                     MAIN 346
 11   FORMAT (1H0,80A1)                                                 MAIN 347
   12  FORMAT(24I3)                                                     MAIN 348
   13  FORMAT(1X,24I3)                                                  MAIN 349
  15  FORMAT(20F4.0)                                                    MAIN 350
   16  FORMAT(118H THE NUMBER OF INTERNAL REPLICATES, OR THE NUMBER OF EMAIN 351
     1XTERNAL REPLICATES HAS NOT BEEN SPECIFIED. IT IS TAKEN TO BE 1 . )MAIN 352
C                                                                       MAIN 353
C      COMPUTATION                 *************************************MAIN 354
 2000 CONTINUE                                                          MAIN 355
C                                                                       MAIN 356
       FN=FLOAT (N)                                                     MAIN 357
       SQRTN=SQRT (FN)                                                  MAIN 358
       FNGRPS = FLOAT(NOGRPS)                                           MAIN 359
       LDIM=LDIMX                                                       MAIN 360
      RR=R                                                              MAIN 361
      RTYPE=3                                                           MAIN 362
      IF(R.EQ.1.0) RTYPE=1                                              MAIN 363
      IF(R    .EQ.2.0) RTYPE=2                                          MAIN 364
      RM1=R-1.0                                                         MAIN 365
      RECR=1.0/R                                                        MAIN 366
       IF(LDAPRT.EQ.2)                                                  MAIN 367
     1     CALL DATAPR(GRNO,MM,NOGRPS,1)                                MAIN 368
       WRITE  (LPRINT,23)                                               MAIN 369
C                                                                       MAIN 370
C        FINISH STARTING CONFIGURATION. LCONSW =   2,3,4.               MAIN 371
C      4=USE INICON, 3=SAVED, 2=WAS READ IN.                            MAIN 372
C      TO FILL IN ADDITIONAL POINTS, IF NEEDED                          MAIN 373
C        LRANDC<0, USE ARBITRARY, ELSE USE RANDOM.                      MAIN 374
C                                                                       MAIN 375
 2100  CONTINUE                                                         MAIN 376
      WRITE(LPRINT,11) TITLE                                            MAIN 377
      IF((LCONSW.EQ.4) .AND. (LRANDC.EQ.(-101)))GO TO 2140              MAIN 378
       ISTCON = 1                                                       MAIN 379
       IF( LCONSW .EQ. 2 ) ISTCON = NCON + 1                            MAIN 380
       IF( LCONSW .EQ. 3 ) ISTCON = NPREVZ + 1                          MAIN 381
       IF( ISTCON .GT. N ) GO TO 2200                                   MAIN 382
       IF(LRANDC.LE.0) GO TO 2110                                       MAIN 383
       DO 2105  K=1,LRANDC                                              MAIN 384
 2105  TEMP=RUNIFV(1)                                                   MAIN 385
 2110  TEMP1=0.01                                                       MAIN 386
       DO 2130 I=ISTCON,N                                               MAIN 387
       DO 2120 L=1,LDIMX                                                MAIN 388
 2120  X(I,L)=0.0                                                       MAIN 389
       K= MOD (I-1,LDIM) +1                                             MAIN 390
       X(I,K)=TEMP1                                                     MAIN 391
      IF(LRANDC.LT.0) GO TO 2130                                        MAIN 392
       DO 2125 L=1,LDIMX                                                MAIN 393
       TEMP = RUNIFV(1)                                                 MAIN 394
 2125  X(I,L) = ALOG( TEMP/(1.0-TEMP))                                  MAIN 395
 2130  TEMP1=TEMP1+0.01                                                 MAIN 396
      GO TO 2200                                                        MAIN 397
 2140  SDSWIT=MOD(GRSDIS(1),100)                                        MAIN 398
      IF(SDSWIT.EQ.11) SDSIN=-1.0                                       MAIN 399
      IF(N.GT.60) NOITIN=0                                              MAIN 400
      CALL INICON(N,LDIM,MMIN,NOITIN,LHIPRT,SDSIN,LSCH)                 MAIN 401
C                                                                       MAIN 402
C      SORT DATA AND IJ AND WW.                                         MAIN 403
C      ALSO RECORD BLOCKS OF EQUAL DATA VALUES.                         MAIN 404
C                                                                       MAIN 405
 2200  CONTINUE                                                         MAIN 406
       NPREVZ = N                                                       MAIN 407
       DO 2250 NG = 1,NOGRPS                                            MAIN 408
       MX = GRNO(NG)                                                    MAIN 409
       MY = GRNO(NG+1) - 1                                              MAIN 410
       MZ = MY - MX + 1                                                 MAIN 411
       SDSWIT = GRSDIS(NG)                                              MAIN 412
      SDSWIT = MOD(SDSWIT,100)                                          MAIN 413
      IF(SDSWIT.EQ.11)SDSWIT = -10                                      MAIN 414
       CALL SORT( DATA(MX),MZ,IJ(MX),WW(MX),DUMMY,2,SDSWIT)             MAIN 415
C                                                                       MAIN 416
C               SORT  WILL SORT THE MM ELEMNTS OF  DATA  IN ALGEBRAIC   MAIN 417
C              ORDER, ASCENDING OR DESCENDING ACCORDING TO WHETHER      MAIN 418
C              SDSWIT IS + OR -.                                        MAIN 419
C              AT THE SAME TIME, THE ELEMENTS IN  IJ  AND IN  WW  WILL BMAIN 420
C              REARRANGED IN EXACTLY THE SAME ORDER. THUS THE           MAIN 421
C              CORRESPONDENCE BETWEEN THE ELEMENTS OF  DATA  AND  IJ    MAIN 422
C              AND  WW  IS PRESERVED.                                   MAIN 423
C                                                                       MAIN 424
       DO 2240 MB = MX,MY                                               MAIN 425
      IF((DATA(MB+1).EQ.DATA(MB)).AND.(MB.NE.MY)) GO TO 2240            MAIN 426
      IJ(MB)=MOD(IJ(MB),KPACK)                                          MAIN 427
      IJ(MB)=IJ(MB)+KPACK                                               MAIN 428
 2240  CONTINUE                                                         MAIN 429
 2250  CONTINUE                                                         MAIN 430
C                                                                       MAIN 431
C      START COMPUTATION IN CURRENT DIMENSION                           MAIN 432
C                                                                       MAIN 433
 2300  FLDIM=FLOAT (LDIM)                                               MAIN 434
       ITNO=0                                                           MAIN 435
       COSAV=0.0                                                        MAIN 436
       SRATAV=0.8                                                       MAIN 437
       ACSAV=0.0                                                        MAIN 438
       STEP = 0.0                                                       MAIN 439
       NBAKUP = 0                                                       MAIN 440
C                                                                       MAIN 441
C      INITIALIZE GRADIENT                                              MAIN 442
C                                                                       MAIN 443
 2400  TEMP1=SQRT (1.0/FLDIM)                                           MAIN 444
       DO 2410 I=1,N                                                    MAIN 445
       DO 2410 L=1,LDIM                                                 MAIN 446
 2410  GR(I,L)=TEMP1                                                    MAIN 447
       SFGR=SQRTN                                                       MAIN 448
C                                                                       MAIN 449
C      PRINT HEADING FOR INFORMATION ABOUT SCALING IN CURRENT DIMENSION MAIN 450
C                                                                       MAIN 451
 2500  WRITE  (LPRINT, 20) N, MM, NOGRPS, LDIM                          MAIN 452
       IF(LHIPRT.EQ.2)  WRITE  (LPRINT,21)                              MAIN 453
       WRITE            (LPRINT, 22)                                    MAIN 454
 20    FORMAT(28H0HISTORY OF COMPUTATION. N= , I4,                      MAIN 455
     1 15H.     THERE ARE  , I6,                                        MAIN 456
     2 26H   DATA VALUES, SPLIT INTO  , I4,                             MAIN 457
     3 27H    LISTS.     DIMENSION =  , I4     )                        MAIN 458
 21    FORMAT(52H0ITERATION STRESS   SRAT SRATAV CAGRGL  COSAV  ACSAV,  MAIN 459
     1 16H    SFGR    STEP )                                            MAIN 460
   22  FORMAT(1X)                                                       MAIN 461
 23    FORMAT(1H1)                                                      MAIN 462
C                                                                       MAIN 463
C      START CURRENT ITERATION     *************************************MAIN 464
C                                                                       MAIN 465
C      NORMALIZE CONFIGURATION. MOVE AND CLEAR GRADIENT.                MAIN 466
C                                                                       MAIN 467
 3000 TEMP1=1.0                                                         MAIN 468
      IF(LCOORS.EQ.0) GO TO 3035                                        MAIN 469
      TEMP1=0.0                                                         MAIN 470
       DO 3030 L=1,LDIM                                                 MAIN 471
       TEMP2=0.0                                                        MAIN 472
       DO 3010 I=1,N                                                    MAIN 473
 3010  TEMP2=TEMP2+X(I,L)                                               MAIN 474
       TEMP2=TEMP2/FN                                                   MAIN 475
       DO 3020 I=1,N                                                    MAIN 476
       X(I,L)=X(I,L)-TEMP2                                              MAIN 477
 3020  TEMP1=TEMP1+X(I,L)**2                                            MAIN 478
 3030  CONTINUE                                                         MAIN 479
       TEMP1=SQRT (FN/TEMP1)                                            MAIN 480
 3035  DO 3040 L=1,LDIM                                                 MAIN 481
       DO 3040 I=1,N                                                    MAIN 482
       X(I,L)=TEMP1*X(I,L)                                              MAIN 483
       GL(I,L)=TEMP1*GR(I,L)                                            MAIN 484
 3040  GR(I,L)=0.0                                                      MAIN 485
       SFGL=TEMP1*SFGR                                                  MAIN 486
C                                                                       MAIN 487
       STBAMU = TEMP1                                                   MAIN 488
C      LOOP THROUGH THE DATA GROUPS                                     MAIN 489
C                                                                       MAIN 490
       STRLST = STRESS                                                  MAIN 491
       STRESS = 0.0                                                     MAIN 492
       DO 3340 NG = 1,NOGRPS                                            MAIN 493
       MX = GRNO(NG)                                                    MAIN 494
       MY = GRNO(NG+1) - 1                                              MAIN 495
       MZ = MY - MX + 1                                                 MAIN 496
       SDSWIT = GRSDIS(NG)                                              MAIN 497
      LFITRM = SDSWIT/100                                               MAIN 498
      SDSWIT = MOD(SDSWIT,100)                                          MAIN 499
      IF(SDSWIT.EQ.11)SDSWIT = -10                                      MAIN 500
C                                                                       MAIN 501
C      COMPUTE DISTANCES AND FIND BEST-FITTING MONOTONE PSEUDO-DISTANCESMAIN 502
C                                                                       MAIN 503
       SUMW = 0.0                                                       MAIN 504
       DBAR = 0.0                                                       MAIN 505
       DO 3120 M=MX,MY                                                  MAIN 506
      LTEMP1=MOD(IJ(M),KPACK)                                           MAIN 507
      I=LTEMP1/IJPACK+1                                                 MAIN 508
      J=MOD(LTEMP1,IJPACK)+1                                            MAIN 509
       TEMP1=0.0                                                        MAIN 510
       DO 3110 L=1,LDIM                                                 MAIN 511
 3110  TEMP1=TEMP1+RPOWER (X(I,L)-X(J,L))                               MAIN 512
       DIST(M)=RROOT (TEMP1)                                            MAIN 513
      DBAR=DBAR+DIST(M)*WW(M)                                           MAIN 514
       SUMW = SUMW + WW(M)                                              MAIN 515
 3120  CONTINUE                                                         MAIN 516
      DBAR=DBAR/SUMW                                                    MAIN 517
C              DBAS IS EITHER  0  OR  DBAR  ACCORDING TO WHETHER        MAIN 518
C              STRESS FORMULA 1 OR 2 IS BEING USED.                     MAIN 519
       DBAS = 0.0                                                       MAIN 520
       IF(LSCH.EQ.2) DBAS = DBAR                                        MAIN 521
       IF(IABS(SDSWIT).GE.10)                                           MAIN 522
     1   CALL FITM(MX,MY,LFITSW)                                        MAIN 523
       IF(IABS(SDSWIT).LT.10)                                           MAIN 524
     1     CALL FITP(MX,MY,SDSWIT,LFITRM)                               MAIN 525
C                                                                       MAIN 526
C      CALCULATE U, T, AND STRESS                                       MAIN 527
C                                                                       MAIN 528
 3200  U=0.0                                                            MAIN 529
       T=0.0                                                            MAIN 530
       DO 3210 M=MX,MY                                                  MAIN 531
      U=U+(DIST(M)-DHAT(M))**2*WW(M)                                    MAIN 532
 3210 T=T+(DIST(M)-DBAS)**2*WW(M)                                       MAIN 533
 3215  U=SQRT (U)                                                       MAIN 534
       TEMP1=T                                                          MAIN 535
       T=SQRT (T)                                                       MAIN 536
       GRSTRS(NG) = U/T                                                 MAIN 537
       STRESS = STRESS + GRSTRS(NG)**2                                  MAIN 538
       IF(U.EQ.0.0)   GO TO 3340                                        MAIN 539
 3220  RUT=1.0/(U*T)                                                    MAIN 540
       UOT3=U/(TEMP1*T)                                                 MAIN 541
C                                                                       MAIN 542
C      CALCULATE THE (NEGATIVE) GRADIENT                                MAIN 543
C                                                                       MAIN 544
 3300  DO 3330 M = MX,MY                                                MAIN 545
       IF(DIST(M).EQ.0.0) GO TO 3330                                    MAIN 546
      LTEMP1=MOD(IJ(M),KPACK)                                           MAIN 547
      I=LTEMP1/IJPACK+1                                                 MAIN 548
      J=MOD(LTEMP1,IJPACK)+1                                            MAIN 549
       FACTOR=UOT3*(DIST(M)-DBAS)-RUT*(DIST(M)-DHAT(M))                 MAIN 550
       FACTOR = (FACTOR/RM1POW(DIST(M)) ) * WW(M)                       MAIN 551
       FACTOR = FACTOR * GRSTRS(NG)                                     MAIN 552
       DO 3320 L=1,LDIM                                                 MAIN 553
       TEMP1 = FACTOR * RM1POW(X(I,L)-X(J,L))                           MAIN 554
       GR(I,L)=GR(I,L)+TEMP1                                            MAIN 555
 3320  GR(J,L)=GR(J,L)-TEMP1                                            MAIN 556
 3330  CONTINUE                                                         MAIN 557
 3340  CONTINUE                                                         MAIN 558
       IF(STRESS .EQ. 0.0 ) GO TO 3700                                  MAIN 559
       STRESS = SQRT( STRESS / FNGRPS )                                 MAIN 560
C                                                                       MAIN 561
C             AVOID MOVING FIXED POINTS                                 MAIN 562
C                                                                       MAIN 563
       IF( LNFIXZ .EQ. 0) GO TO 3400                                    MAIN 564
       DO 3360 I=1,LNFIXZ                                               MAIN 565
       DO 3360 L=1,LDIM                                                 MAIN 566
 3360  GR(I,L) = 0.0                                                    MAIN 567
C                                                                       MAIN 568
C      FIND SCALE FACTOR OF GRADIENT, ANGLE-COSINE BETWEEN GRADIENT     MAIN 569
C      AND PREVIOUS GRADIENT.                                           MAIN 570
C                                                                       MAIN 571
 3400  SFGR=0.0                                                         MAIN 572
       CAGRGL=0.0                                                       MAIN 573
       DO 3410 I=1,N                                                    MAIN 574
       DO 3410 L=1,LDIM                                                 MAIN 575
       SFGR=SFGR+GR(I,L)**2                                             MAIN 576
 3410  CAGRGL=CAGRGL+GR(I,L)*GL(I,L)                                    MAIN 577
       SFGR=SQRT (SFGR/FN)                                              MAIN 578
C              IF GRADIENT   0.0, SKIP AHEAD.                           MAIN 579
       IF(SFGR) 3420,3700,3420                                          MAIN 580
 3420  TEMP1=SFGR*SFGL*FN                                               MAIN 581
       CAGRGL=CAGRGL/TEMP1                                              MAIN 582
C                                                                       MAIN 583
       IF(ITNO.EQ.0  .OR.  NBAKUP.GE.4) GO TO 3500                      MAIN 584
       IF(CAGRGL.GT.(-0.95)  .AND.  STRESS/STRLST.LT. 1.2001 ) GOTO 3500MAIN 585
C                                                                       MAIN 586
C      BACK UP ALONG LAST GRADIENT                                      MAIN 587
C                                                                       MAIN 588
       NBAKUP = NBAKUP + 1                                              MAIN 589
      IF(NBAKUP.EQ.1) STBASC=1.0                                        MAIN 590
      STBASC=STBASC*STBAMU                                              MAIN 591
       TEMP1 = STEP                                                     MAIN 592
       STEP = STEP / 10.0                                               MAIN 593
       WRITE  (LPRINT,38)   STRESS, CAGRGL,  STEP                       MAIN 594
 38    FORMAT(10X, F7.3, 14X, F7.3, 22X, F8.4 )                         MAIN 595
       TEMP1 = (TEMP1 - STEP) / SFGL                                    MAIN 596
      TEMP1=STBASC*TEMP1                                                MAIN 597
       DO 3430 I = 1,N                                                  MAIN 598
       DO 3430 L = 1,LDIM                                               MAIN 599
       X(I,L) = X(I,L) - TEMP1*GL(I,L)                                  MAIN 600
 3430  GR(I,L) = GL(I,L)                                                MAIN 601
       SFGR = SFGL                                                      MAIN 602
       STRESS = STRLST                                                  MAIN 603
       GO TO 3000                                                       MAIN 604
C                                                                       MAIN 605
C      STEP SIZE CALCULATIONS                                           MAIN 606
C                                                                       MAIN 607
3500   IF(ITNO) 9999, 3510, 3520                                        MAIN 608
 3510  SRAT=0.8                                                         MAIN 609
       GO TO 3530                                                       MAIN 610
 3520  SRAT=STRESS/STRLST                                               MAIN 611
 3530  CALL       NEWSTP( STEP, ITNO, SFGR, STRESS,                     MAIN 612
     1 CAGRGL, COSAV, ACSAV, COSAVW, ACSAVW, SRAT, SRATAV )             MAIN 613
       NBAKUP = 0                                                       MAIN 614
C                                                                       MAIN 615
C      PRINT CURRENT STATUS OF COMPUTATION                              MAIN 616
C                                                                       MAIN 617
 3700  IF(LHIPRT.EQ.2) WRITE (LPRINT,30) ITNO,STRESS,SRAT,SRATAV,       MAIN 618
     1 CAGRGL,COSAV,ACSAV,SFGR,STEP                                     MAIN 619
   30  FORMAT(I10,F7.3,F7.3,F7.3,F7.3,F7.3,F7.3,F8.4,   F8.4)           MAIN 620
C                                                                       MAIN 621
C      DECIDE WHETHER TO CONTINUE ITERATING                             MAIN 622
C                                                                       MAIN 623
 3800  IF(STRESS) 9999, 3840, 3810                                      MAIN 624
 3810  IF(SFGR-SFGRMN) 3850, 3850, 3815                                 MAIN 625
 3815  TEMP1 = 0.5 * (1.0+SRATST)                                       MAIN 626
       TEMP2 = 1.0 - TEMP1                                              MAIN 627
       IF( ABS (SRAT-TEMP1) - TEMP2 ) 3816, 3816, 3820                  MAIN 628
 3816  IF( ABS (SRATAV-TEMP1) - TEMP2 ) 3850, 3850, 3820                MAIN 629
 3820  IF(STRESS-STRMIN) 3860, 3860, 3830                               MAIN 630
 3830  IF(ITNO-NOIT) 3900, 3870, 9999                                   MAIN 631
 3840  CONTINUE                                                         MAIN 632
       WRITE            (LPRINT, 31)                                    MAIN 633
 31    FORMAT(24H0ZERO STRESS WAS REACHED )                             MAIN 634
       GO TO 4000                                                       MAIN 635
 3850  CONTINUE                                                         MAIN 636
       WRITE            (LPRINT, 32)                                    MAIN 637
 32    FORMAT(21H0MINIMUM WAS ACHIEVED )                                MAIN 638
       GO TO 4000                                                       MAIN 639
 3860  CONTINUE                                                         MAIN 640
       WRITE            (LPRINT, 33)                                    MAIN 641
 33    FORMAT(32H0SATISFACTORY STRESS WAS REACHED )                     MAIN 642
       GO TO 4000                                                       MAIN 643
 3870  CONTINUE                                                         MAIN 644
       WRITE            (LPRINT, 34)                                    MAIN 645
 34    FORMAT(39H0MAXIMUM NUMBER OF ITERATIONS WERE USED )              MAIN 646
       GO TO 4000                                                       MAIN 647
C                                                                       MAIN 648
C      CONTINUE ITERATING                                               MAIN 649
C                                                                       MAIN 650
 3900  ITNO=ITNO+1                                                      MAIN 651
       TEMP1=STEP/SFGR                                                  MAIN 652
       DO 3910 I=1,N                                                    MAIN 653
       DO 3910 L=1,LDIM                                                 MAIN 654
 3910  X(I,L)=X(I,L)+TEMP1*GR(I,L)                                      MAIN 655
       GO TO 3000                                                       MAIN 656
C                                                                       MAIN 657
C      STOP ITERATING              *************************************MAIN 658
C                                                                       MAIN 659
C      ROTATE FINAL CONFIGURATION TO PRINCIPAL COMPONENTS               MAIN 660
C                                                                       MAIN 661
 4000 IF(LCOORS.LT.2) GO TO 4002                                        MAIN 662
C      COMPUTE X TRANSPOSE TIMES X AND STORE UPPER HALF IN RVEC         MAIN 663
      KK=0                                                              MAIN 664
      DO 4410 K=1,LDIM                                                  MAIN 665
      DO 4410 J=1,K                                                     MAIN 666
      SUM=0.0                                                           MAIN 667
      DO 4405 I=1,N                                                     MAIN 668
 4405 SUM=SUM+X(I,J)*X(I,K)                                             MAIN 669
      KK=KK+1                                                           MAIN 670
 4410 RVEC(KK)=SUM                                                      MAIN 671
C     COMPUTE MATRIX OF EIGENVECTORS, GL                                MAIN 672
      CALL SGEV(LDIM,LDIM,GL)                                           MAIN 673
C      COMPUTE X TIMES GL                                               MAIN 674
      DO 4420 K=1,LDIM                                                  MAIN 675
      DO 4420 J=1,N                                                     MAIN 676
      SUM=0.0                                                           MAIN 677
      DO 4415 I=1,LDIM                                                  MAIN 678
 4415 SUM=SUM+X(J,I)*GL(I,K)                                            MAIN 679
 4420 GR(J,K)=SUM                                                       MAIN 680
      DO 4425 I=1,N                                                     MAIN 681
      DO 4425 J=1,LDIM                                                  MAIN 682
 4425 X(I,J)=GR(I,J)                                                    MAIN 683
      WRITE(LPRINT,152)                                                 MAIN 684
 152  FORMAT(66H0THE FINAL CONFIGURATION HAS BEEN ROTATED TO PRINCIPAL CMAIN 685
     1OMPONENTS.)                                                       MAIN 686
      GO TO 4004                                                        MAIN 687
C                                                                       MAIN 688
 4002 IF(LCOORS.EQ.1) WRITE(LPRINT,153)                                 MAIN 689
 153  FORMAT(110H0THE CONFIGURATION HAS BEEN NORMALIZED DURING THE ITERAMAIN 690
     1TIONS BUT HAS NOT BEEN ROTATED TO PRINCIPAL COMPONENTS.)          MAIN 691
      IF(LCOORS.EQ.0) WRITE(LPRINT,154)                                 MAIN 692
 154  FORMAT(70H0NO NORMALIZATION WAS DONE TO THE CONFIGURATION DURING TMAIN 693
     1HE ITERATIONS.)                                                   MAIN 694
 4004  WRITE            (LPRINT, 40)N,LDIM,STRESS,LSCH,(L,L=1,LDIM)     MAIN 695
 40    FORMAT(27H0THE FINAL CONFIGURATION OF,I4,                        MAIN 696
     1 10H POINTS IN,I3, 22H DIMENSIONS HAS STRESS,F7.3,8H FORMULA ,I2  MAIN 697
     2 /1X/29HLABEL FOR CONFIGURATION PLOTS,10X,19HFINAL CONFIGURATION, MAIN 698
     3     /,39X,10I7)                                                  MAIN 699
       IF(LPUNSW.EQ.2) GO TO 4005                                       MAIN 700
      WRITE    (LPUNCH,41) (TITLE(I),I=1,80),N,LDIM                     MAIN 701
 41   FORMAT (14H CONFIGURATION/80A1/2I3)                               MAIN 702
      WRITE     (LPUNCH,52)                                             MAIN 703
 52   FORMAT (12H (2X,10F7.3))                                          MAIN 704
 4005  DO 4010 I=1,N                                                    MAIN 705
      WRITE(LPRINT,42) PTID(I),I,(X(I,L),L=1,LDIM)                      MAIN 706
       IF(LPUNSW.EQ.2) GO TO 4010                                       MAIN 707
       WRITE            (LPUNCH, 43)I,(X(I,L),L=1,LDIM)                 MAIN 708
 4010  CONTINUE                                                         MAIN 709
 42   FORMAT(20X,A1,18X,I2,10F7.3)                                      MAIN 710
   43  FORMAT(I2,10F7.3)                                                MAIN 711
       WRITE  (LPRINT,46)                                               MAIN 712
       DO 4020 NG=1,NOGRPS                                              MAIN 713
       MZ = GRNO(NG+1) - GRNO(NG)                                       MAIN 714
      SDSWT1=MOD(GRSDIS(NG),100) +1                                     MAIN 715
      IF(SDSWT1-11) 4016,4013,4014                                      MAIN 716
 4013 WRITE(LPRINT,60) NG,MZ,GRSTRS(NG)                                 MAIN 717
 60   FORMAT(1X,I5,2X,I5,F7.3,11H  ASCENDING)                           MAIN 718
      GO TO 4020                                                        MAIN 719
 4014 WRITE(LPRINT,62) NG,MZ,GRSTRS(NG)                                 MAIN 720
 62   FORMAT(1X,I5,2X,I5,F7.3,11H DESCENDING)                           MAIN 721
      GO TO 4020                                                        MAIN 722
 4016 WRITE(LPRINT,61) NG,MZ,GRSTRS(NG),(PCOEFF(I),I=1,SDSWT1)          MAIN 723
 61   FORMAT(1X,I5,2X,I5,F7.3,12H  POLYNOMIAL,5F15.7)                   MAIN 724
 4020  CONTINUE                                                         MAIN 725
 46   FORMAT(14H0DATA GROUP(S) ,/73H0SERIAL  COUNT STRESS REGRESSION COEMAIN 726
     .FFICIENTS (FROM DEGREE 0 TO MAX OF 4))                            MAIN 727
 4030  CONTINUE                                                         MAIN 728
       IF(LDIPRT.EQ.2)                                                  MAIN 729
     1       CALL DATAPR(GRNO,MM,NOGRPS,2)                              MAIN 730
C                                                                       MAIN 731
C        PLOTTING SECTION                                               MAIN 732
C                                                                       MAIN 733
C        PLOT DIST AND DHAT VS. DATA                                    MAIN 734
      IF(LPLSCT.EQ.0) GO TO 5100                                        MAIN 735
      WRITE(LPRINT,44) LDIM,LSCH,STRESS,TITLE                           MAIN 736
 44   FORMAT(52H1DIST(D) AND DHAT(-) (Y-AXIS) VS. DATA (X-AXIS), FOR,I3,MAIN 737
     .   27H DIMENSIONS. STRESS,FORMULA,I2,2H,=F8.4,/,30X,80A1)         MAIN 738
      PTID(1)=DEE                                                       MAIN 739
      PTID(2)=AIE                                                       MAIN 740
      CALL PLOT(DATA,PMAT,0.,0.,0.,0.,MM,2,-1,1800,1)                   MAIN 741
C        PLOT ADDITIONAL SCATTER DIAGRAMS BY GROUPS--FIRST  5 GROUPS    MAIN 742
C         IF LPLSCT=1,   ALL GROUPS IF LPLSCT=2                         MAIN 743
      IF(NOGRPS.EQ.1) GO TO 5090                                        MAIN 744
      KPLT=NOGRPS                                                       MAIN 745
      IF(LPLSCT.EQ.1) KPLT=MIN0(KPLT,5)                                 MAIN 746
      DO 5060 NG=1,KPLT                                                 MAIN 747
      MX=GRNO(NG)                                                       MAIN 748
      MY=GRNO(NG+1)-1                                                   MAIN 749
      MZ=MY-MX+1                                                        MAIN 750
      WRITE(LPRINT,4062) NG,LSCH,GRSTRS(NG),LDIM,TITLE                  MAIN 751
 4062 FORMAT(43H1DIST(D) AND DHAT(-) VS. DATA FOR GROUP NO.,I3,17H.  STRMAIN 752
     .ESS,FORMULA,I2,2H,=,F8.4,12H  DIMENSION=,I3,/,30X,80A1)           MAIN 753
 5060 CALL PLOT(DATA,PMAT,0.,0.,0.,0.,MZ,2,-1,1800,MX)                  MAIN 754
 5090 CONTINUE                                                          MAIN 755
      PTID(1)=AAA                                                       MAIN 756
      PTID(2)=BEE                                                       MAIN 757
C                                                                       MAIN 758
C        PLOT CONFIGURATION                                             MAIN 759
 5100 IF(LPLCON.EQ.0) GO TO 5200                                        MAIN 760
      IF(LDIM-1) 5190,5190,5125                                         MAIN 761
 5125 IF(LPLCON.EQ.1) GO TO 5160                                        MAIN 762
C        PLOT ALL PAIRS OF DIMENSIONS                                   MAIN 763
      DO 5145 J=2,LDIM                                                  MAIN 764
      JMIN1=J-1                                                         MAIN 765
      DO 5145 I=1,JMIN1                                                 MAIN 766
      WRITE(LPRINT,5015) J,I,TITLE                                      MAIN 767
 5015 FORMAT(31H1CONFIGURATION PLOT   DIMENSION,I2,23H (Y-AXIS) VS. DIMEMAIN 768
     .NSION,I3,9H (X-AXIS),/,30X,80A1)                                  MAIN 769
 5145 CALL PLOT(X(1,I),X(1,J),0.,0.,0.,0.,N,1,2,100,1)                  MAIN 770
      GO TO 5200                                                        MAIN 771
C        PLOT PAIRS OF DIMENSIONS SO THAT ALL DIMENSIONS ARE PLOTTED ONCMAIN 772
 5160 DO 5180 J=2,LDIM,2                                                MAIN 773
      IF((LDIM.EQ.5).AND.(J.EQ.4)) GO TO 5185                           MAIN 774
      JMIN1=J-1                                                         MAIN 775
      WRITE(LPRINT,5015) J,JMIN1,TITLE                                  MAIN 776
      CALL PLOT(X(1,JMIN1),X(1,J),0.,0.,0.,0.,N,1,2,100,1)              MAIN 777
      IF((J.NE.2) .OR. (MOD(LDIM,2).EQ.0)) GO TO 5180                   MAIN 778
      JP1=J+1                                                           MAIN 779
      WRITE(LPRINT,5015) JP1,JMIN1,TITLE                                MAIN 780
      CALL PLOT(X(1,JMIN1),X(1,JP1),0.,0.,0.,0.,N,1,2,100,1)            MAIN 781
 5180 CONTINUE                                                          MAIN 782
      GO TO 5200                                                        MAIN 783
 5185 WRITE(LPRINT,5015) LDIM,J,TITLE                                   MAIN 784
      CALL PLOT(X(1,J),X(1,LDIM),0.,0.,0.,0.,N,1,2,100,1)               MAIN 785
      GO TO 5200                                                        MAIN 786
C        ONE DIMENSIONAL CONFIGURATION PLOT                             MAIN 787
 5190 WRITE(LPRINT,5195) TITLE                                          MAIN 788
 5195 FORMAT(1H1,20X,34HPLOT OF CONFIGURATION, DIMENSION 1,/,30X,80A1)  MAIN 789
      CALL PLOT(X(1,1),X(1,1),0.,0.,0.,0.,N,1,3,100,1)                  MAIN 790
C                                                                       MAIN 791
C      CHANGE DIMENSION                                                 MAIN 792
C                                                                       MAIN 793
 5200 ISTC=ISTC+1                                                       MAIN 794
      STPL(ISTC)= STRESS                                                MAIN 795
      DIMSV(ISTC)=LDIM                                                  MAIN 796
 4100  LDIM=LDIM-LDIMD                                                  MAIN 797
       IF(LDIM-LDIMN)4110,2300,2300                                     MAIN 798
C                                                                       MAIN 799
C     PLOT STRESS VERSUS DIMENSION IF MORE THAN TWO POINTS              MAIN 800
C                                                                       MAIN 801
 4110 CONTINUE                                                          MAIN 802
      IF(ISTC.LT.3) GO TO 100                                           MAIN 803
      WRITE(LPRINT, 50)  TITLE                                          MAIN 804
 50   FORMAT (1H1,10X,31HPLOT OF STRESS VERSUS DIMENSION, /,30X,80A1)   MAIN 805
C      INVERT TO IMPROVE READING OF PLOT                                MAIN 806
       NVT=ISTC/2                                                       MAIN 807
       DO 8050 I=1,NVT                                                  MAIN 808
       IC=ISTC+1-I                                                      MAIN 809
       TEMP = DIMSV(I)                                                  MAIN 810
       SEMP = STPL(I)                                                   MAIN 811
       DIMSV(I) = DIMSV(IC)                                             MAIN 812
      STPL(I) =STPL(IC)                                                 MAIN 813
       DIMSV(IC) = TEMP                                                 MAIN 814
       STPL(IC) = SEMP                                                  MAIN 815
 8050  CONTINUE                                                         MAIN 816
C      FIND MAX STRESS FOR Y-AXIS , AND FILL IN PLOTTING CHARACTERS     MAIN 817
      YMA=0.0                                                           MAIN 818
      DO 8231 I=1,ISTC                                                  MAIN 819
      IDIM=DIMSV(I)                                                     MAIN 820
      PTID(I)=XNUM(IDIM)                                                MAIN 821
 8231 IF(YMA .LE. STPL(I)) YMA=STPL(I)                                  MAIN 822
      CALL PLOT(DIMSV,STPL,10.,YMA,0.,0.,ISTC,1,-2,10,1)                MAIN 823
      PTID(1)=AAA                                                       MAIN 824
      PTID(2)=BEE                                                       MAIN 825
      PTID(3)=CEE                                                       MAIN 826
      PTID(4)=DEE                                                       MAIN 827
      PTID(5)=EEE                                                       MAIN 828
      PTID(6)=FFF                                                       MAIN 829
C                                                                       MAIN 830
C              REINITIALIZE, AND RETURN FOR MORE INPUT.                 MAIN 831
       GO TO 100                                                        MAIN 832
C                                                                       MAIN 833
C      NORMAL TERMINATION, AFTER READING  STOP  ONCONTROL CARD          MAIN 834
 9000  STOP                                                             MAIN 835
C                                                                       MAIN 836
C      TROUBLE EXIT                                                     MAIN 837
C                                                                       MAIN 838
 9999  WRITE            (LPRINT, 99)                                    MAIN 839
 99    FORMAT(52H0KRUSKAL. IMPOSSIBLE BRANCH TAKEN FROM IF STATEMENT. ) MAIN 840
       STOP                                                             MAIN 841
C                                                                       MAIN 842
      END                                                               MAIN 843
$      FORTRAN
CBLOCK DATA                                                             BLOCK  1
       BLOCK   DATA                                                     BLOCK  2
C BLDA             FOR KYST                  JANUARY,1973               BLOCK  3
C WRITTEN BY J.KRUSKAL                                                  BLOCK  4
C MODIFIED FOR KYST BY J.KRUSKAL AND J.SEERY JANUARY,1973               BLOCK  5
C                                                                       BLOCK  6
C      THIS PROGRAM HOLDS THE  TABLE OF WORDS WHICH CCACT CONSULTS      BLOCK  7
C      LOGICALLY IT CONSISTS OF SEVERAL TABLES                          BLOCK  8
C      THE VALUES IN MTAB INDICATE WHICH ROWS START NEW TABLES          BLOCK  9
C      EACH ENTRY IN THE TABLE CONTAINS 4 ITEMS                         BLOCK 10
C      THE FIRST IS THE CODED ALPHABETIC WORD                           BLOCK 11
C             (ONE WORD MAY BE ENTERED SEVERAL TIMES )                  BLOCK 12
C      THE SECOND INDICATES THE NATURE OF THIS ENTRY                    BLOCK 13
C              1 INDICATES INTEGER PARAMETERS WHICH MUST BE  EXPLIC9TLY BLOCK 14
C                  READ IN                                              BLOCK 15
C              2 INDICATES REAL PARAMETERS WHICH MUST BE EXPLICITLY     BLOCK 16
C                  READ IN                                              BLOCK 17
C              3 INDICATES AN INTEGER PARAMETER WHICH IS IMPLICITLY     BLOCK 18
C                  DEFINED                                              BLOCK 19
C              4 INDICATES A REAL     PARAMETER WHICH IS IMPLICITLY     BLOCK 20
C                  DEFINED                                              BLOCK 21
C              5 INDICATES A PARAMETER WHICH BELONGS TO CCACT ONLY      BLOCK 22
C              6 INDICATES AN IMPLICITLY DEFINED INTEGER PARAMETER      BLOCK 23
C                EFFECTING PROGRAM EXECUTION UPON RETURN FROM CCACT     BLOCK 24
C      THE THIRD INDICATES WHICH PARAMTER IS INVOLVED                   BLOCK 25
C              SEPARATE NUMBERING FOR MAIN PROGRAM PARAMETERS AND CCACT BLOCK 26
C      THE FOURTH GIVES THE VALUE FOR ANY INTERNALLY DEFINED PARAMETERS BLOCK 27
C      IF THE PARAMETER IS OF TYPE 3,4 OR 6.IF OF TYPE 1 OR 2, THIS     BLOCK 28
C      ENTRY GIVES THE NUMBER OF ITEMS TO BE READ IN.                   BLOCK 29
C                                                                       BLOCK 30
       INTEGER                         LPAR(   42 )                     BLOCK 31
       REAL                             PAR(   42 )                     BLOCK 32
       EQUIVALENCE           (LPAR,PAR)                                 BLOCK 33
C                                                                       BLOCK 34
      INTEGER    MTAB(13)                                               BLOCK 35
      INTEGER TAB1(4,16),TAB2(4,17),TAB3(4,16),TAB4(4,19),TAB5(4,18)    BLOCK 36
      REAL PTID(100),ITEM(101)                                          BLOCK 37
C                                                                       BLOCK 38
       COMMON /ACCUR/  PRECSN,XMAG,XMAXN                                BLOCK 39
       COMMON  /MDSCL1/  LREAD, LPRINT, LPUNCH, LSCRAT                  BLOCK 40
       COMMON  /MDSCL2/  LPAR                                           BLOCK 41
       COMMON  /MDSCL3/  MTAB,TAB1,TAB2,TAB3,TAB4,TAB5                  BLOCK 42
      COMMON /PLTCHR/ PTID,ITEM                                         BLOCK 43
C                                                                       BLOCK 44
      DATA  PRECSN,XMAG,XMAXN /1.5E-8,1.0E-38,1.0E38/    ,              BLOCK 45
     .      LREAD,LPRINT,LPUNCH,LSCRAT /5,6,7,8/                        BLOCK 46
      DATA MTAB /1,50,52,57,63,65,74,76,78,81,84,87,87/                 BLOCK 47
C                                                                       BLOCK 48
C     TAB1 ENTRIES PERTAIN SUCCESSIVELY TO THESE CONTROL WORDS          BLOCK 49
C          DIMMAX,    DIMMIN,    DIMDIF,    CUTOFF,    STRMIN,          BLOCK 50
C          SFGRMN,    COSAVW,    ACSAVW,    DIAGON,    MATRIX,          BLOCK 51
C          HALFMA,    LOWERH,    UPPERH,    LOWERC,    ARBITR,          BLOCK 52
C          TORSCA                                                       BLOCK 53
       DATA  TAB1/                                                      BLOCK 54
     1  6989,    1,    1,    1,                                         BLOCK 55
     2  5375,    1,    2,    1,                                         BLOCK 56
     3  6923,    1,    3,    1,                                         BLOCK 57
     4 13520,    2,    4,    1,                                         BLOCK 58
     5   390,    2,    5,    1,                                         BLOCK 59
     6  2553,    2,    6,    1,                                         BLOCK 60
     7  7303,    2,    7,    1,                                         BLOCK 61
     8 11226,    2,    8,    1,                                         BLOCK 62
     9 11136,    5,    1,    2,                                         BLOCK 63
     A  9277,    3,   10,    1,                                         BLOCK 64
     B  3527,    3,   10,    2,                                         BLOCK 65
     C  6069,    3,   10,    2,                                         BLOCK 66
     D  8938,    3,   10,    3,                                         BLOCK 67
     E  1754,    3,   10,    4,                                         BLOCK 68
     F 10499,    3,   20,  -99,                                         BLOCK 69
     G  2048,    3,   27,    4/                                         BLOCK 70
C                                                                       BLOCK 71
C     TAB2 ENTRIES PERTAIN SUCCESSIVELY TO THESE CONTROL WORDS          BLOCK 72
C          UPPERC,    BLOCKD,    DISSIM,    DISSIM,    SIMILA,          BLOCK 73
C          SIMILA,    CONFIG,    COMPUT,    PRIMAR,    SECOND,          BLOCK 74
C          R     ,    ITERAT,    FIX   ,    SRATST,    STOP  ,          BLOCK 75
C          PRE-IT,    COORDI                                            BLOCK 76
      DATA  TAB2/                                                       BLOCK 77
     1  4623,    3,   10,    5,                                         BLOCK 78
     2  8683,    5,    1,    5,                                         BLOCK 79
     3 15096,    6,   12,    2,                                         BLOCK 80
     4 15096,    3,   11,   10,                                         BLOCK 81
     5  6868,    6,   12,    2,                                         BLOCK 82
     6  6868,    3,   11,   11,                                         BLOCK 83
     7 14846,    6,   12,    5,                                         BLOCK 84
     8 11754,    6,   12,    6,                                         BLOCK 85
     9   765,    3,   13,    1,                                         BLOCK 86
     A  4119,    3,   13,    2,                                         BLOCK 87
     B 16326,    2,   14,    1,                                         BLOCK 88
     C 15170,    1,   15,    1,                                         BLOCK 89
     D  1855,    1,   28,    1,                                         BLOCK 90
     E  3720,    2,   16,    1,                                         BLOCK 91
     F 13171,    6,   12,    7,                                         BLOCK 92
     G 10738,    1,   39,    1,                                         BLOCK 93
     H  7261,    5,    1,   11/                                         BLOCK 94
C                                                                       BLOCK 95
C     TAB3 ENTRIES PERTAIN SUCCESSIVELY TO THESE CONTROL WORDS          BLOCK 96
C          WEIGHT,     WFUNCT,   SFORM1,    SFORM2,    CARDS,           BLOCK 97
C          NOCARD,    SPLIT ,    DATA  ,    SAVE  ,    RANDOM,          BLOCK 98
C          PRINT ,    DFUNCT,    REGRES,    DCONST,    WCONST,          BLOCK 99
C          PLOT                                                         BLOCK100
      DATA  TAB3/                                                       BLOCK101
     1 14543,    6,   12,    3,                                         BLOCK102
     2 11427,    6,   12,    4,                                         BLOCK103
     3  5318,    3,   17,    1,                                         BLOCK104
     4  6181,    3,   17,    2,                                         BLOCK105
     5  6927,    3,   18,    1,                                         BLOCK106
     6 16009,    3,   18,    2,                                         BLOCK107
     7  1966,    5,    1,    3,                                         BLOCK108
     8  6675,    6,   12,    2,                                         BLOCK109
     9  9189,    5,    1,    7,                                         BLOCK110
     A  8330,    1,   20,    1,                                         BLOCK111
     B  2775,    5,    1,    4,                                         BLOCK112
     C 10575,    6,   12,    8,                                         BLOCK113
     D 14439,    5,    1,    6,                                         BLOCK114
     E   267,    2,   29,    5,                                         BLOCK115
     F  1119,    2,   34,    5,                                         BLOCK116
     G  6878,    5,    1,    8/                                         BLOCK117
C                                                                       BLOCK118
C     TAB4 ENTRIES PERTAIN SUCCESSIVELY TO THESE CONTROL WORDS          BLOCK119
C          ABSENT,    PRESEN,    BYROWS,    BYGROU,    BYDECK,          BLOCK120
C          NOMORE,    NOMORE,    DATA  ,    DISTAN,    HISTOR,          BLOCK121
C          NODATA,    NODIST,    NOHIST,    NO    ,    YES   ,          BLOCK122
C          MONOTO,    ASCEND,    DESCEN,    POLYNO                      BLOCK123
      DATA  TAB4/                                                       BLOCK124
     1  4258,    3,    9,    2,                                         BLOCK125
     2  2575,    3,    9,    1,                                         BLOCK126
     3  7761,    3,   19,  100,                                         BLOCK127
     4 11782,    3,   19,  200,                                         BLOCK128
     5 11288,    3,   19,  300,                                         BLOCK129
     6  5274,    3,   19,  400,                                         BLOCK130
     7  5274,    3,   19,    2,                                         BLOCK131
     8  6675,    3,   21,    2,                                         BLOCK132
     9  9824,    3,   22,    2,                                         BLOCK133
     A 12801,    3,   24,    2,                                         BLOCK134
     B 16057,    3,   21,    1,                                         BLOCK135
     C  5863,    3,   22,    1,                                         BLOCK136
     D  9395,    3,   24,    1,                                         BLOCK137
     E  9622,    3,   10,  100,                                         BLOCK138
     F 11125,    3,   10,  200,                                         BLOCK139
     G 15634,    3,   23,    0,                                         BLOCK140
     H  7782,    3,   23,   10,                                         BLOCK141
     I 11188,    3,   23,   11,                                         BLOCK142
     J  3756,    3,   26,    1/                                         BLOCK143
C                                                                       BLOCK144
C     TAB5 ENTRIES PERTAIN SUCCESSIVELY TO THESE CONTROL WORDS          BLOCK145
C          POLYNO,    CONSTA,    NOCONS,    MULTIV,    MULTIV,          BLOCK146
C          DATA  ,    CONFIG,    CONFIG,    SCATTE,    ALL   ,          BLOCK147
C          SOME  ,    NONE  ,    ALL   ,    SOME  ,    NONE  ,          BLOCK148
C          ROTATE,    STANDA,    AS-IS                                  BLOCK149
      DATA  TAB5/                                                       BLOCK150
     1  3756,    1,   23,    1,                                         BLOCK151
     2 14387,    3,   26,  100,                                         BLOCK152
     3  5018,    3,   26,  200,                                         BLOCK153
     4  3608,    3,   26,    0,                                         BLOCK154
     5  3608,    1,   23,    1,                                         BLOCK155
     6  6675,    3,   25,    2,                                         BLOCK156
     7 14846,    3,   27,    3,                                         BLOCK157
     8 14846,    5,    1,    9,                                         BLOCK158
     9 11109,    5,    1,   10,                                         BLOCK159
     A  5766,    3,   40,    2,                                         BLOCK160
     B 13660,    3,   40,    1,                                         BLOCK161
     C 10008,    3,   40,    0,                                         BLOCK162
     D  5766,    3,   41,    2,                                         BLOCK163
     E 13660,    3,   41,    1,                                         BLOCK164
     F 10008,    3,   41,    0,                                         BLOCK165
     G  4503,    3,   42,    2,                                         BLOCK166
     H  3418,    3,   42,    1,                                         BLOCK167
     I  9499,    3,   42,    0/                                         BLOCK168
C                                                                       BLOCK169
      DATA PTID /1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,1HI,1HJ,1HK,1HL,1HM,   BLOCK170
     1           1HN,1HO,1HP,1HQ,1HR,1HS,1HT,1HU,1HV,1HW,1HX,1HY,       BLOCK171
     2           1HZ,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,1HA,1HB,1HC,   BLOCK172
     3           1HD,1HE,1HF,1HG,1HH,1HI,1HJ,1HK,1HL,1HM,1HN,1HO,       BLOCK173
     4           1HP,1HQ,1HR,1HS,1HT,1HU,1HV,1HW,1HX,1HY,1HZ,1H1,1H2,   BLOCK174
     5           1H3,1H4,1H5,1H6,1H7,1H8,1H9,1HA,1HB,1HC,1HD,1HE,       BLOCK175
     6           1HF,1HG,1HH,1HI,1HJ,1HK,1HL,1HM,1HN,1HO,1HP,1HQ,       BLOCK176
     7           1HR,1HS,1HT,1HU,1HV,1HW,1HX,1HY,1HZ,1H1,1H2,1H3,1H4/   BLOCK177
      DATA ITEM/101*1H /                                                BLOCK178
       END                                                              BLOCK179
$      FORTREX STAB
CBSEC1                                                                  BSEC1  1
      SUBROUTINE BSEC1(IB,N,ETA,A,B,W,E,S)                              BSEC1  2
C BSEC1            FOR KYST                  JANUARY,1973               BSEC1  3
C WRITTEN AND TECHNIQUES ADVOCATED IN                                   BSEC1  4
C REFERENCE IMPLEMENTED FOR SGEV                                        BSEC1  5
C BY P.BUSINGER                              JANUARY,1972               BSEC1  6
      DIMENSION A(1),B(1),W(1)                                          BSEC1  7
      REAL ETA,E,S                                                      BSEC1  8
C                                                                       BSEC1  9
C W.KAHAN, ACCURATE EIGENVALUES OF A SYM-                               BSEC1 10
C METRIC TRIDIAGONAL MATRIX. TECH.REPORT                                BSEC1 11
C NO.CS41, JULY 22,1966, COMP.SC.DEPT.,                                 BSEC1 12
C STANFORD UNIVERSITY.                                                  BSEC1 13
C                                                                       BSEC1 14
C BSEC1 IS CALLED BY SGEV, AND                                          BSEC1 15
C USES IB BISECTION STEPS TO FIND THE ALGEBRAICALLY                     BSEC1 16
C GREATEST EIGENVALUE E OF THE SYMMETRIC                                BSEC1 17
C TRIDIAGONAL MATRIX WHOSE DIAGONAL ELEMENTS ARE                        BSEC1 18
C         A(1) THROUGH A(N)                                             BSEC1 19
C AND WHOSE SUB-AND-SUPERDIAGONAL ELEMENTS ARE                          BSEC1 20
C         B(1) THROUGH B(N-1)                                           BSEC1 21
C FURTHERMORE,                                                          BSEC1 22
C         W(1) THROUGH W(N)                                             BSEC1 23
C ARE SET EQUAL TO THE DIAGONAL ELEMENTS OF THE                         BSEC1 24
C MATRIX U IN THE LU-DECOMPOSITION REQUIRED                             BSEC1 25
C FOR SUBSEQUENT INVERSE ITERATION. A AND B                             BSEC1 26
C ARE RESCALED, THE CHANGE IF SCALE BEING                               BSEC1 27
C REFLECTED IN S. ETA IS LEAST POSITIVE MA-                             BSEC1 28
C CHINE NUMBER.                                                         BSEC1 29
C                                                                       BSEC1 30
      REAL BB,D,T,X                                                     BSEC1 31
      IF(N.GT.1)GOTO 5                                                  BSEC1 32
      E=A(1)                                                            BSEC1 33
      W(1)=ETA                                                          BSEC1 34
      GOTO 70                                                           BSEC1 35
    5 B(N)=0.E0                                                         BSEC1 36
      T=0.E0                                                            BSEC1 37
      DO 10 I=1,N                                                       BSEC1 38
         T=AMAX1(T,ABS(A(I)))                                           BSEC1 39
   10    T=AMAX1(T,ABS(B(I)))                                           BSEC1 40
      IF(T.NE.0.E0)GOTO 20                                              BSEC1 41
      E=0.E0                                                            BSEC1 42
      DO 19 I=1,N                                                       BSEC1 43
   19    W(I)=ETA                                                       BSEC1 44
      GOTO 70                                                           BSEC1 45
   20 S=S*T                                                             BSEC1 46
      DO 30 I=1,N                                                       BSEC1 47
         A(I)=A(I)/T                                                    BSEC1 48
   30    B(I)=B(I)/T                                                    BSEC1 49
      E=3.E0                                                            BSEC1 50
      D=-E                                                              BSEC1 51
      DO 50 K=1,IB                                                      BSEC1 52
         X=(D+E)/2.E0                                                   BSEC1 53
         U=1.E0                                                         BSEC1 54
         NU=0                                                           BSEC1 55
         BB=0.E0                                                        BSEC1 56
         DO 40 I=1,N                                                    BSEC1 57
            U=(A(I)-BB/U)-X                                             BSEC1 58
            IF(U.GT.ETA)GOTO 35                                         BSEC1 59
            U=AMIN1(U,-ETA)                                             BSEC1 60
            NU=NU+1                                                     BSEC1 61
   35       BB=B(I)**2                                                  BSEC1 62
   40       W(I)=U                                                      BSEC1 63
         IF(NU.LT.N)D=X                                                 BSEC1 64
         IF(NU.EQ.N)E=X                                                 BSEC1 65
   50    CONTINUE                                                       BSEC1 66
      IF(NU.EQ.N)GOTO 70                                                BSEC1 67
      U=1.E0                                                            BSEC1 68
      BB=0.E0                                                           BSEC1 69
      DO 60 I=1,N                                                       BSEC1 70
         U=(A(I)-BB/U)-E                                                BSEC1 71
         IF(U.GT.ETA)GOTO 55                                            BSEC1 72
         U=AMIN1(U,-ETA)                                                BSEC1 73
   55    BB=B(I)**2                                                     BSEC1 74
   60    W(I)=U                                                         BSEC1 75
   70 RETURN                                                            BSEC1 76
      END                                                               BSEC1 77
$      FORTREX STAB
CCACT                                                                   CACT   1
       SUBROUTINE CCACT                                                 CACT   2
C CCACT            FOR KYST                  JANUARY,1973               CACT   3
C WRITTEN BY J.KRUSKAL                                                  CACT   4
C MODIFIED FOR KYST BY J.KRUSKAL AND J.SEERY JANUARY,1973               CACT   5
C      CCACT--CONTROL CARD ACTIVITY. READS AND INTERPRETS CONTROL CARDS.CACT   6
C                                                                       CACT   7
      REAL MAP(40),DUMMY(27),ATAB(4,86)                                 CACT   8
      INTEGER TAB(4,86),IFAC(6)                                         CACT   9
      INTEGER    MTAB(13)                                               CACT  10
       EQUIVALENCE                         (ATAB,TAB)                   CACT  11
C                                                                       CACT  12
       INTEGER                         LPAR(   42 )                     CACT  13
       REAL                             PAR(   42 )                     CACT  14
      EQUIVALENCE (MAP(1),DUMMY(1)), (MAP(28),CHTAB(1))                 CACT  15
       EQUIVALENCE           (LPAR,PAR)                                 CACT  16
C                                                                       CACT  17
       INTEGER  IPARAM(5)                                               CACT  18
       REAL    C(81)                                                    CACT  19
      REAL WORD(18)                                                     CACT  20
       REAL CHTAB(13)                                                   CACT  21
C                                                                       CACT  22
      INTEGER DFLTSW                                                    CACT  23
       REAL BLANK, DOT                                                  CACT  24
       INTEGER BLSW,NUMSW,DECSW,TYPE,XTYPE,T,TA, PARNO, TABNO           CACT  25
C                                                                       CACT  26
       COMMON  /MDSCL1/  LREAD, LPRINT, LPUNCH, LSCRAT                  CACT  27
       COMMON  /MDSCL2/  LPAR                                           CACT  28
       COMMON  /MDSCL3/  MTAB,  TAB                                     CACT  29
C                                                                       CACT  30
        DATA    BLANK, DOT, EQUALS, COMMA, C(81) /1H ,1H.,1H=,1H,,1H  / CACT  31
      DATA DOLLAR/1H$/                                                  CACT  32
      DATA MODP,IFAC/ 16381, 907, 887, 883, 881, 877, 863/              CACT  33
C                                                                       CACT  34
      DATA  MAP/1H ,1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,1HI,1HJ,1HK,1HL,    CACT  35
     .          1HM,1HN,1HO,1HP,1HQ,1HR,1HS,1HT,1HU,1HV,1HW,1HX,1HY,    CACT  36
     .          1HZ,1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,1H-,1H+,1H./CACT  37
C                                                                       CACT  38
C                                                                       CACT  39
C                                                                       CACT  40
 1     FORMAT(80A1)                                                     CACT  41
 11    FORMAT(1H0,80A1)                                                 CACT  42
 2     FORMAT(1X,I1)                                                    CACT  43
 3    FORMAT(1X,80A1)                                                   CACT  44
 4     FORMAT(1X,I18)                                                   CACT  45
 5     FORMAT(1X,F18.3)                                                 CACT  46
 9     FORMAT(51H  CCACT DIAGNOSTIC. CONTROL CARD ABOVE IS IMPROPER.)   CACT  47
  99    FORMAT( 12H ITEM NUMBER, I5, 20H   HAS EXPECTED TYPE ,  I5,     CACT  48
     1 18H  AND ACTUAL TYPE , I5, 16H . THIS ITEM IS , 18A1)            CACT  49
 98   FORMAT(12H ITEM NUMBER,I5,39H  HAS ILLEGAL CHARACTER. THIS ITEM ISCACT  50
     .  ,18A1)                                                          CACT  51
C                                                                       CACT  52
C      READ AND PRINT CONTROL CARD                                      CACT  53
C                                                                       CACT  54
 100   READ  (LREAD,1)  (C(I),I=1,80)                                   CACT  55
       WRITE(LPRINT,11) (C(I),I=1,80)                                   CACT  56
C      ANALYZE CONTROL CARD INTO TOKENS                                 CACT  57
C                                                                       CACT  58
C              EACH  TOKEN  IS DELIMITED BY BLANKS                      CACT  59
C              THERE ARE  FOUR TYPES OF TOKENS.                         CACT  60
C              ALPHABETIC, INTEGER, DECIMAL, AND DEFAULT                CACT  61
C              ALPHABETIC UNLESS FIRST CHARACTER IS DIGIT OR DEC POINT  CACT  62
C                      OR PLUS OR MINUS  OR DOLLAR SIGN                 CACT  63
C              DECIMALS       DISTINGUISHED BY DECIMAL POINT            CACT  64
C              DEFAULT = $                                              CACT  65
C                                                                       CACT  66
       REWIND LSCRAT                                                    CACT  67
       BLSW=1                                                           CACT  68
       NUMSW=1                                                          CACT  69
       DECSW=1                                                          CACT  70
      DFLTSW=1                                                          CACT  71
       K=0                                                              CACT  72
C                                                                       CACT  73
       DO 300 I=1,81                                                    CACT  74
C                                                                       CACT  75
       X=C(I)                                                           CACT  76
       GO TO (110,120), BLSW                                            CACT  77
C                                                                       CACT  78
 110   IF(X.EQ.BLANK  .OR.   X.EQ.EQUALS  .OR.   X.EQ. COMMA ) GO TO 300CACT  79
       BLSW=2                                                           CACT  80
       JA=I                                                             CACT  81
       DO 115 KX=1,13                                                   CACT  82
 115   IF(X.EQ.CHTAB(KX))NUMSW=2                                        CACT  83
      IF(X.EQ.DOLLAR) DFLTSW=2                                          CACT  84
       IF(X.EQ.DOT) DECSW=2                                             CACT  85
       GO TO 300                                                        CACT  86
C                                                                       CACT  87
 120   IF(X.EQ.BLANK  .OR.   X.EQ.EQUALS  .OR.   X.EQ. COMMA ) GO TO 130CACT  88
       IF(X.EQ.DOT) DECSW=2                                             CACT  89
       GO TO 300                                                        CACT  90
C                                                                       CACT  91
 130   K=K+1                                                            CACT  92
       JB = MIN0 (I-1,JA+16)                                            CACT  93
       JC=18-(JB-JA+1)                                                  CACT  94
       TYPE=1                                                           CACT  95
       IF(NUMSW.EQ.2) TYPE=NUMSW+DECSW-1                                CACT  96
      IF(DFLTSW.EQ.2) TYPE=4                                            CACT  97
C                                                                       CACT  98
       WRITE (LSCRAT,2) TYPE                                            CACT  99
      IF(TYPE.EQ.4) GO TO 160                                           CACT 100
       GO TO (140,150),NUMSW                                            CACT 101
 140  WRITE (LSCRAT,3) (C(J),J=JA,JB), (BLANK,J=1,JC)                   CACT 102
       GO TO 160                                                        CACT 103
 150  WRITE (LSCRAT,3) (BLANK,J=1,JC), (C(J),J=JA,JB)                   CACT 104
 160   BLSW=1                                                           CACT 105
       NUMSW=1                                                          CACT 106
       DECSW=1                                                          CACT 107
      DFLTSW=1                                                          CACT 108
       GO TO 300                                                        CACT 109
C                                                                       CACT 110
 300   CONTINUE                                                         CACT 111
C                                                                       CACT 112
C      ANALYZE TOKENS AND SET PARAMETER VALUES ACCORDINGLY              CACT 113
C                                                                       CACT 114
       KB=K                                                             CACT 115
       IF(KB.EQ.0) RETURN                                               CACT 116
       REWIND LSCRAT                                                    CACT 117
       XTYPE=1                                                          CACT 118
       IPARAM(1) = 1                                                    CACT 119
      NT=0                                                              CACT 120
      NTOKNS=0                                                          CACT 121
      NOPCC=0                                                           CACT 122
C                                                                       CACT 123
       DO 1000 K=1,KB                                                   CACT 124
C                                                                       CACT 125
       READ  (LSCRAT,2) TYPE                                            CACT 126
      IF( (XTYPE.EQ.1) .AND. (TYPE.NE.1) ) GO TO 995                    CACT 127
      IF( (XTYPE.NE.TYPE). AND. (TYPE.NE.4) ) GO TO  995                CACT 128
 350  GO TO (400,410,420,430), TYPE                                     CACT 129
C                                                                       CACT 130
 400  READ (LSCRAT,3) (WORD(L),L=1,18)                                  CACT 131
C     CONVERT FIRST SIX LETTERS OF ALPHABETIC WORD INTO CODE NUMBER.    CACT 132
C     LARGEST POSSIBLE CODE NUMBER IS LESS THAN 2**15.                  CACT 133
      ICODE=0                                                           CACT 134
      DO 408 M1=1,6                                                     CACT 135
      X=WORD(M1)                                                        CACT 136
      DO 402 M2=1,38                                                    CACT 137
      IF( X.EQ.MAP(M2) ) GO TO 404                                      CACT 138
 402  CONTINUE                                                          CACT 139
C     ILLEGAL CHARACTER                                                 CACT 140
      WRITE (LPRINT,9)                                                  CACT 141
      WRITE (LPRINT,98) K,(WORD(L),L=1,18)                              CACT 142
      GO TO 999                                                         CACT 143
 404  M3=(M2-1)*IFAC(M1)                                                CACT 144
      IF( M3.GE.MODP ) M3=M3-MODP                                       CACT 145
      ICODE=ICODE+M3                                                    CACT 146
      IF( ICODE.GE.MODP ) ICODE=ICODE-MODP                              CACT 147
 408  CONTINUE                                                          CACT 148
       GO TO 510                                                        CACT 149
C                                                                       CACT 150
 410   READ  (LSCRAT,4) INTPAR                                          CACT 151
      ISUB=PARNO+NT                                                     CACT 152
      LPAR(ISUB)=INTPAR                                                 CACT 153
       GO TO 430                                                        CACT 154
C                                                                       CACT 155
 420   READ  (LSCRAT,5) DECPAR                                          CACT 156
      ISUB=PARNO+NT                                                     CACT 157
      PAR(ISUB)=DECPAR                                                  CACT 158
C                                                                       CACT 159
 430  NT=NT+1                                                           CACT 160
      IF(NT.EQ.NTOKNS) XTYPE=1                                          CACT 161
       GO TO 1000                                                       CACT 162
C                                                                       CACT 163
 510   TABNO = IPARAM(1)                                                CACT 164
       MA = MTAB(TABNO)                                                 CACT 165
       MB = MTAB(TABNO+1) - 1                                           CACT 166
C                                                                       CACT 167
       LMISSW = 0                                                       CACT 168
       DO 700 M=MA,MB                                                   CACT 169
      IF( ICODE.NE.TAB(1,M) ) GO TO 700                                 CACT 170
       LMISSW = 1                                                       CACT 171
 600   XTYPE=1                                                          CACT 172
      NT=0                                                              CACT 173
       IPARAM(1) = 1                                                    CACT 174
      PARNO=TAB(3,M)                                                    CACT 175
      LTEMP=TAB(2,M)                                                    CACT 176
      NTOKNS=TAB(4,M)                                                   CACT 177
      IF(LTEMP.GT.2) NTOKNS=0                                           CACT 178
      GO TO(610,620,630,640,650,660), LTEMP                             CACT 179
C                                                                       CACT 180
C      NAME OF INTEGER PARAMETER                                        CACT 181
 610   XTYPE=2                                                          CACT 182
       GO TO 700                                                        CACT 183
C                                                                       CACT 184
C      NAME OF DECIMAL PARAMETER                                        CACT 185
 620   XTYPE=3                                                          CACT 186
       GO TO 700                                                        CACT 187
C                                                                       CACT 188
C      IMPLICITLY SPECIFIED INTEGER PARAMETER                           CACT 189
C      A SINGLE IMPLICITLY SPECIFIED PARAMTER CAN IN PRINCIPLE HOLD     CACT 190
C      AS MANY AS THREE PARAMTERS WHICH THE MAIN PROGRAM CONSIDERS      CACT 191
C      CONCEPTUALLY DISTINCT.  ONE IS IN THE UNITS POSITION, ANOTHER IN CACT 192
C      IN THE HUNDREDS POSITION, AND ANOTHER IN THE TEN-THOUSANDS       CACT 193
C      POSITION.                                                        CACT 194
 630   LP = 1                                                           CACT 195
      LT=TAB(4,M)                                                       CACT 196
       IF(LT.GE.100) LP = 100                                           CACT 197
       IF(LT.GE.10000) LP = 10000                                       CACT 198
       LQ = 100*LP                                                      CACT 199
       LA = LPAR(PARNO)                                                 CACT 200
      LPAR(PARNO)= LA-(MOD(LA,LQ)/LP)*LP +TAB(4,M)                      CACT 201
       GO TO 700                                                        CACT 202
C                                                                       CACT 203
C      IMPLICITLY SPECIFIED DECIMAL PARAMETER                           CACT 204
 640  PAR(PARNO)=ATAB(4,M)                                              CACT 205
       GO TO 700                                                        CACT 206
C                                                                       CACT 207
C      INTERNAL PARAMETER OF CCACT PROGRAM                              CACT 208
 650  IPARAM(PARNO)=TAB(4,M)                                            CACT 209
       GO TO 700                                                        CACT 210
C                                                                       CACT 211
C     IMPLICITLY SPECIFIED INTEGER PARAMETER                            CACT 212
C     DETERMINES PROGRAM FLOW UPON RETURN FROM CCACT--SO ONLY ONE OF    CACT 213
C     THIS TYPE PER CONTROL CARD ALLOWED.                               CACT 214
 660  IF(NOPCC)9999,665,996                                             CACT 215
 665  NOPCC=1                                                           CACT 216
      LPAR(PARNO)=TAB(4,M)                                              CACT 217
      GO TO 700                                                         CACT 218
C                                                                       CACT 219
 700   CONTINUE                                                         CACT 220
      IF(LMISSW.EQ.0) GO TO 994                                         CACT 221
 1000  CONTINUE                                                         CACT 222
      IF(NT.LT.NTOKNS) GO TO 997                                        CACT 223
C                                                                       CACT 224
 1001  RETURN                                                           CACT 225
C                                                                       CACT 226
 994  WRITE(LPRINT,9)                                                   CACT 227
      WRITE(LPRINT,998) (WORD(L),L=1,18)                                CACT 228
 998  FORMAT(27H UNDEFINED CONTROL PHRASE  ,18A1)                       CACT 229
      GO TO 999                                                         CACT 230
 995  WRITE (LPRINT,9)                                                  CACT 231
      WRITE (LPRINT,99) K,XTYPE,TYPE, (WORD(L),L=1,18)                  CACT 232
      GO TO 999                                                         CACT 233
 996  WRITE(LPRINT,9)                                                   CACT 234
      WRITE(LPRINT,6)                                                   CACT 235
 6    FORMAT(69H0NO MORE THAN ONE CONTROL PHRASE OF TYPE 6 ALLOWED ON A CACT 236
     .CONTROL CARD.)                                                    CACT 237
      GO TO 999                                                         CACT 238
 997  WRITE(LPRINT,9)                                                   CACT 239
      WRITE(LPRINT,97)                                                  CACT 240
 97   FORMAT(47H CONTROL PHRASE NOT COMPLETED ON A SINGLE CARD.)        CACT 241
 999  CALL EXIT                                                         CACT 242
C                                                                       CACT 243
 9999 WRITE(LPRINT,7)                                                   CACT 244
 7    FORMAT(23H0ERROR EXIT FROM CCACT.)                                CACT 245
      GO TO 999                                                         CACT 246
       END                                                              CACT 247
$      FORTREX STAB
CCONFIG                                                                 CONFIG 1
      SUBROUTINE CONFIG(N,ND)                                           CONFIG 2
C CONFIG           FOR KYST                  JANUARY,1973               CONFIG 3
C WRITTEN BY F.YOUNG                                                    CONFIG 4
C MODIFIED FOR KYST BY J.KRUSKAL AND J.SEERY JANUARY,1973               CONFIG 5
C ROUTINE TO OBTAIN  A METRIC MULTIDIMENSIONAL SCALING SOLUTION         CONFIG 6
C                                                                       CONFIG 7
      DIMENSION DATA(1800),WW(1800),IJ(1800),X(100,6)                   CONFIG 8
      DIMENSION STORE(494),RWMEAN(100),ROOTS(6),XBEST(100,6),           CONFIG 9
     .     DATAIN(5430)                                                 CONFIG10
      REAL MEAN                                                         CONFIG11
      COMMON /KYST1/ DATA,WW,IJ,X                                       CONFIG12
      COMMON /KYST2/ STORE,RWMEAN,ROOTS,XBEST,DATAIN                    CONFIG13
C                                                                       CONFIG14
C PARAMETERS                                                            CONFIG15
C X      - CONFIGURATION AT OUTPUT                                      CONFIG16
C N      - NUMBER OF POINTS IN SCALING                                  CONFIG17
C ND     - NUMBER OF DIMENSIONS OF SOLUTION                             CONFIG18
C DATAIN - ON INPUT, DATA FROM INICON.  DESTROYED BY CONFIG             CONFIG19
C RWMEAN - ROW AVERAGES                                                 CONFIG20
C ROOTS - ON RETURN FROM SGEV, THE FIRST ND EIGENVALUES                 CONFIG21
C                                                                       CONFIG22
      ISUB(I1,I2)=((I2-1)*(I2-2))/2+I1                                  CONFIG23
      MEAN=0.0                                                          CONFIG24
      FN=N                                                              CONFIG25
C                                                                       CONFIG26
C    COMPUTE ROW,COLUMN, AND OVERALL AVERAGES OF DATIN VALUES SQUARED   CONFIG27
C                                                                       CONFIG28
      DO 130 J=1,N                                                      CONFIG29
      RWMEAN(J)=0.0                                                     CONFIG30
      DO 129 I=1,N                                                      CONFIG31
      IF(I-J) 126,129,127                                               CONFIG32
 126  IJC=ISUB(I,J)                                                     CONFIG33
      GO TO 128                                                         CONFIG34
 127  IJC=ISUB(J,I)                                                     CONFIG35
 128  RWMEAN(J)=DATAIN(IJC)**2+RWMEAN(J)                                CONFIG36
 129  CONTINUE                                                          CONFIG37
 130  MEAN=RWMEAN(J)+MEAN                                               CONFIG38
      MEAN=MEAN/FN**2                                                   CONFIG39
C                                                                       CONFIG40
C    CONVERT DATAIN VALUES TO SCALAR-PRODUCT-LIKE VALUES BY THE         CONFIG41
C    HOUSEHOLDER-TORGERSON PROCEDURE.  WORK FROM LAST ROW AND           CONFIG42
C    COLUMN BACKWARDS.                                                  CONFIG43
C                                                                       CONFIG44
      DO 135 JJ=1,N                                                     CONFIG45
      J=N-JJ+1                                                          CONFIG46
      JM1=J-1                                                           CONFIG47
      J1=(JM1*(J-2))/2                                                  CONFIG48
      J2=(J*JM1)/2                                                      CONFIG49
      DO 135 II=1,J                                                     CONFIG50
      I=J-II+1                                                          CONFIG51
      IJ1=J1+I                                                          CONFIG52
      IJ2=J2+I                                                          CONFIG53
      IF(I.EQ.J) GO TO 133                                              CONFIG54
C    HOUSEHOLDER-TORGERSON FORMULAS                                     CONFIG55
      DATAIN(IJ2)=((RWMEAN(I)+RWMEAN(J))/FN-MEAN-DATAIN(IJ1)**2)/2.0    CONFIG56
      GO TO 135                                                         CONFIG57
 133  DATAIN(IJ2)=((2.0*RWMEAN(I))/FN-MEAN)/2.0                         CONFIG58
 135  CONTINUE                                                          CONFIG59
C                                                                       CONFIG60
C OBTAIN EIGENROOTS AND EIGENVECTORS                                    CONFIG61
      CALL SGEV(N,ND,X)                                                 CONFIG62
C                                                                       CONFIG63
C COMPUTE CONFIGURATION                                                 CONFIG64
C                                                                       CONFIG65
      DO 136 I=1,ND                                                     CONFIG66
 136  ROOTS(I)=SQRT(ABS(ROOTS(I)))                                      CONFIG67
      DO 137 I=1,N                                                      CONFIG68
      DO 137 J=1,ND                                                     CONFIG69
 137  X(I,J)=X(I,J)*ROOTS(J)                                            CONFIG70
      RETURN                                                            CONFIG71
      END                                                               CONFIG72
$      FORTRAN
CDATAPR                                                                 DATAPR 1
      SUBROUTINE DATAPR(GRNO,MM,NOGRPS,LCODE)                           DATAPR 2
C DATAPR           FOR KYST                  JANUARY,1973               DATAPR 3
C WRITTEN BY J.KRUSKAL                                                  DATAPR 4
C                                                                       DATAPR 5
C           DATA PRINT                                                  DATAPR 6
C                                                                       DATAPR 7
      DIMENSION DATA(1800),WW(1800),IJ(1800),X(100,6)                   DATAPR 8
      DIMENSION GR(100,6),GL(100,6),STORE(1830),DIST(1800),DHAT(1800),  DATAPR 9
     .     P(50)                                                        DATAPR10
      INTEGER Q(50),GRNO(1)                                             DATAPR11
      EQUIVALENCE(P(1),Q(1),STORE(1))                                   DATAPR12
      COMMON /KYST1/ DATA,WW,IJ,X                                       DATAPR13
      COMMON /KYST2/ GR,GL,STORE,DIST,DHAT                              DATAPR14
       COMMON   /MDSCL1/ LREAD, LPRINT, LPUNCH, LSCRAT                  DATAPR15
      DATA KPACK,IJPACK /10000,100/                                     DATAPR16
C                                                                       DATAPR17
 1     FORMAT( 1H1, 4(30H   I   J   DATA   WGHT GP  NO  ) )             DATAPR18
 2     FORMAT(1H1, 3(40H   I   J   DATA   DIST DHAT WGHT GP  NO ) )     DATAPR19
 3     FORMAT( 1X, 5( 2I4, F9.3,        F5.2, I3, I4, 1H, ) )           DATAPR20
 4     FORMAT( 1X, 5( 2I4, F9.3, 2F5.2, F5.2, I3, I4, 1H, ) )           DATAPR21
C                                                                       DATAPR22
 100   NCOPA=3                                                          DATAPR23
       IF(LCODE.EQ.1) NCOPA=4                                           DATAPR24
C                                                                       DATAPR25
       MMLEFT=MM                                                        DATAPR26
 200   MMNOW=MIN0( MMLEFT, 50*NCOPA )                                   DATAPR27
       MMUSED=MM-MMLEFT                                                 DATAPR28
       MMLEFT=MMLEFT-MMNOW                                              DATAPR29
       IF(LCODE.EQ.1) WRITE (LPRINT,1)                                  DATAPR30
       IF(LCODE.EQ.2) WRITE (LPRINT,2)                                  DATAPR31
C                                                                       DATAPR32
       DO 300 LINENO=1,50                                               DATAPR33
       NCONOW=(MMNOW-LINENO+50)/50                                      DATAPR34
       IF(NCONOW.LE.0) GO TO 310                                        DATAPR35
C                                                                       DATAPR36
       L=0                                                              DATAPR37
       DO 270 LC=1,NCONOW                                               DATAPR38
       M = MMUSED + 50*(LC-1) + LINENO                                  DATAPR39
C                                                                       DATAPR40
      IJM=MOD(IJ(M),KPACK)                                              DATAPR41
      I=IJM/IJPACK+1                                                    DATAPR42
      J=MOD(IJM,IJPACK)+1                                               DATAPR43
C                                                                       DATAPR44
       L=L+1                                                            DATAPR45
       Q(L)=I                                                           DATAPR46
       L=L+1                                                            DATAPR47
       Q(L)=J                                                           DATAPR48
C                                                                       DATAPR49
       L=L+1                                                            DATAPR50
       P(L)=DATA(M)                                                     DATAPR51
C                                                                       DATAPR52
       IF(LCODE.EQ.1) GO TO 250                                         DATAPR53
C                                                                       DATAPR54
       L=L+1                                                            DATAPR55
       P(L)=DIST(M)                                                     DATAPR56
       L=L+1                                                            DATAPR57
       P(L)=DHAT(M)                                                     DATAPR58
C                                                                       DATAPR59
 250   CONTINUE                                                         DATAPR60
C                                                                       DATAPR61
       L=L+1                                                            DATAPR62
       P(L)=WW(M)                                                       DATAPR63
C                                                                       DATAPR64
       DO 220 NG=1,NOGRPS                                               DATAPR65
       IF(M.GE.GRNO(NG)) NGM=NG                                         DATAPR66
 220   CONTINUE                                                         DATAPR67
       L=L+1                                                            DATAPR68
       Q(L)=NGM                                                         DATAPR69
C                                                                       DATAPR70
       L=L+1                                                            DATAPR71
       Q(L)=M-GRNO(NGM)+1                                               DATAPR72
C                                                                       DATAPR73
 270   CONTINUE                                                         DATAPR74
C                                                                       DATAPR75
       LL=L                                                             DATAPR76
       IF(LCODE.EQ.1) WRITE (LPRINT,3) (P(L),L=1,LL)                    DATAPR77
       IF(LCODE.EQ.2) WRITE (LPRINT,4) (P(L),L=1,LL)                    DATAPR78
 300   CONTINUE                                                         DATAPR79
C                                                                       DATAPR80
 310   IF(MMLEFT.GT.0) GO TO 200                                        DATAPR81
C                                                                       DATAPR82
       RETURN                                                           DATAPR83
C                                                                       DATAPR84
       END                                                              DATAPR85
$      FORTRAN
CDFLT                                                                   DFLT   1
      SUBROUTINE DFLT(N,A,B,E,X)                                        DFLT   2
C DFLT             FOR KYST                  JANUARY,1973               DFLT   3
C WRITTEN AND DEFLATION TECHNIQUE                                       DFLT   4
C ADVOCATED IN REFERENCE IMPLEMENTED                                    DFLT   5
C FOR SGEV BY P.BUSINGER                     JANUARY,1972               DFLT   6
C                                                                       DFLT   7
C G.GOLUB AND W.KAHAN, CALCULATING THE                                  DFLT   8
C SINGULAR VALUES AND PSEUDO-INVERSE OF A                               DFLT   9
C MATRIX.  J.SIAM NUMER. ANAL. 2,205-224(1965).                         DFLT  10
C DFLT IS CALLED BY SGEV, AND                                           DFLT  11
      REAL A(1),B(1),E,X(1)                                             DFLT  12
C                                                                       DFLT  13
C DEFLATES THE SYMMETRIC TRIDIAGONAL MATRIX                             DFLT  14
C WITH DIAGONAL ELEMENTS A(1) THROUGH A(N)                              DFLT  15
C AND SUB- AND SUPER-DIAGONAL ELEMENTS B(1)                             DFLT  16
C THROUGH B(N-1) INTO A SYMMETRIC TRIDIAGONAL                           DFLT  17
C MATRIX OF ORDER N-1. E IS THE EIGENVALUE                              DFLT  18
C BEING REMOVED AND X AN ASSOCIATED EIGEN-                              DFLT  19
C VECTOR. THE TANGENTS OF HALF THE ROTATION                             DFLT  20
C ANGLES OVERWRITE THE FIRST (N-1) COMPONENTS                           DFLT  21
C OF X.                                                                 DFLT  22
C                                                                       DFLT  23
      REAL H,W,F,V,T,Q,S,C,AINEW                                        DFLT  24
      H=A(1)-E                                                          DFLT  25
      W=B(1)                                                            DFLT  26
      N1=N-1                                                            DFLT  27
      IF(N1.LT.1)GOTO 30                                                DFLT  28
      DO 20 I=1,N1                                                      DFLT  29
         F=X(I)                                                         DFLT  30
         V=X(I+1)                                                       DFLT  31
         IF(AMAX1(ABS(H),ABS(W)).LE.AMAX1(ABS(F),ABS(V)))GOTO 10        DFLT  32
         F=W                                                            DFLT  33
         V=-H                                                           DFLT  34
   10    T=0.E0                                                         DFLT  35
         IF(F*F.EQ.0.E0)GOTO 20                                         DFLT  36
         T=F                                                            DFLT  37
         IF(V.GE.0.E0)T=-T                                              DFLT  38
         T=T/(SQRT(F*F+V*V)+ABS(V))                                     DFLT  39
         S=(T+T)/(1.E0+T*T)                                             DFLT  40
         C=(1.E0-T)*(1.E0+T)/(1.E0+T*T)                                 DFLT  41
         IF(I.GT.1)B(I-1)=C*H+S*W                                       DFLT  42
         F=C*A(I)+S*B(I)                                                DFLT  43
         Q=C*B(I)+S*A(I+1)                                              DFLT  44
         W=S*B(I+1)                                                     DFLT  45
         B(I+1)=C*B(I+1)                                                DFLT  46
         AINEW=F*C+Q*S                                                  DFLT  47
         H=Q*C-F*S                                                      DFLT  48
         A(I+1)=A(I)+A(I+1)-AINEW                                       DFLT  49
         A(I)=AINEW                                                     DFLT  50
         X(I+1)=C*X(I+1)-S*X(I)                                         DFLT  51
   20    X(I)=T                                                         DFLT  52
      B(N-1)=0.E0                                                       DFLT  53
   30 RETURN                                                            DFLT  54
      END                                                               DFLT  55
$      FORTRAN
CDTRAN                                                                  DTRAN  1
      FUNCTION DTRAN(XX)                                                DTRAN  2
C DTRAN            FOR KYST                  JANUARY,1973               DTRAN  3
C WRITTEN FOR KYST BY J.SEERY                JANUARY,1973               DTRAN  4
C                                                                       DTRAN  5
C      DTRAN--DATA TRANSFORMATION.                                      DTRAN  6
      DIMENSION DUMMY(28),DUM(4)                                        DTRAN  7
      COMMON /MDSCL2/ DUMMY,A,B,P,S,T,WCON1,WCON2,WCON3,WCON4,WCON5,DUM DTRAN  8
      DTRAN=S+T*((A+B*XX)**P)                                           DTRAN  9
      RETURN                                                            DTRAN 10
      END                                                               DTRAN 11
$      FORTRAN
CEQSOLV                                                                 EQSOLV 1
       SUBROUTINE EQSOLV(A,B,C,NT,NF)                                   EQSOLV 2
C EQSOLV           FOR KYST                  JANUARY,1973               EQSOLV 3
C WRITTEN BY J.KRUSKAL                                                  EQSOLV 4
CEQSOLV        SOLVES SIMULTANEOUS LINEAR EQUATIONS                     EQSOLV 5
       REAL    A(5,5), B(5), C(5)                                       EQSOLV 6
C                                                                       EQSOLV 7
       NTM1=NT-1                                                        EQSOLV 8
       IF(NT.EQ.NF) GO TO 60                                            EQSOLV 9
       DO 50 II=NF,NTM1                                                 EQSOLV10
       IIP1=II+1                                                        EQSOLV11
       DO 40 I=IIP1,NT                                                  EQSOLV12
       E = A(I,II)/A(II,II)                                             EQSOLV13
       DO 30 J=NF,NT                                                    EQSOLV14
       A(I,J)=A(I,J)-E*A(II,J)                                          EQSOLV15
 30    CONTINUE                                                         EQSOLV16
       B(I)=B(I)-E*B(II)                                                EQSOLV17
 40    CONTINUE                                                         EQSOLV18
 50    CONTINUE                                                         EQSOLV19
C                                                                       EQSOLV20
 60    DO 90 IC=NF,NT                                                   EQSOLV21
       I=NT+NF-IC                                                       EQSOLV22
       E=B(I)                                                           EQSOLV23
       IF(I.EQ.NT) GO TO 80                                             EQSOLV24
       IP1=I+1                                                          EQSOLV25
       DO 70 J=IP1,NT                                                   EQSOLV26
       E=E-C(J)*A(I,J)                                                  EQSOLV27
 70    CONTINUE                                                         EQSOLV28
 80    C(I) = E / A(I,I)                                                EQSOLV29
 90    CONTINUE                                                         EQSOLV30
       RETURN                                                           EQSOLV31
       END                                                              EQSOLV32
$      FORTRAN
CFITM                                                                   FITM   1
      SUBROUTINE FITM(ISTRT,MM,LFITSW)                                  FITM   2
C FITM             FOR KYST                  JANUARY,1973               FITM   3
C WRITTEN BY J.KRUSKAL                                                  FITM   4
C MODIFIED FOR KYST BY J.KRUSKAL AND J.SEERY JANUARY,1973               FITM   5
C                                                                       FITM   6
C      FITM    PERFORMS WEIGHTED-LEAST-SQUARES MONOTONE REGRESSION      FITM   7
C      THIS ROUTINE FINDS THE VALUES FOR  DHAT  WHICH ARE MONOTONIC     FITM   8
C      AND WHICH BEST APPROXIMATE  THE VALUES OF  DIST,                 FITM   9
C      IN THE SENSE THAT THE SUM OF THE SQUARED DEVIATIONS,             FITM  10
C      EACH ONE WEIGHTED BY THE CORRESPONDING VALUE IN  WW ,            FITM  11
C      IS A MINIMUM.                                                    FITM  12
C      IT DIRECTLY USES  SORT.                                          FITM  13
C                                                                       FITM  14
      DIMENSION DATA(1800),IJ(1800),WW(1800),DIST(1800),DHAT(1800),     FITM  15
     .     LBLOCK(1800),X(100,6),DUMMY(30),GR(100,6),GL(100,6)          FITM  16
      INTEGER SDSWIT                                                    FITM  17
       COMMON  /MDSCL1/  LREAD, LPRINT, LPUNCH, LSCRAT                  FITM  18
      COMMON /KYST1/ DATA,WW,IJ,X                                       FITM  19
      COMMON /KYST2/ GR,GL,LBLOCK,DUMMY,DIST,DHAT                       FITM  20
C                                                                       FITM  21
      DATA KPACK/10000/                                                 FITM  22
C                                                                       FITM  23
C      FORM FIRST APPROXIMATION TO CORRECT PARTITON                     FITM  24
C                                                                       FITM  25
C              IF LFITSW 1,  USE PRIMARY APPROACH.  THUS WE SORT EACH   FITM  26
C              BLOCK OF EQUAL VALUES OF DATA ACCORDING TO DIST VALUES.  FITM  27
C              THEN WE CREATE PARTITON INTO BLOCKS OF SIZE 1 TO START.  FITM  28
C                                                                       FITM  29
C              IF LFITSW 2,   USE SECONDARY APPROACH. THUS WE START WITHFITM  30
C              PARTITION INTO BLOCKS OF EQUAL DATA VALUES.              FITM  31
C                                                                       FITM  32
C              WITHIN EACH BLOCK  OF WHATEVER SIZE, THE FIRST DHAT VALUEFITM  33
C      GIVES THE WEIGHTED TOTAL OF THE DIST VALUES IN THAT BLOCK, WHILE FITM  34
C              THE LAST DHAT VALUE GIVES THE SUM OF THE WEIGHTS FOR THATFITM  35
C              BLOCK.  SIMILARLY, WITHIN EACH BLOCK, THE FIRST LBLOCK   FITM  36
C              VALUE AND THE LAST LBLOCK VALUE BOTH GIVE THE SIZE OF THEFITM  37
C              BLOCK.                                                   FITM  38
C      BLOCKS OF SIZE 1 FORM A PARTIAL EXCEPTION. EVERYTHIN IS THE SAME FITM  39
C              EXCEPT THAT THE SUM OF THE W VALUES IS NOT FOUND IN THE  FITM  40
C              LAST DHAT VALUE IN THE BLOCK.                            FITM  41
C                                                                       FITM  42
 100  MA=ISTRT                                                          FITM  43
 110  KSUB=MA-1                                                         FITM  44
 115  KSUB=KSUB+1                                                       FITM  45
      KK=IJ(KSUB)-KPACK                                                 FITM  46
      IF(KK) 115,120,120                                                FITM  47
 120  MB=KSUB                                                           FITM  48
      K=MB-MA+1                                                         FITM  49
       GO TO ( 200, 300 ), LFITSW                                       FITM  50
C                                                                       FITM  51
C              PRIMARY APPROACH                                         FITM  52
C                                                                       FITM  53
  200  IF(K-1) 9999, 220, 210                                           FITM  54
C              IF K 1, SAVE SORTING TIME                                FITM  55
C                                                                       FITM  56
 210  IJ(MB)=IJ(MB)-KPACK                                               FITM  57
       CALL SORT( DIST(MA), K, IJ(MA), DATA(MA),WW(MA),3,+1)            FITM  58
      IJ(MB)=IJ(MB)+KPACK                                               FITM  59
C               SORT  WILL SORT THE K ELEMENTS OF  DIST  IN ALGEBRAIC   FITM  60
C                ORDER, ASCENDING                                       FITM  61
C              THE ELEMENTS OF IJ  AND  DATA  WILL BE REARRANGED IN     FITM  62
C              EXACTLY THE SAME ORDER AS THOSE OF  DIST .               FITM  63
C      ALSO THE ELEMENTS OF  WW .                                       FITM  64
C                                                                       FITM  65
  220  DO 230 M=MA,MB                                                   FITM  66
       DHAT(M) = DIST(M) * WW(M)                                        FITM  67
  230  LBLOCK(M)=1                                                      FITM  68
       GO TO 400                                                        FITM  69
C                                                                       FITM  70
C              SECONDARY APPROACH                                       FITM  71
C                                                                       FITM  72
C                                                                       FITM  73
  300  LBLOCK(MA)=K                                                     FITM  74
       LBLOCK(MB)=K                                                     FITM  75
       TEMP1=0.0                                                        FITM  76
       TEMP2 = 0.0                                                      FITM  77
       DO 310 M=MA,MB                                                   FITM  78
       TEMP1 = TEMP1 + DIST(M) * WW(M)                                  FITM  79
       TEMP2 = TEMP2 + WW(M)                                            FITM  80
 310   CONTINUE                                                         FITM  81
       DHAT(MA)=TEMP1                                                   FITM  82
       IF(K-1) 9999, 400, 320                                           FITM  83
 320   DHAT(MB) = TEMP2                                                 FITM  84
C              PROCEED TO NEXT BOCK. QUERY. IS IT THE LAST              FITM  85
C                                                                       FITM  86
  400  MA=MA+K                                                          FITM  87
       IF(MA-MM-1) 110, 500, 9999                                       FITM  88
C                                                                       FITM  89
C      START MAIN COMPUTATION      *************************************FITM  90
C                                                                       FITM  91
 500  MA=ISTRT                                                          FITM  92
  510  LUD=2                                                            FITM  93
       NSATIS=0                                                         FITM  94
  520  K=LBLOCK(MA)                                                     FITM  95
       MB=MA+K-1                                                        FITM  96
       IF(K-1) 9999, 530, 540                                           FITM  97
 530   WT =WW(MB)                                                       FITM  98
       GO TO 550                                                        FITM  99
 540   WT = DHAT(MB)                                                    FITM 100
 550   DAV = DHAT(MA) / WT                                              FITM 101
       GO TO ( 600, 700 ), LUD                                          FITM 102
C                                                                       FITM 103
C              IS BLOCK DOWN-SATISFIED. IF NOT, MERGE                   FITM 104
C                                                                       FITM 105
 600  IF(MA-ISTRT) 9999,630,610                                         FITM 106
  610  MBD=MA-1                                                         FITM 107
       KD=LBLOCK(MBD)                                                   FITM 108
       MAD=MBD-KD+1                                                     FITM 109
       IF(KD-1) 9999, 613, 615                                          FITM 110
 613   WTD =WW(MBD)                                                     FITM 111
       GO TO 617                                                        FITM 112
 615   WTD = DHAT(MBD)                                                  FITM 113
 617   DAVD = DHAT(MAD) / WTD                                           FITM 114
       IF( DAVD-DAV ) 630, 620, 620                                     FITM 115
  620  KNEW=K+KD                                                        FITM 116
       LBLOCK(MAD)=KNEW                                                 FITM 117
       LBLOCK(MB)=KNEW                                                  FITM 118
       DTONEW = DHAT(MAD) + DHAT(MA)                                    FITM 119
       DHAT(MAD) = DTONEW                                               FITM 120
       DHAT(MB) = WT + WTD                                              FITM 121
       NSATIS=0                                                         FITM 122
       MA=MAD                                                           FITM 123
       GO TO 800                                                        FITM 124
  630  NSATIS=NSATIS+1                                                  FITM 125
       GO TO 800                                                        FITM 126
C                                                                       FITM 127
C              IS BLOCK UP-SATISFIED. IF NOT, MERGE                     FITM 128
C                                                                       FITM 129
  700  IF(MB-MM) 710, 730, 9999                                         FITM 130
  710  MAU=MB+1                                                         FITM 131
       KU=LBLOCK(MAU)                                                   FITM 132
       MBU=MAU+KU-1                                                     FITM 133
       IF(KU-1) 9999, 713, 715                                          FITM 134
 713   WTU =WW(MBU)                                                     FITM 135
       GO TO 717                                                        FITM 136
 715   WTU = DHAT(MBU)                                                  FITM 137
 717   DAVU = DHAT(MAU) / WTU                                           FITM 138
       IF( DAV-DAVU ) 730, 720, 720                                     FITM 139
  720  KNEW=K+KU                                                        FITM 140
       LBLOCK(MA)=KNEW                                                  FITM 141
       LBLOCK(MBU)=KNEW                                                 FITM 142
       DTONEW = DHAT(MA) + DHAT(MAU)                                    FITM 143
       DHAT(MA) = DTONEW                                                FITM 144
       DHAT(MBU) = WT + WTU                                             FITM 145
       NSATIS=0                                                         FITM 146
       GO TO 800                                                        FITM 147
  730  NSATIS=NSATIS+1                                                  FITM 148
       GO TO 800                                                        FITM 149
C                                                                       FITM 150
C              PROCEED TO NEXT BLOCK IF READY.                          FITM 151
C                                                                       FITM 152
  800  LUD = 3-LUD                                                      FITM 153
C              QUERY. IS BLOCK BOTH UP AND DOWN SATISFIED. IF NOT, RETURFITM 154
       IF(NSATIS-1) 520, 520, 810                                       FITM 155
C              QUERY. IS THIS LAST BLOCK. IF NOT, GO ON TO NEXT BLOCK.  FITM 156
  810  IF(MB-MM) 820, 900, 9999                                         FITM 157
  820  MA=MB+1                                                          FITM 158
       GO TO 510                                                        FITM 159
C                                                                       FITM 160
C      MAIN COMPUTATION COMPLETE. PLACE ANSWERS IN DHAT.                FITM 161
C                                                                       FITM 162
 900  MA=ISTRT                                                          FITM 163
  910  K=LBLOCK(MA)                                                     FITM 164
       MB=MA+K-1                                                        FITM 165
       IF(K-1) 9999, 940, 920                                           FITM 166
 920   TEMP1 = DHAT(MA) / DHAT(MB)                                      FITM 167
       DO 930 M=MA,MB                                                   FITM 168
  930  DHAT(M)=TEMP1                                                    FITM 169
       GO TO 945                                                        FITM 170
 940   DHAT(MA) = DIST(MA)                                              FITM 171
 945   MA = MB + 1                                                      FITM 172
       IF(MA-MM-1) 910, 950, 9999                                       FITM 173
  950  RETURN                                                           FITM 174
C                                                                       FITM 175
C      TROUBLE EXIT                                                     FITM 176
C                                                                       FITM 177
 9999  WRITE            (LPRINT, 99)                                    FITM 178
   99  FORMAT(50H0KRUSKAL. IMPOSSIBLE BRANCH TAKEN ON IF STATEMENT. )   FITM 179
       STOP                                                             FITM 180
C                                                                       FITM 181
      END                                                               FITM 182
$      FORTRAN
CFITMS                                                                  FITMS  1
      SUBROUTINE FITMS(ISTRT,MM)                                        FITMS  2
C FITMS            FOR KYST                  JANUARY,1973               FITMS  3
C WRITTEN FOR KYST BY J.SEERY                JANUARY,1973               FITMS  4
C                                                                       FITMS  5
C       FITMS PERFORMS LEAST-SQUARES MONOTONE REGRESSION                FITMS  6
C      THIS ROUTINE FINDS THE VALUES FOR DHATIN WHICH ARE MONOTONIC     FITMS  7
C      AND WHICH BEST APPROXIMATE THE VALUES OF DISTIN IN THE SENSE     FITMS  8
C       THAT THE SUM OF THE SQUARED DEVIATIONS IS A MINIMUM.            FITMS  9
C      CALLED ONLY FROM INICON                                          FITMS 10
C      ASSUMES NO TIES.  DISTIN AND DHATIN SHARE THE SAME LOCATIONS.    FITMS 11
C                                                                       FITMS 12
      DIMENSION GR(100,6),GL(100,6),LBLOCK(1800),DUMMY(30),DHATIN(1800),FITMS 13
     .     IJIN(1800),DISTIN(1800)                                      FITMS 14
       COMMON  /MDSCL1/  LREAD, LPRINT, LPUNCH, LSCRAT                  FITMS 15
      COMMON /KYST2/ GR,GL,LBLOCK,DUMMY,DISTIN,IJIN                     FITMS 16
      EQUIVALENCE(DISTIN(1),DHATIN(1))                                  FITMS 17
C                                                                       FITMS 18
C      WITHIN EACH BLOCK, THE FIRST LBLOCK VALUE AND THE LAST LBLOCK    FITMS 19
C                 VALUE BOTH GIVE THE SIZE OF THE BLOCK.                FITMS 20
C                                                                       FITMS 21
C      FORM FIRST APPROXIMATION TO CORRECT PARTITON                     FITMS 22
C      INITIALIZE LBLOCK TO 1                                           FITMS 23
      DO 230 M=ISTRT,MM                                                 FITMS 24
  230  LBLOCK(M)=1                                                      FITMS 25
C                                                                       FITMS 26
C      START MAIN COMPUTATION      *************************************FITMS 27
C                                                                       FITMS 28
      MA=ISTRT                                                          FITMS 29
 510   LUD=2                                                            FITMS 30
       NSATIS=0                                                         FITMS 31
  520  K=LBLOCK(MA)                                                     FITMS 32
       MB=MA+K-1                                                        FITMS 33
      WT=K                                                              FITMS 34
 550  DAV=DHATIN(MA)/WT                                                 FITMS 35
       GO TO ( 600, 700 ), LUD                                          FITMS 36
C                                                                       FITMS 37
C              IS BLOCK DOWN-SATISFIED. IF NOT, MERGE                   FITMS 38
C                                                                       FITMS 39
 600  IF(MA-ISTRT)9999,630,610                                          FITMS 40
  610  MBD=MA-1                                                         FITMS 41
       KD=LBLOCK(MBD)                                                   FITMS 42
       MAD=MBD-KD+1                                                     FITMS 43
      WTD=KD                                                            FITMS 44
 617  DAVD=DHATIN(MAD)/WTD                                              FITMS 45
       IF( DAVD-DAV ) 630, 620, 620                                     FITMS 46
  620  KNEW=K+KD                                                        FITMS 47
       LBLOCK(MAD)=KNEW                                                 FITMS 48
       LBLOCK(MB)=KNEW                                                  FITMS 49
      DTONEW=DHATIN(MAD)+DHATIN(MA)                                     FITMS 50
      DHATIN(MAD)=DTONEW                                                FITMS 51
       NSATIS=0                                                         FITMS 52
       MA=MAD                                                           FITMS 53
       GO TO 800                                                        FITMS 54
  630  NSATIS=NSATIS+1                                                  FITMS 55
       GO TO 800                                                        FITMS 56
C                                                                       FITMS 57
C              IS BLOCK UP-SATISFIED. IF NOT, MERGE                     FITMS 58
C                                                                       FITMS 59
  700  IF(MB-MM) 710, 730, 9999                                         FITMS 60
  710  MAU=MB+1                                                         FITMS 61
       KU=LBLOCK(MAU)                                                   FITMS 62
       MBU=MAU+KU-1                                                     FITMS 63
      WTU=KU                                                            FITMS 64
 717  DAVU=DHATIN(MAU)/WTU                                              FITMS 65
       IF( DAV-DAVU ) 730, 720, 720                                     FITMS 66
  720  KNEW=K+KU                                                        FITMS 67
       LBLOCK(MA)=KNEW                                                  FITMS 68
       LBLOCK(MBU)=KNEW                                                 FITMS 69
      DTONEW=DHATIN(MA)+DHATIN(MAU)                                     FITMS 70
      DHATIN(MA)=DTONEW                                                 FITMS 71
       NSATIS=0                                                         FITMS 72
       GO TO 800                                                        FITMS 73
  730  NSATIS=NSATIS+1                                                  FITMS 74
       GO TO 800                                                        FITMS 75
C                                                                       FITMS 76
C              PROCEED TO NEXT BLOCK IF READY.                          FITMS 77
C                                                                       FITMS 78
  800  LUD = 3-LUD                                                      FITMS 79
C              QUERY. IS BLOCK BOTH UP AND DOWN SATISFIED. IF NOT, RETURFITMS 80
       IF(NSATIS-1) 520, 520, 810                                       FITMS 81
C              QUERY. IS THIS LAST BLOCK. IF NOT, GO ON TO NEXT BLOCK.  FITMS 82
  810  IF(MB-MM) 820, 900, 9999                                         FITMS 83
  820  MA=MB+1                                                          FITMS 84
       GO TO 510                                                        FITMS 85
C                                                                       FITMS 86
C      MAIN COMPUTATION COMPLETE. PLACE ANSWERS IN DHATIN.              FITMS 87
C                                                                       FITMS 88
 900  MA=ISTRT                                                          FITMS 89
  910  K=LBLOCK(MA)                                                     FITMS 90
       MB=MA+K-1                                                        FITMS 91
       IF(K-1) 9999, 940, 920                                           FITMS 92
 920  XK=K                                                              FITMS 93
      TEMP1=DHATIN(MA)/XK                                               FITMS 94
       DO 930 M=MA,MB                                                   FITMS 95
 930  DHATIN(M)=TEMP1                                                   FITMS 96
 940   MA = MB + 1                                                      FITMS 97
       IF(MA-MM-1) 910, 950, 9999                                       FITMS 98
  950  RETURN                                                           FITMS 99
C                                                                       FITMS100
C      TROUBLE EXIT                                                     FITMS101
C                                                                       FITMS102
 9999  WRITE            (LPRINT, 99)                                    FITMS103
 99   FORMAT(34H0ERROR EXIT FROM SUBROUTINE FITMS.)                     FITMS104
      CALL EXIT                                                         FITMS105
      END                                                               FITMS106
$      FORTRAN
CFITP                                                                   FITP   1
      SUBROUTINE FITP(ISTRT,MM,NTI,K)                                   FITP   2
C FITP             FOR KYST                  JANUARY,1973               FITP   3
C WRITTEN BY J.KRUSKAL                                                  FITP   4
C      FITP--LEAST SQUARES POLYNOMIAL REGRESSION                        FITP   5
C                                                                       FITP   6
      DIMENSION DATA(1800),WW(1800),IJ(1800),X(100,6)                   FITP   7
      DIMENSION GR(100,6),GL(100,6),LBLOCK(1800),C(30),DIST(1800),      FITP   8
     .     DHAT(1800),A(5,5),B(5),V(5)                                  FITP   9
      EQUIVALENCE (A(1,1),LBLOCK(1)) ,(B(1),LBLOCK(26)),                FITP  10
     .     (V(1),LBLOCK(31))                                            FITP  11
      COMMON /KYST1/ DATA,WW,IJ,X                                       FITP  12
      COMMON /KYST2/ GR,GL,LBLOCK,C,DIST,DHAT                           FITP  13
C                                                                       FITP  14
       NT = MAX0( 1, MIN0( 5,NTI ) )                                    FITP  15
C                                                                       FITP  16
       SUMW=0.0                                                         FITP  17
       DO 1105 I=K,NT                                                   FITP  18
       B(I) = 0.0                                                       FITP  19
       DO 1105 J=K,I                                                    FITP  20
       A(I,J) = 0.0                                                     FITP  21
 1105  CONTINUE                                                         FITP  22
C                                                                       FITP  23
      DO 1150 M=ISTRT,MM                                                FITP  24
       DA=DATA(M)                                                       FITP  25
       SUMW=SUMW+WW(M)                                                  FITP  26
C                                                                       FITP  27
       DO 1110 I=K,NT                                                   FITP  28
       V(I)=REGR(DA,I)                                                  FITP  29
 1110  CONTINUE                                                         FITP  30
C                                                                       FITP  31
       DO 1140 I=K,NT                                                   FITP  32
       DO 1130 J=K,I                                                    FITP  33
       A(I,J)=A(I,J)+V(I)*V(J)*WW(M)                                    FITP  34
 1130  CONTINUE                                                         FITP  35
       B(I)=B(I)+V(I)*DIST(M)*WW(M)                                     FITP  36
 1140  CONTINUE                                                         FITP  37
C                                                                       FITP  38
 1150  CONTINUE                                                         FITP  39
C                                                                       FITP  40
C                                                                       FITP  41
       DO 1155 I = K,NT                                                 FITP  42
       DO 1155 J= K,I                                                   FITP  43
       A(J,I) = A(I,J)                                                  FITP  44
 1155  CONTINUE                                                         FITP  45
C                                                                       FITP  46
       CALL EQSOLV(A,B,C,NT,K)                                          FITP  47
C                                                                       FITP  48
      DO 1170 M=ISTRT,MM                                                FITP  49
       DA = DATA(M)                                                     FITP  50
       TEMP=0.0                                                         FITP  51
       DO 1160 I=K,NT                                                   FITP  52
       TEMP=TEMP+C(I)*REGR(DA,I)                                        FITP  53
 1160  CONTINUE                                                         FITP  54
       DHAT(M)=TEMP                                                     FITP  55
 1170  CONTINUE                                                         FITP  56
C                                                                       FITP  57
       RETURN                                                           FITP  58
       END                                                              FITP  59
$      FORTRAN STAB
CINICON                                                                 INICON 1
      SUBROUTINE INICON(N,ND,M1,NITER,LHIPRT,SDS,LSCH)                  INICON 2
C INICON           FOR KYST                  JANUARY,1973               INICON 3
C WRITTEN BY F.YOUNG                                                    INICON 4
C MODIFIED FOR KYST BY J.KRUSKAL AND J.SEERY JANUARY,1973               INICON 5
C                                                                       INICON 6
C WRITTEN BY FORREST W. YOUNG JUNE 1971                                 INICON 7
C ROUTINE TO COMPUTE INITIAL CONFIGURATION                              INICON 8
C        USES THE TORGERSON-YOUNG PROCEDURE TO DERIVE THE BEST FITTING  INICON 9
C        QUASI-NONMETRIC EUCLIDIAN CONFIGURATION                        INICON10
C                                                                       INICON11
      DIMENSION DATA(1800),WW(1800),IJ(1800),X(100,6)                   INICON12
      DIMENSION GR(100,6),XBEST(100,6),DATAIN(5430),IJIN(1800),         INICON13
     .     DHATIN(1800),DUMMY(1)                                        INICON14
      COMMON /ACCUR/ PRECSN,XMAG,XMAXN                                  INICON15
      COMMON /MDSCL1/ LREAD,LPRINT,LPUNCH,LSCRAT                        INICON16
      COMMON /KYST1/ DATA,WW,IJ,X                                       INICON17
      COMMON /KYST2/ GR,XBEST,DATAIN                                    INICON18
      EQUIVALENCE (DHATIN(1),DATAIN(1831)) , (IJIN(1),DATAIN(3631))     INICON19
      DATA NCOLS,KPACK,IJPACK /6,10000,100/                             INICON20
      ISUB(I1,I2)=((I2-1)*(I2-2))/2+I1                                  INICON21
C                                                                       INICON22
C PARAMETERS                                                            INICON23
C                                                                       INICON24
C CONSTANTS                                                             INICON25
C                                                                       INICON26
      NM1=N-1                                                           INICON27
      NM2=N-2                                                           INICON28
      N2=(N*NM1)/2                                                      INICON29
      FN2=N2                                                            INICON30
      MTIMES=NM1/NCOLS+1                                                INICON31
      XMAGM=-XMAXN                                                      INICON32
      STRBST = 1.0                                                      INICON33
C                                                                       INICON34
C      CONVERT RAW DATA VALUES IN DATA TO AVERAGED VALUES, AND STORE    INICON35
C       THEM IN CORRESPONDING POSITIONS IN DATAIN.  WORK WITH NCOLS (=6)INICON36
C        ROWS AND COLUMNS AT A TIME.  FLAG MISSING DATA VALUES.         INICON37
C        AT THE SAME TIME COMPUTE MAX,MIN, AND AVERAGE OF NON-MISSING   INICON38
C        DATAIN VALUES                                                  INICON39
C                                                                       INICON40
      AVG=0.0                                                           INICON41
      NUM=0                                                             INICON42
      DINMAX=XMAGM                                                      INICON43
      DINMIN=1.0/XMAG                                                   INICON44
      INSUB=0                                                           INICON45
      MEND=0                                                            INICON46
      DO 30 MT=1,MTIMES                                                 INICON47
      MSTART=MEND+1                                                     INICON48
      MEND=MEND+NCOLS                                                   INICON49
      MEND=MIN0(MEND,N)                                                 INICON50
      MNDM1=MEND-1                                                      INICON51
      DO 5 L1=1,MNDM1                                                   INICON52
      DO 5 L2=1,NCOLS                                                   INICON53
      X(L1,L2)=0.                                                       INICON54
 5    XBEST(L1,L2)=0.                                                   INICON55
      DO 15 L=1,M1                                                      INICON56
      IF(WW(L)) 15,15,8                                                 INICON57
 8    KK=MOD(IJ(L),KPACK)                                               INICON58
      I=KK/IJPACK+1                                                     INICON59
      J=MOD(KK,IJPACK)+1                                                INICON60
      K=MIN0(I,J)                                                       INICON61
      KK=MAX0(I,J)                                                      INICON62
      IF((KK.LT.MSTART).OR.(KK.GT.MEND)) GO TO 15                       INICON63
      IF(K.EQ.KK) GO TO 15                                              INICON64
      KK=MOD(KK-1,NCOLS)+1                                              INICON65
      X(K,KK)=X(K,KK)+DATA(L)*WW(L)                                     INICON66
      XBEST(K,KK)=XBEST(K,KK)+WW(L)                                     INICON67
 15   CONTINUE                                                          INICON68
      KK=0                                                              INICON69
      IF(MSTART.NE.1) GO TO 18                                          INICON70
      MSTART=2                                                          INICON71
      KK=1                                                              INICON72
 18   DO 25 M=MSTART,MEND                                               INICON73
      MM1=M-1                                                           INICON74
      KK=KK+1                                                           INICON75
      DO 25 K=1,MM1                                                     INICON76
      INSUB=INSUB+1                                                     INICON77
      IF(XBEST(K,KK)) 20,20,22                                          INICON78
 20   DATAIN(INSUB)=XMAGM                                               INICON79
      GO TO 25                                                          INICON80
 22   NUM=NUM+1                                                         INICON81
      DIN=X(K,KK)/XBEST(K,KK)                                           INICON82
      AVG=AVG+DIN                                                       INICON83
      DATAIN(INSUB)=DIN                                                 INICON84
      DINMAX=AMAX1(DINMAX,DIN)                                          INICON85
      DINMIN=AMIN1(DINMIN,DIN)                                          INICON86
 25   CONTINUE                                                          INICON87
 30   CONTINUE                                                          INICON88
      FNUM=NUM                                                          INICON89
      AVG=AVG/FNUM                                                      INICON90
C                                                                       INICON91
C        REPLACE MISSING DATAIN VALUES WITH AVERAGE VALUE, AND          INICON92
C        CONVERT DATAIN TO DISTANCE-LIKE VALUES                         INICON93
C                                                                       INICON94
      IF(SDS) 35,99,38                                                  INICON95
C        HERE IF SIMILARITIES                                           INICON96
 35   XFAC=(10.0*DINMAX-DINMIN)/9.0                                     INICON97
      GO TO 40                                                          INICON98
C        HERE IF DISSIMILARITIES                                        INICON99
 38   XFAC=(DINMAX-10.0*DINMIN)/9.0                                     INICON00
 40   DO 60 L=1,N2                                                      INICON01
      IF(DATAIN(L).GT.XMAGM) GO TO 55                                   INICON02
      DATAIN(L)=XFAC+SDS*AVG                                            INICON03
      GO TO 60                                                          INICON04
 55   DATAIN(L)=XFAC+SDS*DATAIN(L)                                      INICON05
 60   CONTINUE                                                          INICON06
C                                                                       INICON07
C NORMALIZE THE CURRENT DATA                                            INICON08
C                                                                       INICON09
      SSQ=0.0                                                           INICON10
      DO 70  I=1,N2                                                     INICON11
 70   SSQ=SSQ+DATAIN(I)*DATAIN(I)                                       INICON12
      SSQ=SQRT(SSQ/FN2)                                                 INICON13
      DO 75 L=1,N2                                                      INICON14
 75   DATAIN(L)=DATAIN(L)/SSQ                                           INICON15
C                                                                       INICON16
      WRITE(LPRINT,1000) N,ND                                           INICON17
 1000 FORMAT(46H0INITIAL CONFIGURATION COMPUTATION   NO. PTS.=,I4,      INICON18
     .    8H    DIM=,I4)                                                INICON19
C                                                                       INICON20
C COMPUTE THE BEST FITTING METRIC EUCLIDIAN CONFIGURATION FOR THIS DATA INICON21
C                                                                       INICON22
      CALL CONFIG(N,ND)                                                 INICON23
C                                                                       INICON24
C     RECONVERT RAW DATA VALUES IN DATA TO AVERAGED VALUES   STORE THEM INICON25
C     IN DHATIN, AND THEIR POSITIONS IN IJIN.  IF PRE-ITERATIONS ARE    INICON26
C     TO BE PERFORMED, FLAG THE MISSING DATA AND POSITIONS.  OTHERWISE, INICON27
C     PACK THE AVERAGED VALUES TOGETHER.                                INICON28
C                                                                       INICON29
      ISTRT=1                                                           INICON30
      INSUB=0                                                           INICON31
      MEND=0                                                            INICON32
      DO 160 MT=1,MTIMES                                                INICON33
      MSTART=MEND+1                                                     INICON34
      MEND=MEND+NCOLS                                                   INICON35
      MEND=MIN0(MEND,N)                                                 INICON36
      MNDM1=MEND-1                                                      INICON37
      DO 130 L1=1,MNDM1                                                 INICON38
      DO 130 L2=1,NCOLS                                                 INICON39
      GR(L1,L2)=0.                                                      INICON40
 130  XBEST(L1,L2)=0.                                                   INICON41
      DO 140 L=1,M1                                                     INICON42
      IF(WW(L)) 140,140,134                                             INICON43
 134  KK=MOD(IJ(L),KPACK)                                               INICON44
      I=KK/IJPACK+1                                                     INICON45
      J=MOD(KK,IJPACK)+1                                                INICON46
      K=MIN0(I,J)                                                       INICON47
      KK=MAX0(I,J)                                                      INICON48
      IF((KK.LT.MSTART).OR.(KK.GT.MEND)) GO TO 140                      INICON49
      IF(K.EQ.KK) GO TO 140                                             INICON50
      KK=MOD(KK-1,NCOLS)+1                                              INICON51
      GR(K,KK)=GR(K,KK)+DATA(L)*WW(L)                                   INICON52
      XBEST(K,KK)=XBEST(K,KK)+WW(L)                                     INICON53
 140  CONTINUE                                                          INICON54
      KK=0                                                              INICON55
      IF(MSTART.NE.1) GO TO 144                                         INICON56
      MSTART=2                                                          INICON57
      KK=1                                                              INICON58
 144  DO 155 M=MSTART,MEND                                              INICON59
      MM1=M-1                                                           INICON60
      KK=KK+1                                                           INICON61
      DO 155 K=1,MM1                                                    INICON62
      IF(XBEST(K,KK)) 146,146,152                                       INICON63
 146  IF(NITER.EQ.0) GO TO 155                                          INICON64
      ISTRT=ISTRT+1                                                     INICON65
      INSUB=INSUB+1                                                     INICON66
      DHATIN(INSUB)=XMAGM                                               INICON67
      IJIN(INSUB)=IJPACK*(K-1)+MM1                                      INICON68
      GO TO 155                                                         INICON69
 152  INSUB=INSUB+1                                                     INICON70
C     NEGATE DATA VALUES IF SIMILARITIES                                INICON71
      DHATIN(INSUB)=SDS*GR(K,KK)/XBEST(K,KK)                            INICON72
      IJIN(INSUB)=IJPACK*(K-1)+MM1                                      INICON73
 155  CONTINUE                                                          INICON74
 160  CONTINUE                                                          INICON75
      NELE=INSUB                                                        INICON76
      NUM=INSUB-ISTRT+1                                                 INICON77
      FNUM=NUM                                                          INICON78
C                                                                       INICON79
C    SORT DHATIN INTO ASCENDING ORDER, AND REARRANGE IJIN.              INICON80
 200  CALL SORT(DHATIN,NELE,IJIN,DUMMY,DUMMY,1,1)                       INICON81
C                                                                       INICON82
      IF(NITER) 99,228,210                                              INICON83
C                                                                       INICON84
C    COMPUTE ITERATIVE BEST-FITTING EUCLIDIAN INITIAL CONFIGURATION     INICON85
C                                                                       INICON86
 210  IF(LHIPRT.EQ.2) WRITE(LPRINT,1010)                                INICON87
 1010 FORMAT(30H0     PRE-ITERATION     STRESS,/)                       INICON88
      IT=-1                                                             INICON89
C                                                                       INICON90
C    START OF ITERATIVE COMPUTATION                                     INICON91
C                                                                       INICON92
 220  IT=IT+1                                                           INICON93
C                                                                       INICON94
C COMPUTE MINKOWSKI DISTANCES                                           INICON95
C                                                                       INICON96
 228  DBAR=0.0                                                          INICON97
      DO 235 L=1,NELE                                                   INICON98
      K=IJIN(L)                                                         INICON99
      I=K/IJPACK+1                                                      INICON00
      J=MOD(K,IJPACK)+1                                                 INICON01
      SUM=0.0                                                           INICON02
      DO 230 K=1,ND                                                     INICON03
 230  SUM=SUM+RPOWER(X(I,K)-X(J,K))                                     INICON04
      DHATIN(L)=RROOT(SUM)                                              INICON05
      IF(L.LT.ISTRT) GO TO 235                                          INICON06
      DBAR=DBAR+DHATIN(L)                                               INICON07
 235  CONTINUE                                                          INICON08
C                                                                       INICON09
C COMPUTE THE BEST-FITTING PSEUDO-DISTANCES USING MONOTONE REGRESSION   INICON10
C                                                                       INICON11
      CALL FITMS(ISTRT,NELE)                                            INICON12
C                                                                       INICON13
C COMPUTE STRESS OF THE CURRENT CONFIGURATION                           INICON14
C                                                                       INICON15
      U=0.0                                                             INICON16
      T=0.0                                                             INICON17
C                                                                       INICON18
C IF STRESS FORMULA 1 OPTED BY USER THEN DBAR=0                         INICON19
C                                                                       INICON20
      DBAR=DBAR/FNUM                                                    INICON21
      IF(LSCH.EQ.1) DBAR=0.0                                            INICON22
      DO 240 L=ISTRT,NELE                                               INICON23
      K=IJIN(L)                                                         INICON24
      I=K/IJPACK+1                                                      INICON25
      J=MOD(K,IJPACK)+1                                                 INICON26
      SUM=0.0                                                           INICON27
      DO 238 K=1,ND                                                     INICON28
 238  SUM=SUM+RPOWER(X(I,K)-X(J,K))                                     INICON29
      DISTL=RROOT(SUM)                                                  INICON30
      U=U+(DISTL-DHATIN(L))**2                                          INICON31
 240  T=T+(DISTL-DBAR)**2                                               INICON32
      STRESS=SQRT(U/T)                                                  INICON33
C                                                                       INICON34
      IF(NITER) 99,875,245                                              INICON35
C                                                                       INICON36
C WRITE RESULTS OF CURRENT ITERATION                                    INICON37
C                                                                       INICON38
 245  IF(LHIPRT.EQ.2) WRITE(LPRINT,1015) IT,STRESS                      INICON39
 1015 FORMAT(1H ,I15,F13.4)                                             INICON40
C DETERMINE IF THIS IS THE BEST ITERATION SO FAR.                       INICON41
C                                                                       INICON42
      IF(STRESS) 99,247,248                                             INICON43
 247  STRBST=STRESS                                                     INICON44
      GO TO 850                                                         INICON45
 248  IF(STRESS-STRBST) 250,250,330                                     INICON46
C                                                                       INICON47
C    THIS IS THE BEST PRE-ITERATION SO FAR.  IF THERE ARE MORE PRE-     INICON48
C     ITERATIONS TO GO, GO ON TO SEE IF A BETTER CONFIGURATION CAN BE   INICON49
C     FOUND.                                                            INICON50
C                                                                       INICON51
 250  STRBST=STRESS                                                     INICON52
      IF(IT-NITER) 255,290,290                                          INICON53
C                                                                       INICON54
C    PUT DHATIN ENTRIES INTO CORRECT DATAIN POSITIONS                   INICON55
C    FILL MISSING ENTRIES OF DATAIN WITH THE CORRESPONDING DISTANCES--  INICON56
C     WHICH ARE STORED IN THE FIRST POSITIONS OF DHATIN                 INICON57
C                                                                       INICON58
 255  DO 260 L=1,NELE                                                   INICON59
      K=IJIN(L)                                                         INICON60
      I=K/IJPACK+1                                                      INICON61
      J=MOD(K,IJPACK)+1                                                 INICON62
      M=ISUB(I,J)                                                       INICON63
 260  DATAIN(M)=DHATIN(L)                                               INICON64
C                                                                       INICON65
C      SAVE THIS CONFIGURATION AS BEST SO FAR                           INICON66
C                                                                       INICON67
      DO 275 I=1,N                                                      INICON68
      DO 275 J=1,ND                                                     INICON69
 275  XBEST(I,J)=X(I,J)                                                 INICON70
C                                                                       INICON71
      CALL CONFIG(N,ND)                                                 INICON72
C                                                                       INICON73
C    GO TO TO NEXT PRE-ITERATION                                        INICON74
      GO TO 220                                                         INICON75
C                                                                       INICON76
 290  WRITE(LPRINT,1020) NITER                                          INICON77
 1020 FORMAT(33H0MAXIMUM NUMBER OF PRE-ITERATIONS,I3,10H, REACHED.)     INICON78
      GO TO 850                                                         INICON79
C                                                                       INICON80
C IF THIS SOLUTION FITS NOT AS WELL AS THE PREVIOUS ONE, COME HERE      INICON81
C TO RECOVER THE BEST SOLUTION.                                         INICON82
C                                                                       INICON83
 330  IT=IT-1                                                           INICON84
      WRITE(LPRINT,1025) IT                                             INICON85
 1025 FORMAT(73H0STRESS STARTING TO INCREASE  BEST VALUE ACHIEVED ON PREINICON86
     .-ITERATION NUMBER, I3)                                            INICON87
      DO 31 I=1,N                                                       INICON88
      DO 31 J=1,ND                                                      INICON89
   31 X(I,J)=XBEST(I,J)                                                 INICON90
C                                                                       INICON91
 850  WRITE(LPRINT,1030) N,ND,STRBST,LSCH                               INICON92
 1030 FORMAT(34H0THE BEST INITIAL CONFIGURATION OF,I4,10H POINTS IN,I4, INICON93
     .   22H DIMENSIONS HAS STRESS,F7.3,8H FORMULA,I2)                  INICON94
      GO TO 900                                                         INICON95
C                                                                       INICON96
C      NO PRE-ITERATIONS PERFORMED.                                     INICON97
1005  FORMAT(38H0NO PRE-ITERATIONS PERFORMED.  STRESS=,F8.4,            INICON98
     .        9H ,FORMULA,I2)                                           INICON99
 875  WRITE(LPRINT,1005) STRESS ,LSCH                                   INICON00
 900  RETURN                                                            INICON01
C                                                                       INICON02
   99 WRITE(LPRINT,98)                                                  INICON03
 98   FORMAT(22H0ERROR EXIT FOR INICON)                                 INICON04
      CALL EXIT                                                         INICON05
      END                                                               INICON06
$      FORTRAN
CNRMLZ                                                                  NRMLZ  1
      SUBROUTINE NRMLZ(N,V)                                             NRMLZ  2
C NRMLZ            FOR KYST                  JANUARY,1973               NRMLZ  3
C WRITTEN FOR SGEV BY P.BUSINGER             JANUARY,1972               NRMLZ  4
C CALLED BY NVIT1 TO NORMALIZE A VECTOR OF LENGTH N                     NRMLZ  5
      DIMENSION V(1)                                                    NRMLZ  6
C                                                                       NRMLZ  7
      REAL S                                                            NRMLZ  8
      S=0.E0                                                            NRMLZ  9
      DO 10 K=1,N                                                       NRMLZ 10
   10    S=AMAX1(S,ABS(V(K)))                                           NRMLZ 11
      S=SQRT(S)                                                         NRMLZ 12
      DO 20 K=1,N                                                       NRMLZ 13
   20    V(K)=V(K)/S                                                    NRMLZ 14
      RETURN                                                            NRMLZ 15
      END                                                               NRMLZ 16
$      FORTRAN
CNEWSTP                                                                 NEWSTP 1
       SUBROUTINE NEWSTP( STEP, ITNO, SFGR, STRESS,                     NEWSTP 2
     1 CAGRGL, COSAV, ACSAV, COSAVW, ACSAVW, SRAT, SRATAV )             NEWSTP 3
C NEWSTP           FOR KYST                  JANUARY,1973               NEWSTP 4
C WRITTEN BY J.KRUSKAL                                                  NEWSTP 5
C                                                                       NEWSTP 6
C      NEWSTP  THIS SUBROUTINE COMPUTES THE STEP SIZE.                  NEWSTP 7
C                                                                       NEWSTP 8
C              THE MAIN PURPOSE OF THIS ROUTINE IS TO COMPUTE THE NEW   NEWSTP 9
C              VALUE OF  STEP .                                         NEWSTP10
C              INCIDENTALLY, IT UPDATES  COSAV ,  ACSAV , AND  SRATAV . NEWSTP11
C                                                                       NEWSTP12
C      UPDATE THREE AVERAGE QUANTITIES                                  NEWSTP13
C                                                                       NEWSTP14
       COSAV = CAGRGL*COSAVW  +  COSAV*(1.0-COSAVW)                     NEWSTP15
       ACSAV = ABS (CAGRGL)*ACSAVW  +  ACSAV*(1.0-ACSAVW)               NEWSTP16
       SRATAV = (SRAT**0.33334)  *  (SRATAV**0.66666)                   NEWSTP17
       IF(ITNO) 100, 100, 200                                           NEWSTP18
C                                                                       NEWSTP19
C              GUESS INITIAL STEP SIZE                                  NEWSTP20
C                                                                       NEWSTP21
 100   STEP= (25.0*STRESS) * SFGR                                       NEWSTP22
       RETURN                                                           NEWSTP23
C                                                                       NEWSTP24
C              FIND NEW STEP SIZE                                       NEWSTP25
C                                                                       NEWSTP26
 200   ANG=4.0**COSAV                                                   NEWSTP27
       TEMP1 = 1.0 + (AMIN1 (1.0,SRATAV) ) ** 5                         NEWSTP28
       TEMP2 = 1.0 + ( ACSAV - ABS (COSAV) )                            NEWSTP29
       BIAS = 1.6 / (TEMP1*TEMP2)                                       NEWSTP30
       GOODLK = SQRT (AMIN1 (1.0,SRAT) )                                NEWSTP31
       STEP = STEP * ANG * BIAS * GOODLK                                NEWSTP32
       RETURN                                                           NEWSTP33
       END                                                              NEWSTP34
$      FORTRAN
CNVIT1                                                                  NVIT1  1
      SUBROUTINE NVIT1(N,OMEGA,B,U,V)                                   NVIT1  2
C NVIT1            FOR KYST                  JANUARY,1973               NVIT1  3
C WRITTEN AND TECHNIQUES ADVOCATED IN                                   NVIT1  4
C REFERENCE IMPLEMENTED FOR SGEV                                        NVIT1  5
C BY P.BUSINGER                              JANUARY,1972               NVIT1  6
      DIMENSION B(1),U(1),V(1)                                          NVIT1  7
C                                                                       NVIT1  8
C G.GOLUB AND W.KAHAN, CALCULATING THE SINGULAR                         NVIT1  9
C VALUES AND PSEUDO-INVERSE OF A MATRIX. J.SIAM                         NVIT1 10
C NUMER.ANAL., SER.B, VOL.2 (1965), 205-224.                            NVIT1 11
C                                                                       NVIT1 12
C NVIT1 IS CALLED BY SGEV, AND ASSUME THAT                              NVIT1 13
C GIVEN ARE THE DIAGONAL ELEMENTS                                       NVIT1 14
C         U(1) THROUGH U(N)                                             NVIT1 15
C (ALL OF THE SAME SIGN) OF THE MATRIX                                  NVIT1 16
C U IN THE LU-DECOMPOSITION OF A NEARLY                                 NVIT1 17
C SINGULAR SYMMETRIC TRIDIAGONAL MATRIX                                 NVIT1 18
C A. ALSO GIVEN ARE THE OFF-DIAGONAL                                    NVIT1 19
C ELEMENTS                                                              NVIT1 20
C         B(1) THROUGH B(N-1)                                           NVIT1 21
C OF A AND OMEGA = LARGEST MACHINE NUMBER.                              NVIT1 22
C COMPUTED ARE THE COMPONENTS                                           NVIT1 23
C         V(1) THROUGH V(N)                                             NVIT1 24
C OF THE EIGENVECTOR ASSOCIATED WITH THE                                NVIT1 25
C APPROXIMATE ZERO EIGENVALUE OF A.                                     NVIT1 26
C                                                                       NVIT1 27
C CALLS NRMLZ TO NORMALIZE A VECTOR                                     NVIT1 28
      REAL W,S,OMGA                                                     NVIT1 29
      OMGA=1.E-1*OMEGA                                                  NVIT1 30
      V(1)=1.E0                                                         NVIT1 31
      IF(N.LT.2)GOTO 40                                                 NVIT1 32
      DO 10 I=2,N                                                       NVIT1 33
         I1=I-1                                                         NVIT1 34
         IF(ABS(B(I1)*V(I1)).GE.ABS(OMGA*U(I1)))CALL NRMLZ(I1,V)        NVIT1 35
         W=-(B(I1)*V(I1))/U(I1)                                         NVIT1 36
   10    V(I)=W+SIGN(V(1),W)                                            NVIT1 37
      CALL NRMLZ(N,V)                                                   NVIT1 38
      B(N)=0.E0                                                         NVIT1 39
      S=0.E0                                                            NVIT1 40
      DO 30 II=1,N                                                      NVIT1 41
         I=N+1-II                                                       NVIT1 42
         IF(ABS(V(I)-B(I)*S).GE.ABS(OMGA*U(I)))CALL NRMLZ(I1,V)         NVIT1 43
         W=(V(I)-B(I)*S)/U(I)                                           NVIT1 44
   20    V(I)=W                                                         NVIT1 45
   30    S=W                                                            NVIT1 46
      S=0.E0                                                            NVIT1 47
      DO 35 I=1,N                                                       NVIT1 48
   35    S=S+V(I)**2                                                    NVIT1 49
      S=SQRT(S)                                                         NVIT1 50
      DO 39 I=1,N                                                       NVIT1 51
   39    V(I)=V(I)/S                                                    NVIT1 52
   40 RETURN                                                            NVIT1 53
      END                                                               NVIT1 54
$      FORTRAN
CPLOT                                                                   PLOT   1
      SUBROUTINE PLOT(X,Y,XA,YA,XI,YI,NPOI,NY,ID,LONG,MX)               PLOT   2
C PLOT             FOR KYST                  JANUARY,1973               PLOT   3
C WRITTEN BY FORREST W. YOUNG, NOVEMBER, 1965                           PLOT   4
C     UPDATED OCTOBER 11,1971  JAY R LEVINSOHN                          PLOT   5
C ALTERED FOR KYST BY J.KRUSKAL AND J.SEERY  JANUARY,1973               PLOT   6
C                                                                       PLOT   7
C ROUTINE TO GENERATE A ONE PAGE PLOT OF A MATRIX Y VS SOME VECTOR X    PLOT   8
C                                                                       PLOT   9
C Y MAY BE A MATRIX AND IT IS ASSUMED THAT X IS A VECTOR                PLOT  10
C                                                                       PLOT  11
C THE PARAMETERS XA,YA,XI,YI INDICATE THE UPPER                         PLOT  12
C AND LOWER BOUNDS FOR EACH AXIS.                                       PLOT  13
C THE USER MAY PASS THE MIN AND MAX FOR X AND Y OR HE MAY LET THE       PLOT  14
C PROGRAM FIND THEM.  IF XA=XI THEN THE PROGRAM WILL FIND THE MAX AND MIPLOT  15
C FOR THE X VECTOR,  IF YA=YI THE PROGRAM WILL FIND THE MAX AND MIN     PLOT  16
C FOR THE Y MATRIX.                                                     PLOT  17
C                                                                       PLOT  18
C IF  -ID-  IS POSITIVE, AXES WILL BE INCLUDED ON THE PLOT,             PLOT  19
C IF  -ID-  IS NEGATIVE, NO AXES WILL APPEAR.                           PLOT  20
C IF ID IS PLUS OR MINUS 2, THE POINTS WILL BE IDENTIFIED BY ROWS       PLOT  21
C IF ID IS NOT EQUAL TO 2 THEN POINTS WILL BE LABELED BY COLS.          PLOT  22
C                                                                       PLOT  23
C IF ID IS PLUS OR MINUS THREE THE PLOT IS ASSUMED TO BE 1 DIMENSIONAL. PLOT  24
C ONE DIMENSIONAL PLOTS ARE PLOTTED WITH NO AXES, AND ALL Y VALUES      PLOT  25
C FIXED AT MAXY-MINY/2                                                  PLOT  26
C POINTS IN THE 1 DIM. PLOTS ARE IDENDTIFIED BY ROW.                    PLOT  27
C                                                                       PLOT  28
C DEFINITIONS                                                           PLOT  29
C  PARAMETERS                                                           PLOT  30
C    Y--A MATRIX (LONG,NX) IT CONTAINS VALUES TO BE PLOTTED AGAINST X   PLOT  31
C    X--A VECTOR (LONG); IT CONTAINS VALUES TO BE PLOTTED AGAINST Y     PLOT  32
C    XA-- X AXIS UPPER BOUND   (MAY BE PASSED AS ZERO)                  PLOT  33
C    YA-- Y AXIS  "     "     "      "     "                            PLOT  34
C     XI-- X AXIS LOWER BOUND  (MAY BE PASSED AS ZERO)                  PLOT  35
C    YI-- Y AXIS LOWER BOUND(MAY BE PASSED AS ZERO )                    PLOT  36
C    NPOI-- NO. OF ELEMENTS IN X, NO. OF ROWS IN Y                      PLOT  37
C    NY--NUMBER OF COLUMNS IN Y MATRIX                                  PLOT  38
C    ID-- PRINT CONTROL PARAMETER, CONTROLS LABELING, AND AXIS GEN.     PLOT  39
C     LONG-- MAXIMUM  NUMBER OF ELEMENTS IN X, MAX. NO. ROWS IN Y       PLOT  40
C  MX  -- WORK WITH ELEMENTS MX THROUGH MX+NPOI-1                       PLOT  41
C  VARIABLES                                                            PLOT  42
C    DELX-- X AXIS INCREMENT                                            PLOT  43
C    DELY-- Y AXIS INCREMENT                                            PLOT  44
C    FINC-- FLOATING INCREMENT FOR XLABEL                               PLOT  45
C    IDAT-- COUNTER OVER ROWS OF X AND Y                                PLOT  46
C    ILAB-- INTEGER LABEL FLAG                                          PLOT  47
C    ITEM(101)--ITEMS TO BE PRINTED INITIALLY BLANK                     PLOT  48
C    LINE-- COUNTER OF CURRENT LINE                                     PLOT  49
C     NLPP-- NUMBER OF LINES PLOTTED PER PAGE                           PLOT  50
C    SPANX-- RANGE OF X                                                 PLOT  51
C    SPANY-- RANGE OF Y                                                 PLOT  52
C    VALUE-- DATUM PRINTED ON Y AXIS                                    PLOT  53
C    XLABEL(11)--X-AXIS LABELS                                          PLOT  54
C     XNOW-- CURRENT VALUE OF X VECTOR NOW BEING WORKED ON              PLOT  55
C     YNOW-- CURRENT VALUE OF Y MATRIX NOW BEING WORKED ON              PLOT  56
C                                                                       PLOT  57
      REAL ITEM(101),MINUS,PTID(100),YMINV(2),YMAXV(2)                  PLOT  58
      DIMENSION XLABL(11),X(1),Y(LONG,NY),XS(10),YS(10),IND(2)          PLOT  59
      COMMON /MDSCL1/ LREAD,LPRINT,LPUNCH,LSCRAT                        PLOT  60
      COMMON /ACCUR/ PRECSN,XMAG,XMAXN                                  PLOT  61
      COMMON /PLTCHR/ PTID,ITEM                                         PLOT  62
      DATA MINUS,DD,STAR,BLANK,AIE /1H-,1H.,1H*,1H ,1HI/                PLOT  63
      DATA XS/0.,0.,0.,0.,0.,4.,2.,4.,1.,0./                            PLOT  64
      DATA YS/0.,0.,0.,0.,2.,4.,3.,4.,7.,2./                            PLOT  65
C                                                                       PLOT  66
      JY=1                                                              PLOT  67
      KODE=1                                                            PLOT  68
      INDEX=0                                                           PLOT  69
      LINE=1                                                            PLOT  70
      IY0L=0                                                            PLOT  71
      IX0L=0                                                            PLOT  72
      IDAT=1                                                            PLOT  73
      NLPP=53                                                           PLOT  74
      XMAX=XA                                                           PLOT  75
      YMAX=YA                                                           PLOT  76
      YMIN=YI                                                           PLOT  77
      XMIN=XI                                                           PLOT  78
      NPOI2=NY*NPOI                                                     PLOT  79
      NPOIRL=MX+NPOI-1                                                  PLOT  80
      LXFLG=0                                                           PLOT  81
      LYFLG=0                                                           PLOT  82
C                                                                       PLOT  83
C DETERMINE MINIMUM AND MAXIMUM OF X AXIS, IF NECESSARY                 PLOT  84
C                                                                       PLOT  85
      IF (XMAX-XMIN)122,119,122                                         PLOT  86
 119  XMAX=-XMAXN                                                       PLOT  87
      XMIN=XMAXN                                                        PLOT  88
      DO 121 I=MX,NPOIRL                                                PLOT  89
      IF(X(I)-XMAX) 125,123,123                                         PLOT  90
 123  XMAX=X(I)                                                         PLOT  91
 125  IF(X(I)-XMIN) 124,124,121                                         PLOT  92
 124  XMIN=X(I)                                                         PLOT  93
  121 CONTINUE                                                          PLOT  94
C                                                                       PLOT  95
C     CALL SERCH TO INITIALIZE; IF DESIRED SET YMAX AND YMIN            PLOT  96
C                                                                       PLOT  97
C                                                                       PLOT  98
C     IF WE HAVE A 1D PLOT NO NEED FOR SERCH                            PLOT  99
C                                                                       PLOT 100
      IF (IABS(ID)-3)122,305,122                                        PLOT 101
 122  DO 295 I=1,NY                                                     PLOT 102
      IY=I                                                              PLOT 103
      CALL SERCH(Y(1,IY),NPOIRL,INDEX,KODE,IY,MX)                       PLOT 104
      YMAXV(I)=Y(INDEX,I)                                               PLOT 105
      YMINV(I)=Y(KODE,I)                                                PLOT 106
      IND(I)=INDEX                                                      PLOT 107
      IF(NY.EQ.1) GO TO 300                                             PLOT 108
 295  KODE=1                                                            PLOT 109
      YMAX=AMAX1(YMAXV(1),YMAXV(2))                                     PLOT 110
      YMIN=AMIN1(YMINV(1),YMINV(2))                                     PLOT 111
      IF(YMAX.NE.YMAXV(1)) JY=2                                         PLOT 112
      YNOW=YMAX                                                         PLOT 113
      YNOW1=AMIN1(YMAXV(1),YMAXV(2))                                    PLOT 114
      GO TO 305                                                         PLOT 115
 300   YNOW=Y(INDEX,1)                                                  PLOT 116
      YNOW1=-XMAXN                                                      PLOT 117
      IF(YMAX-YMIN) 305,302,305                                         PLOT 118
 302  YMAX=Y(INDEX,1)                                                   PLOT 119
      YMIN=Y(KODE,1)                                                    PLOT 120
C                                                                       PLOT 121
C DETERMINE RANGE AND INCREMENT OF BOTH AXES                            PLOT 122
 305  P=ALOG10(XMAX-XMIN+.0001)                                         PLOT 123
      IP=IFIX(P)                                                        PLOT 124
      IF(P.LT.0.) IP=IP-1                                               PLOT 125
      SPAN=(XMAX-XMIN)/(10.**IP)                                        PLOT 126
      IF(SPAN.LT.5.0) GO TO 310                                         PLOT 127
      DEL=10.**IP                                                       PLOT 128
      GO TO 330                                                         PLOT 129
 310  IF(SPAN.LT.2.0) GO TO 320                                         PLOT 130
      DEL=(10.**IP)/2.0                                                 PLOT 131
      GO TO 330                                                         PLOT 132
 320  DEL=(10.**IP)/5.0                                                 PLOT 133
 330  IF(XMAX.LT.0.) GO TO 336                                          PLOT 134
      XMAX=FLOAT(IFIX(XMAX/DEL+.9999))*DEL                              PLOT 135
      IF(XMIN) 334,332,332                                              PLOT 136
 332  XMIN=FLOAT(IFIX(XMIN/DEL))*DEL                                    PLOT 137
      GO TO 338                                                         PLOT 138
 334  XMIN=FLOAT(IFIX(XMIN/DEL-.9999))*DEL                              PLOT 139
      GO TO 338                                                         PLOT 140
 336  XMAX=FLOAT(IFIX(XMAX/DEL))*DEL                                    PLOT 141
      XMIN=FLOAT(IFIX(XMIN/DEL+.9999))*DEL                              PLOT 142
 338  SPANX=XMAX-XMIN                                                   PLOT 143
      NXUNIT=SPANX/DEL+.00001                                           PLOT 144
      IF(NXUNIT.LE.10) GO TO 3338                                       PLOT 145
      LXFLG=LXFLG+1                                                     PLOT 146
      IF(LXFLG-1) 9999,305,9999                                         PLOT 147
 3338 NXP1=NXUNIT+1                                                     PLOT 148
      FINC=DEL                                                          PLOT 149
      DELX=SPANX/(100.-XS(NXUNIT))                                      PLOT 150
      IXS=IFIX(XS(NXUNIT)/2.+.1)                                        PLOT 151
      XSF=IXS                                                           PLOT 152
      IF(ID.LT.0) GO TO 339                                             PLOT 153
      IF((XMIN.GE.0.).OR.(XMAX.LE.0.) ) GO TO 339                       PLOT 154
      IX0L=-XMIN/DELX+1.5                                               PLOT 155
      IX0L=IX0L+IXS                                                     PLOT 156
C                                                                       PLOT 157
 339  IF(ID.EQ.3) GO TO 345                                             PLOT 158
      P=ALOG10(YMAX-YMIN+.0001)                                         PLOT 159
      IP=IFIX(P)                                                        PLOT 160
      IF(P.LT.0.) IP=IP-1                                               PLOT 161
      SPAN=(YMAX-YMIN)/(10.**IP)                                        PLOT 162
      IF(SPAN.LT.5.0) GO TO 340                                         PLOT 163
      DEL=10.**IP                                                       PLOT 164
      GO TO 360                                                         PLOT 165
 340  IF(SPAN.LT.2.0) GO TO 350                                         PLOT 166
      DEL=(10.**IP)/2.0                                                 PLOT 167
      GO TO 360                                                         PLOT 168
 350  DEL=(10.**IP)/5.0                                                 PLOT 169
 360  IF(YMAX.LT.0.) GO TO 366                                          PLOT 170
      YMAX=FLOAT(IFIX(YMAX/DEL+.9999))*DEL                              PLOT 171
      IF(YMIN) 364,362,362                                              PLOT 172
 362  YMIN=FLOAT(IFIX(YMIN/DEL))*DEL                                    PLOT 173
      GO TO 368                                                         PLOT 174
 364  YMIN=FLOAT(IFIX(YMIN/DEL-.9999))*DEL                              PLOT 175
      GO TO 368                                                         PLOT 176
 366  YMAX=FLOAT(IFIX(YMAX/DEL))*DEL                                    PLOT 177
      YMIN=FLOAT(IFIX(YMIN/DEL+.9999))*DEL                              PLOT 178
 368  SPANY = YMAX-YMIN                                                 PLOT 179
      NYUNIT=SPANY/DEL+.00001                                           PLOT 180
      IF(NYUNIT.LE.10) GO TO 3668                                       PLOT 181
      LYFLG=LYFLG+1                                                     PLOT 182
      IF(LYFLG-1)9999,339,9999                                          PLOT 183
 3668 DELY=SPANY/(FLOAT(NLPP-1)-YS(NYUNIT))                             PLOT 184
      DELTY=DEL                                                         PLOT 185
      VALUE=YMAX                                                        PLOT 186
      YININC=(NLPP-1)/NYUNIT                                            PLOT 187
      LININC=IFIX(YININC)                                               PLOT 188
      IYLAB=LININC                                                      PLOT 189
      IF(ID.LT.0) GO TO 369                                             PLOT 190
      IF((YMAX.LE.0.).OR.(YMIN.GE.0.)) GO TO 369                        PLOT 191
      IY0L=YMAX/DELY+1.5                                                PLOT 192
      GO TO 369                                                         PLOT 193
 345  VALUE=0.0                                                         PLOT 194
      DELTY=0.0                                                         PLOT 195
      LININC=27                                                         PLOT 196
      IYLAB=1                                                           PLOT 197
C                                                                       PLOT 198
C LABEL PLOT AT TOP                                                     PLOT 199
C                                                                       PLOT 200
 369  XLABL(1)=XMIN                                                     PLOT 201
      DO 120 I=1,NXUNIT                                                 PLOT 202
  120 XLABL(I+1)=XLABL(I)+FINC                                          PLOT 203
      GO TO (9999,9999,9999,365,370,375,380,385,390,395), NXUNIT        PLOT 204
 365  WRITE(LPRINT,3304)  (XLABL(I),I=1,NXP1)                           PLOT 205
      WRITE(LPRINT,3314)                                                PLOT 206
 3304 FORMAT(1H0,F20.4,4F25.4)                                          PLOT 207
 3314 FORMAT(1H ,14X,103H*.************************.********************PLOT 208
     1****.************************.************************.*)         PLOT 209
      GO TO 400                                                         PLOT 210
 370  WRITE(LPRINT,3305)  (XLABL(I),I=1,NXP1)                           PLOT 211
      WRITE(LPRINT,3315)                                                PLOT 212
 3305 FORMAT(1H0,6F20.4)                                                PLOT 213
 3315 FORMAT(1H ,14X,103H*.*******************.*******************.*****PLOT 214
     2**************.*******************.*******************.*)         PLOT 215
      GO TO 400                                                         PLOT 216
 375  WRITE(LPRINT,3306)  (XLABL(I),I=1,NXP1)                           PLOT 217
      WRITE(LPRINT,3316)                                                PLOT 218
 3306 FORMAT(1H0,6X,7F16.4)                                             PLOT 219
 3316 FORMAT(1H ,14X,103H***.***************.***************.***********PLOT 220
     1****.***************.***************.***************.***)         PLOT 221
      GO TO 400                                                         PLOT 222
 380  WRITE(LPRINT,3307)  (XLABL(I),I=1,NXP1)                           PLOT 223
      WRITE(LPRINT,3317)                                                PLOT 224
 3307 FORMAT(1H0,7X,8F14.4)                                             PLOT 225
 3317 FORMAT(1H ,14X,103H**.*************.*************.*************.**PLOT 226
     1***********.*************.*************.*************.**)         PLOT 227
      GO TO 400                                                         PLOT 228
 385  WRITE(LPRINT,3308)  (XLABL(I),I=1,NXP1)                           PLOT 229
      WRITE(LPRINT,3318)                                                PLOT 230
 3308 FORMAT(1H0,10X,9F12.4)                                            PLOT 231
 3318 FORMAT(1H ,14X,103H***.***********.***********.***********.*******PLOT 232
     1****.***********.***********.***********.***********.***)         PLOT 233
      GO TO 400                                                         PLOT 234
 390  WRITE(LPRINT,3309)  (XLABL(I),I=1,NXP1)                           PLOT 235
      WRITE(LPRINT,3319)                                                PLOT 236
 3309 FORMAT(1H0,8X,10F11.3)                                            PLOT 237
 3319 FORMAT(1H ,14X,103H*.**********.**********.**********.**********.*PLOT 238
     1*********.**********.**********.**********.**********.**)         PLOT 239
      GO TO 400                                                         PLOT 240
 395  WRITE(LPRINT,3310)  (XLABL(I),I=1,NXP1)                           PLOT 241
      WRITE(LPRINT,3320)                                                PLOT 242
 3310 FORMAT(1H0,9X,11F10.3)                                            PLOT 243
 3320 FORMAT(1H ,14X 103H*.*********.*********.*********.*********.*****PLOT 244
     1****.*********.*********.*********.*********.*********.*)         PLOT 245
 400  CONTINUE                                                          PLOT 246
C                                                                       PLOT 247
C     PREPARE PLOT                                                      PLOT 248
C                                                                       PLOT 249
      KODE=0                                                            PLOT 250
 4000 IF (LINE-NLPP)130,130,2000                                        PLOT 251
  130 IF (IDAT-1)700,701,700                                            PLOT 252
  700 IF(IABS(ID)-3)706,701,706                                         PLOT 253
  706 IF (LFLAG)702,702,701                                             PLOT 254
 702  IF(IDAT-NPOI2) 703,703,701                                        PLOT 255
 703  CALL SERCH(Y(1,JY),NPOIRL,INDEX,KODE,JY,MX)                       PLOT 256
      IF(INDEX.EQ.(-1)) GO TO 7704                                      PLOT 257
      IND(JY)=INDEX                                                     PLOT 258
      YNOW=Y(INDEX,JY)                                                  PLOT 259
      IF(YNOW-YNOW1) 7702,701,701                                       PLOT 260
 7704 JY=NY-JY+1                                                        PLOT 261
      YNOW=YNOW1                                                        PLOT 262
      YNOW1=-XMAXN                                                      PLOT 263
      GO TO 701                                                         PLOT 264
 7702 JY=NY-JY+1                                                        PLOT 265
      TEMP=YNOW                                                         PLOT 266
      YNOW=YNOW1                                                        PLOT 267
      YNOW1=TEMP                                                        PLOT 268
  701 CONTINUE                                                          PLOT 269
C                                                                       PLOT 270
C     CHECK TO SEE IF ANY POINTS IN Y REMAIN                            PLOT 271
C                                                                       PLOT 272
      IF(IDAT-NPOI2) 806,806,1000                                       PLOT 273
  806 IF (IABS(ID)-3)805,808,805                                        PLOT 274
 805  IF (YMAX-YNOW)811,801,801                                         PLOT 275
  801 IF (YNOW-YMIN)811,802,802                                         PLOT 276
  802 I=(YMAX-YNOW)/DELY+1.50001                                        PLOT 277
      GOTO 807                                                          PLOT 278
  808 IF (LINE-27)807,809,807                                           PLOT 279
 809  DO 10 II=MX,NPOIRL                                                PLOT 280
      I=II                                                              PLOT 281
 10   CALL XITEM(I,XMAX,XMIN,ID,DELX,X,JY,XSF)                          PLOT 282
      GOTO 606                                                          PLOT 283
  807 IF (I-LINE)820,820,810                                            PLOT 284
  810 LFLAG=1                                                           PLOT 285
      GOTO 1000                                                         PLOT 286
  811 LFLAG=0                                                           PLOT 287
      GOTO 606                                                          PLOT 288
  820 LFLAG=0                                                           PLOT 289
      INDEX=IND(JY)                                                     PLOT 290
      CALL XITEM(INDEX,XMAX,XMIN,ID,DELX,X,JY,XSF)                      PLOT 291
  606 IDAT=IDAT+1                                                       PLOT 292
      GOTO 4000                                                         PLOT 293
C                                                                       PLOT 294
C PRINT ONE LINE OF PLOT                                                PLOT 295
C                                                                       PLOT 296
C     CHECK ID -- IF MINUS NO AXES,IF + AXES                            PLOT 297
C                IF ID=3 1D PLOT                                        PLOT 298
C                                                                       PLOT 299
 1000 IF (ID)133,901,901                                                PLOT 300
 901  IF(IABS(ID)-3) 902,133,902                                        PLOT 301
 902  IF(IX0L) 907,907,903                                              PLOT 302
 903  IF(ITEM(IX0L)-BLANK) 907,906,907                                  PLOT 303
 906  ITEM(IX0L)=AIE                                                    PLOT 304
 907  IF(IY0L) 133,133,909                                              PLOT 305
 909  IF(LINE-IY0L) 133,904,133                                         PLOT 306
  904 DO 905 I=1,101                                                    PLOT 307
      IF (ITEM(I)-BLANK)905,908,905                                     PLOT 308
  908 ITEM(I)=MINUS                                                     PLOT 309
  905 CONTINUE                                                          PLOT 310
  133 CONTINUE                                                          PLOT 311
C                                                                       PLOT 312
      IF(IYLAB.NE.LININC) GO TO 915                                     PLOT 313
      WRITE(LPRINT,3300) VALUE,DD,(ITEM(J),J=1,101),DD,VALUE            PLOT 314
 3300  FORMAT(1H ,F12.4,2X,103A1,F13.4)                                 PLOT 315
      VALUE=VALUE-DELTY                                                 PLOT 316
      IYLAB=0                                                           PLOT 317
      GO TO 920                                                         PLOT 318
 915  WRITE(LPRINT,3410) STAR, (ITEM(J),J=1,101),STAR                   PLOT 319
 3410 FORMAT(1H ,14X,103A1)                                             PLOT 320
 920  IYLAB=IYLAB+1                                                     PLOT 321
      DO 3350 I=1,101                                                   PLOT 322
 3350 ITEM(I)=BLANK                                                     PLOT 323
      LINE=LINE+1                                                       PLOT 324
      GOTO 4000                                                         PLOT 325
C                                                                       PLOT 326
C LABEL PLOT ALONG THE BOTTOM                                           PLOT 327
C                                                                       PLOT 328
 2000 GO TO (9999,9999,9999,965,970,975,980,985,990,995), NXUNIT        PLOT 329
 965  WRITE(LPRINT,3314)                                                PLOT 330
      WRITE(LPRINT,3304) (XLABL(I),I=1,NXP1)                            PLOT 331
      GO TO 998                                                         PLOT 332
 970  WRITE(LPRINT,3315)                                                PLOT 333
      WRITE(LPRINT,3305) (XLABL(I),I=1,NXP1)                            PLOT 334
      GO TO 998                                                         PLOT 335
 975  WRITE(LPRINT,3316)                                                PLOT 336
      WRITE(LPRINT,3306) (XLABL(I),I=1,NXP1)                            PLOT 337
      GO TO 998                                                         PLOT 338
 980  WRITE(LPRINT,3317)                                                PLOT 339
      WRITE(LPRINT,3307) (XLABL(I),I=1,NXP1)                            PLOT 340
      GO TO 998                                                         PLOT 341
 985  WRITE(LPRINT,3318)                                                PLOT 342
      WRITE(LPRINT,3308) (XLABL(I),I=1,NXP1)                            PLOT 343
      GO TO 998                                                         PLOT 344
 990  WRITE(LPRINT,3319)                                                PLOT 345
      WRITE(LPRINT,3309) (XLABL(I),I=1,NXP1)                            PLOT 346
      GO TO 998                                                         PLOT 347
 995  WRITE(LPRINT,3320)                                                PLOT 348
      WRITE(LPRINT,3310) (XLABL(I),I=1,NXP1)                            PLOT 349
 998  CONTINUE                                                          PLOT 350
C                                                                       PLOT 351
      RETURN                                                            PLOT 352
C                                                                       PLOT 353
 9999 WRITE(LPRINT,9995)                                                PLOT 354
 9995 FORMAT(22H0ERROR EXIT FROM PLOT.)                                 PLOT 355
      STOP                                                              PLOT 356
C                                                                       PLOT 357
      END                                                               PLOT 358
$      FORTRAN
CREGR                                                                   REGR   1
       FUNCTION REGR(DA,I)                                              REGR   2
C REGR             FOR KYST                  JANUARY,1973               REGR   3
C WRITTEN BY J.KRUSKAL                                                  REGR   4
CREGR          CALCULATES VALUES OF FUNCTIONS REGRESSED OVER            REGR   5
C                                                                       REGR   6
       J=MIN0(I,4)                                                      REGR   7
       GO TO (10,20,30,40), J                                           REGR   8
C                                                                       REGR   9
 10    REGR=1.0                                                         REGR  10
       RETURN                                                           REGR  11
C                                                                       REGR  12
 20    REGR=DA                                                          REGR  13
       RETURN                                                           REGR  14
C                                                                       REGR  15
 30    REGR=DA*DA                                                       REGR  16
       RETURN                                                           REGR  17
C                                                                       REGR  18
 40    REGR=DA**(I-1)                                                   REGR  19
       RETURN                                                           REGR  20
C                                                                       REGR  21
       END                                                              REGR  22
$      FORTRAN
CRM1POW                                                                 RM1POW 1
      FUNCTION RM1POW(YY)                                               RM1POW 2
C RM1POW           FOR KYST                  JANUARY,1973               RM1POW 3
C WRITTEN BY J.KRUSKAL                                                  RM1POW 4
C MODIFIED FOR KYST BY J.SEERY               JANUARY,1973               RM1POW 5
C      RM1POW--CALCULATES (R-1)-ST POWER                                RM1POW 6
C                                                                       RM1POW 7
      INTEGER RTYPE                                                     RM1POW 8
      COMMON /METRIC/ RTYPE,RM1,RECR,R                                  RM1POW 9
      Y=YY                                                              RM1POW10
      GO TO (210,220,230), RTYPE                                        RM1POW11
 210  RM1POW=0.0                                                        RM1POW12
      IF(Y.NE.0.0) RM1POW=SIGN(1.0,Y)                                   RM1POW13
      RETURN                                                            RM1POW14
 220  RM1POW=Y                                                          RM1POW15
      RETURN                                                            RM1POW16
 230  RM1POW=SIGN(ABS(Y)**RM1,Y)                                        RM1POW17
      RETURN                                                            RM1POW18
      END                                                               RM1POW19
$      FORTRAN
CRPOWER                                                                 RPOWER 1
      FUNCTION RPOWER(XX)                                               RPOWER 2
C RPOWER           FOR KYST                  JANUARY,1973               RPOWER 3
C WRITTEN BY J.KRUSKAL                                                  RPOWER 4
C MODIFIED FOR KYST BY J.SEERY               JANUARY,1973               RPOWER 5
C      RPOWER--CALCULATES R-TH POWER                                    RPOWER 6
C                                                                       RPOWER 7
      INTEGER RTYPE                                                     RPOWER 8
      COMMON /METRIC/ RTYPE,RM1,RECR,R                                  RPOWER 9
      X=XX                                                              RPOWER10
      GO TO (110,120,130), RTYPE                                        RPOWER11
 110  RPOWER=ABS(X)                                                     RPOWER12
      RETURN                                                            RPOWER13
 120  RPOWER=X*X                                                        RPOWER14
      RETURN                                                            RPOWER15
 130  RPOWER=ABS(X)**R                                                  RPOWER16
      RETURN                                                            RPOWER17
      END                                                               RPOWER18
$      FORTRAN
CRROOT                                                                  RROOT  1
      FUNCTION RROOT(ZZ)                                                RROOT  2
C RROOT            FOR KYST                  JANUARY,1973               RROOT  3
C WRITTEN BY J.KRUSKAL                                                  RROOT  4
C MODIFIED FOR KYST BY J.SEERY               JANUARY,1973               RROOT  5
C      RROOT--CALCULATES R-TH ROOT                                      RROOT  6
C                                                                       RROOT  7
      INTEGER RTYPE                                                     RROOT  8
      COMMON /METRIC/ RTYPE,RM1,RECR,R                                  RROOT  9
      Z=ZZ                                                              RROOT 10
      GO TO (310,320,330), RTYPE                                        RROOT 11
 310  RROOT=Z                                                           RROOT 12
      RETURN                                                            RROOT 13
 320  RROOT=SQRT(Z)                                                     RROOT 14
      RETURN                                                            RROOT 15
 330  RROOT=Z**RECR                                                     RROOT 16
      RETURN                                                            RROOT 17
      END                                                               RROOT 18
$      FORTRAN
CRUNIFV                                                                 RUNIFV 1
       FUNCTION RUNIFV(A)                                               RUNIFV 2
C RUNIFV           FOR KYST                  JANUARY,1973               RUNIFV 3
C WRITTEN BY J.KRUSKAL                                                  RUNIFV 4
C                                                                       RUNIFV 5
C      THIS IS A SIMPLE ROUTINE FOR GENERATING RANDOM NUMBERS.          RUNIFV 6
C      THEY ARE UNIFORMLY DISTRIBUTED BETWEEN 0 AND 1.                  RUNIFV 7
C      IT IS NOT FAST, NOR DOES IT PRODUCE VERY HIGH QUALITY NUMBERS.   RUNIFV 8
C      ITS MAIN VIRTUE IS THAT IT IS IN FORTRAN AND WILL WORK ON        RUNIFV 9
C      ALMOST ANY MACHINE ON WHICH FORTRAN IS AVAILABLE.                RUNIFV10
C                                                                       RUNIFV11
       DATA MODULO,FLMOD,K/8192,8192.0,1/                               RUNIFV12
C      MODULO   2**13                                                   RUNIFV13
C                                                                       RUNIFV14
       DO 10 I=1,15                                                     RUNIFV15
 10    K = MOD(5*K, MODULO)                                             RUNIFV16
       Z = FLOAT(K)/FLMOD                                               RUNIFV17
       RUNIFV = Z                                                       RUNIFV18
       RETURN                                                           RUNIFV19
C                                                                       RUNIFV20
       END                                                              RUNIFV21
$      FORTRAN
CSBK                                                                    SBK    1
      SUBROUTINE SBK(N,K,A,E,Z)                                         SBK    2
C SBK              FOR KYST                  JANUARY,1973               SBK    3
C TRANSLATED FROM ALGOL AND ADAPTED                                     SBK    4
C FOR SGEV BY P.BUSINGER                     JANUARY,1972               SBK    5
      REAL A(1),E(N),Z(100,K)                                           SBK    6
C                                                                       SBK    7
C R.S.MARTIN, C.REINSCH, AND J.H.WILKINSON,                             SBK    8
C HOUSEHOLDER-S TRIDIAGONALIZATION OF A SYM-                            SBK    9
C METRIC MATRIX. NUMER.MATH.11,181-195(1968).                           SBK   10
C                                                                       SBK   11
C SBK IS CALLED BY SGEV, AND                                            SBK   12
C USES THE INFORMATION PROVIDED BY THE SUB-                             SBK   13
C ROUTINE STO3 IN A AND E TO BACKTRANSFORM                              SBK   14
C THE EIGENVECTOR MATRIX Z.                                             SBK   15
C                                                                       SBK   16
      REAL H,S                                                          SBK   17
      IF(N.LT.2)GOTO 50                                                 SBK   18
      DO 40 I=2,N                                                       SBK   19
         IF(E(I).EQ.0.E0)GOTO 40                                        SBK   20
         L=I-1                                                          SBK   21
         MM=I*(I-1)/2+I-1                                               SBK   22
         H=E(I)*A(MM)                                                   SBK   23
         DO 30 J=1,K                                                    SBK   24
            S=0.E0                                                      SBK   25
            DO 10 M=1,L                                                 SBK   26
                  MM=I*(I-1)/2+M                                        SBK   27
   10             S=S+A(MM)*Z(M,J)                                      SBK   28
            S=S/H                                                       SBK   29
            DO 20 M=1,L                                                 SBK   30
               MM=I*(I-1)/2+M                                           SBK   31
   20          Z(M,J)=Z(M,J)+S*A(MM)                                    SBK   32
   30       CONTINUE                                                    SBK   33
   40    CONTINUE                                                       SBK   34
   50 RETURN                                                            SBK   35
      END                                                               SBK   36
$      FORTRAN
CSERCH                                                                  SERCH  1
      SUBROUTINE  SERCH(VEC,LEN,INDEX,KODE,JY,MX)                       SERCH  2
C SERCH            FOR KYST                  JANUARY,1973               SERCH  3
C     OCTOBER 10,1971     JAY R LEVINSOHN                               SERCH  4
C MODIFIED FOR KYST BY J.SEERY               JANUARY,1973               SERCH  5
C                                                                       SERCH  6
C     SEARCH SERVES TWO FUNCTIONS                                       SERCH  7
C        1 FIND THE MIN AND MAX OF SOME VECTOR VEC                      SERCH  8
C        2 FIND THE NEXT BIGGEST ITEM IN SOME VECTOR VEC                SERCH  9
C      CAN KEEP TRACK OF TWO VECTORS (OF THE SAME LENGTH) AT A TIME     SERCH 10
C     PARAMETERS                                                        SERCH 11
C        VEC-- VECTOR TO BE SEARCHED                                    SERCH 12
C        LEN-- NUMBER OF ELEMENTS IN VEC                                SERCH 13
C        MX--      WORK WITH ELEMENTS MX THROUGH LEN                    SERCH 14
C        INDEX-- POINTER TO NEXT BIGGEST ITEM IN VEC                    SERCH 15
C        KODE-- OP. CODE 1=FIND MIN AND MAX,0= FIND NEXT BIGGEST ITEM   SERCH 16
C     IF KODE=1 THEN POINTER TO MIN IS RETURNED IN KODE WHILE POINTER   SERCH 17
C     TO MAX IS RETURNED IN INDEX                                       SERCH 18
C     ROUTINE ASSUMES IT WILL BE CALLED TO FIRST FIND MIN AND MAX OF    SERCH 19
C     VECTOR AND THEN CALLED ITERATIVELY TO FIND THE SECOND LARGEST     SERCH 20
C     ELEMENT, THIRD LARGEST ELEMENT,, ETC.... UNTIL  FINAL ELEMENT HAS SERCH 21
C     HAS BEEN EXTRACTED.  ONE CAN TEST THE FINAL ELEMENT AGAINST THE   SERCH 22
C     MIN TO CHECK FOR ACCURATE COMPLETION.                             SERCH 23
C                                                                       SERCH 24
      DIMENSION VEC(1)                                                  SERCH 25
      DIMENSION FMIN(2),IMIN(2),INDB4(2),KNT(2)                         SERCH 26
      COMMON /MDSCL1/ LREAD,LPRINT,LPUNCH,LSCRAT                        SERCH 27
      COMMON /ACCUR/ PRECSN,XMAG,XMAXN                                  SERCH 28
C                                                                       SERCH 29
      IF  (KODE)500,200,100                                             SERCH 30
C                                                                       SERCH 31
C     OP CODE =1 THEN FIND MIN AND MAX                                  SERCH 32
C                                                                       SERCH 33
 100  FMAX=-XMAXN                                                       SERCH 34
      FMIN(JY)=XMAXN                                                    SERCH 35
      KNT(JY)=MX                                                        SERCH 36
      DO 150 I=MX,LEN                                                   SERCH 37
      IF(VEC(I)-FMAX) 160,155,155                                       SERCH 38
 155  FMAX=VEC(I)                                                       SERCH 39
         INDEX=I                                                        SERCH 40
 160   IF(VEC(I)-FMIN(JY)) 165,165,150                                  SERCH 41
 165  FMIN(JY)=VEC(I)                                                   SERCH 42
         KODE=I                                                         SERCH 43
  150 CONTINUE                                                          SERCH 44
C                                                                       SERCH 45
C     SAVE MIN IN SAVE,MIN POINTER IN IMIN,RESET FMIN, SAVE MAX         SERCH 46
C                                                                       SERCH 47
      IMIN(JY)=KODE                                                     SERCH 48
      FMIN(JY)=-XMAXN                                                   SERCH 49
      INDB4(JY)=INDEX                                                   SERCH 50
      RETURN                                                            SERCH 51
C                                                                       SERCH 52
C     OP. CODE NOT EQUAL 1 FIND NEXT LARGEST ELE IN VEC                 SERCH 53
C                                                                       SERCH 54
 200  IF(KNT(JY).LT.LEN) GO TO 210                                      SERCH 55
      INDEX=-1                                                          SERCH 56
      GO TO 350                                                         SERCH 57
 210  INDEX=INDB4(JY)                                                   SERCH 58
      FMB4=VEC(INDEX)                                                   SERCH 59
      VEC(INDEX)=FMB4+1.0                                               SERCH 60
      KNT(JY)=KNT(JY)+1                                                 SERCH 61
      DO 250 I=MX,LEN                                                   SERCH 62
      IF(VEC(I)-FMIN(JY)) 250,255,255                                   SERCH 63
 255  IF(VEC(I)-FMB4) 260,260,250                                       SERCH 64
 260  IF(I-INDB4(JY)) 265,250,265                                       SERCH 65
 265  FMIN(JY)=VEC(I)                                                   SERCH 66
         INDEX=I                                                        SERCH 67
  250 CONTINUE                                                          SERCH 68
C                                                                       SERCH 69
C     SUBTRACT 1 FROM EACH ELEMENT OF VEC TO GET IT BACK AS BEFORE      SERCH 70
C     DON'T DO SUBTRACTION UNLESS AT LAST ITEM OF VEC                   SERCH 71
C                                                                       SERCH 72
      IF(KNT(JY)-LEN) 325,320,325                                       SERCH 73
 320  DO 330 I=MX,LEN                                                   SERCH 74
         IF (I-INDEX)335,330,335                                        SERCH 75
  335    VEC(I)=VEC(I)-1.0                                              SERCH 76
  330 CONTINUE                                                          SERCH 77
C                                                                       SERCH 78
C     SAVE INDEX AND RESET FMIN                                         SERCH 79
C                                                                       SERCH 80
 325  INDB4(JY)=INDEX                                                   SERCH 81
      FMIN(JY)=-XMAXN                                                   SERCH 82
 350  RETURN                                                            SERCH 83
C                                                                       SERCH 84
C     ERROR ROUTINE ENTER HERE IF KODE<0                                SERCH 85
C                                                                       SERCH 86
 500  WRITE(LPRINT,505) KODE                                            SERCH 87
  505 FORMAT(24H ERROR IN SERCH, KODE= ,I4)                             SERCH 88
      RETURN                                                            SERCH 89
      END                                                               SERCH 90
$      FORTRAN
CSGEV                                                                   SGEV   1
      SUBROUTINE SGEV(N,K,Z)                                            SGEV   2
C SGEV             FOR KYST                  JANUARY,1973               SGEV   3
C PROGRAMMED BY P.BUSINGER                   JANUARY,1972               SGEV   4
C CALLS BSEC1,NVIT1,STO3,DFLT,SBK                                       SGEV   5
      DIMENSION B(100),C(100),D(100),U(100),STORE(194),E(6),            SGEV   6
     .     XBEST(100,6),A(5430),Z(100,6)                                SGEV   7
      COMMON /KYST2/ B,C,D,U,STORE,E,XBEST,A                            SGEV   8
      COMMON /ACCUR/ EPS,ETA,OMEGA                                      SGEV   9
C                                                                       SGEV  10
C FINDS THE K ALGEBRAICALLY GREATEST EIGENVALUES                        SGEV  11
C     E(1) .GE. E(2) .GE. ... .GE. E(K)                                 SGEV  12
C OF THE N BY N REAL SYMMETRIC MATRIX WHOSE LOWER                       SGEV  13
C TRIANGULAR PART IS GIVEN ROW-COMPACTED IN A.                          SGEV  14
C ALSO PRODUCED IS                                                      SGEV  15
C AN ASSOCIATED ORTHOGONAL SET Z OF EIGENVECTORS.                       SGEV  16
C                                                                       SGEV  17
C                - CAUTION -                                            SGEV  18
C NO A PRIORI ERROR BOUNDS ARE AVAILABLE AND NO                         SGEV  19
C A POSTERIORI CHECKS ARE PERFORMED ON THE RESULTS                      SGEV  20
C COMPUTED BY THIS SUBROUTINE. ONE CHECK WHICH A                        SGEV  21
C USER CAN PERFORM INEXPENSIVELY CONSISTS IN COMPU-                     SGEV  22
C TING THE RESIDUAL MATRICES AZ-ZE AND (ZT)Z-I                          SGEV  23
C WHOSE ELEMENTS SHOULD BE NEGLIGIBLE.                                  SGEV  24
C                                                                       SGEV  25
      REAL S,SN,CS,T,Z1,Z2                                              SGEV  26
      IF(K.LT.1)GOTO 99                                                 SGEV  27
      IF(N.LT.1)GOTO 99                                                 SGEV  28
      E(1)=A(1)                                                         SGEV  29
      Z(1,1)=1.E0                                                       SGEV  30
      IF(N.EQ.1)GOTO 99                                                 SGEV  31
      CALL STO3(N,ETA/EPS,A,D,B)                                        SGEV  32
      DO 10 I=2,N                                                       SGEV  33
   10    C(I-1)=B(I)                                                    SGEV  34
      IB=5.E0-ALOG(EPS)/ALOG(2.E0)                                      SGEV  35
      S=1.E0                                                            SGEV  36
      DO 30 J=1,K                                                       SGEV  37
         M=N+1-J                                                        SGEV  38
         CALL BSEC1(IB,M,ETA,D,C,U,E(J),S)                              SGEV  39
         DO 20 I=1,N                                                    SGEV  40
   20       Z(I,J)=0.E0                                                 SGEV  41
         CALL NVIT1(M,OMEGA,C,U,Z(1,J))                                 SGEV  42
         IF(J.NE.K)CALL DFLT(M,D,C,E(J),Z(1,J))                         SGEV  43
   30    E(J)=S*E(J)                                                    SGEV  44
      IF(K.EQ.1)GOTO 80                                                 SGEV  45
      DO 70 JJ=2,K                                                      SGEV  46
         J=K+1-JJ                                                       SGEV  47
         NJ=N-K+JJ                                                      SGEV  48
         Z(NJ,J)=1.E0                                                   SGEV  49
         NJ=N-J                                                         SGEV  50
         DO 60 II=1,NJ                                                  SGEV  51
            I=NJ+1-II                                                   SGEV  52
            T=Z(I,J)                                                    SGEV  53
            IF(T.EQ.0.E0)GOTO 60                                        SGEV  54
            Z(I,J)=0.E0                                                 SGEV  55
            SN=(T+T)/(1.E0+T*T)                                         SGEV  56
            CS=(1.E0-T)*(1.E0+T)/(1.E0+T*T)                             SGEV  57
            DO 50 M=J,K                                                 SGEV  58
               Z1=Z(I,M)                                                SGEV  59
               Z2=Z(I+1,M)                                              SGEV  60
               Z(I,M)=CS*Z1-SN*Z2                                       SGEV  61
   50          Z(I+1,M)=CS*Z2+SN*Z1                                     SGEV  62
   60       CONTINUE                                                    SGEV  63
   70    CONTINUE                                                       SGEV  64
   80 CALL SBK(N,K,A,B,Z)                                               SGEV  65
   99 RETURN                                                            SGEV  66
      END                                                               SGEV  67
$      FORTRAN
CSORT                                                                   SORT   1
       SUBROUTINE SORT (A, N, B, C, D, K, SWITCH )                      SORT   2
C SORT             FOR KYST                  JANUARY,1973               SORT   3
C WRITTEN BY J.KRUSKAL                                                  SORT   4
C                                                                       SORT   5
C     THIS ROUTINE SORTS INPUT ARRAY 'A' AND REARRANGES, OPTIONALLY,    SORT   6
C     ARRAYS 'B', 'C', AND 'D', IN ORDER CORRESPONDING TO 'A'.          SORT   7
C     N = NUMBER OF ITEMS IN 'A' (AND 'B', 'C', 'D', IF USED)           SORT   8
C     K = 0--SORT 'A' ONLY, 1--REARRANGE 'B', 2--REARRANGE 'B' AND 'C', SORT   9
C         3--REARRANGE 'B', 'C', AND 'D'.                               SORT  10
C     IF 'SWITCH' IS POSITIVE, SORT WILL BE IN ASCENDING ORDER,         SORT  11
C                 IF ZERO OR NEGATIVE, IN DESCENDING ORDER.             SORT  12
C     ALGORITHM FROM CACM, JULY 1959, PAGE 30 BY D. L. SHELL            SORT  13
C                                                                       SORT  14
      DIMENSION A(1),B(1),C(1),D(1)                                     SORT  15
       INTEGER  SWITCH                                                  SORT  16
 105   KP1=K+1                                                          SORT  17
       IF(N.LE.1) GO TO 999                                             SORT  18
      M = 1                                                             SORT  19
  106  M = M + M                                                        SORT  20
       IF( M .LE. N ) GO TO 106                                         SORT  21
       M = M - 1                                                        SORT  22
 994       M = M/2                                                      SORT  23
       IF(M.EQ.0) GO TO 999                                             SORT  24
       KK = N-M                                                         SORT  25
       J = 1                                                            SORT  26
992     I = J                                                           SORT  27
996     IM = I + M                                                      SORT  28
       IF(SWITCH)  810,810,800                                          SORT  29
800     IF(A(I).GT.A(IM)) GO TO 110                                     SORT  30
       GO TO 995                                                        SORT  31
810    IF(A(I).LT.A(IM))  GO TO 110                                     SORT  32
995     J = J+1                                                         SORT  33
       IF(J.GT.KK)  GO TO 994                                           SORT  34
       GO TO 992                                                        SORT  35
 110   TEMP=A(I)                                                        SORT  36
       A(I) = A(IM)                                                     SORT  37
       A(IM) = TEMP                                                     SORT  38
       GO TO ( 140, 130, 120, 115), KP1                                 SORT  39
 115   TEMP = D(I)                                                      SORT  40
       D(I) = D(IM)                                                     SORT  41
        D(IM) = TEMP                                                    SORT  42
 120   TEMP=C(I)                                                        SORT  43
       C(I) = C(IM)                                                     SORT  44
       C(IM) = TEMP                                                     SORT  45
 130   TEMP=B(I)                                                        SORT  46
       B(I) = B(IM)                                                     SORT  47
       B(IM) = TEMP                                                     SORT  48
140       I = I-M                                                       SORT  49
       IF(I.LT.1)  GO TO 995                                            SORT  50
        GO TO 996                                                       SORT  51
999      RETURN                                                         SORT  52
      END                                                               SORT  53
$      FORTRAN
CST03                                                                   ST03   1
      SUBROUTINE STO3(N,TOL,A,D,E)                                      ST03   2
C STO3             FOR KYST                  JANUARY,1973               ST03   3
C TRANSLATED FROM ALGOL AND ADAPTED                                     ST03   4
C FOR SGEV BY P.BUSINGER                     JANUARY,1972               ST03   5
      REAL TOL,A(1),D(N),E(N)                                           ST03   6
C                                                                       ST03   7
C R.S.MARTIN, C.REINSCH, AND J.H.WILKINSON,                             ST03   8
C HOUSEHOLDER-S TRIDIAGONALIZATION OF A SYM-                            ST03   9
C METRIC MATRIX. NUMER.MATH.11,181-195(1968).                           ST03  10
C                                                                       ST03  11
C STO3 IS CALLED BY SGEV, AND                                           ST03  12
C USES HOUSEHOLDER TRANSFORMATIONS TO REDUCE                            ST03  13
C THE REAL SYMMETRIC N BY N MATRIX WHOSE                                ST03  14
C LOWER TRIANGULAR PART IS GIVEN ROW-COMPACTED                          ST03  15
C IN A TO TRIDIAGONAL FORM. THE DIAGONAL ELEMENTS                       ST03  16
C OF THE REDUCED MATRIX ARE                                             ST03  17
C        D(1) THROUGH D(N)                                              ST03  18
C AND THE OFF-DIAGONAL ELEMENTS ARE                                     ST03  19
C        E(2) THROUGH E(N).                                             ST03  20
C TOGETHER WITH THE E-S, THE INFORMATION                                ST03  21
C LEFT IN A ALLOWS THE SUBROUTINE SBAK TO                               ST03  22
C APPLY THE REDUCING TRANSFORMATION TO THE                              ST03  23
C EIGENVECTORS OF THE TRIDIAGONAL MATRIX.                               ST03  24
C TOL MUST BE GIVEN AS (LEAST POSITIVE MA-                              ST03  25
C CHINE NUMBER)/(RELATIVE MACHINE PRECISION).                           ST03  26
C                                                                       ST03  27
      REAL F,G,H                                                        ST03  28
      DO 10 I=1,N                                                       ST03  29
         M=I*(I-1)/2+I                                                  ST03  30
   10    D(I)=A(M)                                                      ST03  31
      DO 80 II=1,N                                                      ST03  32
         I=N+1-II                                                       ST03  33
         L=I-1                                                          ST03  34
         H=0.E0                                                         ST03  35
         IF(L.LT.1)GOTO 21                                              ST03  36
         DO 20 K=1,L                                                    ST03  37
            M=I*(I-1)/2+K                                               ST03  38
   20       H=H+A(M)**2                                                 ST03  39
   21    IF(H.GT.TOL)GOTO 30                                            ST03  40
         E(I)=0.E0                                                      ST03  41
         GOTO 72                                                        ST03  42
   30    M=I*(I-1)/2+I-1                                                ST03  43
         F=A(M)                                                         ST03  44
         G=SQRT(H)                                                      ST03  45
         IF(F.GE.0.E0)G=-G                                              ST03  46
         E(I)=G                                                         ST03  47
         H=H-F*G                                                        ST03  48
         A(M)=F-G                                                       ST03  49
         F=0.E0                                                         ST03  50
         IF(L.LT.1)GOTO 61                                              ST03  51
         DO 60 J=1,L                                                    ST03  52
            G=0.E0                                                      ST03  53
            IF(J.LT.1)GOTO 41                                           ST03  54
            DO 40 K=1,J                                                 ST03  55
               M1=J*(J-1)/2+K                                           ST03  56
               M2=I*(I-1)/2+K                                           ST03  57
   40          G=G+A(M1)*A(M2)                                          ST03  58
   41       J1=J+1                                                      ST03  59
            IF(L.LT.J1)GOTO 51                                          ST03  60
            DO 50 K=J1,L                                                ST03  61
               M1=K*(K-1)/2+J                                           ST03  62
               M2=I*(I-1)/2+K                                           ST03  63
   50          G=G+A(M1)*A(M2)                                          ST03  64
   51       G=G/H                                                       ST03  65
            E(J)=G                                                      ST03  66
            M=I*(I-1)/2+J                                               ST03  67
   60       F=F+G*A(M)                                                  ST03  68
   61    H=F/(H+H)                                                      ST03  69
         IF(L.LT.1)GOTO 72                                              ST03  70
         DO 71 J=1,L                                                    ST03  71
            M=I*(I-1)/2+J                                               ST03  72
            F=A(M)                                                      ST03  73
            G=E(J)-H*F                                                  ST03  74
            E(J)=G                                                      ST03  75
            IF(J.LT.1)GOTO 71                                           ST03  76
            DO 70 K=1,J                                                 ST03  77
               M1=J*(J-1)/2+K                                           ST03  78
               M2=I*(I-1)/2+K                                           ST03  79
   70          A(M1)=A(M1)-F*E(K)-G*A(M2)                               ST03  80
   71       CONTINUE                                                    ST03  81
   72    H=D(I)                                                         ST03  82
         M=I*(I-1)/2+I                                                  ST03  83
         D(I)=A(M)                                                      ST03  84
   80    A(M)=H                                                         ST03  85
      RETURN                                                            ST03  86
      END                                                               ST03  87
$      FORTRAN
CWTRAN                                                                  WTRAN  1
      FUNCTION WTRAN(XX)                                                WTRAN  2
C WTRAN            FOR KYST                  JANUARY,1973               WTRAN  3
C WRITTEN FOR KYST BY J.SEERY                JANUARY,1973               WTRAN  4
C      WTRAN--WEIGHTS TRANSFORMATION.                                   WTRAN  5
C                                                                       WTRAN  6
      DIMENSION DUMMY(28),DUM(4)                                        WTRAN  7
      COMMON /MDSCL2/ DUMMY,DCON1,DCON2,DCON3,DCON4,DCON5,A,B,P,S,T,DUM WTRAN  8
      WTRAN=S+T*((A+B*XX)**P)                                           WTRAN  9
      RETURN                                                            WTRAN 10
      END                                                               WTRAN 11
$      FORTRAN
CXITEM                                                                  XITEM  1
      SUBROUTINE XITEM(INDEX,XMAX,XMIN,ID,DELX,X,JX,XSF)                XITEM  2
C XITEM            FOR KYST                  JANUARY,1973               XITEM  3
C     JAY R LEVINSOHN     OCTOBER 31,1971                               XITEM  4
C ADAPTED FOR KYST BY J.SEERY                JANUARY,1973               XITEM  5
C                                                                       XITEM  6
C     ITS FUNCTION IS TO PLACE  ELEMENTS OF THE X VECTOR IN THE         XITEM  7
C     PRINT LINE.                                                       XITEM  8
C                                                                       XITEM  9
C                                                                       XITEM 10
      REAL ITEM(101),PTID(100),X(1)                                     XITEM 11
      COMMON /PLTCHR/ PTID,ITEM                                         XITEM 12
C                                                                       XITEM 13
      XNOW=X(INDEX)                                                     XITEM 14
         IF (XNOW-XMAX)610,610,605                                      XITEM 15
  610    IF (XNOW-XMIN)605,615,615                                      XITEM 16
  615    J=(XNOW-XMIN)/DELX + 1.5+XSF                                   XITEM 17
C                                                                       XITEM 18
C     DETERMINE LABEL FOR OBSERVATION                                   XITEM 19
C                                                                       XITEM 20
      IF(IABS(ID)-3)503,501,503                                         XITEM 21
  503 IF(IABS(ID)-2)500,501,500                                         XITEM 22
 500  IF(ITEM(J).NE.PTID(2)) ITEM(J)=PTID(JX)                           XITEM 23
      GOTO 605                                                          XITEM 24
  501 ITEM(J)=PTID(INDEX)                                               XITEM 25
  605 CONTINUE                                                          XITEM 26
      RETURN                                                            XITEM 27
      END                                                               XITEM 28