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